5. User-supplied POTENL routines

As described above, many potentials can be specified either in terms of exponentials and inverse powers or using a user-supplied VRTP or VSTAR routine. It is best to use one of these mechanisms if it can be adapted to your needs. However, in cases where this is not possible or not desirable, you will need to write a complete POTENL routine. This section describes the specification that is required of subroutine POTENL, in order to allow users to write their own potential routines for complicated potentials.

In each run, POTENL is called once for initialisation purposes, and may on this call read any data necessary to specify the potential. It returns information about the terms present in the potential expansion for processing by MOLSCAT. Subsequently, POTENL is called many times to evaluate the potential coefficients for particular intermolecular distances.

The syntax of a call to POTENL is

CALL POTENL (IC,MXLAM,NPOTL,LAMBDA,RR,P,ITYP)

RR and P are DOUBLE PRECISION arguments; the rest are INTEGER. LAMBDA and P are arrays which should be dimensioned at LAMBDA(1), P(1) to switch off FORTRAN array bound checking.

There are three basic types of call to POTENL;

IC = -1

Initialisation call. POTENL is called once with IC = -1, before any other calls to it, to allow it to read any necessary data and set up various parameters. On this call it must return information about the symmetries in the potential.

IC = 0

At subsequent calls to POTENL, IC is 0, and the routine must evaluate the intermolecular potential

IC = 1,2

If POTENL is called with IC = 1 or 2, it must evaluate the IC'th radial derivative of each potential coefficient.

The VIVS propagator is most efficient if radial derivatives of the potential coefficients are available. The WKB integrator also requires first derivatives to start the integration (in locating turning points). Even then, there is an option, controlled by the &INPUT logical variable NUMDER, which, if set to .TRUE., allows the derivatives to be evaluated numerically without making an IC = 1 or 2 call to POTENL. For most other propagators, IC = 1 and 2 are not used, and it is adequate to supply a POTENL routine that prints an error message and exits if it is called with IC = 1 or 2.

The specification of POTENL for initialisation and evaluation calls will be described separately.


5.1 Initialisation

IC

On entry, IC = -1 is the flag for an initialisation call.

MXLAM

On entry, MXLAM specifies the size of the LAMBDA array which has been provided externally. This value is used to check for array bound errors.

