SUBROUTINE POTENL(ICNTRL, MXLMB, MPOTL, LAM, R, P, ITYP) C C ----------------------------------------------------------------- C * MOLSCAT GENERAL POTENL ROUTINE; DESCRIPTION OF FUNCTIONS: C ----------------------------------------------------------------- C VERSION 14 IMPLEMENTS THREE OPTIONS FOR DESCRIBING THE POT'L C 1. POT'L EXPANDED IN ANGULAR FUNCTIONS - SYMMETRIES DESCRIBED C BY MXLAM,LAMBDA INPUT, RADIAL COEFFS DESCRIBED BY INPUT C POWERS AND EXPONENTIALS *OR* VSTAR MECHANISM. C THIS IS THE ORIGINAL MOLSCAT OPTION. C MXLAM.GT.0 MUST BE INPUT AND LVRTP MUST BE .FALSE. (DEFAULT) C 2. POT'L EXPANDED IN ANGULAR FUNCTIONS - PROJECTED VIA VRTP C MECHANISM. SYMMETRIES MAY BE DESCRIBED *EITHER* BY C A.) SYMMETRY DESCRIPTIONS INPUT VIA LAMBDA ARRAY C MXLAM.GT.0 AND LVRTP=.TRUE. *MUST* BE INPUT C B.) SYMMETRY DESCRIPTIONS GENERATED FROM LMAX (MMAX) C MXLAM.LE.0 (DEFAULT) AND LMAX.GE.0 MUST BE INPUT C IF BOTH ARE SPECIFIED (MXLAM.GT.0.AND.LMAX.GE.0) LMAX IS C IGNORED, I.E., CASE (2-A) TAKES PRECEDENCE C ALLOWED ONLY FOR NON-IOS CASES (ITYPE.LT.100) C 3. POT'L IS NOT EXPANDED IN ANGULAR FUNCTIONS (SUITABLE FOR C IOS CALCULATIONS ONLY) AND IS OBTAINED VIA THE VRTP MECHANISM C MXLAM.LE.0 MUST BE SPECIFIED (AND ITYPE.GT.100 IN &BASIS) C C ----------------------------------------------------------------- C * NOTES ON HISTORY OF ROUTINE: C ----------------------------------------------------------------- C C *INTRODUCES XPT(MXPT,MXDIM), XWT(MXPT,MXDIM), INX(MXDIM),* C * NPTS(MXDIM); NPTS IS IN NAMELIST /BASIN/ * C * TO ALLOW GENERAL, MULTI-DIMENSIONAL PROJECTIONS * C ********************************************************** C CORRECTIONS 19 OCT 95 AFFECTING PROJECTION OF ITYPE=3 (SG) C PROJECTION FOR ITYPE=3 ADDED 20 JUL 94 C CODE FOR ITYPE=4 ADDED BY SG 30 JUN 94 (FOLLOWING TRP CODE) C ITYPE=9 INTERFACE ADDED BY JMH 15 AUG 94 C PREVIOUS REVISION DATES 1 FEB 1994 (SG); 3 JAN 1994 (JMH). C ----------------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE LVRTP,ITSAVE,NTERM,NPOWER,A,E,NPT,NPS,NPTS,XPT,XWT, 1 IXFAC,NDIM c fl save Rbig,coeff,nbpt C C ----------------------------------------------------------------- C * NOTES FOR PROGRAMS OTHER THAN MOLSCAT (STORAGE CONSIDERATIONS): C ----------------------------------------------------------------- C THE X() ARRAY CAN BE DEFINED INTERNALLY OR, IF THIS ROUTINE C IS USED WITH THE MOLSCAT CODE, IT IS TAKEN FROM THE C /MEMORY/...,X() STORAGE MECHANISM IN MOLSCAT. C THIS DECK SHOULD BE MODIFIED ACCORDINGLY IN THE STATEMENTS BELOW C AND THE STATEMENTS WHICH FOLLOW STATEMENT NUMBER 2000 C ----------------------------------------------------------------- C ----- NEXT TWO STATEMENTS ARE USED FOR INTERNAL X() STORAGE ----- C ----- ALSO, "X" MUST BE ADDED TO THE "SAVE" STATEMENT ABOVE ----- C PARAMETER (MXX=30000) C DIMENSION X(MXX) C ----- NEXT TWO STATEMENTS ARE FOR MOLSCAT /MEMORY/ MECHANISM----- DIMENSION X(1) COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X C C ----------------------------------------------------------------- C * SPECIFICATION STATEMENTS: C ----------------------------------------------------------------- C MXDIM IS MAX NUMBER OF DIMENSIONS FOR PROJECTION C MXPT LIMITS POINTS PER DIMENSION FOR PROJECTION C MXHERM LIMITS HERMITE POLYNOMIALS FOR VIBRATIONAL PROJECTION PARAMETER (MXPT=96, MXDIM=3, MXHERM=20) C DIMENSIONS FOR LAMBDA AND POWER/EXPONENTIAL TERMS PARAMETER (MXL=360, IXMX=200, IEXMX=200, NPXMX=20) INTEGER CFLAG LOGICAL QOUT,LVRTP, XLAM,LCALC CHARACTER*8 QNAME(10), QTYPE(10) C CHARACTER*6 PNAMES C DIMENSION PNAMES(25),LOCN(25),INDX(25) DIMENSION P(MXLMB), LAM(MXLMB) DIMENSION NTERM(MXL),NPOWER(IXMX),NPUNI(NPXMX),LAMBDA(MXL) DIMENSION A(IXMX), E(IEXMX) DIMENSION H(MXHERM) DIMENSION XPT(MXPT,MXDIM),XWT(MXPT,MXDIM),INX(MXDIM),NPTS(MXDIM) EQUIVALENCE (NPT,NPTS(1)), (NPS,NPTS(2)) C EQUIVALENCE (MXLAM,MXSYM),(LMAX,L1MAX) COMMON/NPOT/NPTL COMMON/ANGLES/COSANG(7),FACTOR,IHOMO,ICNSYM,IHOMO2,ICNSY2 c COMMON/FL/Rbig(200),coeff(200,200) c fl dimension Rbig(500),coeff(500,500) c fl dimension coef(500),b0(500),c0(500),d0(500) C C ----------------------------------------------------------------- C * NAMELIST SPECIFICATION (AND DESCRIPTION OF PARAMETERS): C ----------------------------------------------------------------- NAMELIST/POTL/ RM,EPSIL,MXLAM,MXSYM,LAMBDA,NPOTL,NTERM, 1 A,NPOWER,E,CFLAG,LVRTP,NPT,NPS, 2 IHOMO,ICNSYM,LMAX,L1MAX,L2MAX,MMAX,IVMIN,IVMAX, 3 NPTS,IHOMO2,ICNSY2 C DATA PNAMES/'RM','EPSIL','MXLAM','MXSYM','LAMBDA','NPOTL','NTERM', C 1 'A','NPOWER','E','CFLAG','LVRTP','NPT','NPS', C 2 'IHOMO','ICNSYM','LMAX','L1MAX','L2MAX','MMAX','IVMIN','IVMAX', C 3 'NPTS','IHOMO2','ICNSY2'/ C DATA INDX/25*0/ C C RM - LENGTH SCALING FACTOR; VALUE IN ANGSTROMS C EPSIL - ENERGY SCALING FACTOR; VALUE IN WAVENUMBERS (1/CM) C MXLAM - NUMBER OF POTENTIAL TERMS RETURNED C MXSYM - A SYNONYM FOR MXLAM, RETAINED FOR COMPATIBILITY C LAMBDA - SYMMETRY INDICES FOR POTENTIAL C NPOTL - NO LONGER A RELEVANT INPUT PARAMETER C -------- BELOW DESCRIBE TERMS AS EXPONENETIALS/INVERSE POWERS C NTERM - ARRAY: NTERM(I) IS NUMBER OF TERMS CONTRIBUTING TO P(I) C NTERM(I) .LT. 0 CALLS VINIT/VSTAR FOR POTENTIAL TERM I C A - ARRAY OF PRE-EXPONENTIAL (OR PRE-POWER) FACTORS C FIRST NTERM(1) ELEMENTS REFER TO P(1), C NEXT NTERM(2) ELEMENTS REFER TO P(2) ETC. C NPOWER - ARRAY OF POWERS FOR POTENTIAL TERMS C NPOWER HAS SAME ORDERING AS A C NPOWER(J) .EQ. 0 INDICATES EXPONENTIAL C E - ARRAY OF EXPONENTS: EACH ELEMENT OF THIS ARRAY C CORRESPONDS TO A ZERO IN THE NPOWER ARRAY, C IE E(1) CORRESPONDS TO FIRST ZERO, E(2) TO SECOND ETC. C CFLAG - FLAG FOR SCALING POTENTIAL FOR ITYPE = 5 OR 6: C SET CFLAG=1 IF INPUT A COEFFICIENTS ARE FOR AN C EXPANSION IN C_LM INSTEAD OF Y_LM C -------- BELOW ARE FOR POTENTIALS PROJECTION VIA VRTP MECHANISM C LVRTP - LOGICAL FLAG FOR NON-EXPANDED POTENTIAL: C MXLAM.LE.0 (DEFAULT) FORCES LVRTP=.TRUE. C NPTS - NUMBERS OF GAUSS POINTS FOR PROJECTING POTENTIAL C NPT - EQUIVALENT TO NPTS(1) C NPS - EQUIVALENT TO NPTS(2) C IHOMO - 2 IF POTENTIAL IS SYMMETRIC ABOUT THETA=90, 1 OTHERWISE C ICNSYM - ORDER OF ROTATIONAL SYMMETRY ABOUT PRINCIPAL AXIS C ALSO USED FOR 2ND MOLECULE (I.E., IHOMO2) IN ITYPE=3 C (NOTE: IHOMO & ICNSYM ARE NORMALLY COMPUTED C AUTOMATICALLY OR SET BY THE SUPPLIED VRTP ROUTINE) C -------- BELOW ARE FOR AUTOMATIC GENERATION OF LAMBDA ARRAY; C ONLY FOR PROJECTED POT'LS (LVRTP = TRUE) IF LMAX.GE.0 C LMAX - INCLUDE ALL TERMS FROM 0 TO LMAX IN STEPS OF IHOMO C L1MAX - MAX L1 VALUE FOR MOLECULE-1 (ITYPE=3) C L2MAX - MAX L2 VALUE FOR MOLECULE-2 (ITYPE=3) C MMAX - FOR ITYPE = 5 OR 6, EXCLUDE TERMS WITH M.GT.MMAX C IVMIN, IVMAX - FOR ITYPE = 2, V LOOPS FROM IVMIN TO IVMAX C DATA QTYPE/'LAMBDA =','ABS(MU)=',' MU = ',' L1 = ', 1 ' L2 = ',' L = ',' V = ','V-PRIME=', 2 ' J = ','J-PRIME='/ C C STATEMENT FUNCTION ... F(I)=DBLE(I+I+1) C c define R grid for dvpt of radial coeff do iabc=1,56 Rbig(iabc)=4.6d0+dble(iabc-1)*0.1d0 enddo do iabc=57,118 Rbig(iabc)=10.2d0+dble(iabc-57)*0.25d0 enddo do iabc=119,195 Rbig(iabc)=25.7d0+dble(iabc-119)*1d0 nbpt=iabc enddo IF (ICNTRL.GE.0) GOTO 1000 IF (ICNTRL.EQ.-1) GOTO 2000 WRITE(6,633) ICNTRL,R STOP 1000 continue c write(*,*) 'bonjour12' do i=1,MXLMB do ii=1,nbpt coef(ii)=coeff(ii,i) c write(*,*) coef(ii) enddo call spline(nbpt, Rbig, coef, b0, c0, d0) P(i)=seval(nbpt, R, Rbig, coef, b0, c0, d0) enddo c write(*,*) R,MXLMB,P(1),P(2),P(3) RETURN C C C ****************************************************************** C ** ** C ** CODE BELOW IS FOR AN INITIALIZATION CALL ** C ** ** C ****************************************************************** C 2000 PI=ACOS(-1.D0) C C ------ NEXT TWO STATEMENTS ARE NEEDED IF USING AN INTERNAL ----- C ------ X() ARRAY; I.E. NOT USING MOLSCAT /MEMORY/ MECHANISM C MX=MXX C IXNEXT=1 C ------------------------------------------------------------------ C WRITE(6,634) C INITIALIZE NAMELIST VARIABLES BEFORE READ RM=1.D0 EPSIL=1.D0 MXLAM=0 LMAX=-1 L2MAX=-1 MMAX=-1 IVMIN=-1 IVMAX=-1 IHOMO=2 ICNSYM=1 ICNSY2=1 IHOMO2=1 CFLAG=0 DO 2999 ID=1,MXDIM 2999 NPTS(ID)=0 LVRTP=.FALSE. NPOTL=-1 C C NAMELIST/POTL/ RM,EPSIL,MXLAM,MXSYM,LAMBDA,NPOTL,NTERM, C 1 A,NPOWER,E,CFLAG,LVRTP,NPT,NPS, C 2 IHOMO,ICNSYM,LMAX,L1MAX,L2MAX,MMAX,IVMIN,IVMAX, C 3 NPTS,IHOMO2,ICNSY2 C------------------------------------------------------------------- C ARRAYS FOR NAMELIST SIMULATOR C LOCN(1)=LOC(RM) C LOCN(2)=LOC(EPSIL) C LOCN(3)=LOC(MXLAM) C LOCN(4)=LOC(MXSYM) C LOCN(5)=LOC(LAMBDA) C LOCN(6)=LOC(NPOTL) C LOCN(7)=LOC(NTERM) C LOCN(8)=LOC(A) C LOCN(9)=LOC(NPOWER) C LOCN(10)=LOC(E) C LOCN(11)=LOC(CFLAG) C INDX(11)=4 C LOCN(12)=LOC(LVRTP) C INDX(12)=3 C LOCN(13)=LOC(NPT) C LOCN(14)=LOC(NPS) C LOCN(15)=LOC(IHOMO) C LOCN(16)=LOC(ICNSYM) C LOCN(17)=LOC(LMAX) C LOCN(18)=LOC(L1MAX) C LOCN(19)=LOC(L2MAX) C LOCN(20)=LOC(MMAX) C LOCN(21)=LOC(IVMIN) C LOCN(22)=LOC(IVMAX) C LOCN(23)=LOC(NPTS) C LOCN(24)=LOC(IHOMO2) C LOCN(25)=LOC(ICNSY2) C CALL NAMLIS('&POTL ',PNAMES,LOCN,INDX,25,IEOF) C------------------------------------------------------------------- READ(5,POTL) C IF (NPOTL.NE.-1) WRITE(6,602) NPOTL ITYPE=ITYP-10*(ITYP/10) LCALC=.FALSE. C XLAM=.FALSE. C C CHECK FOR LVRTP OR MXLAM.LE.0, ("UNEXPANDED" POTENTIAL CASE). C IF(MXLAM.LE.0) LVRTP=.TRUE. C WRITE(6,636) IF (ITYP.GT.100 .OR. 1 ITYPE.EQ.1 .OR. ITYPE.EQ.2 .OR. ITYPE.EQ.3. OR. 2 ITYPE.EQ.5 .OR. ITYPE.EQ.6) 2 GO TO 3010 WRITE(6,637) ITYP STOP C 3010 FACTOR=1.D0 ITSAVE=ITYP C CHECK/PROCESS IHOMO/ICNSYM INPUT IF (IHOMO.NE.1.OR.ICNSYM.NE.1.OR.IHOMO2.NE.1.OR.ICNSY2.NE.1) THEN WRITE(6,681) IHOMO,ICNSYM,IHOMO2,ICNSY2 C FOR COMPATIBILITY WITH OLD CODE WHICH USES ICNSYM FOR IHOMO2 ... IF (ITYPE.EQ.3.AND.IHOMO2.NE.1.AND.ICNSYM.EQ.1) ICNSYM=IHOMO2 ENDIF C INITIALIZATION CALL TO VRTP. C IT MAY USE OR SET RM, EPSIL, IHOMO, ICNSYM CALL VRTP(ICNTRL,RM,EPSIL) C C CHECK FOR VALID IHOMO, ICNSYM AFTER CALL TO VRTP IF (IHOMO.NE.1) THEN IF (IHOMO.NE.2) THEN WRITE(6,606) IHOMO STOP ELSE WRITE(6,*) ' &POTL INPUT OR VRTP SPECIFIES', 1 ' HOMONUCLEAR SYMMETRY.' ENDIF ENDIF IF (ICNSYM.LE.0) THEN WRITE(6,*) ' *** POTENL. ILLEGAL &POTL ICNSYM =',ICNSYM STOP ELSEIF (ICNSYM.GT.1) THEN WRITE(6,*) ' ICNSYM INPUT OR FROM VRTP SPECIFIES', 1 ' AXIAL SYMMETRY, ICNSYM =',ICNSYM ENDIF C C C ****************************************************************** C CODE BELOW IS CASE 2 -- PROJECTED EXPANSION USING VRTP C SYMMETRIES INPUT VIA *EITHER* LMAX *OR* MXLAM,LAMBDA C LATTER TAKES PRECEDENCE IF BOTH ARE SPECIFIED. C ****************************************************************** C NPOTL FOR THIS CASE CALCULATED IN CODE BEGINNING AT 3999 C IF (LMAX.LT.0.AND.MXLAM.LE.0) THEN WRITE(6,*) ' *** POTENL. LMAX.LT.0.AND.MXLAM.LE.0' WRITE(6,*) ' YOU MUST SPECIFY SYMMETRIES VIA ONE', 1 ' OR THE OTHER IN &POTL' STOP ELSEIF (LMAX.GE.0.AND.MXLAM.GT.0) THEN WRITE(6,607) LMAX LMAX=-1 ENDIF C C FOR THE COUPLING CASES ALLOWED HERE (1,2,3,5,6) PROJECTING OUT C THE POTENTIAL COEFFICIENTS INVOLVES A QUADRATURE OVER THETA. C GET GAUSS-LEGENDRE QUADRATURE POINTS AND WEIGHTS, C AND CHECK THAT NUMBER OF POINTS IS SENSIBLE. C WRITE(6,635) NPT IF(LMAX.GE.0) THEN MXLM=LMAX+1 ELSE IF(ITYPE.EQ.3) THEN IADD=3 ELSE WRITE(6,638) ITYPE STOP ENDIF MXLM=0 IND=1 DO 3011 I=1,MXLAM MXLM=MAX(MXLM,LAMBDA(IND)) 3011 IND=IND+IADD MXLM=MXLM+1 ENDIF C IF(MXLM.GT.NPT) WRITE(6,648) NPT,MXLM NPT=MAX(NPT,MXLM) C NB ABOVE CODE GUARANTEES THAT NPT.GE.1 IF(NPT.GT.MXPT) GO TO 9400 C NEXP=NPT CALL GAUSSP(-1.D0,1.D0,NPT,XPT(1,1),XWT(1,1)) IF(NPT.NE.NEXP) THEN C NB THIS BRANCH SHOULD NOT OCCUR WITH CURRENT GAUSSP WRITE(6,653) NEXP STOP ENDIF C IF(IHOMO.EQ.2) THEN DO 3014 IPT=1,NPT/2 3014 XWT(IPT,1)=2.D0*XWT(IPT,1) NPT=(NPT+1)/2 WRITE(6,*) ' HOMONUCLEAR SYMMETRY: ONLY HALF OF', 1 ' THE THETA-1 POINTS WILL BE USED' ENDIF C C SET UP OTHER QUADRATURE POINTS AND WEIGHTS FOR PROJECTING C POTENTIAL COMPONENTS FOR ITYPE=1,2,3,5,6 C IF (ITYPE.EQ.3) GO TO 3003 WRITE(6,638) ITYPE STOP C C ITYPE = 3 LINEAR ROTOR - LINEAR ROTOR C 3003 IF (L1MAX.GE.0) THEN LCALC=.TRUE. IF (L2MAX.LT.0) L2MAX=L1MAX MXLAM=0 DO 3034 L1=0,L1MAX,IHOMO DO 3034 L2=0,L2MAX,ICNSYM LMIN=ABS(L1-L2) C 19 OCT 95: LMAX->LTOP BELOW LTOP=L1+L2 DO 3034 LL=LMIN,LTOP,2 MXLAM=MXLAM+1 IF (3*MXLAM.GT.MXLMB) GO TO 9500 LAM(3*MXLAM-2)=L1 LAM(3*MXLAM-1)=L2 3034 LAM(3*MXLAM)=LL ELSE IF (3*MXLAM.GT.MXLMB) GO TO 9500 L1MAX=0 L2MAX=0 DO 3134 I=1,MXLAM LAM(3*I-2)=LAMBDA(3*I-2) L1MAX=MAX(L1MAX,LAMBDA(3*I-2)) LAM(3*I-1)=LAMBDA(3*I-1) L2MAX=MAX(L2MAX,LAMBDA(3*I-1)) 3134 LAM(3*I)=LAMBDA(3*I) ENDIF XLAM=.TRUE. MAXV=MIN(L1MAX,L2MAX) IF (NPS.LE.L2MAX) WRITE(6,648) NPS,L2MAX NPS=MAX(NPS,L2MAX) NEXP=NPS WRITE(6,685) NPS CALL GAUSSP(-1.D0,1.D0,NPS,XPT(1,2),XWT(1,2)) IF(NPS.NE.NEXP) THEN C NB THIS BRANCH SHOULD NOT OCCUR WITH CURRENT GAUSSP WRITE(6,653) NEXP STOP ENDIF C IF SYMMETRIC 2ND MOLECULE, REDUCE POINTS/INCREASE WEIGHTS IF(ICNSYM.EQ.2) THEN DO 3314 IPT=1,NPS/2 3314 XWT(IPT,2)=2.D0*XWT(IPT,2) NPS=(NPS+1)/2 WRITE(6,*) ' HOMONUCLEAR MOLECULE 2:', 1 ' ONLY HALF THE POINTS WILL BE USED' ENDIF C SET UP GAUSS-MEHLER INTEGRATION FOR PHI ON (0,PI) WRITE(6,675) NPTS(3) IF (NPTS(3).LT.MAXV) THEN WRITE(6,*) ' *** POTENL. INSUFFICIENT NUMBER OF POINTS', 1 ' REQUESTED FOR PHI',NPTS(3) WRITE(6,*) ' INCREASED TO',MAXV NPTS(3)=MAXV ENDIF IF (NPTS(3).GT.MXPT) GO TO 9400 FACTL=PI/DBLE(NPTS(3)) TH=-FACTL/2.D0 DO 3342 IX=1,NPTS(3) TH=TH+FACTL XWT(IX,3)=(2.D0*FACTL) 3342 XPT(IX,3)=TH NDIM=3 IF (NDIM.GT.MXDIM) GO TO 9300 NTOT=NPTS(1)*NPTS(2)*NPTS(3) IXFAC=MX-MXLAM*NTOT MX=IXFAC IF (MX+1.LT.IXNEXT) GO TO 9600 C IX=IXFAC C N.B. USE OF YRR MIGHT BE EXPENSIVE; CODE COULD BE MODIFIED C SIMILARLY TO THAT IN IOSB1 PI8=8.D0*PI*PI DO 3329 IP3=1,NPTS(3) DO 3329 IPX=1,NPTS(2) DO 3329 IPT=1,NPTS(1) DO 3329 IL=1,MXLAM L1=LAM(3*IL-2) L2=LAM(3*IL-1) LL=LAM(3*IL) IX=IX+1 3329 X(IX)=YRR(L1,L2,LL,XPT(IPT,1),XPT(IPX,2),XPT(IP3,3))*PI8/F(LL) GOTO 3999 C C ATTEMPT TO PROCESS ITYPE AND POTENTIAL DESCRIPTION NUMBERS 3999 QOUT=.TRUE. NPOTL=MXLAM NQPL=1 WRITE(6,639) ITYPE=ITYP-10*(ITYP/10) IF(ITYPE.EQ.3) GOTO 2003 C WRITE(6,640) ITYPE QOUT=.FALSE. GOTO 2100 2003 NQPL=3 QNAME(1)=QTYPE(4) QNAME(2)=QTYPE(5) QNAME(3)=QTYPE(6) WRITE(6,644) IF(LCALC) THEN WRITE(6,*) ' FOR MOLECULE - 1' WRITE(6,615) L1MAX,IHOMO WRITE(6,*) ' FOR MOLECULE - 2' WRITE(6,615) L2MAX,ICNSYM ENDIF GOTO 2100 C C ******************************************************************* C CODE BELOW IS MAINLY FOR CASE 1 - EXPANDED POTENTIAL C USING NTERM,NPOWER,A,E *OR* VSTAR MECHANISM C ******************************************************************* C HOWEVER, CASE 2 - EXPANDED POT'L PROJECTED FROM VRTP C ALSO RUNS THROUGH THIS CODE, BUT DOES LITTLE C ******************************************************************* C 2100 IF (.NOT.XLAM.AND.NQPL*MXLAM.GT.MXL) WRITE(6,650) MXLAM,NQPL,MXL IX=0 IEX=0 IQ=0 NPX=0 DO 9000 I=1,MXLAM C OUTPUT SYMMETRY DESCRIPTION ONLY IF MXLAM,LAMBDA WERE USED C I.E., LCALC=.FALSE. IF(.NOT.LCALC) THEN WRITE(6,651) I IF(QOUT) WRITE(6,652) (QNAME(J),LAMBDA(IQ+J),J=1,NQPL) WRITE(6,654) ENDIF IQ=IQ+NQPL NT=NTERM(I) C FOR CASE 2, LVRTP=.TRUE. AND WE SKIP PROCESSING IF(LVRTP .OR. NT.EQ.0) GOTO 9000 9000 CONTINUE WRITE(6,663) NPX C IF(NPX.EQ.0) GOTO 9020 -- NOT REQUIRED IN FORTRAN 77 DO 9010 I=1,NPX 9010 WRITE(6,664) I, NPUNI(I) 9020 CONTINUE C C IF LAM HAS NOT YET BEEN FILLED, GET FROM LAMBDA IF (XLAM) GO TO 9050 IF (MXLAM*NQPL.GT.MXLMB) GO TO 9500 DO 9030 I=1,MXLAM*NQPL 9030 LAM(I)=LAMBDA(I) XLAM=.TRUE. C C COMMON RETURN POINT FOR ALL INITIALIZATIONS. C SET VALUES BACK IN CALLING PARAMETERS. C 9050 WRITE(6,665) EPSIL,RM,MXLAM,NPOTL C R=RM P(1)=EPSIL MPOTL=NPOTL MXLMB=MXLAM c fl save radial coefficient do iabc=1,nbpt DO 1700 I=1,MXLMB 1700 coeff(iabc,I)=0.D0 NTOT=1 DO 1710 ID=1,NDIM INX(ID)=1 c BELOW COULD BE ELIMINATED BY 'SAVE NTOT' 1710 NTOT=NTOT*NPTS(ID) c START YPT() INDEX = IX; c IX COUNTS DOWN LAMBDA THEN 1ST DIMENSION, 2ND DIMENSION, ... IX=IXFAC DO 1800 I=1,NTOT WEIGHT=1.D0 DO 1810 ID=1,NDIM COSANG(ID)=XPT(INX(ID),ID) 1810 WEIGHT=WEIGHT*XWT(INX(ID),ID) CALL VRTP(0,Rbig(iabc),SUM) SUM=SUM*WEIGHT c ACCUMULATE CONTRIBUTIONS TO EACH P() DO 1820 IL=1,MXLMB IX=IX+1 1820 coeff(iabc,IL)=coeff(iabc,IL)+SUM*X(IX) c INCREMENT THE INDICES FOR EACH DIMENSION, INX(ID), STARTING W/ 1ST ID=1 1830 INX(ID)=INX(ID)+1 IF (INX(ID).LE.NPTS(ID)) GO TO 1800 c WE REACH HERE IF WE'VE HIT MAX FOR THIS DIMENSION; START NEXT, INX(ID)=1 ID=ID+1 IF (ID.LE.NDIM) GO TO 1830 c IF WE REACH HERE WE SHOULD HAVE COUNTED ALL NTOT ELEMENTS IF (I.EQ.NTOT) GO TO 1800 WRITE(6,*) ' POTENL. ERROR IN PROJECTION. NO. TERMS',I STOP 1800 CONTINUE enddo RETURN C C ********** ERROR CONDITIONS ********** C 9300 WRITE(6,9301) NDIM,MXDIM 9306 FORMAT(/' *** POTENL. PROJECTED POTENTIAL HAS',I3, 1 ' DIMENSIONS, BUT MXDIM=',I3) STOP 9400 WRITE(6,9401) NPT,NPS,MXPT 9401 FORMAT(/' *** POTENL. EITHER NPT OR NPS EXCEEDS MXPT' 2 /' NPT =',I6,' NPS =',I6/' MXPT=',I7) C WRITE(6,649) NPT,MXPT STOP 9500 WRITE(6,9501) MXLMB,MXLAM 9501 FORMAT(/' *** POTENL. DIMENSION OF EXTERNAL LAM ARRAY EXCEEDED'/ 1 ' SIZE PASSED FROM CALLING PROGRAM (MXLMB) =',I8/ 2 ' OFFENDING VALUE OF MXLAM =',I8) STOP C C BELOW IS REACHED IF THERE WAS NOT ENOUGH ROOM IN THE X ARRAY TO C STORE THE PROJECTION COEFFS. IF USING /MEMORY/...X, IT IS C POSSIBLE FOR THE CODE HERE TO OVERWRITE THE LAM ARRAY WITH C COEFFS. HOWEVER, THE PROGRAM SHOULD THEN TERMINATE WHEN CHKSTR C IS CALLED FROM DRIVER AFTER RETURN FROM POTENL INITIALIZATION. 9600 NREQ=MXLAM*(NPT+NPS) MXSTRT=MX+NREQ WRITE(6,9601) NPT,NPS,MXLAM,NREQ,MXSTRT,MXSTRT-IXNEXT+1 9601 FORMAT(' *** POTENL. NOT ENOUGH ROOM FOR PROJECTION COEFFICIENTS'/ 1 ' REQUIRES (',I4,' +',I4,') * ',I4,' =',I8/ 2 ' OF',I8,' ORIGINALLY SUPPLIED IN X(), ONLY',I8, 3 ' WERE AVAILABLE.') STOP c format list 9301 FORMAT(/' *** POTENL. DIMENSION NPXMX EXCEEDED',I6) 633 FORMAT(/' *** ERROR IN POTENL, ICNTRL =',I6,' R =',E16.8) 634 FORMAT(/' STANDARD MOLSCAT POTENL ROUTINE (AUG 94) ', 1 'CALLED FOR POTENTIAL.'// 2 ' /POTL/ DATA ARE --') 602 FORMAT(/' *** POTENL. CURRENT CODE IGNORES INPUT &POTL NPOTL =', 1 I6) 636 FORMAT(/' POTENTIAL IS **NOT** EXPANDED IN ANGULAR FUNCTIONS.'// 1 ' A SUITABLE VRTP ROUTINE MUST BE SUPPLIED.') 637 FORMAT(' *** POTENL. ERROR VRTP NOT SUPPORTED FOR ITYPE =',I6) 681 FORMAT(' INPUT VALUES ARE IHOMO =',I2,', ICNSYM =',I2, 1 ', IHOMO2 =',I2,', ICNSY2 =',I2/ 2 ' THESE MAY BE OVERRIDDEN BY VRTP') 606 FORMAT(/' *** POTENL. ILLEGAL IHOMO =',I6,' FROM &POTL INPUT' 1 ,' OR VRTP') 607 FORMAT(' *** POTENL. IGNORING INPUT &POTL LMAX =',I4, 1 ' IN FAVOR OF MXLAM, LAMBDA() VALUES') 635 FORMAT(I4,'-POINT GAUSSIAN QUADRATURE REQUESTED TO PROJECT', 1 ' COMMON THETA-1 COMPONENT') 638 FORMAT(/' *** POTENL. ILLEGAL LOGICAL PATH. ITYPE =',I6) 648 FORMAT(I4,'-POINT QUADRATURE IS INSUFFICIENT TO PROJECT OUT', 1 ' LEGENDRE COMPONENTS REQUESTED'/' INCREASED TO ',I3, 2 ' ACCORDINGLY') 653 FORMAT(I4,'-POINT GAUSS-LEGENDRE QUADRATURE IS NOT AVAILABLE') 685 FORMAT(I4,'-POINT GAUSSIAN QUADRATURE REQUESTED TO PROJECT', 1 ' LEGENDRE COMPONENTS - MOLECULE 2') 675 FORMAT(I4,'-POINT QUADRATURE REQUESTED TO PROJECT OUT', 1 ' PHI COMPONENTS') 639 FORMAT(/' ANGULAR DEPENDENCE OF POTENTIAL EXPANDED IN TERMS OF') 640 FORMAT(/' *** POTENL. ITYPE =',I4,' CANNOT BE PROCESSED TO', 1 ' DETERMINE THE POTENTIAL SYMMETRY LABLES') 644 FORMAT(' CONTRACTED NORMALISED SPHERICAL HARMONICS, SUM', 1 '(M1,M2,M) C(L1,M1,L2,M2,L,M) Y(L1,M1) Y(L2,M2) Y(L,M)'/ 2 ' SEE RABITZ, J. CHEM. PHYS. 57, 1718 (1972)') 615 FORMAT(' POTENTIAL SYMMETRIES GENERATED FROM LMAX =',I3, 1 ' AND IHOMO =',I2) 650 FORMAT(/' *** POTENL. MXLAM =',I4,' AND NQPL =',I2, 1 ' APPEAR TO EXCEED INTERNAL STORAGE IN LAMBDA(',I5,')'/ 2 ' WILL ATTEMPT TO PROCEED.') 651 FORMAT(/' INTERACTION POTENTIAL FOR SYMMETRY TYPE NUMBER',I4) 652 FORMAT(' WHICH HAS ',6(A8,I3,3X)) 654 FORMAT(1X) 663 FORMAT(/' NUMBER OF UNIQUE POWERS =',I4) 664 FORMAT(' POWER',I3,' =',I4) 665 FORMAT(/' POTENL PROCESSING FINISHED.'// 1 ' ENERGY IN UNITS OF EPSILON =',F15.5,' CM-1'/ 2 ' R IN UNITS OF RM =',F15.5,' ANGSTROMS'// 3 ' MXLAM =',I5/' NPOTL =',I5) END c routine for spline double precision function seval(n, u, x, y, b, c, d) integer n double precision u, x(n), y(n), b(n), c(n), d(n) c c this subroutine evaluates the cubic spline function c c seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3 c c where x(i) .lt. u .lt. x(i+1), using horner's rule c c if u .lt. x(1) then i = 1 is used. c if u .ge. x(n) then i = n is used. c c input.. c c n = the number of data points c u = the abscissa at which the spline is to be evaluated c x,y = the arrays of data abscissas and ordinates c b,c,d = arrays of spline coefficients computed by spline c c if u is not in the same interval as the previous call, then a c binary search is performed to determine the proper interval. c integer i, j, k double precision dx data i/1/ if ( i .ge. n ) i = 1 if ( u .lt. x(i) ) go to 10 if ( u .le. x(i+1) ) go to 30 c c binary search c 10 i = 1 j = n+1 20 k = (i+j)/2 if ( u .lt. x(k) ) j = k if ( u .ge. x(k) ) i = k if ( j .gt. i+1 ) go to 20 c c evaluate spline c 30 dx = u - x(i) seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i))) return end subroutine spline (n, x, y, b, c, d) integer n double precision x(n), y(n), b(n), c(n), d(n) c c the coefficients b(i), c(i), and d(i), i=1,2,...,n are computed c for a cubic interpolating spline c c s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2 + d(i)*(x-x(i))**3 c c for x(i) .le. x .le. x(i+1) c c input.. c c n = the number of data points or knots (n.ge.2) c x = the abscissas of the knots in strictly increasing order c y = the ordinates of the knots c c output.. c c b, c, d = arrays of spline coefficients as defined above. c c using p to denote differentiation, c c y(i) = s(x(i)) c b(i) = sp(x(i)) c c(i) = spp(x(i))/2 c d(i) = sppp(x(i))/6 (derivative from the right) c c the accompanying function subprogram seval can be used c to evaluate the spline. c c integer nm1, ib, i double precision t c nm1 = n-1 if ( n .lt. 2 ) return if ( n .lt. 3 ) go to 50 c c set up tridiagonal system c c b = diagonal, d = offdiagonal, c = right hand side. c d(1) = x(2) - x(1) c(2) = (y(2) - y(1))/d(1) do 10 i = 2, nm1 d(i) = x(i+1) - x(i) b(i) = 2.*(d(i-1) + d(i)) c(i+1) = (y(i+1) - y(i))/d(i) c(i) = c(i+1) - c(i) 10 continue c c end conditions. third derivatives at x(1) and x(n) c obtained from divided differences c b(1) = -d(1) b(n) = -d(n-1) c(1) = 0. c(n) = 0. if ( n .eq. 3 ) go to 15 c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1)) c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3)) c(1) = c(1)*d(1)**2/(x(4)-x(1)) c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3)) c c forward elimination c 15 do 20 i = 2, n t = d(i-1)/b(i-1) b(i) = b(i) - t*d(i-1) c(i) = c(i) - t*c(i-1) 20 continue c c back substitution c c(n) = c(n)/b(n) do 30 ib = 1, nm1 i = n-ib c(i) = (c(i) - d(i)*c(i+1))/b(i) 30 continue c c c(i) is now the sigma(i) of the text c c compute polynomial coefficients c b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2.*c(n)) do 40 i = 1, nm1 b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.*c(i)) d(i) = (c(i+1) - c(i))/d(i) c(i) = 3.*c(i) 40 continue c(n) = 3.*c(n) d(n) = d(n-1) return c 50 b(1) = (y(2)-y(1))/(x(2)-x(1)) c(1) = 0. d(1) = 0. b(2) = b(1) c(2) = 0. d(2) = 0. return end