815 lines
28 KiB
Fortran
815 lines
28 KiB
Fortran
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
|
|
|
|
|