The standard POTENL routine currently implements three options:
The parameters that may be input in &POTL NAMELIST are as follows:
RM is the length scaling factor that will be used by MOLSCAT. RM must be supplied in Angstroms (default 1.0), and defines the units of RMIN, RMAX, etc., which are input in &INPUT. Subsequent calls to POTENL will also handle distances in these units.
EPSIL is the energy scaling factor for the potential (default 1.0). EPSIL must be supplied in 1/cm, and subsequent calls to POTENL must return potential coefficients in these energy units. The value given to EPSIL does not affect MOLSCAT's interpretation of energies input in &INPUT and &BASIS.
MXLAM is the number of terms in the potential expansion.
MXSYM is a synonym for MXLAM.
LAMBDA is an array specifying the symmetries of the MXLAM different terms. It is structured differently for the different values of ITYP; see the documentation for an initialisation call to POTENL in Section 5.1 for more details. Note that potential symmetries specified in LAMBDA may be in any order and may be repeated, except for ITYP=2 and 7 for which each symmetry (including vibration-rotation labels) must be unique.
LVRTP is a logical variable. If LVRTP = .FALSE. (the default), the potential is specified in terms of its expansion in angular functions. If LVRTP = .TRUE., subroutine VRTP is called as described in Section 4.3 to evaluate the potential at specified collision coordinates. Note that if MXLAM .le. 0 (which is the default value) LVRTP is forced to .TRUE..
When LVRTP=.TRUE. and MXLAM .le. 0, a nonnegative LMAX (default=-1) must be input to generate the LAMBDA array internally (but only for non-IOS cases, i.e., ITYPE .lt. 100). LMAX specifies the highest Legendre term to include (ITYP=1,2). Similarly, for ITYP=5,6 LMAX and MMAX indicate the highest spherical harmonic to include; if MMAX is not specified (default = -1) MMAX is set equal to LMAX. For ITYP=3, L1MAX and L2MAX are the highest Legendre terms for the two molecules (L1MAX is a synonym for LMAX). Note that IHOMO and ICNSYM (described below) may be used to exclude values not allowed by symmetry.
For ITYP=2, if LMAX is used with LVRTP=.TRUE. to generate the LAMBDA array internally, IVMIN must be given a nonnegative value (default=-1) to indicate the lowest vibrational state. IVMAX, if greater than IVMIN, indicates the highest vibrational state; otherwise IVMAX is set equal to IVMIN.
If angular components of the potential are obtained using the VRTP mechanism, the number of Gauss integration points used in the projection is obtained from the NPTS array. NPTS(1) is the number of points for the theta integration common to all supported ITYPE. NPTS(2) is for the integration over (harmonic oscillator) vibrational functions for ITYP=2; for the integration over phi for ITYP=5,6, and for the integration over theta-2 for ITYP = 3. NPTS(3) is for the phi integration forITYP = 3.NPT is a synonym for NPTS(1) and NPS is a synonym for NPTS(2). NPT should be somewhat larger than the largest value of lambda in the potential expansion. If IHOMO = 2 (see below) only half the points will be used. For ITYP = 5 or 6, NPS should be somewhat larger than the largest value of m in the potential expansion; for ITYP = 2 NPS should be somewhat larger than the largest value of v. If ICNSYM .gt. 1 (see below), then NPS points are used on the range 0 .lt. phi .lt. 2*pi/ICNSYM.
If set to 2, indicates homonuclear symmetry (reflection about theta=pi/2). Default value of 1 indicates heteronuclear symmetry. For potentials which are not expanded in angular functions this value must be set either here or internally in VRTP; for potentials which are expanded in angular functions, it is determined from the LAMBDA array.
For nonlinear molecules denotes symmetry about z-axis; e.g., for ammonia ICNSYM = 3 indicates three-fold symmetry. Default is 1. For potentials which are not expanded in angular functions this value must be set here or internally in VRTP; for potentials which are expanded in angular functions, it is determined from the LAMBDA array.IHOMO2 and ICNSY2 have been added to support a 2nd molecule in version 14. For historical reasons ICNSYM is used to describe homonuclear symmetry of the second molecule in ITYP = 3; however, input of IHOMO2 will produce the same effect. Note: values of IHOMO and ICNSYM in &POTL may be overridden by values set in the initialization call to a user supplied VRTP.
NTERM is an array of MXLAM integers, describing how to evaluate the potential if LVRTP .eq. .FALSE.. Each element of the NTERM array corresponds to one element in the expansion of the potential.If NTERM(I) is less than zero, this element of the potential array will be evaluated using the VSTAR mechanism described below.
If NTERM(I) is positive, this element of the potential array will be evaluated as a sum of NTERM exponential or inverse power terms, specified by A, E and NPOWER.
For each positive value in the NTERM array, there must be NTERM values of NPOWER specified. If an element of NPOWER is 0, that potential term is an exponential described by A and E. If an element of NPOWER is negative, that potential term is an inverse power term described by A and NPOWER.
E is an array of exponential factors. One value must be provided for each zero in the NPOWER array. N.b. values should be negative.
A is an array of pre-exponential and pre-power factors. One value of A must be supplied for each number in the NPOWER array. If NPOWER(I) is 0, the corresponding potential term isA(I) * EXP(E(IE)*R)while if NPOWER(I) is negative, it isA(I) * R**NPOWER(I)
For ITYP=5 and 6, there are two common definitions of the potential expansion: the potential may be expanded either in terms of normalised spherical harmonics Ylm or in terms of renormalised spherical harmonics Clm. Internally, the program uses the former definition. If the coefficients are input in terms of renormalised functions, CFLAG should be set to 1 to instruct the program to multiply each input coefficient by an appropriate factor.
VINIT will be called once for each negative value in the NTERM array during the initialisation call to POTENL. The syntax of the call is
CALL VINIT( I, RM, EPSIL )On entry, I is the index of the potential element in the NTERM (and LAMBDA) arrays, and RM and EPSIL have the values read in in &POTL. None of these arguments should be altered by the routine, but it may be desirable to store their values in local variables so that they are available on subsequent calls. VINIT should also read any data required for the potential evaluation.
VSTAR is called many times during evaluation calls to POTENL. The syntax of the call is
CALL VSTAR( I, RR, SUM )On entry, I is the index of the potential element required, and RR is the intermolecular distance (in units of RM). VSTAR must evaluate the I'th potential element and return its value in SUM (in units of EPSIL).
The specifications for VSTAR1 and VSTAR2 are identical to that for VSTAR, except that they must return the first and second radial derivatives of the I'th potential element respectively. For most propagators, the radial derivatives are not used: in these cases, it is adequate to supply VSTAR1 and VSTAR2 routines that simply print an error message and exit.
The calling sequence is
CALL VRTP(IDERIV,R,V)On an initialization call, IDERIV=-1, and RM(=R) and EPSIL(=V) can be set by the routine or, if read in &POTL, used by the routine. The initialization call may also set the variables IHOMO and ICNSYM in the common block ANGLES:
COMMON /ANGLES/ COSANG(7),FACTOR,IHOMO,ICNSYM,IHOMO2,ICNSY2NOTE: dimensions in this common block were changed in version 14 to anticipate more complex collision cases; earlier versions of VRTP must be modified accordingly when run with version 14 of MOLSCAT.
VRTP should set IHOMO = 1 if the potential has heteronuclear symmetry or IHOMO = 2 if it has homonuclear symmetry (with respect to theta), and ICNSYM to specify the order of any rotational symmetry about the z axis. For ITYPE = 103, ICNSYM is used to specify heteronuclear or homonuclear symmetry for the second linear molecule. Note that, for IOS calculations, IHOMO and ICNSYM would normally be determined automatically by examining the LAMBDA array, but this can be done only if the potential is expanded in angular functions.
Subsequent calls (IDERIV = 0, 1, 2 for the potential and its first and second derivatives respectively) require the routine to return in V the potential at distance R and "angles" specified in the COSANG array.
Some propagators do not use POTENL/VRTP calls with IDERIV .gt. 0, and for these it is adequate to supply a VRTP routine that traps calls with IDERIV .gt. 0, prints an error message, and exits.
For ITYP = 1, 2, 3, 5, and 6, COSANG(1) is cos(theta).
For ITYP = 3 COSANG(2) is cos(theta-2) and COSANG(3) is phi.
For ITYP = 5 and 6, COSANG(2) is phi.
For non-IOS ITYP = 2 cases, COSANG(2) is q, the reduced harmonic oscillator coordinate, such that the ground-state vibrational wavefunction of the diatom is exp(-q**2/2); the code ASSUMES that the molecule is a harmonic oscillator.
For ITYPE = 102, VRTP does not use COSANG(2), and must return a vector of matrix elements <v1|V(R,theta)|v2> in the order expected by POTENL. Appropriate communication is provided by a negative MXLAM in &POTL input, with v1,v2 values supplied as LAMBDA(3*I-1),LAMBDA(3*I),I=1,abs(MXLAM); LAMBDA(3*I-2), which would normally have the order of the Legendre polynomial, is not relevant and is ignored.
A sample VRTP routine is supplied with MOLSCAT which evaluates a simple exponential/inverse power potential for a mock asymmetric top-atom case and which is used in tests of the ITYP = 6 code. This routine must be replaced in calculations using the VRTP mechanism.