MOLSCAT, Version 14

Tutorial on Input Data

Last update 13 March 1995

This tutorial is designed to help new users of the MOLSCAT code learn about the input data required for a calculation. It walks through some standard or typical cases and provides examples of input decks which can actually be run; the output from such runs should illustrate the purpose and effect of various options. A more complete description of all the input options is contained in the full documentation which should be consulted as necessary. Literature references are also given at the end of the full documentation.

Introduction

MOLSCAT calculates the outcome of nonreactive collisions of a molecule with an atom or with another molecule. A typical application is the calculation of state-to-state cross sections or rate constants for rotational (and possibly vibrational) excitation of the colliding species. MOLSCAT does this by using quantum coupled channel methods. These solve the time-independent Schrodinger equation to obtain a wavefunction for the whole system. The wavefunction is expanded as sums of products of the (asymptotic) rotation and/or vibration wavefunctions of the two colliding species, a partial wave (spherical harmonic) expansion of the angle dependence of the collision coordinate (distance between the two species), and functions of the collision distance. The latter are determined by solving coupled second-order differential equations. The coupling arises from the angle (and vibrational) dependence of the forces between the two species, i.e., the forces which cause rotational and vibrational excitation. Information about the outcome of the collisions is contained in the large-distance behavior of the wavefunction and this is conveniently summarized in terms of collisional S-matrices. Many phenomenological cross sections can be described in terms of the amplitudes and phases of the S-matrix elements and these are the main result of a MOLSCAT calculation.

A difficulty of this approach is the fact that the colliding species have an infinite set of rotational states, and it is necessary to truncate the expansion of the total wavefunction to some finite number of these. In general, the wavefunction will converge by including higher energy basis functions. It is typically necessary to include all levels which are energetically accessible at the collision energy of interest (open channels) plus some energetically inaccessible levels (closed channels). Convergence is slower for more anisotropic interaction forces and for interactions with stronger attractive forces. Note that S-matrix elements are only defined between open channels.

This approach, which is exact except for the truncation of the basis set expansion, is called the close coupling method and is computationally feasible only for systems which have a rather small number of rotational and vibrational levels accessible at collision energies of interest. By making some approximations to the coupling terms it is possible to decouple the problem into smaller blocks and MOLSCAT is equipped to do this for several decoupling schemes. Of these, the coupled states approximation has been found to be reasonably accurate, especially for systems dominated by short-range forces and at higher collision energies. The infinite order sudden (IOS) approximation has been found to be useful in cases where the rotational energy spacings are small compared with the collision energy. For systems where coupled channel methods are not feasible other techniques such as classical trajectory calculations or semi-classical methods should be used, but MOLSCAT does not handle such calculations.

Format of input data

For each calculation it is necessary to provide the program with information about the rotational and/or vibrational wavefunctions which should be included and to specify the intermolecular forces as a function of collision distance and relative orientations. The type of basis functions and the coordinate system needed to describe the interaction potential depend on the kinds of colliding species. Several possible combinations are supported by MOLSCAT, and, before continuing, it is important to ascertain whether the collision system of interest is one of these. The collision types are described by an internal variable ITYP which has the following allowed values.

ITYP = 1
Collision of a rigid linear rotor with a structureless atom
ITYP = 2
Collision of a vibrating diatomic molecule with a structureless atom
ITYP = 3
Collision of two linear rigid rotors
ITYP = 4
Collision of an (a)symmetric top rigid rotor with a linear rigid rotor
ITYP = 5
Collision of a symmetric top rigid rotor with a structureless atom
ITYP = 6
Collision of an asymmetric top rigid rotor with a structureless atom
ITYP = 7
Same as ITYP=2, but the interaction potential can depend on rotational as well as vibrational quantum numbers
ITYP = 8
A structureless atom with a rigid corrugated surface
ITYP = 9
User defined collision type. This requires that the user supply routines to handle the required basis functions, calculate the required matrix elements, etc.

If the collision system of interest is described by ITYP = 1 - 6 the present tutorial should give a reasonable introduction to the required input data. For other collision types the full documentation should be consulted.

Besides expansion basis functions and an interaction potential, it is also necessary to provide input data which specify the collision energies, the method to use for integrating the coupled equations, approximate coupling scheme, optional processing, etc.

