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`

;

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.

At subsequent calls to `POTENL`

, `IC`

is 0, and the routine must evaluate the intermolecular potential

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.

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

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`

.

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.

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.

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)

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**.

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**.

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

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).

Specification of `POTENL`

is the same as for
`ITYP=5`

.

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**.

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.

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.

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.

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`

.

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.

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

- IC = 0 evaluate the potential V
- IC = 1 evaluate dV/dR
- IC = 2 evaluate d2V/dR2

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.

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`

).

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).

Forward to Section 6

Back to Section 4

Back to the Table of Contents