MOLSCAT User's Manual


MOLSCAT obtains the intermolecular potential via calls to subroutine POTENL which returns values at input distance R as an expansion in appropriate angular functions. The general-purpose version of POTENL supplied with MOLSCAT is described in this section. Simple potentials, whose coefficients can be expressed as a sum of exponentials and inverse powers, can be input directly as described in Section 4.1. For more complicated potential forms users must supply their own potential routine(s). This may be done by (1) using the VSTAR mechanism described in Section 4.2 to describe (some of) the radial expansion coefficients, (2) using the VRTP mechanism described in Section 4.3 to provide an angle-dependent potential rather than radial expansion coefficients, or (3) replacing the POTENL routine completely. The last method is sometimes the most efficient as it may avoid repeating computations required if the different potential terms are evaluated separately. The specification of the POTENL subroutine is described in detail in Section 5. Versions of POTENL that handle numerical tabulations of the angular coefficients (which often result from ab initio surfaces) by polynomial or spline interpolation are available from SG.

The standard POTENL routine currently implements three options:

  1. The potential is expanded in angular functions appropriate to the collision type. These are specified by a nonnegative input value of MXLAM (the number of angular terms) and ITYPE-dependent input values for LAMBDA to describe the symmetry of each term. The radial coefficients may be described by powers and exponentials or by the VSTAR mechanism (Section 4.2) for each term.

  2. The potential is expanded in angular functions as above, but the radial coefficients are projected via the VRTP mechanism (Section 4.3). The angular terms may be described by either In case of conflict (MXLAM .gt. 0, LMAX .ge. 0) LMAX is ignored, i.e., case (a) takes precedence.

  3. The potential is NOT expanded in angular functions. This is suitable only for IOS calculations. The potential as a function of collision coordinates must be supplied by a VRTP routine (Section 4.3). MXLAM .le. 0 (default) and ITYPE .gt. 100 must be specified.

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.

4.1 Exponentials and inverse powers

The following variables describe the input required when an element of the potential array is to be represented as a sum of exponentials and inverse powers.


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 is
         A(I) * EXP(E(IE)*R)
while if NPOWER(I) is negative, it is
         A(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.

4.2 The VSTAR mechanism

If any element of the NTERM array is negative, subroutines VINIT, VSTAR, VSTAR1 and VSTAR2 must be supplied by the user. However, VSTAR1 and VSTAR2 will be called only if the VIVS or WKB propagator is being used in a mode requiring explicit derivatives (see Section 2.12). The dummy routines supplied with MOLSCAT, if called, will print an error message and terminate the program.

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

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

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.

4.3 The VRTP mechanism

It is often more convenient to specify V(R,angles) explicitly than to expand the potential in angular functions. The general-purpose version of POTENL allows this for IOS calculations for ITYPE = 101, 102, 103, 105, and 106 and for other calculations for ITYP = 1, 2, 3, 5 and 6. The option is invoked by specifying LVRTP = .TRUE. or MXSYM .le. 0; POTENL then expects the potential to be supplied by a routine VRTP. (N.b. specifying MXSYM .gt. 0 and LVRTP=.TRUE. for an IOS case will cause Legendre terms to be projected via the VRTP mechanism, but this is generally a very inefficient way to handle IOS cases.)

The calling sequence is

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:
NOTE: 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.

Forward to Section 5

Back to Section 3

Back to the Table of Contents