In MOLSCAT, the input data are divided into three sets of NAMELIST input. NAMELIST input is not standard FORTRAN but is implemented on most platforms and is too convenient to forego. In general it consists of data cards of the form, &<name> data1=value1, data2=value2, ... &END where <name> is the name associated with the input set; data1, data2, etc. are names of allowed variables in that set; and &END specifies the end of data for this set. &<name> generally must begin in column 2 and all other input must begin in column 2 or beyond, but, being nonstandard FORTRAN, the implementation specifics for a given platform may vary and local documentation should be consulted. One of the advantages of NAMELIST input is the ability to have default values for variables and many MOLSCAT input variables do, in fact, have sensible defaults and need not be included in the input data. The three sets of required data are &INPUT, &BASIS, and &POTL, and they are expected in this order. The &INPUT contains general control variables and &BASIS and &POTL describe the expansion basis set and interaction potential, respectively. Appropriate input for the latter two depends strongly on the kinds of collision species.

The NAMELIST data format is described more completely in Section 1.1 of the complete documentation.

A simple example: rotational excitation of CO by He

The basic input required for nearly all runs will be illustrated here by generating an input deck that might be used to calculate cross sections for rotational excitation of carbon monoxide by helium atoms. This case falls into MOLSCAT's ITYP = 1 category. A low collision energy will be chosen so calculations are relatively cheap.

Describing the interaction potential (&POTL)

In many ways the determination and description of the interaction potential is one of the more difficult parts, and so it will be considered first. For a linear rigid rotor and a structureless atom the interaction depends on the collision distance, R, and on the angle between the collision coordinate and the linear molecule axis, theta. R is measured from the CO center of mass to the He atom. For coupled channel calculations it is generally advantageous to expand the angle dependence of the interaction in terms of some set of angular functions. For ITYP = 1 MOLSCAT uses the complete set of Legendre polynomials, P_L(cos(theta)), where L is the order of the Legendre polynomial. Thus,

```   V(R,theta) = sum-over-L  V_L(R) * P_L(cos(theta)).
```

(For more complicated collision systems there may be more than one sensible angular expansion available, but it is important to use the one expected by the MOLSCAT code; these are described in Section 5.1 of the full documention.)

Consider a simple "asymmetric Lennard-Jones" potential which was used in early work on the CO-He system (S. Green and P. Thaddues, Astrophys. J. 205, 766 (1976)) and which included only P_0, P_1, and P_2 Legendre polynomials. The spherical term was described by a Lennard-Jones function, and the anisotropic terms were described by similar functional forms:

```   V_0(R)/EPS = (RM/R)**12 - 2 * (RM/R)**6

V_1(R)/EPS = -0.03 * (RM/R)**12 + 0.0073 * (RM/R)**7

V_2(R)/EPS = 0.2 * (RM/R)**12 - 0.34 * (RM/R)**6
```

Here `RM` is the position of the minimum, approximately 3.5 Angstroms, and `EPS` is the well depth. approximately 21 1/cm.

The following &POTL parameters can be used to describe this potential. First, MOLSCAT has to know how many angular terms there are; this is specified by `MXLAM` (`MXSYM` is a synonym and may be used interchangeably):

```   MXLAM=3,
```

It needs to know the indices for each of the terms, i.e. the degree of each Legendre polynomials. Note that the Legendre polynomials are each described by a single index; more complex systems may require several indices to describe each angular term. The input variable for these indices is `LAMBDA` and it should contain values for the number of terms specified by `MXLAM`, in this case three symmetries times one term per symmetry:

```   LAMBDA=0,1,2,
```

Next, the radial coefficients, V_L(R) must be specified for each of these angular terms. MOLSCAT has a built in mechanism for potentials which can be described as sums of exponentials and inverse powers of the collision distance, which is the case for this simple potential. The required input variable is `NTERM` which must be given a value for each of the angular symmetries specified by `MXLAM`. The values specify how many exponential and/or inverse power terms are in each of the angular functions. In this case there are two inverse powers in each:

```   NTERM=2,2,2,
```

So far this has specified that there are six terms, two for each of the three angular functions. The inverse power for each of these terms is specified, sequentially, in the variable `NPOWER`:

```   NPOWER=-12,-6,-12,-7,-12,-6,
```