On exit, MXLAM must specify the number of distinct terms in the expansion of the potential (i.e., the size of the P array that will be returned by subsequent calls to POTENL.

NPOTL

On entry, NPOTL is not set.

On exit, NPOTL must contain the appropriate value for the collision partners involved, specified by ITYP = MOD(ITYPE,10). The required values of NPOTL are

ITYP = 1, 3, 4, 5, 6: NPOTL=MXLAM

ITYP = 2 and 7: NPOTL=(KMAX+1), where KMAX is the order of the largest Legendre component in the potential. See code in COUPLE for precise use.

ITYP = 8: This is complicated, and depends on the particular version of SURBAS in use. Basically, NPOTL = 2 for the usual version of SURBAS, and NPOTL = 6 or 8 for special normal incidence versions (depending on symmetry of lattice). Existing versions of SURBAS place the necessary value in COMMON/NPOT/, and POTENL picks it up from there.

LAMBDA

Not set on entry.

On exit, the LAMBDA array must contain indices specifying the symmetry of the potential terms that will be used. The LAMBDA array is used differently for each collision type (ITYP). Although it is externally a one-dimensional array, it is conceptually two-dimensional for some collision types, and may be handled explicitly as a two-dimensional array in POTENL if desired. Each element (or column) of LAMBDA corresponds to an element of the P array returned by subsequent calls to POTENL. MOLSCAT does not require that the symmetry terms be supplied in any particular order, but just that the Ith column of the LAMBDA array should correspond to the Ith element of the P array.

Detailed specifications of LAMBDA for the different ITYP values follow. This should be read in conjunction with the documentation for the P array for an evaluation call.

ITYP = 1

Potential is expanded in Legendre polynomials (these are unnormalized, such that P(2)=3*x*x-1).

LAMBDA is a one-dimensional array.
LAMBDA(I) contains the Legendre polynomial order of P(I)

ITYP = 2

Potential is expanded in Legendre polynomials for each (v,v') pair.

LAMBDA is a two-dimensional array LAMBDA(3, ).
LAMBDA(1,I) contains the Legendre polynomial order.
LAMBDA(2,I) contains the vibrational quantum number v.
LAMBDA(3,I) contains the vibrational quantum number v'.

Note that the matrix elements are invariant to exchange of v and v', and only one of these should be supplied.

ITYP = 3

Potential is expanded as contracted spherical harmonics as described by Green, J. Chem. Phys. 62, 2271 (1975).

LAMBDA is a two-dimensional array LAMBDA(3, )
LAMBDA(1,I) is the index corresponding to l1
LAMBDA(2,I) is the index corresponding to l2
LAMBDA(3,I) is the index corresponding to the vector sum l = (l1+l2)

Note that, for identical molecules, the potential is invariant to exchanging l1 and l2; nonetheless, for l1 .ne. l2 both terms MUST be supplied.

ITYP = 4

Potential is expanded in functions that are rotationally invariant contractions of rotation matrices and spherical harmonics as described by Phillips, Maluendes, McLean, and Green, J. Chem. Phys. 101, 5824 (1994).

LAMBDA is a two-dimensional array LAMBDA(4, )
LAMBDA(1,I) is p1, the tensor order for the asymmetric rotor
LAMBDA(2,I) is q1, projection of p1 on the asymmetric rotor axis
LAMBDA(3,I) is p2, the tensor order for the linear rotor
LAMBDA(4,I) is p, the contraction of p1 and p2

ITYP = 5

Potential is expanded as a sum over angular functions

[(Y(l,m) + (-)**m Y(l,-m))] / [(1 + delta(m,0))].

The Y(l,m) are normalized spherical harmonics; note that V(0,0) then differs from the 'spherical average' by a factor. The angle theta is measured from the z-axis, which is the symmetry axis of the top. The angle phi is measured from the x-z plane, which MUST be a reflection plane in the current implementation.

LAMBDA is a two-dimensional array LAMBDA(2, ).
LAMBDA(1,I) is l(I).
LAMBDA(2,I) is m(I).

ITYP = 6

Specification of POTENL is the same as for ITYP = 5.

ITYP = 7

LAMBDA is a two-dimensional array LAMBDA(5,)
LAMBDA(1,I) is the Legendre polynomial order
LAMBDA(2,I) is the vibrational quantum number, v
LAMBDA(3,I) is the rotational quantum number, j
LAMBDA(4,I) is the vibrational quantum number, v'
LAMBDA(5,I) is the rotational quantum number, j'

Note that the matrix elements are invariant to exchange of vj and v'j', and only one of these should be supplied.

ITYP = 8

LAMBDA is a two-dimensional array LAMBDA(2, ).

LAMBDA(1,I) and LAMBDA(2,I) are the reciprocal lattice indices corresponding to potential element P(I). They must be specified in the unique form returned by subroutine ORDER. This last requirement makes sure that the lattice symmetry is properly accounted for.

RR

Not set on entry.

On exit, RR must contain a length scaling factor which is applied to (virtually) all the distances used internally in MOLSCAT. It is specified in Angstroms. It is often convenient to set RR = 1.0 in the initialisation call to POTENL, and to handle everything in Angstroms thereafter.

P

Not set on entry.

On exit, P(1) must contain the energy scaling factor EPSIL to be used internally in MOLSCAT, and subsequent calls to POTENL must return energies in units of EPSIL. EPSIL is specified in 1/cm. It may be convenient to set EPSIL = 1.0 in the initialisation call to POTENL, and to handle everything in 1/cm thereafter. The value given to EPSIL here does NOT affect the interpretation of energy parameters input in &INPUT and &BASIS.

ITYP

On entry, ITYP is the collision type.

If POTENL is coded specifically for a particular value of ITYP, it should check that the correct value is passed, as a precaution against the accidental use of the wrong executable version of MOLSCAT. The value of this parameter should not be changed by POTENL.


5.2 Evaluation

For an evaluation call to POTENL, only the IC, MXLAM, NPOTL, RR and P arguments are passed. LAMBDA and ITYP do NOT contain the values they were given in the initialisation call, so these must be stored internally in POTENL if you need them for an evaluation call.

IC

On entry, IC = 0, 1 or 2 is the flag for an evaluation call:

Note that calls to POTENL with IC = 1 or 2 only occur for the VIVS propagator if IVP, IVPP, ISHIFT or IDIAG is set and for the WKB integrator. Even then, there is an option (controlled by the &INPUT logical variable NUMDER) which, if set to .TRUE., allows the derivatives to be evaluated numerically without making an IC = 1 or 2 call to POTENL. It is thus not altogether necessary for POTENL to cope with IC = 1 and 2 calls, but it should at least trap an attempt to call it this way and print an error message.

RR

On entry this is the intermolecular distance at which the potential is to be evaluated, in units of RM (see the RR argument for an initialisation call to POTENL).

P

Not set on entry.

On exit, P(I),I=1,MXLAM must contain the potential array in the order specified earlier by the LAMBDA array returned by the initialisation call to POTENL. The P array is required in units of EPSIL (see discussion of the initialisation call).


ONWARDOn to Section 6

TOPBack to the Table of Contents


Maintained by Sheldon Green, agxsg@giss.nasa.gov