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

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.

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.

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.

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.

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

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

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

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`

.

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.

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

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.

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.

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.

Back to Documentation Home