Note that the first two of these correspond to the P_0 term, the next two to the P_1 term, etc. (To specify an exponential, rather than an inverse power, set `NPOWER = 0` for that term.) The coefficient that multiplies each inverse power must also be given, and these are input in the same order in the variable `A`. For the potential described above:

```   A=1.,-2.,-0.03,0.0073,0.2,-0.34,
```

A description of capabilities for potentials expressed as sums of exponential and inverse power terms is given in Section 4.1 of the full documentation.

Finally, note that there are scaling factors for both distance and energy in the above potential and these must be given in the variables `RM` (in units of Angstroms) and `EPSIL` (in units of 1/cm):

```   RM=3.5, EPSIL=21.,
```

In cases where the potential terms do not have scaling factors like this, it is still necessary to describe the units of distance (in terms of Angstroms) and energy (in terms of 1/cm) in the variables `RM` and `EPSIL`. Thus, if distances are in Angstroms and energies are in 1/cm:

```   RM=1.0, EPSIL=1.0,
```

and, if distances are in Bohr radii (atomics units) and energies are in electron volts:

```   RM=0.529177, EPSIL=8065.7,
```

Finally, for the above example, the complete required &POTL input is:

``` &POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2,
NTERM=2,2,2, NPOWER=-12,-6,-12,-7,-12,-6,
A=1., -2., -0.03, 0.0073, 0.2, -0.34,       &END
```

In general NAMELIST variables can be given in any order and arrays can be filled with sequential values as shown. An equivalent input for LAMBDA would be

```   LAMBDA(1)=0, LAMBDA(2)=1, LAMBDA(3)=2,
```

Because NAMELIST input is not part of standard FORTRAN, however, one should check local documentation for specific rules.

It is worth remarking that the angular symmetry terms may be listed in any order and (for most ITYP) they may be repeated. For example, suppose it were convenient for some reason to separate the repulsive and attractive terms in the above potential. One would have three repulsive terms (P_0, P_1, and P_2) and, similarly, three attractive terms, and each would consist of a single inverse power. One would then have the following input:

``` &POTL RM=3.5, EPSIL=21., MXLAM=6, LAMBDA=0,1,2,0,1,2,
NTERM=1,1,1,1,1,1, NPOWER=-12,-12,-12, -6,-7,-6,
A=1., -0.03, 0.2, -2., 0.0073, -0.34  &END
```

which gives exactly the same interaction potential as the previous set. Note that the `NPOWER` and `A` values have been rearrange to correspond to the order in the `LAMBDA` array. It is generally less efficient to do the calculation with repeated symmetries, and this discussion should be considered as an illustration of the flexibility of the program and not as an indication of recommended procedure. On the other hand, there is no penalty in efficiency for specifying the symmetries in a different order.

Describing the rotor basis set (&BASIS)

Consider next the description of the expansion basis set. It is always necessary to specify the collision type. The required variable, `ITYPE`, is obtained by adding to ITYP a number, IADD, which is a multiple of ten and which specifies an approximate method. In particular, for full close coupling IADD = 0, for coupled states, IADD = 20, and for the infinite order sudden approximation IADD = 100. For full close coupling for the current case:

```
ITYPE=1,
```

The rotational basis set for this case requires specifying the rigid rotor quantum numbers, J, and this can be most easily done by setting `JMIN`, `JMAX`, and `JSTEP`, which have rather obvious meanings. Since the default values give `JMIN = 0` and `JSTEP = 1` it is generally only necessary to specify `JMAX`. The rotational energies can be specified by giving the usual spectroscopic rotation constant (in 1/cm). For CO:

```   BE=1.92265,
```

Since rotational energies for a linear rigid rotor are given by `BE*J*(J+1)` it is readily verified that J = 5 is a closed channel at a collision energy of 50 1/cm, which is the energy that will be chosen below for this calculation. Therefore, to include all open channels and one closed channel:

```   JMAX=5,
```

The following &BASIS data provide a complete description

``` &BASIS ITYPE=1, JMAX=5, BE=1.92265,   &END
```

Controlling the calculation (&INPUT)

Consider finally the &INPUT section. It is always necessary to specify the collisional reduced mass, in atomic mass units, in the variable `URED`. If the masses of the molecule and the atom are M1 and M2, respectively, the reduced mass is M1*M2/(M1+M2). For CO-He, with masses of 12, 16, and 4, respectively:

