MOLSCAT, Version 14

A Tutorial on input data

HTML version; 3 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. 

See also 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 Angstoms, 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,

Note that the symmetries can be listed in any order and, for most ITYP,
they may be repeated. 

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,

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. For example,
if distances are in Bohr radii (atomics units) and energies in electron volts:

   RM=0.529177, EPSIL=8065.7,

Thus, 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. 

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 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. For
CO-He:

   URED = 3.503,

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.503, 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.503, 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 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; 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 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. 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. The following simple FORTRAN
routine will provide the same potential as specified initially for the CO-He
problem. 


      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     COMMUNICATION WITH POTENL IS PROVIDED BY (VERSION 14 DIMENSIONS) 
      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.503, 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 or
until NGMP(3) is reached. Typical effective values are: 

   NGMP=20, 3, 29,

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

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