```   URED = 3.5,
```

It is necessary to specify the collision energies (these are total energy, kinetic plus internal rotor energy); generally `NNRG` specifies the number of energies and the values (in 1/cm) are given in the array `ENERGY`. To do calculations at a single, relatively low collision energy:

```   NNRG=1, ENERGY=50.,
```

Values for the starting and ending radial distances for the radial integration are generally required (in units of `RM`, which is specified in &POTL), although the default values are reasonable if `RM` is approximately the distance of the minimum, as it is here. Note that these values are normally used as guesses and the program may adjust them, so it is not crucial to get precise values. `RMIN` should be well into the classically forbidden region, however, and `RMAX` should be far into the asymptotic region. Reasonable values for the current system are:

```   RMIN=0.7, RMAX=10.,
```

It is necessary to choose a numerical method for solving the coupled equations. For many systems a good choice is the modified log-derivative method of Manolopoulos, J. Chem. Phys. 85, 6425 (1986), since it is a good compromise between speed and accuracy and is very easy to control. It is specified by:

```
INTFLG=6,
```

Accuracy of this propagator is controlled by a single input variable, `STEPS`, which is the number of steps per asymptotic wavelength of the radial wavefunction. Good accuracy is achieved with values between 10. and 20. unless the well depth is particularly large compared with the collision energy.

Some other values which should generally be set because the defaults probably do not provide what you desire include the following. `PRNTLV` default is 0 which provides virtually no output; sensible values are probably from 3 to 5; higher values give more detailed partial results. `ISIGPR` default of 0 must be changed in order to print cross sections. `LABEL` can be set to a character string which will be printed with the output to give some identification to the run. Thus:

```  LABEL='CO-He L-J potential of Green and Thaddeus',
PRNTLV=3, ISIGPR=1,
```

It is generally necessary to do calculations for a range of partial waves; these are the orbital angular momenta and correspond classically to an impact parameter. They are labeled by the variable JTOT. JTOT = 0 corresponds to head-on collisions, and large JTOT correspond to glancing collisions. For interactions which decrease sufficiently rapidly with distance, contributions to the cross sections become negligible for sufficiently large JTOT. It is often desirable to specify the range of partial waves by specifying `JTOTL`, `JTOTU`, and `JSTEP` as the lowest and highest JTOT and the increment, respectively. However, the default values will cause the program to start at zero and increment JTOT by steps of one until state-to-state cross sections have achieved reasonable convergence.

Using automatic cross section convergence, the complete data set would be:

``` &INPUT URED = 3.5, NNRG=1, ENERGY=50.,
INTFLG=6, STEPS=10., RMIN=0.7, RMAX=10.,
LABEL='CO-He L-J potential of Green and Thaddeus',
PRNTLV=3, ISIGPR=1,
&END
```

A complete input deck

The &INPUT, &BASIS, and &POTL data described above and in this order, which is just the reverse of the order in which they were discussed, are enough for a complete calculation which can be ascertained by using this as input for a MOLSCAT run.

``` &INPUT URED = 3.5, NNRG=1, ENERGY=50.,
INTFLG=6, STEPS=10., RMIN=0.7, RMAX=10.,
LABEL='CO-He L-J potential of Green and Thaddeus',
PRNTLV=3, ISIGPR=1,
&END
&BASIS ITYPE=1, JMAX=5, BE=1.92265,   &END
&POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2,
NTERM=2,2,2, NPOWER=-12,-6,-12,-7,-12,-6,
A=1., -2., -0.03, 0.0073, 0.2, -0.34,       &END
```

Expanding on the example

It is often desirable to calculate cross sections at several collision energies in order to obtain thermal averages. MOLSCAT can do calculations at up to 100 energies in a single run. The number of desired energies is specified in `NNRG` and the input for `ENERGY` is expanded accordingly. For example,

```   NNRG=4, ENERGY = 120., 110., 100., 90.,
```

It is generally best to put the highest energy first since the step size for some propagators is calculated from the first energy and an input variable, `STEPS` (as in the `INTFLG = 6` example above) and the highest energy will give the most conservative estimate. When calculating multiple energies it must be recalled that the expansion basis set should also be chosen to be adequate for the highest energy as it will then likely be adequate for the lower energies as well.

Most of the propagators can make use of a scratch data set, if supplied, to reduce the work at energies subsequent to the first, and it is generally advantageous to supply such a data set be using the variable `ISCRU`. For example, to request unit(2) as the scratch data set:

```   ISCRU=2,
```

MOLSCAT automatically accumulates state-to-state integral cross sections. However, for many other purposes it is necessary to save the S-matrices on a file for further processing. This is requested by setting the variable `ISAVEU` to the desired unit number. To request unit(13) as the save file:

```   ISAVEU=13,
```

Both `ISCRU` and `ISAVEU` are part of the &INPUT data.

It is generally desirable to check that the parameters chosen for a propagator give adequate accuracy. This can be done easily by running the program again with a different value for `STEPS` (in the case of `INTFLG = 6`); larger values should give higher accuracy. It also is easy to change to a different propagator by simply changing the value of `INTFLG`. The propagators which appear to be most generally useful are `INTFLG = 6`, discussed above, `INTFLG = 8`, and `INTFLG = 3`.

`INTFLG = 8` is particularly useful for systems which have strong long range interactions. It uses the same method as `INTFLG = 6` at short-range, and requires the `STEPS` variable as above. It also requires additional control variables, for switching to an Airy propagator at long-range and for control of the latter. `RMID` is the distance at which it switches methods (as usual, in units of `RM`); values somewhat beyond the potential minimum are recommended. A value for `TOLHI`, used to increase step size in the long-range part, is also required, and values around `TOLHI = 1.05` are useful.

`INTFLG = 3`, the Walker-Light R-matrix propagator is another generally useful method which also takes bigger steps at long-range. It is controlled by `DR`, which is the initial step size in units of `RM`. Values on the order of 0.01 to 0.1 are generally useful. To control the increase of step size at long-range a second parameter is required, `RVFAC`. Step size is increased at distances greater than `RVFAC` times the classical turning distance; values of `RVFAC` around 1.5 are generally useful.

If many calculations, or large calculations will be run, it is useful to try different propagators and vary the tolerance parameters to reach a compromise between accuracy and speed. See Section 2.12 of the full documentation for a more complete description of the different propagators and their required and allowed control parameters. See also Section 6 which describes MOLSCAT facilities for automating convergence testing.

Alternate mechanisms exist for specifying the rotational levels in the &BASIS data. Rather than specifying a range with `JMIN`, `JMAX`, and `JSTEP`, one can specifically list the levels using input parameters `NLEVEL` to specify the number of functions and `JLEVEL` to list the rotational quantum numbers. Thus, the following two data sets are equivalent:

``` &BASIS ITYPE=1, JMAX=5, BE=1.92265,   &END
```

and

``` &BASIS ITYPE=1, NLEVEL=6, JLEVEL=0,1,2,3,4,5,
BE=1.92265,   &END
```

This form can be used if one wanted, for some reason, to specify the basis in a different order. It also has the advantage of allowing the rotor energies to be specified in an input variable, `ELEVEL`, which is not allowed if the basis set list is generated from limits like `JMAX`. For example, one could specify

``` &BASIS ITYPE=1, NLEVEL=6, JLEVEL=0,1,2,3,4,5,
BE=1.92265,
ELEVEL=0., 3.845, 11.5, 22.8, 38.4, 57.0,
&END
```

Values given by `ELEVEL` will override any value supplied by `BE` in this case. This mechanism can be useful, for example, if the rotational energies are perturbed for some reason from rigid rotor values.

For linear rotors excited by atoms (`ITYPE=1`) use of `NLEVEL`, etc. is not generally helpful, but for more complex collision systems, it often provides a desirable option. In general, specifying a value for `NLEVEL` will override values for `JMAX`, etc. and require input values for `JLEVEL`.

More complex potentials - the VSTAR mechanism

For realistic collision systems the interaction potential generally cannot be described by simple exponential and inverse power functions as were used above. MOLSCAT provides two alternative mechanisms for more general interactions. The VSTAR mechanism is useful if the potential has been expanded into the appropriate angular terms and if each of the V_L(R) can be readily calculated. In this case one can request that MOLSCAT call a user supplied VSTAR routine for some or all of the radial functions.

To illustrate the use of this VSTAR mechanism, consider using it for the V_1(R) and V_2(R) terms in the CO-He example. If `NTERM` is set to a negative value for a given angular term, the program will attempt to obtain that coefficient from a routine called VSTAR (with entry VINIT for initialization). The modified &POTL input would be

``` &POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2,
NTERM=2,-1,-1, NPOWER=-12,-6,
A=1., -2.,      &END
```

Note that we still specify three symmetries, and must still give the Legendre function indices. Also, we specify that V_0(R) is still obtained from two inverse power terms. However, the 2nd and 3rd symmetry terms will be obtained by the VSTAR mechanism and the `NPOWER` and `A` values for these have been eliminated from the &POTL input.

The VSTAR routine supplied with MOLSCAT is a dummy. If called it will print an error message and terminate the program. Therefore it must be replaced with a user written routine. It is necessary to relink or reload the program so that the user supplied routines replace those supplied with MOLSCAT.

The following simple FORTRAN routine will provide the same potential as specified for the second and third angular terms of the CO-He potential described above.

```      SUBROUTINE VSTAR(I,R,V)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IF (I.EQ.2) THEN
V= -0.03*(1/R)**12 + 0.0073*(1/R)**7
ELSEIF (I.EQ.3) THEN
V= 0.2*(1/R)**12 - 0.34*(1/R)**6
ELSE
WRITE(6,*) ' VSTAR CALLED FOR ILLEGAL SYMMETRY',I
ENDIF
RETURN
ENTRY VINIT(I,RM,EPSIL)
IF (I.EQ.2.OR.I.EQ.3) THEN
WRITE(6,*) ' VSTAR INITILIZED FOR SYMMETRY',I
RETURN
ELSE
WRITE(6,*) ' VINIT CALLED FOR ILLEGAL SYMMETRY',I
STOP
ENDIF
ENTRY VSTAR1
ENTRY VSTAR2
WRITE(6,*) ' VSTAR1 OR VSTAR2 CALLED; NOT SUPPORTED'
STOP
END
```

Several things should be noted. Calls to VINIT and VSTAR will specify the symmetry number in the first argument and it should be checked to see if calls for this symmetry are expected. There will always be a VINIT initialization call for each symmetry for which the VSTAR mechanism is being used and it will be called before any calls to VSTAR. The present routine simply prints an informative initialization message. One could do other initialization processing, read input data, etc.; the values of `RM` and `EPSIL` from &POTL are made available during the initialization call. The value of R passed to VSTAR will be in units of `RM` and the value of V returned will be multiplied by `EPSIL`; note that these scaling factors are not applied by the VSTAR routine, itself. For some propagators first and/or second derivatives are required and these are provided by entries VSTAR1 and VSTAR2; it is generally not necessary to deal with these cases, but such calls should be trapped, as is done here.

More complex potentials - the VRTP mechanism

The VRTP mechanism is useful if the potential is more readily available as a function of distance and angles, that is, it has not already been expanded in terms of angular functions. In this case the user can supply a VRTP routine. The VRTP routine supplied with MOLSCAT provides a test case for ITYP = 6 and must be replaced by an appropriate routine.

For coupled channel calculations MOLSCAT actually requires that the potential be expanded in terms of angular functions, but can do this expansion automatically via an appropriate Gauss quadrature using results from the VRTP routine. For some IOS calculations, the angle dependent potential may be used directly.

Details for specifying the VRTP routine are given in Section 4.3 of the full documentation. As an example, the following simple FORTRAN code will generate the CO-He potential discussed above.

```      SUBROUTINE VRTP(IDERIV,R,P)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION  P(1)
C
C     THE FOLLOWING NAMED COMMON BLOCK MUST BE PRESENT
C     IT PROVIDES COMMUNICATION WITH POTENL
C     THE FOLLOWING HAS DIMENSIONS FOR VERSION 14
COMMON /ANGLES/COSANG(7),FACTOR,IHOMO,ICNSYM,IHOMO2,ICNSY2
C
C     LEGENDRE POLYNOMIAL FUNCTIONS
P0(CT)=1.D0
P1(CT)=CT
P2(CT)=0.5D0*(3.D0*CT*CT-1.D0)
C
IF (IDERIV.LT.0) THEN
C       THIS IS THE 'INITIALIZATION' CALL
WRITE(6,*) ' VRTP INITIALIZED FOR CO-HE TUTORIAL POTENTIAL'
RETURN
ELSEIF (IDERIV.EQ.0) THEN
C       THIS CALL REQUESTS CALCULATION OF POTL AT R, COS(THETA)
C       GET INVERSE POWERS
RM1=1.D0/R
RM6=RM1**6
RM12=RM6*RM6
C       ASSEMBLE V0, V1, V2
V0=RM12-2.D0*RM6
V1=-0.03*RM12+0.0073*RM6*RM1
V2=0.2*RM12-0.34*RM6
C       COMBINE WITH LEGENDRE POLYNOMIALS
C       FOR ITYP=1 COSANG(1) CONTAINS COS(THETA)
CT=COSANG(1)
P(1)=(V0*P0(CT)+V1*P1(CT)+V2*P2(CT))*FACTOR
C       FACTOR IS AN ITYP DEPENDENT FACTOR SUPPLIED BY MOLSCAT IN /ANGLES/
ELSE
WRITE(6,*) ' VRTP NOT SUPPORTED FOR IDERIV',IDERIV
STOP
ENDIF
RETURN
END
```

Several things should be noted. During the initialization call to VRTP (signalled by IDERIV < 0) R may contain `RM` and P(1) may contain `EPSIL` from the &POTL input; alternatively, values for `RM` and `EPSIL` may be set here and returned during this initialization call. Note that the latter will override the former. The initialization call may do other processing, for example, reading input data or scaling parameters to &POTL input values for `RM` or `EPSIL`.

In subsequent calls to VRTP, IDERIV = 0 requests the potential and IDERIV = 1 or 2 requests the first or second derivative; the latter two cases generally do not have to be supported but such calls should be trapped as indicated. For an evaluation call (IDERIV = 0) R contains the distance, in units of `RM`, and the values for the appropriate angles are in COSANG; use of COSANG by different ITYP is detailed in Section 4.3 of the full documentation. Finally, the calculated value of the potential should be multiplied by FACTOR, an ITYP-dependent geometrical factor which is passed in COMMON/ANGLES/, and the resulting value should be returned in P(1).

To invoke the VRTP mechanism, appropriate values for &POTL data must be specified. In particular `LVRTP` must have the value .TRUE.; it will be set to this value if `MXLAM` is less than or equal to 0 (default is 0). However, to specify projection of appropriate angular terms (Legendre polynomials for ITYP = 1), `MXLAM` and `LAMBDA` generally will be specified as above, so `LVRTP` must be specifically included in the &POTL data:

```   LVRTP=.TRUE.,
```

Also, the number of Gauss points for the numerical integration to project the Legendre components must be specified. A conservative value for this case is:

```
NPTS=7,
```

Thus, the revised &POTL data to invoke the VRTP mechanism for this case is:

``` &POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2,
LVRTP=.TRUE., NPTS=7,  &END
```

Note that `NTERM`, `NPOWER`, and `A` are no longer relevant.

For cases like this an alternative form for specifying the angular terms is also available. One generally desires all the allowed symmetries up to some maximum and it is possible to specify this using `LMAX` instead of `MXLAM` and `LAMBDA`; then the default value of `MXLAM = 0` automatically sets `LVRTP` to .TRUE. so that input equivalent to the above is given by

``` &POTL RM=3.5, EPSIL=21., LMAX=2, NPTS=7, &END
```

Approximate scattering methods

Close coupling calculations rapidly become too expensive with increasing collision energies, especially for nonlinear rotors or for collisions of two rotors. MOLSCAT supports several approximate decoupling methods and two of these, coupled states (CS) and the infinite order sudden (IOS) approximation, have been found to be particularly useful; it is always wise, of course, to test accuracy of an approximate method by comparison with some converged close coupling results.

In general, the only change required to request a decoupling approximation is to modify ITYPE in &BASIS. For the CS approximation one adds 20 and for IOS one adds 100. However, these options do allow certain additional parameters which can be useful in a real calculation and some of these are discussed here.

Coupled States

The coupled states approximation is best viewed in a body-fixed coordinate system in which basis functions with different projections, m, of the rotor momentum on the body-fixed z-axis are decoupled. Separate calculations are thus done for different m blocks at each partial wave. Rotor basis functions must have j>m to be included in a block. As a consequence, contributions to a cross section from j to j' come only from blocks with m < min(j,j') and it is sensible to skip m blocks which do not contribute to cross sections of interest.

The input parameter `JZCSMX` can be used to limit calculations to m blocks less than or equal to `JZCSMX`. As a simple example, in the CO-He calculations at 50 1/cm collision energy and a basis designated by JMAX = 5, the j=5 level is closed; it is sensible to then specify

``` &BASIS ITYPE=21, JMAX=5, JZCSMX=4, BE=1.92265, &END
```

Sometimes only cross sections out of the lowest rotational levels are required and `JZCSMX` can be used to avoid unnecessary calculations.

IOS approximation

The IOS approximation ignores rotational energy spacings compared with the collision energy. It effectively sets all rotational energies to zero. It is therefore not necessary to provide information about rotational energies to MOLSCAT for IOS cases; any information provided will be ignored. Further, state-to-state cross sections in the IOS approximation can be expressed in terms of "generalized IOS cross sections", Q(L), which, for linear rotors, are just the j=0 to j=L cross sections. Because the program calculates these quantities directly and only as a final step uses them to generate cross sections among the levels specified in &BASIS, it is really not necessary to provide a rotational basis set at all for IOS cases. The "generalized IOS cross sections" themselves are often the desired result.

The number of generalized cross sections to compute is determined by input parameters in the &INPUT set. In particular, for a linear rotor excited by an atom, `LMAX` determines the highest Q(L) value which will be computed. (If `LMAX` is not specified, the program will try to chose a value consistent with the rotational levels specified in the &BASIS input, but this is a less desirable method for controlling this parameter.)

The IOS approximation does not actually expand a system wavefunction in terms of rotational basis functions; it is, in fact, equivalent to expanding in the complete, infinite set of basis functions. Rather, it calculates S-matrices as a function of rotor orientations and it is the orientation dependence of the IOS S-matrices which determines the generalized IOS cross sections. This orientation dependence is handled with Gauss quadratures and, rather than specify a basis set, it is necessary to specify the number of points to use in the Gauss quadratures (in fact, the number of points for each dimension of the orientation dependence must be specified for more complex collision systems). This value is specified by the input paramter `IOSNGP` in the &BASIS set. It should generally be somewhat larger than the value of `LMAX`.

For the CO-He example discussed initially, the following input would be appropriate for an IOS calculation:

``` &INPUT URED = 3.5, NNRG=1, ENERGY=50.,
INTFLG=6, STEPS=10., RMIN=0.7, RMAX=10.,
LABEL='CO-He L-J potential of Green and Thaddeus',
PRNTLV=3, ISIGPR=1,
LMAX=10,
&END
&BASIS ITYPE=101, JMAX=5,  IOSNGP=16, &END
&POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2,
NTERM=2,2,2, NPOWER=-12,-6,-12,-7,-12,-6,
A=1., -2., -0.03, 0.0073, 0.2, -0.34,       &END
```

Points to note: `LMAX` was chosen to provide all the Q(L) that will be required to calculate state-to-state cross sections among the levels specified by `JMAX`, and `IOSNGP` was chosen to give reasonable accuracy for this `LMAX`. The rotation constant has been removed from the &BASIS data; this is not necessary, but this value will be ignored in an IOS calculation.

Finally, it should be noted that the WKB propagator `INTFLG=-1` is often sufficiently accurate for IOS calculations and is generally very fast. It essentially approximates a radial integral from the classical turning point to infinity using Gauss quadrature. The number of quadrature points is specified as a triplet of values in the input parameter `NGMP` in the &INPUT data set. The first value is an initial number of points, the second value is an increment, and the final value the maximum number of points. Calculations are done first with NGMP(1) points and then with NGMP(1)+NGMP(2) points, and continue by incrementing with NGMP(2) until two successive S-matrices are converged to a specified tolerance, STEST, or until NGMP(3) is reached. Typical effective values are:

```   NGMP=20, 3, 29,  STEST=0.002,
```

It should be recalled that the WBK propagator is only applicable to problems with a single classical turning point and should therefore be used only when the collision energy is large compared with the potential well depth.