CCP6 MOLSCAT User manual
Version 14
Jeremy M. Hutson
Department of Chemistry, University of Durham,
Science Laboratories, South Road, Durham, DH1 3LE, England
Electronic mail: J.M.Hutson@durham.ac.uk
and Sheldon Green
NASA Goddard Space Flight Center, Institute for Space Studies,
2880 Broadway, New York, NY 10025, U.S.A.
Electronic mail: agxsg@nasagiss.giss.nasa.gov
MOLSCAT is a general-purpose package for performing quantum nonreactive
molecular scattering calculations. It is written in near-standard FORTRAN 77,
and runs on (among others) DEC VAX, DEC Ultrix, IBM mainframes, IBM RS/6000,
Sun Sparc, Intel iPSC and CRAY computers. This manual describes the
capabilities and modes of use of the program, and explains the input data
required.
JMH and SG would be grateful to be informed of any bugs you encounter,
either in MOLSCAT or in this documentation.
Version of 19 September 1994
Table of Contents
1. Introduction
1.1 Input data format
1.2 Collision types
1.3 Propagators
1.4 Important changes in version 14
1.5 Acknowledging use of MOLSCAT
2. NAMELIST block &INPUT
2.1 General purpose variables
2.2 Print level control
2.3 Angular momentum
2.4 Scattering energies
2.5 S-matrix output file
2.6 Partial cross section output file
2.7 Total cross section output file
2.8 Propagator scratch file
2.9 Pressure broadening calculations
2.10 Characterising scattering resonances
2.11 Surface scattering calculations
2.12 Control of propagators
- Units of length
- Range of integration
- Step size
- DeVogelaere propagator
- R-matrix propagator
- Hybrid log-derivative/VIVS propagator
- WKB propagator
3. NAMELIST block &BASIS
3.1 ITYP = 1
3.2 ITYP = 2
3.3 ITYP = 3
3.4 ITYP = 4
3.5 ITYP = 5
3.6 ITYP = 6
3.7 ITYP = 7
3.8 ITYP = 8
3.9 ITYP = 9
3.10 Decoupling approximations
4. NAMELIST block &POTL
4.1 Exponentials and inverse powers
4.2 The VSTAR mechanism
4.3 The VRTP mechanism
5. User-supplied POTENL routines
5.1 Initialisation
5.2 Evaluation
6. NAMELIST block &CONV
7. Array dimensioning
8. COMMON blocks
9. Subroutines of the MOLSCAT system
10. Machine-dependent features
10.1 Main program
10.2 NAMELIST
10.3 Integer length
10.4 Date and time routines
10.5 CPU time routines
10.6 Floating point underflow
10.7 Linear algebra
10.8 File handling
11. References
1. Introduction
_______________
The outcome of a molecular collision process is described quantum-mechanically
by the S-matrix, which contains information on the probability amplitudes and
phases for the various possible outcomes. Experimental observables, such as
differential and integral cross sections, can be written in terms of S-matrix
elements. For nonreactive collisions a standard computational technique for
solving the time-independent Schrodinger equation involves expanding the
total wavefunction in the (rotation-vibration) basis sets of the colliding
species and a partial wave (spherical harmonic) expansion for the angular
part of the collision coordinate. This results in coupled second-order
differential equations for radial functions which are labeled by the quantum
numbers of the asymptotic basis and the partial waves. The coupling, which
vanishes asymptotically, is due to the intermolecular potential. Truncation
of the infinite asymptotic basis sets (generally on an energy criterion) leads
to the close-coupling method. The scattering S-matrix is obtained by matching
the resulting radial functions at large distances to those which would have
been obtained in the absence of an interaction potential. It is generally
advantageous to transform the asymptotic basis set to a total angular
momentum representation. MOLSCAT constructs the coupled equations appropriate
to the collision problem of interest, and then solves these equations to
obtain S matrices. By default the S-matrices are accumulated to obtain
state-to-state integral cross sections. They may optionally be used to
compute line shape cross sections. They may also optionally be saved to a
data set for processing to other kinds of collision properties.
It is necessary to provide the program with information about the collision
system of interest. In particular, one must specify the types of collision
partners -- e.g., linear rigid rotors, vibrating diatomic molecules,
asymmetric rotors -- and specify which rotational (and possibly
vibrational) levels of each to include in the expansion of the total
wavefunction. It is necessary to describe the intermolecular forces between
the two collision partners. It is also necessary to describe the collision
energies of interest. It may be desirable to specify a less computationally
intensive approximate method. From these data the program can set up the
required coupled equations. MOLSCAT implements several numerical methods for
solving the coupled equations and it is necessary to specify which one of
these to use and to provide some parameters to control numerical accuracy.
Other input parameters may be used to provide finer control of the
calculation, specify optional processing such as the calculation of line
shape cross sections, or options such as saving the S-matrices to a data
set so they may be used to calculate further scattering properties.
These parameters are input to the program in NAMELIST format as described
below. The supported collision types are summarized in Section 1.2 and the
implemented numerical methods for solving the coupled differential equations
are summarized in Section 1.3.
1.1 Input data format
MOLSCAT's main data file is read on channel 5, and consists of 3 (or 4)
blocks of NAMELIST data in the following order:
&INPUT, read in subroutine DRIVER, provides overall control of the
calculation (collision energies, partial wave limits), print options,
optional calculations (line shape cross sections, resonance searches,
output data sets), and choice of and tolerances for the numerical
propagator. Allowed parameters are described in detail in Section 2.
&BASIS, read in entry BASIN of subroutine BASE, describes the collision
type, the rotation-vibration basis set, and dynamical approximations.
Input parameters are described in Section 3.
&POTL, read in the initialisation call of subroutine POTENL if the
general-purpose version of POTENL is used, describes the intermolecular
potential. Input parameters are described in Section 4. For some cases it
may be desirable to substitute a user-supplied POTENL routine (which does
not necessarily read &POTL input data); specifications for this routine are
given in Section 5.
&CONV, read in subroutine CONVRG if the automatic convergence testing option
is selected (&INPUT ICONV = 1). Input parameters are described in Section 6.
Each NAMELIST block starts with &"name", where "name" is INPUT, BASIS, POTL
or CONV as appropriate, and is terminated by &END. Between these delimiters,
the NAMELIST input data consist of "KEYWORD=value" entries in the data
file, where KEYWORD is the variable name. Note that all lines of data must
start in column 2 or later; characters in column 1 will cause errors. In some
implementations of NAMELIST, every value must be terminated by a comma. Many
of the variables have sensible default values, which are used if the variable
is not listed in the data file.
1.2 Collision types
MOLSCAT can perform close-coupling calculations for the following collision
types, controlled by the &BASIS ITYPE input parameter (see section 3):
ITYPE = 1 Atom - linear rigid rotor scattering
ITYPE = 2 Atom - vibrating diatom scattering (rotationally and/or
vibrationally inelastic)
ITYPE = 3 Linear rigid rotor - linear rigid rotor scattering
ITYPE = 4 Asymmetric top - linear molecule scattering
ITYPE = 5 Atom - symmetric top scattering (also handles near-symmetric tops)
ITYPE = 6 Atom - asymmetric top scattering (also handles spherical tops)
ITYPE = 7 A generalised form of ITYPE = 2, allowing the intermolecular
potential matrix elements to depend on the diatom j quantum number
ITYPE = 8 Atom - corrugated surface diffractive (elastic) scattering.
At present, the code is restricted to centrosymmetric lattices,
for which the potential matrices are real.
ITYPE = 9 MOLSCAT calls user-supplied routines to define the coupling case.
The computer time required to solve a set of N coupled equations is nominally
proportional to N**3, and in practice the limit on N is from 30 to 500,
depending on the speed of the computer and the amount of memory available. The
basis sets necessary for converged close-coupling calculations quickly exceed
this limit as scattering energies increase or rotational constants decrease,
particularly for collision types 2 to 7. However, MOLSCAT also provides for
various approximate (decoupling) methods which reduce the dimensionality of
the coupled equations. The methods supported at present are
IADD = 10 Effective potential approximation
IADD = 20 Coupled states (centrifugal sudden) approximation
IADD = 30 Decoupled l-dominant approximation.
IADD = 100 Infinite order sudden approximation.
These approximate methods are invoked by adding the appropriate value of
IADD to the ITYPE values listed above: thus, for example, ITYPE = 22 invokes
the coupled states approximation for atom - vibrating diatom scattering.
Not all these approximations are supported for all collision types, though
themost common ones are. The program will print a warning message and exit
if you request an option that is not supported.
1.3 Propagators
The coupled equations may be solved using any one of several methods,
controlled by the input parameter &INPUT INTFLG (see Section 2.1).
In the absence of special circumstances, INTFLG = 8 is usually a good choice.
INTFLG = 2 DeVogelaere's method. This is a reliable and accurate method,
but is rather slow, especially for large reduced masses or high
scattering energies. It is not well suited to problems with very
long-range potentials.
INTFLG = 3 R-matrix propagator. This is a very stable method. It is very
simple to use, but has relatively poor step size convergence
properties, so that it is not generally suitable for obtaining
very high accuracy. It has largely been superseded by the
modified log-derivative propagators available as INTFLG = 6 to 8
(see below).
INTFLG = 4 The VIVAS hybrid propagator, using the log-derivative method at
small R and the variable interval variable step (VIVS) method at
large R. This is sometimes the most efficient propagator if many
different energies are required, but control of it is rather more
complicated than for other propagators. In particular, it requires
the first and second derivatives of the potential for maximum
efficiency. It sometimes presents stability problems, and is not
generally recommended for inexperienced users.
With appropriate choices of the switchover distance, either the
log-derivative method or VIVS can be used over the whole range of R
if desired. This will usually be useful only in special cases.
INTFLG = 5 Log-derivative propagator of Johnson. This is a very stable and
reliable propagator, and is particularly efficient on vector
processors. However, it has now largely been superseded by the
diabatic modified log-derivative propagator available as
INTFLG = 6 (see below).
INTFLG = 6 Diabatic modified log-derivative method of Manolopoulos. This is
a very efficient and stable method, especially at short range.
INTFLG = 6 automatically detects single-channel cases (including
IOS cases) and uses a more efficient implementation of the
propagator.
INTFLG = 7 Quasiadiabatic modified log-derivative method of Manolopoulos.
This method offers better accuracy than INTFLG = 6 for very
strongly coupled problems, but is relatively expensive (even at
subsequent energies). It is recommended for very strongly coupled
problems only.
INTFLG = 8 Hybrid modified log-derivative Airy propagator of Manolopoulos and
Alexander. Uses the same method as INTFLG = 6 at short range, but
changes to the Airy propagator at long range. This is the
recommended general-purpose propagator for cross section
calculations.
INTFLG = -1 WKB semi-classical solution using Gauss-Mehler quadrature. This is
suitable only for 1-channel problems and hence is usually useful
only for IOS (ITYPE.gt.100) cases.
Note that Gordon's propagator, which was available as INTFLG = 1 in early
versions of MOLSCAT, is not implemented in this version.
All these propagators have options which allow them to use potential matrices
stored at the first total energy when doing calculations at subsequent
energies. For the R-matrix propagator, VIVS, and the modified log-derivative
propagators, some of the remaining work is also avoided at subsequent
energies, so that they may cost only 30% as much as the first. Parameters
which control accuracy of the propagators are in put in &INPUT; details are
given in Section 2.12.
1.4 Important changes in version 14
____________________________________
In modifying MOLSCAT efforts are made to provide compatibility with earlier
versions: ideally, earlier input decks and ancillary routines (e.g., POTENL
routines and routines which process the S-matrix save tape) will produce
identical results. This is not always possible, and changes in version
14 which may require modifications are noted here.
The value of JTOTU required to request automatic checking of partial wave
convergence for state-to-state cross sections has been increased from 999
to 999999 and the default value has also been increased accordingly; see
Section 2.3.
The format of the S-matrix save tape (ISAVEU) has been changed. The number
of open channels for each JTOT, parity case, energy set is now in the record
which describes the set, rather than in the next record which lists levels,
partial waves, and wavevectors; see Section 2.5.
Interpretation of IFEGEN, the flag to specify using the ENERGY values as
kinetic energies for pressure broadening calculations, has been changed; see
Section 2.9. Beginning in version 14, for IFEGEN .gt. 0, values in the
ENERGY list will be retained in the expanded list of total energies ONLY if
they match a total energy required for a requested pressure broadening
calculation. Further, for IFEGEN .gt. 1, calculations for some JTOT, parity
case, energy sets may be skipped if they are not required for the requested
pressure broadening. THIS MAY RESULT IN STATE-TO-STATE CROSS SECTIONS WHICH
ARE INCOMPLETE FOR SOME CASES.
Dimensions in common block CMBASE which contains information about the basis
set (see Section 3) have been significantly increased to accomodate larger
calculations; this may necessitate changes to user-supplied BASE9 (and
related) routines.
Dimensions in the common block ANGLES, which provides communication with
VRTP routines have been increased; see Section 4.3. Existing VRTP routines,
which use this common block for communication, must be modified accordingly.
FINDRM has been completely rewritten to provide a more robust algorithm for
choosing a starting point for the propagation; however, the new routine does
not generally obtain the same value as the earlier versions and this may
cause results to vary slightly from those in earlier versions.
The "FILE=filename" parameters have been removed from OPEN statements. See
Section 10.8
1.5 Acknowledging use of MOLSCAT
Publications that present results obtained with MOLSCAT should cite it as:
J. M. Hutson and S. Green, MOLSCAT computer code, version 14 (1994),
distributed by Collaborative Computational Project No. 6 of the
Science and Engineering Research Council (UK).
The date that the program was last updated is usually held in the variable
PDATE in subroutine DRIVER, and is printed out in each run. PDATE is used to
keep track of minor modifications, but the version number is only incremented
when major changes are made.
2. NAMELIST block &INPUT
________________________
2.1 General purpose variables
LABEL Title for run (up to 80 characters).
Note that this must appear in the data file as a literal string
enclosed in quotation marks, e.g. LABEL='TEST RUN'
URED Reduced mass for collision, m1*m2/(m1+m2), in atomic mass units
(mass of carbon-12 is 12.0)
INTFLG Selects propagator to be used (see Section 1.3 above)
LASTIN (default 1) If LASTIN is set to 1, the run will terminate after the
current scattering calculation is complete. If LASTIN is 0, another
complete set of data will be read once this set has been processed.
ICONV (default 0). If ICONV is set to 1 on input, NAMELIST block &CONV
(see Section 6) will be read and the program will perform automatic
convergence testing with respect to step size, RMIN, RMID or RMAX
(see Section 2.12). If ICONV = 0 (the default), NAMELIST block
&CONV is not read.
MXSIG (default 0). By default, MOLSCAT will accumulate total cross sections
between every pair of levels in the basis set, including closed
channels. This uses a fairly large amount of storage, and is often not
required. If MXSIG is set to a positive integer, only cross sections
between the first MXSIG levels will be calculated.
LMAX (default 0). For IOS calculations, these are the highest L and M
MMAX values for which generalized IOS cross sections, Q(L,M,M') will be
accumulated. For ITYPE=103 (see Section 3.3), LMAX and MMAX identify
the maximum L for rotors 1 and 2 respectively.
IRSTRT (default 0) If not zero requests that a calculation be restarted from
saved S-matrices on unit ISAVEU (see Section 2.5). IRSTRT=3 restarts
after the last complete (JTOT,M,INRG)-set found on ISAVEU. IRSTRT=2
restarts after the last COMPLETED symmetry block (M) found on
ISAVEU; similarly, IRSTRT=1 restarts after a completed JTOT block.
If IRSTRT=-1, the program tries to restart at JTOT=JTOTL taken from
the current &INPUT data. This option is useful if a job terminates
abnormally; it can also be used to change tolerances in the
propagator or to switch to a different propagator.
2.2 Print control
The amount of printed output produced by MOLSCAT is controlled by the integer
variable PRNTLV (default 0); sensible values of PRNTLV vary from 1 (when only
total cross sections are required) to 35 (when debugging). PRNTLV = 3 prints
out some information about each set of coupled equations solved, and complete
S-matrices are printed out if PRNTLV = 11 or more. Voluminous debugging output
starts appearing at PRNTLV = 15.
The number of page throws generated is controlled by the parameter ITHROW
(default 0). If ITHROW = 1, each new S-matrix calculation starts on a new
page.
The printing of total and partial cross sections is controlled by the
parameter ISIGPR (default 0). This must be set to 1 if printing of cross
sections is required (2 to get coupled states cross sections which are
incomplete owing to missing m-values).
2.3 Angular momentum
MOLSCAT has built-in loops over total angular momentum JTOT and "symmetry"
type M. The loop over JTOT is controlled by the three input variables JTOTL,
JTOTU and JSTEP; the program will loop from JTOTL to JTOTU in steps of JSTEP.
Note that the orbital angular momentum label appearing in various decoupling
schemes is treated as a total angular momentum for this purpose.
If total cross sections are required, there is an alternative form of input
flagged by the value JTOTU .GE. 999999 (or JTOTU .lt. JTOTL); this alternative
is obtained by the default values. The end of the loop over JTOT is then
controlled by the variables DTOL, OTOL and NCAC: JTOT will start at JTOTL and
be incremented by JSTEP until NCAC (default 4) successive JTOT values each
contribute less that DTOL (default 0.3) to any diagonal cross section matrix
element and less than OTOL (default 0.005) to any off-diagonal cross section.
(Note that MOLSCAT always gives cross sections in units of square Angstroms.)
This option should be used with care, since you must be sure that the
calculation will converge before the job runs out of time. Note that elastic
cross sections usually converge very much more slowly than inelastic ones.
The loop over symmetry types ("parity cases") M is used differently for
different decoupling schemes. For close-coupling, M takes values 1 and 2 for
odd and even parity label respectively. (For ITYPE = 3, if the two linear
molecules are identical M is further subdivided into exchange symmetries.)
For coupled states calculations, M-1 is the body-fixed projection quantum
number. Under normal circumstances, MOLSCAT generally loops over all possible
values of M for the coupling case concerned (see, however, Section 3.10).
If the input variable MSET is non-zero (default 0), MOLSCAT will skip all
calculations except for M=MSET. This option is particularly useful for
resonance searches (see Section 2.10) and convergence checking (see
Section 6). If MHI is also non-zero (default 0), calculations are performed
for the range of M values from MSET to MHI.
2.4 Scattering energies
The total energies at which scattering calculations are performed are
controlled by the array ENERGY and the variables NNRG, NTEMP, TEMP, DNRG and
EUNITS. If DNRG = 0.0 (the default), calculations will be performed for the
NNRG energies listed in the ENERGY array. If DNRG is not 0.0, calculations
will be performed at NNRG equally spaced energies, DNRG apart, starting at
ENERGY(1). The maximum allowed value of NNRG is currently 100.
The units of ENERGY and DNRG are determined by the integer variable EUNITS
(default 1) as follows:
EUNITS = 1 1/cm
EUNITS = 2 Kelvin
EUNITS = 3 MHz
EUNITS = 4 GHz
EUNITS = 5 eV
EUNITS = 6 erg
EUNITS = 7 Hartrees
EUNITS = 8 kJ/mol
EUNITS = 9 kcal/mol
If a character string (length 4 or less) is input to EUNITS (not allowed
on some machines), an attempt is made to match one of the allowed units
and set the conversion factor.
The input energies are converted immediately into 1/cm using the appropriate
conversion factor, and all energies subsequently printed by MOLSCAT are in
1/cm.
There are special interpretations of the NNRG and ENERGY parameters which may
be invoked when
1) Searching for resonances
2) Calculating line broadening cross sections.
These are described separately in the Sections 2.9 and 2.10.
An alternative form of input is available to facilitate Boltzmann averaging of
cross sections using Gaussian quadrature. This is controlled by the array TEMP
and the variables NTEMP and NGAUSS. If NTEMP is greater than 0 (default 0),
scattering calculations are performed at energies corresponding to NGAUSS
point Gaussian quadratures for the NTEMP different temperatures (in Kelvin)
in the TEMP array. The maximum allowed values are NTEMP = 5 and NGAUSS = 6.
Note that Gaussian quadrature is not a very reliable way of thermally
averaging some types of cross sections, and use of this option is generally
not recommended.
2.5 S-matrix output file
In addition to its main printed output on channel 6, MOLSCAT can also write
out a file containing S-matrices for subsequent processing by another
program (for example, in the calculation of differential cross sections or
transport cross sections). This option is invoked if ISAVEU is not
0 (default 0), and the S-matrices are then written to channel ISAVEU.
Early versions of the program saved this in formatted card image files with
format statements as indicated below. From MOLSCAT version 11 onwards (June
92), the data are saved in unformatted (binary) files which provide better
precision, use less space, and reduce the conversion overhead. In the binary
files the data enumerated below are written in the same order, as single
(logical) records (i.e. single unformatted WRITE statements), except for (7),
which is described more fully below. Beginning with version 14 (Aug 94) NOPEN
is found in record (5) rather than in (6).
Contents of the ISAVEU file are as follows:
1.) LABEL,ITYPE,NLEV,NQN,URED,IP (A80/3I4,F8.4,I4)
LABEL = title of run; character length 80.
ITYPE specifies the collision type.
NLEV is the number of levels in the angular basis set.
NQN is the number of (quantum) descriptors per level.
URED is the reduced mass in atomic mass units.
IP is the MOLSCAT version number (14 in this version).
2.) ((JLEV(I,J),I=1,NLEV),J=1,NQN) (20I4)
JLEV(ILEV,J) are the quantum numbers of level ILEV.
The meaning depends on ITYPE; see Section 3.
3.) NLEVEL,(ELEVEL(I),I=1,NLEVEL) (I4/(5E16.8))
Number and values of the asymptotic level energies (1/cm).
4.) NNRG,(ENERGY(I),I=1,NNRG) (I4/(5E16.8))
Number and values of the scattering energies (1/cm).
(NOTE: beginning in version 14 (Aug 94) NOPEN has been moved from
record (6) to record (5), and programs which process saved
S-matrices should check the value of IP in record (1) to
determine the expected format.)
5.) JTOT,INRG,EN,IEXCH,WT,M,NOPEN (2I4,E16.8,I4,E16.8,I4)
These describe a single scattering calculation:
JTOT is the total angular momentum.
INRG is the index of the energy in the list in 4.).
EN is the total energy (1/cm);it should equal ENERGY(INRG).
IEXCH is the exchange parity for identical molecules
IEXCH = 0: no exchange symmetry
IEXCH = 1: odd exchange symmetry
IEXCH = 2: even exchange symmetry
WT (if nonzero) is a statistical weight of this JTOT, M.
M is the serial index of the symmetry block.
NOPEN is the number of open channels in the S-matrix.
6.) (LEV(I),L(I),WV(I),I=1,NOPEN) (I4/(2I4,E16.8))
LEV(I) is a pointer to JLEV (see above) for the basis set in channel I.
L(I) is the orbital angular momentum for channel I.
WV(I) is the wavevector of channel I (1/Angstroms).
7.) (SREAL(I),I=1,NSQ),(SIMAG(I),I=1,NSQ) (5E16.8)
NSQ=NOPEN*NOPEN
SREAL(I) and SIMAG(I) are the real and imaginary parts of the S-matrix
elements, listed in FORTRAN storage order for the NOPEN*NOPEN matrices
SREAL and SIMAG. In unformatted files (version 11 and later), SREAL and
SIMAG are each written as a single record, listing only the lower
triangle, i.e.
((S(I,J),J=1,I),I=1,NOPEN)
These arrays are written by calls to subroutine SWRITE, and it is
recommended that entry SREAD in that routine can be used to read the
records for subsequent processing.
5.), 6.) and 7.) are repeated for each S-matrix calculated, looping
over INRG (innermost), M and JTOT (outermost):
JTOT = JTOTL,JTOTU,JSTEP
M = 1,MXPAR
INRG = 1,NNRG
5.), 6.), 7.)
END LOOP
MXPAR depends on ITYPE. Note that not every S-matrix will necessarily
exist. S-matrices may be missing from the file either because there are
no open channels for that energy, or because there was an error or
convergence failure in the calculation.
The S-matrix output file is not supported for IOS calculations, and the
output that can be generated for some cases is probably not generally
useful.
2.6 Partial cross section output file
The option to write out partial cross sections to a separate file is not
implemented in this version of MOLSCAT, but could be resuscitated if
necessary.
2.7 Total cross section output file
If ISIGU is greater than 0 (default 0), MOLSCAT will maintain a (direct
access) file containing accumulated cross sections on channel ISIGU. This
file is updated for each energy after each complete JTOT, so it contains
valid information about the run so far even if the program terminates
abnormally.
2.8 Propagator scratch file
All the propagators have options to save energy-independent matrices on a
scratch file to save computation at subsequent energies. This is particularly
advantageous for the R-matrix and VIVS propagators. If ISCRU .gt. 0
(default 0), this save file is created on channel ISCRU and used for any
subsequent energies. It can be large.
Note that this option saves CPU time at the expense of disc I/O. It is
usually advantageous for the R-matrix, VIVS, quasiadiabatic and hybrid Airy
propagators if NNRG is not 1, but for the DeVogelaere and basic
log-derivative propagators it may not save resources overall unless the
potential itself is very expensive to calculate. This is particularly true
on fast machines such as CRAYs. Note, however, that ISCRU .gt. 0 saves
considerable time for single-channel IOS cases with INTFLG=6, since the
computer time for these is usually dominated by potential evaluations.
When searching for and characterising resonances, a propagator scratch file
created in one run may be used in a subsequent run to save work. This is
described in Section 2.10 below.
The propagator scratch file option is not implemented in the parallel version
of MOLSCAT.
2.9 Pressure broadening cross sections
Calculations of pressure broadening and shifting cross sections are supported
for ITYPE = 1, 2, 3, 5, 6, 7, 11, 12, 15, 16, 17, 21, 22, 25, 26, 27, 31, 101,
and 105. The pressure broadening option is invoked if NLPRBR .gt. 0
on input (default 0); the value of NLPRBR specifies the number of
spectroscopic lines for which pressure broadening calculations are to
be carried out. The lines themselves are specified by the array LINE, of
4*NLPRBR elements. Each successive quartet of elements of LINE indexes the
rotor levels involved in the transition, according to the elements of the
JLEVEL array described in Section 3. The 4 values indicate JA, JB, JA', JB',
where JA-JB is the spectroscopic transition and the primes indicate post-
collision values. Note that the "diagonal" values JA=JA' and JB=JB' describe
widths and shifts of isolated lines, and off-diagonal values describe
collisional transfer of intensity.
Pressure broadening calculations require S-matrix elements involving the
initial and final (spectroscopic) levels at the same kinetic (not total)
energy. Generation of the necessary total energies is facilitated by the input
parameter IFEGEN. If IFEGEN .gt. 0 (default 0), the program will treat the
input ENERGY values as kinetic energies, and generate the necessary total
energies for the lines requested. If IFEGEN is 0, only those requested
lines which can be constructed from the total energies actually specified by
NNRG,ENERGY will be calculated. Through version 13 of MOLSCAT, input values
of relative kinetic energies in ENERGY were retained in the the list of
generated total energies whether they were needed for a requested pressure
broadening calculation or not; beginning with version 14 only total energies
needed for the requested pressure broadening calculations are retained.
Further, beginning with version 14, specifying IFEGEN .gt. 1 will suppress
calculations for individual JTOT, INRG, M ("parity case") which do not
contribute S-matrices needed for the requested pressure broadening. WARNING:
some subset of the state-to-state cross sections may be incorrect (missing
some partial wave/parity contributions) if this last option is used.
The tensor order of the spectroscopic transition (i.e., 1 for dipole
transitions, 0 for isotropic Raman scattering) is specified in the input array
LTYPE. If the default value (-1) is found, LTYPE is calculated from ABS(JA-JB).
2.10 Characterising scattering resonances
Scattering resonances and predissociating states of Van der Waals molecules
appear as characteristic features in the energy-dependence of S-matrices. In
MOLSCAT, a run to find or characterise a resonance is flagged by setting the
parameter KSAVE .gt. 0 (default 0). This option is NOT supported for IOS
calculations (ITYPE .gt. 100). Setting KSAVE .gt. 0 has principal effects:
1) The S-matrix eigenphase sum is calculated at each energy, and a short
summary of the eigenphase sums is written to channel KSAVE.
2) The output to channel ISAVEU is modified to be suitable for use with the
program RESFIT. S-matrices, K-matrices and eigenphase sums are written to
channel ISAVEU unformatted.
A MOLSCAT run to characterise a resonance must deal with only one total
angular momentum JTOT and symmetry type M in a given run. Thus, JTOTL and
JTOTU must be the same, and the parameter MSET must be set greater than 0
to select the particular symmetry type required.
A special option is available to assist with searching for narrow resonances.
If NNRG is given as a negative integer of the form -5*N, with N integral and
the energies themselves specified by ENERGY(1) and DNRG, MOLSCAT will perform N
groups of 5 equally spaced calculations. After each group, the program will
try to interpret the 5 eigenphase sums as the "tail" of a resonance, and
estimate the position of the resonance centre. This estimate is then used as
the starting point for the next group of 5 energies. This option can be useful
in searching for a resonance if a reasonably good estimate of its energy is
already available, and may succeed in converging on a narrow resonance from as
far as 10**5 widths away. However, convergence is NOT guaranteed, and it is
not usually useful to do more than 3 sets of 5 energies in a single run.
Resonance search calculations should usually be done with either the R-matrix
propagator or one of the modified log-derivative propagators, since many
different energies are always necessary. It is thus always desirable to use
the ISCRU option to save energy-independent matrices on a scratch file. For
the special case of resonance searches, with JTOTU = JTOTL and MSET .gt. 0,
the ISCRU file from one run may be used as input for the next run at a
different set of energies. This option is invoked by setting ISCRU negative
for the subsequent run(s); the program then expects to find a suitable input
file on channel |ISCRU|. It reads the header on this file to check that it
contains valid information, and then proceeds with "subsequent energy"
calculations for all the energies requested.
2.11 Surface scattering calculations
Atom - surface scattering calculations are not plagued by complications
arising from total angular momentum and parity. MOLSCAT uses the internal
loops over JTOT and M to loop over the polar angle theta (measured from the
surface normal) and the azimuthal angle phi (measured relative to the surface
reciprocal lattice vector g1). The loop over theta is controlled by the input
parameters JTOTL, JTOTU, JSTEP, THETLW and THETST, while that over phi is
controlled by MXPHI, PHILW and PHIST. The logic used is equivalent to:
DO 840 JTOT = JTOTL,JTOTU,JSTEP
THETA = THETLW + THETST*JTOT
DO 830 M = 1,MXPHI
PHI = PHILW + PHIST*(M-1)
..
..
Scattering calculation for angles THETA, PHI
..
..
830 CONTINUE
840 CONTINUE
Note that, for surface scattering, subsequent energies have a rather
non-intuitive meaning, because the threshold energies (K+G)**2 depend on the
parallel component of the momentum K. MOLSCAT interprets subsequent energies
as having the SAME parallel momentum as the first energy, but a different
perpendicular momentum. This corresponds to a change in the polar angle as
well as the scattering energy.
2.12 Control of propagators
Convergence of coupled channel calculations with respect to integration range
and step size is very important; lack of convergence can give very poor
results, whereas unnecessarily conservative tolerances can waste large
amounts of computer time. It is always advisable to conduct careful step
size convergence tests on a small basis set before embarking on any major
set of calculations with MOLSCAT.
Units of length
_______________
MOLSCAT operates with units of lengths specified by RM, which is set in
subroutine POTENL (see Section 5) and which may be input via the &POTL data
set (see Section 4). Variables which control the integration, such as RMIN,
RMAX, etc. (described below) are in units of RM. Note that RM itself is
specified in units of Angstroms, and users might find it convenient to have
all distances in Angstroms by setting RM=1.0. As another example, to specify
all distances in bohr radii (atomic units) one should specify RM=0.529177.
Range of integration
____________________
The mechanisms used to choose the starting and finishing points for
integrating the coupled equations are the same for all propagators. In the
simplest case, the input variables RMIN (default 0.8) and RMAX (default 10.0)
are supplied and used exactly as input. (Note that the default values are
appropriate when RM is approximately the distance of the potential minimum.)
However, in certain circumstances RMIN and RMAX are modified from their input
values:
1) If IRMSET .gt. 0 on input (default 9), the program will estimate a suitable
distance to start propagating, using the criterion that the wavefunction
amplitude in all channels should be less than 10**(-IRMSET) at the starting
point. A crude semiclassical estimate of the wavefunction is used to make
this estimate, which is calculated separately for each JTOT, M (at the
first energy only).
2) If there is a centrifugal barrier present, the program checks that all open
channels are classically accessible at the RMAX requested, for all energies
in the ENERGY list. If necessary, RMAX is increased to the value of R at
the furthest classical turning point found in the centrifugal potential
for any energy. However, RMAX is never decreased from the input value.
For special purposes, this action can be suppressed by setting IRXSET = 0
on input (default 1). However, this is not recommended for general cross
section calculations.
Thus, the input value of RMAX is the smallest distance at which the
propagator may terminate. The DeVogelaere and R-matrix propagators also
check some aspects of convergence internally, and may propagate beyond
RMAX under certain circumstances.
It should be noted that propagating out to the centrifugal turning point is
NOT necessarily adequate, particularly when calculating elastic integral
cross sections or line shift cross sections. It is also important to
realise that a value of RMAX which is adequate for low JTOT may not be
adequate for high JTOT, and that calculations using too small a value of
RMAX may appear to be converged with respect to the partial wave sum when
they are not actually converged.
Step size
_________
The DeVogelaere and the basic log-derivative propagators are
wavefunction-following methods, and use a constant step size controlled by the
parameter STEPS (default 10.0). This is interpreted as the number of steps per
half-wavelength for the open channel of highest kinetic energy in the
asymptotic region. A value between 10 and 20 is usually adequate, unless the
depth of the potential well is large compared to the scattering energy.
The two modified log-derivative propagators are actually potential-following
methods, but they nevertheless use the same mechanism (STEPS) to determine
step size, for compatibility with the basic log-derivative code. For these
propagators, values of STEPS around 5.0 are usually adequate.
The R-matrix and VIVS propagators are potential-following methods, so the de
Broglie wavelength is less important. Instead, the required initial step size
is input explicitly in the variable DR (default 0.02). However, in both
propagators the step size may subsequently be modified, as described for the
individual propagators below.
DeVogelaere propagator (INTFLG = 2)
______________________
STABIL The DeVogelaere method is potentially unstable in a classically
forbidden region, because the exponential growth of closed channel
wavefunctions can lead to a loss of linear independence of the
solutions. To avoid this, the DeVogelaere propagator re-imposes
linear independence every STABIL steps. The default (5.0) is usually
adequate.
STEST is an integration tolerance parameter used for two purposes:
1) The K-matrix obtained by matching to the asymptotic boundary
conditions should be exactly symmetric, but this is never the case
because of numerical rounding errors. The program counts and prints
out the number of pairs of symmetry-related off-diagonal elements
which differ by more than STEST.
2) The S-matrix is calculated at successive values of R, starting at
RMAX and incrementing by approximately half a wavelength each time.
This is repeated until two successive S-matrices differ by no more
than STEST in any element, up to a maximum of 20 attempts.
R-matrix propagator (INTFLG = 3)
___________________
RMID (default 9999.0). For R .lt. RMID, the R-matrix propagator uses a
constant step size of DR. For R .gt. RMID, the step size is DR*R/RMID,
so that larger steps are taken in the long-range region. If RVFAC is set
greater than the default value of 0.0, RMID is set automatically to
RVFAC*RTURN, where RTURN is an estimate of the position of the
classical turning point, calculated for each JTOT and parity case.
(RVIVAS is a synonym for RMID.)
SHRINK If the integer variable SHRINK is 1 on input (the default), the
R-matrix propagator will monitor closed channels in the long range
region and accumulate a phase integral for each. If a semiclassical
estimate based on this phase integral indicates that the wavefunction
amplitude in the channel concerned has dropped below 10**(-IRMSET),
the channel will be removed from the basis set for the remainder of
the propagation.
VTOL controls a convergence criterion applied to the eigenvectors of the
potential matrix (default VTOL = 1.D-6). The R-matrix propagator will
not attempt to apply boundary conditions until the eigenvector is a
(permuted) unit matrix to within tolerance VTOL for two successive
steps. This test is applied only to eigenvectors which are
asymptotically non-degenerate.
MAXSTP is the largest number of steps allowed for the R-matrix propagator
(default 10000). If convergence has not been achieved after this many
steps, the program prints an error message and exits.
Hybrid log-derivative/VIVS (VIVAS) propagator (INTFLG = 4)
_____________________________________________
RVIVAS (default 9999.0) is the distance at which the switchover from the
log-derivative method to VIVS takes place (but see also RVFAC below).
A value of RVIVAS just inside the potential minimum is usually most
efficient. If RVIVAS is more than RMAX, the entire propagation will
use the log-derivative method; however, the preferred method of
achieving this effect is to use INTFLG = 5, since this uses
considerably less storage. Conversely, if RVIVAS is less that RMIN,
the entire propagation will use VIVS (not recommended). It must be
remembered that MOLSCAT itself may adjust RMIN and RMAX as described
above; RVIVAS can if necessary be set to a large negative or
positive number to ensure the desired effect. (RMID is a synomym
for RVIVAS.)
RVFAC (default 0.0). If RVFAC .gt. 0.0, RVIVAS is adjusted automatically
for each JTOT / parity case, to be RVFAC times an estimate of the
classical turning point. A value of RVFAC around 1.1 is usually
suitable, and is recommended for cross section calculations
using INTFLG = 4.
The log-derivative propagator has no variable parameters except STEPS, which
determines the step size as described above. However, control of the VIVS
propagator is quite complicated. VIVS operates with both a variable interval
length and a variable step length. A single diagonalising transformation is
used over the whole of each interval, which may consist of several steps. The
algorithms used to determine interval and step sizes will be discussed
separately.
The step size and interval size algorithms used by VIVS attempt to take the
longest step which will give the required accuracy at each point in the
integration. However, the optimum step size at one scattering energy is not
necessarily safe at another, and VIVS can sometimes give problems if the groups
of energies in a particular run are not chosen with care. In particular, one
should avoid: 1) A first energy which is close to (or above) a channel
threshold and a subsequent energy which is far below it. 2) A first energy
which is far above a channel threshold and a subsequent energy which is close
to it.
1) Interval sizes
TOLHI VIVS accumulates perturbation corrections to the wavefunction as it
propagates, and uses these to obtain a suitable length for the next
interval. The input parameter DR described above is used as the size
of the first interval, and subsequent interval lengths are obtained
using the input tolerance TOLHI (default 0.001); the criterion is
that some functional t of the perturbation corrections should be not
greater than TOLHI over any interval. Within an interval, t is
checked against TOLHI at each step, and a new interval is started
(with a new diagonalising transformation) if it appears likely to
exceed TOLHI over the next step.
2) Step sizes
Each interval is divided into IALPHA steps (default IALPHA = 6). Within an
interval, the step sizes increase geometrically, with each step being a factor
alpha larger than the previous one. The quantity alpha may be specified in
either of 2 ways:
i) If the logical input parameter IALFP is .FALSE. (the default), alpha
increases linearly from ALPHA1 (default 1.0) at the starting point for
VIVS to ALPHA2 (default 1.5) at RMAX.
ii) If IALFP is .TRUE. on input, the program starts with an initial value of
ALPHA1, and adjusts this as it propagates. ALPHA2 is then not used.
3) Automatic step and interval lengths
A useful special option is obtained by setting IALPHA = 0 on input. Intervals
then consist of a variable number of steps, and the decision to start a new
interval is based solely on the magnitude of the perturbation corrections; a
new interval is started (and a new diagonalising transformation obtained)
whenever the quantity t approaches TOLHI. The initial step size is taken from
DR, and subsequent steps use a criterion based on TOLHI.
The IALPHA = 0 option can be very efficient, and often requires remarkably few
intervals/steps to produce converged results.
4) Perturbation corrections
There are several logical input variables which control the extent to which
VIVS calculates and uses perturbation corrections to the wavefunction. The
three variables IV, IVP and IVPP (defaults .TRUE., .FALSE. and .FALSE.
respectively) control the calculation of perturbation corrections due to the
potential itself (IV) and to its first (IVP) and second (IVPP) derivatives.
The perturbation corrections thus calculated are used in calculating interval
sizes, but are included in the wavefunction only if IPERT is .TRUE. (the
default). If ISHIFT is .TRUE. (default .FALSE.), the second derivative is used
to shift the reference potential to give the best fit to the true potential.
For production runs, IV, IVP, IVPP, IPERT and ISHIFT should usually all be
.TRUE. This may be forced by setting IDIAG = .TRUE., which overrides any
.FALSE. values for the individual variables.
If IVP, IVPP or ISHIFT is set, VIVS will require radial derivatives of the
potential. These are supplied properly for simple potentials by the standard
version of POTENL described below, but for some potentials they can be
difficult to evaluate. In this case, the input variable NUMDER may be set
.TRUE., in which case the necessary derivatives are calculated numerically,
and POTENL is never called with IC .gt. 0 (see Section 5).
5) Other variables
ISYM (default .TRUE.). If ISYM is .TRUE., the R-matrix is forced to be
symmetric at the end of each interval. This is usually advisable for
production runs.
XSQMAX (default 1.D4) controls the application of perturbation corrections
to deeply closed channels. If a channel is locally closed by more
than XSQMAX reduced units, perturbation corrections for it are not
calculated. The default should be adequate.
Log-derivative propagators (INTFLG = 5, 6 and 7)
_________________________
The only relevant input variables are RMIN, RMAX and STEPS as described at
the beginning of this section.
Hybrid log-derivative/Airy propagator (INTFLG = 8)
_____________________________________
The modified log-derivative (INTFLG = 6) propagator is used from RMIN to RMID.
In general, values of RMID somewhat beyond the distance of the minimum are
recommended. If RVFAC .gt. 0.0 on input, RMID is calculated as described
above; values of 1.2 to 1.4 are generally useful. If RMID is less than RMIN,
or if it is determined that RMID-RMIN is smaller than one step, the
log-derivative propagator will be skipped. The step size for this propagator
is controlled by STEPS as described above unless IABSDR (default 0) is set to
1, in which case the step size is taken from the input DR.
The Airy propagator is used to propagate from RMID to RMAX. If RMID is greater
than RMAX this propagator is not called. By default the initial step size is
taken as the value calculated for the modified log-derivative part, but it
can be further controlled by the following variables.
DRAIRY (default -1.0) can be used to input an absolute step size.
If DRAIRY is less than zero, the initial step size is taken
from the log-derivative propagator.
TOLHI (default 0.001). If less than 1, the propagator uses this as a
tolerance to adjust the step size (by a perturbation method) to try
to maintain this accuracy. If greater than 1, the step size is
increased by this factor in each interval. Values around 1.03 to
1.07 are generally useful.
POWRX (default 3). If TOLHI is less than 1, this is the inverse power used
in estimating the step size from perturbation calculations. The
default value is usually adequate.
WKB propagator (INTFLG = -1)
__________________________________________________
The INTFLG=-1 propagator is suitable only for single channel cases and is
used almost exclusively for IOS calculations. It does the WKB integral
for the phase shift by Gauss-Mehler quadrature. It is controlled by
input variables NGMP and TOL
NGMP Dimension 3. N-point Gauss integration is done starting
with N = NGMP(1), incrementing by NGMP(2), until NGMP(3).
Starting with the 2nd pass the phase shift is compared
with the previously calculated value until it has
converged to within a tolerance specified by TOL.
TOL (default = 0.001) is the convergence tolerance for the
WKB phase shift.
3. NAMELIST block &BASIS
________________________
The parameters input in &BASIS specify the collision type, the quantum numbers
and energies of the levels to be used in the basis set, and the dynamical
approximations (if any) to be used in constructing the coupled equations. The
input parameters common to all collision types are as follows:
ITYPE specifies the collision type, as described in Section 1.1. It is
constructed from ITYP, which specifies the nature of the colliding
particles, and IADD, which specifies the dynamical approximations to
be made (if any). The values allowed for ITYP are as follows:
ITYP = 1 Linear rigid rotor + atom.
ITYP = 2 Diatomic vibrotor + atom.
ITYP = 3 Linear rigid rotor + linear rigid rotor.
ITYP = 4 Asymmetric rigid rotor + linear rigid rotor
ITYP = 5 Symmetric top rigid rotor + atom.
ITYP = 6 Asymmetric rigid rotor + atom.
ITYP = 7 Diatomic vibrotor + atom, with different diagonal
potentials for different rotational levels.
ITYP = 8 Atom + rigid corrugated surface.
ITYP = 9 User-supplied code for other collision types.
Approximate methods are specified as follows:
ITYPE = ITYP Full close-coupling
ITYPE = ITYP + 10 Effective potential method
ITYPE = ITYP + 20 Coupled states approximation
ITYPE = ITYP + 30 Decoupled l-dominant approximation
ITYPE = ITYP + 100 Infinite order sudden approximation
NLEVEL If NLEVEL .gt. 0, it generally indicates the number of levels to be
included in the basis set, as specified by the JLEVEL array.
If NLEVEL = 0, the quantum numbers for the levels to be included in
the basis set will be calculated internally. For exceptions and more
details see discussion for individual collision types below.
JLEVEL is an integer array specifying the quantum numbers of the levels to
be included in the basis set. The JLEVEL array is structured
differently for different values of ITYP as described below. JLEVEL
is a one-dimensional array of dimension 4000, although it is
conceptually two-dimensional for ITYP .gt. 1.
ELEVEL is an array of NLEVEL molecular level energies, corresponding to the
levels in the JLEVEL array. If all the elements of ELEVEL are 0.0, or
NLEVEL is 0, the energies are calculated using molecular constants as
described below, and ELEVEL is not used. ELEVEL is currently
dimensioned at 1000.
EUNITS is an integer specifying the energy units for ELEVEL and molecular
constants input in &BASIS. EUNITS is interpreted in the same way as
the EUNITS parameter of the &INPUT NAMELIST block (see Section 2.4).
The parameters used to control the automatic generation of basis sets and
energy levels from quantum numbers and molecular constants will be described
separately for the different values of ITYP below. The methods of specifying
the molecular levels to be included are independent of any dynamical
approximations employed, so that the information given for each ITYP below is
applicable to ITYPE = ITYP, ITYP+10, ITYP+20 and ITYP+30. For IOS cases (ITYPE
= ITYP+100) the required input is generally the same, except that rotational
energies are not required for IOS calculations.
The input data are used to generate a ((JLEV(LEV,I),LEV=1,NLEV),I=1,NQN)
table which contains a description of the quantum numbers and other necessary
data for each 'level'. The number of columns, NQN, for each ITYPE and the
meaning of each column are described below for the different collision types.
Note that the last column, JLEV(LEV,NQN), is always a pointer to the index of
level LEV in the cross section array (SIGMA).
3.1 ITYP = 1
For atom - linear rigid rotor systems, the JLEVEL array is usually generated
from input parameters JMIN, JMAX and JSTEP, which have the obvious meanings.
For special purposes, NLEVEL may be set greater than zero and a list of NLEVEL
j values supplied in the array JLEVEL.
If JLEVEL is calculated from JMIN, JMAX and JSTEP, the energy levels are
generated from the input parameters BE (rotational constant), ALPHAE
(vibrational dependence of BE), and DE (centrifugal distortion constant). The
rotational constant actually used is, of course, simply BE-0.5*ALPHAE, and
this value may be input directly in BE with ALPHAE = 0.0 (the default) if
preferred.
If NLEVEL .gt. 0, the energy levels may be specified in the input array
ELEVEL. If ELEVEL is not supplied, they are generated from BE etc. as above.
NQN = 2 and JLEV(LEV,1) = j, the rotational quantum number.
3.2 ITYP = 2
For atom - vibrating diatom scattering, the JLEVEL array must contain a list
of (j,v) pairs specifying the rotational and vibrational quantum numbers of
the levels to be included in the basis set. There is no option for generating
this list automatically.
If the level energies are not supplied in ELEVEL, they are generated from the
input parameters WE (omega(e)), WEXE (omega(e)x(e)), BE (rotational constant)
and ALPHAE (alpha(e)), using the standard formulae. DE, the centrifugal
distortion constant is assumed to be independent of vibrational level here.
Note that energies calculated from WE, etc. are set with the zero of energy
at v = j = 0.
For IOS calculations the vibrational quantum numbers are selected from the
JLEVEL array, and any rotational quantum numbers supplied are ignored. If
vibrational energies are supplied in ELEVEL, the first one encountered for
each vibrational state is kept. If all ELEVEL are zero, energies are
generated from WE and WEXE. If the input NLEVEL is negative, a user-supplied
routine GET102 is called to set NLEVEL, JLEVEL and ELEVEL values.
NQN = 3; JLEV(LEV,1) = j, the rotational quantum number, and JLEV(LEV,2) = v,
the vibrational quantum number.
3.3 ITYP = 3
For linear rigid rotor - linear rigid rotor scattering, the JLEVEL array must
contain a list of (j1,j2) pairs describing the rotational quantum numbers of
the two molecules. If NLEVEL = 0 on input, JLEVEL is generated from J1MIN,
J1MAX, J1STEP, J2MIN, J2MAX and J2STEP, which have the obvious meanings for
molecules 1 and 2 respectively.
If the two molecules are identical, the basis functions corresponding to
(j1,j2) and (j2,j1) are indistinguishable. In this case, the input variable
IDENT (default 0) should be set to 1. Only distinguishable pairs (i.e. J1 .gt.
J2) should be supplied. It is also necessary to specify the statistical weights
to be applied to the different symmetry combinations when calculating cross
sections. The statistical weights for anti-symmetric and symmetric
combinations of (j1,j2) and (j2,j1) may either be specified explicitly in the
WT array as WT(1) and WT(2) respectively, or (if both WT(1) and WT(2) are
zero) may be calculated from the single nuclear spin input in SPNUC. Note
that MOLSCAT will not perform a scattering calculation for any symmetry type
with a statistical weight of zero.
Energy values may be specified in ELEVEL or (if all ELEVEL are zero) will be
computed from BE(I), ALPHAE(I) and DE(I). If IDENT=1 and the values for I=2
are zero, the program sets them equal to the values for I=1.
NQN = 4; JLEV(LEV,1) = j1, the rotational quantum number of the first
molecule, JLEV(LEV,2) = j2, the rotational quantum number of the second
molecule, and JLEV(LEV,3) = j12, the vector sum of j1 and j2. N.b.
JLEV(LEV,4), the index of the cross section, depends only on j1 and j2,
not on j12.
3.4 ITYP = 4
The basis functions here are combinations of asymmetric top rotor
wavefunctions and linear rotor wavefunctions. In general, input of the
asymmetric top functions follows the capabilities for ITYP = 6 which are
described more fully below. If rotor constants (A(1), B(1), and C(1)) are
specified, the asymmetric top wavefunctions will be calculated, subject to
limitations in JMIN, JMAX, JSTEP (synonyms are J1MIN, J1MAX, and J1STEP),
EMAX, and ISYM. Otherwise, wavefunctions will be input from unit IASYMU;
if NLEVEL .gt. 0, only this many levels will be input, else input to
end-of-file. If functions are generated from rotor constants and if a valid
IASYMU unit number is also specified, the fucntions will be written to this
unit in a format suitable for future input.
The asymmetric top wavefunctions are combined with linear rotor functions
specified by J2MIN, J2MAX, J2STEP; energies of the linear rotor are obtained
from the rotation constant which must be supplied as BE(2). If specified,
ALPHAE(2) and DE(2) are also used as for ITYP = 1. If EMAX .gt. 0 is
specified, it is used to limit basis functions to energies (asymmetric top
plus linear rotor) less than EMAX.
In order to allow for finer control of the basis set, specifying NLEVEL
.lt. 0 will try to obtain from IASYMU the basis functions specified by
(JLEVEL(3*i-2) = J1, JLEVEL(3*i-1) = ITAU, JLEVEL(3*i) = J2), i = 1,
abs(NLEVEL). The asymmetric top functions for all requested J1, ITAU must
be available on unit IASYMU. In this case it is also possible to specify the
energies (asymmetric top plus linear rotor) in the ELEVEL array; otherwise
energies are taken from IASYMU for the asymmetric top functions and computed
from BE(2) for the linear rotor.
NQN = 8; JLEV(LEV,1) = j12, the vector sum of asymmetric top and linear
rotor angular momenta, JLEV(LEV,2) = j2, JLEV(LEV,3) = j1, JLEV(LEV,4) =
ITAU, JLEV(LEV,5) = PRTY (a parity designation discussed under ITYP = 6),
JLEV(LEV,6) and JLEV(LEV,7) are pointers to the asymmetric top wavefunctions.
3.5 ITYP = 5
Rotor basis sets may be taken either from NLEVEL, JLEVEL input or from quantum
number restrictions. In the former case, NLEVEL should be greater than 0 and
rotor levels should be specified in JLEVEL(3*i-2)=j, JLEVEL(3*i-1)=k, and
JLEVEL(3*i)=parity, i=i,NLEVEL. Note that parity is indicated by an even or an
odd integer for even and odd combinations of (+k) and (-k); for k=0, only even
parity is allowed.
Alternatively, the basis set may be specified as including levels for j from
JMIN to JMAX: JSTEP is not used. The additional input variable KSET (default
0) ay be used to limit the k values included: if KSET is negative, only
levels ith k=|J2MAX| are included; if KSET is zero or positive, only levels
with k ess than or equal to KSET are included. If ALL k levels are required,
KSET should be set to at least JMAX (999 recommended).
The corresponding energy levels may either be input in ELEVEL, or they may be
calculated from the zeroth order near-symmetric top equation using rotational
constants supplied in the input variables A, B and C. Note that these MUST
correspond, respectively, to moments of inertia about the x, y, and z axes
used in the description of the potential (see below); they may not be in the
descending order standardly adopted in spectroscopy. In particular, the
program requires that the z axis be the symmetry axis of the top (or the axis
of near symmetry) and that the x-z plane be a reflection plane. For a
symmetric top, then, A=B and C is the unique constant.
NQN = 4; JLEV(LEV,1) = j, the rotational quantum number, JLEV(LEV,2) = k, the
projection of j on the symmetry axis, and JLEV(LEV,3) is a parity designation.
3.6 ITYP = 6
Asymmetric top wavefunctions are expanded in terms of symmetric top functions.
The necessary functions may either be calculated internally or supplied as
data.
If all three rotation constants (A, B and C) are specified, the program
constructs the asymmetric top Hamiltonian and diagonalises it. Note that
(as for ITYP = 5) A, B and C MUST correspond to the x, y and z axes used in
the potential expansion, and are NOT necessarily in descending order of
magnitude. In this case centrifugal distortion constants DJ, DJK, and DK
may also be used; they are input as ROTI(7), ROTI(8), and ROTI(9),
respectively. Rotational functions will be calculated for J = JMIN, JMAX,
JSTEP and kept only if they have energy below EMAX (if EMAX.gt.0 has been
specified). Functions will also be kept only if they meet the symmetry
restrictions imposed by ISYM, which is described below. If a valid unit
number is specified for IASYMU the resulting wavefunctions will be written
in the format described below for input from IASYMU.
If rotation constants are not input, the program (subroutine SET6) will try
to obtain the wavefunctions by reading data from unit IASYMU (specified in
&BASIS, default = 5, the standard input unit). Each level is described on
IASYMU by J,ITAU,ELVL (2I5,F15.10) where J is the rotational quantum number,
ITAU is an index that is not actually used, and ELVL is the rotor energy (in
1/cm unless &BASIS EUNITS has specified otherwise). For each level, NK=2*J+1
expansion coefficients are required (corresponding to k=-J,-J+1,...,+J) and
these must follow the J,ITAU,ELVL card in format (6F12.8) on (NK+5)/6
subsequent cards. Coefficients need not be normalized, but they will be
checked for valid symmetries (by subroutine IPASYM). Note, again, that the
coefficients must correspond to symmetric top functions in the coordinate
system of the potential energy (cf. ITYP = 5 above); in particular, the
rotational constants used to diagonalize the rotor hamiltonian must
correspond to the x, y, and z-axes in which the potential is described.
If IASYMU is the standard input unit, data cards should follow the &BASIS ...
&END set and precede the &POTL set; in this case a positive value of NLEVEL
must also be set in &BASIS, to specify the number of expected rotor functions.
If IASYMU is an auxiliary unit, it is possible to read levels until an
end-of-file condition occurs, and specifying NLEVEL .lt. 0 will cause SET6 to
do this. If NLEVEL = 0, input functions will be skipped if they do not fall
in the range JMIN, JMAX, JSTEP (if JMAX .gt. 0 is specified) or if they are
above energy EMAX (if EMAX .gt. 0 is specified). It is possible to specify a
subset of the levels on IASYMU by setting NLEVEL to a negative value. In
that case, input functions will be skipped unless JI and ITAU match values in
JLEVEL(2*I-1),JLEVEL(2*I),I=1,abs(NLEVEL). In this case indexing order in the
SIGMA (cross section) array will correspond to the ordering in JLEVEL.
Note that for IOS calculations (ITYPE = 106) rotor wavefunctions are used
only to calculate state-to-state cross sections from the "generalized IOS"
cross sections. For this case input of rotor wavefunctions is currently
limited; in particular, generation of rotational wavefunctions from the
rotational constants is not supported and they must be explicitly supplied as
data on IASYMU if state-to-state cross sections are desired.
Asymmetric top functions are characterised by two symmetries:
(1) the k values involved are EITHER even OR odd (corresponding to + or -
symmetry with respect to a C2 rotation about the z axis)
(2) the functions are of the form
sum_k c_k (D_mk + epsilon D_m,-k)
where epsilon is either +1 or -1.
Internally, the program assigns one of four symmetry types to each asymmetric
top function (in subroutine IPASYM):
PRTY even/odd k epsilon
0 even +1
1 even -1
2 odd +1
3 odd -1
If the monomer functions are calculated from rotational constants (but not if
they are input explicitly from IASYMU), the symmetry types that are actually
included in the basis set may be restricted using the input variable ISYM. The
selection criteria that are possible in MOLSCAT version 11 and later are much
more general than in earlier versions.
In the current version, ISYM is interpreted bitwise; a particular asymmetric
rotor function is included if it passes ALL the appropriate tests.
If bit 0 is set, odd k functions are excluded
If bit 1 is set, even k functions are excluded
If bit 2 is set, functions with epsilon*(-1)**j = -1 are excluded
If bit 3 is set, functions with epsilon*(-1)**j = +1 are excluded
If bit 4 is set, functions with degeneracy 1 are excluded
If bit 5 is set, functions with degeneracy 2 are excluded
If bit 6 is set, functions with degeneracy 3 are excluded
If bit 7 is set, functions with degeneracy .gt. 3 are excluded
To construct ISYM, start with 0 (including all functions) and add 2**n for
each class of functions to be excluded.
To decide which states are coupled, consider the symmetries (l,m) present in
the expansion of the intermolecular potential:
(1) If only terms with m even are present, functions with k even are not
coupled to those with k odd.
(2) If only terms with (l+m) even are present, functions with epsilon*(-1)**j
even are not coupled to those with epsilon*(-1)**j odd.
If either or both of these symmetries is present, it is most efficient to do
separate calculations for the (two or four) different sets.
For the special cases of J=0 and coupled states with K=0, there are additional
restrictions on the functions that are coupled: j+j'+l must be even, and
epsilon must be conserved. Do not be misled by these into believing that the
symmetries hold in more general cases.
For asymmetric rotors, all functions are in principle of degeneracy 1.
However, near-degeneracies can occur, and the program interprets these as
degenerate if the levels (from diagonalisation) are within a tolerance coded
in subroutine SET6C (usually 1.D-8). To be safe, do not set any high-bit
flags for asymmetric rotors.
The flags that test degeneracy are intended for use with spherical tops, and
allow A, E and F levels to be included selectively.
The levels to be included may also be restricted using the variable EMAX if
desired.
NQN = 6; JLEV(LEV,1) = j, the rotational quantum number, JLEV(LEV,2) is an
index to distinguish different asymmetric top levels with the same J-value
(e.g., the 'tau' label; this value is actually not used by the program),
JLEV(LEV,3) is the PRTY value discussed above, JLEV(LEV,4) and JLEV(LEV,5)
are pointers to coefficients of the asymmetric top wavefunctions expanded
in terms of symmetric top functions.
3.7 ITYP = 7
The basis set for ITYP = 7 is specified in the same way as that for ITYP = 2.
NQN = 3; JLEV(LEV,1) = j, the rotational quantum number, and JLEV(LEV,2) = v,
the vibrational quantum number.
3.8 ITYP = 8
For atom - corrugated surface scattering, the JLEVEL array must contain a list
of (g1,g2) pairs specifying the surface reciprocal lattice vectors to be
included in the basis set. If NLEVEL = 0 on input, these will be calculated
from J1MAX, J2MAX and EMAX; g1 loops from -J1MAX to +J1MAX, and g2 loops from
-J2MAX to +J2MAX. However, not all these functions will be included in the
basis set, as described below.
The dimensions and symmetry of the surface unit cell are specified in the
input array ROTI. ROTI(1) and ROTI(2) are the lengths of the real space
lattice vectors in Angstroms (irrespective of the units of length used
elsewhere). ROTI(3) is the lattice angle in degrees.
To understand MOLSCAT's methods of generating basis sets for surface
scattering, it is important to realise that the operation is carried out in
two steps. At the time that the &BASIS NAMELIST block is read, MOLSCAT sets
up a master list of basis functions which MAY be included in subsequent
scattering calculations. This takes place before entering the loops over
angles and energies, so must be angle- and energy-independent. For each
basis function G = (g1,g2) covered by the loops over g1 and g2 described
above, MOLSCAT checks that the parallel kinetic energy
(0.5*hbar**2/URED)*|G|**2 is less than the input parameter EMAX, and
includes the basis function in its master list only if it passes this check.
Note that the test involving EMAX is bypassed if NLEVEL .gt. 0 and the JLEVEL
array is read in explicitly.
Subsequently, for particular values of THETA, PHI and ENERGY, MOLSCAT
calculates the parallel component of the incident momentum K, and selects
from the master list those basis functions for which
(0.5*hbar**2/URED)*|K+G|**2 is less than the input parameter EMAXK. This
occurs even if the JLEVEL array is specified explicitly in &BASIS data.
Thus, only those basis functions which satisfy both the EMAX and EMAXK
criteria will ultimately be included in the basis set. Since the EMAXK
criterion is usually the most sensible physically, EMAX serves principally
to keep an automatically generated master list within manageable bounds; the
current dimensioning of MOLSCAT does not allow the master list to contain
more than 200 functions.
For the common case of surface scattering with the incident beam approaching
along a symmetry direction, MOLSCAT automatically constructs the appropriate
symmetrised linear combinations of basis functions. This takes place as part
of the second step described above, since it is not until that stage that the
program knows the incident angle.
For normal incidence scattering from square or hexagonal lattices, enormous
savings in computer time are possible if the full symmetry is accounted for.
Special versions of the routine SURBAS which take advantage of the symmetry
are available from JMH.
3.9 ITYPE=9: general user-supplied routine.
To facilitate writing user-supplied code which handles different collision
types than those described above, skeleton routines are provided, with
numerous comment cards, which provide a template for the required code.
JHALF (default 1) is introduced to facilitate half-odd-integer JTOT
values in BASE9 codes.
3.10 Decoupling approximations
The following variables are specifically relevant to decoupling
approximations.
JZCSMX (default -1). By default, CS/DLD calculations are performed for
all values of the body-fixed projection number m up to the largest
diatom j value in the basis set. However, if cross sections are
required only between the lowest few levels, this is unnecessary.
If JZCSMX is set .gt. -1 on input, m will be limited by JZCSMX
instead of JMAX. Cross sections involving levels with j .gt. JZCSMX
are then not valid.
JZCSFL (default 0). This is a historical remnant and values other than the
default are not recommended. In CS calculations it sets the orbital
quantum number in each channel as follows (with allowed values of
JZCSFL = -1, 0, 1):
L(I) = IABS( JTOT + JZCSFL*J(I) ).
IBOUND (default 0). For spectroscopic applications using the helicity
decoupling (coupled states) approximation, it is useful to be able to
include the diagonal part of the Coriolis coupling exactly. This is
achieved if IBOUND = 1 on input. Note, however, that the asymptotic
matching still uses Bessel functions for integer l, so that the
convergence with respect to RMAX may be slightly irregular when this
option is used. This option is supported only for ITYPE = 21, 22, 23
and 27.
Another effect of IBOUND = 1 is to suppress coupled states
calculations for projection quantum numbers greater than JTOT.
IOSNGP an integer array of dimension 3. It specifies the number of (Gauss)
integration points to use for orientation dependence of the
scattering matrix. IOSNGP(1) is the number of points for theta.
In ITYPE=105 and 106, IOSNGP(2) is the number of points for phi;
in ITYPE=103, the 3 values are for theta(1), theta(2), delta(phi).
If values are not given (defaults are 0) the program will try to
choose the minimum number needed for requested values of LMAX, MMAX
(see &INPUT, Section 2.1) and/or input basis set rotor levels.
IPHIFL (default 0) controls the type of numerical quadrature on phi for
ITYPE = 103, 105 and 106. The default requests Gauss-Mehler
(recommended) while IPHIFL .ne. 0 requests Gauss-Legendre.
4. NAMELIST block &POTL
_______________________
MOLSCAT obtains the intermolecular potential via calls to subroutine POTENL
which returns values at input distance R as an expansion in appropriate angular
functions. The general-purpose version of POTENL supplied with MOLSCAT is
described in this section. Simple potentials, whose coefficients can be
expressed as a sum of exponentials and inverse powers, can be input directly
as described in Section 4.1. For more complicated potential forms users
must supply their own potential routine(s). This may be done by (1) using
the VSTAR mechanism described in Section 4.2 to describe (some of) the
radial expansion coefficients, (2) using the VRTP mechanism described in
Section 4.3 to provide an angle-dependent potential rather than radial
expansion coefficients, or (3) replacing the POTENL routine completely.
The last method is sometimes the most efficient as it may avoid repeating
computations required if the different potential terms are evaluated
separately. The specification of the POTENL subroutine is described in
detail in Section 5. Versions of POTENL that handle numerical tabulations
of the angular coefficients (which often result from ab initio surfaces)
by polynomial or spline interpolation are available from SG.
The standard POTENL routine currently implements three options:
1. The potential is expanded in angular functions appropriate to the
collision type. These are specified by a nonnegative input value of
MXLAM (the number of angular terms) and ITYPE-dependent input values
for LAMBDA to describe the symmetry of each term. The radial
coefficients may be described by powers and exponentials or by the
VSTAR mechanism (Section 4.2) for each term.
2. The potential is expanded in angular functions as above, but the
radial coefficients are projected via the VRTP mechanism (Section 4.3).
The angular terms may be described by either
a.) nonnegative MXLAM must be input along with appropriate values of
LAMBDA to specify the symmetries; and LVRTP= .TRUE. must be input
to specify projection via the VRTP mechanism (Section 4.3).
b.) MXLAM.LE.0 (default); nonnegative LMAX (default=-1) must be input.
This option is allowed only for non-IOS cases (ITYPE.lt.100).
LMAX (along with other input parameters depending on ITYPE: MMAX,
L1MAX, L2MAX, IVMIN, IVMAX) are used to generate the LAMBDA array.
In case of conflict (MXLAM.gt.0, LMAX.GE.0) LMAX is ignored, i.e.,
case (a) takes precedence.
3. The potential is NOT expanded in angular functions. This is suitable
only for IOS calculations. The potential as a function of collision
coordinates must be supplied by a VRTP routine (Section 4.3).
MXLAM.LE.0 (default) and ITYPE.gt.100 must be specified.
The parameters that may be input in &POTL NAMELIST are as follows:
RM is the length scaling factor that will be used by MOLSCAT. RM must be
supplied in Angstroms (default 1.0), and defines the units of RMIN,
RMAX etc. input in &INPUT. Subsequent calls to POTENL will also handle
distances in these units.
EPSIL is the energy scaling factor for the potential (default 1.0). EPSIL
must be supplied in 1/cm, and subsequent calls to POTENL must return
potential coefficients in these energy units. The value given to EPSIL
does not affect MOLSCAT's interpretation of energies input in &INPUT
and &BASIS.
MXLAM is the number of terms in the potential expansion.
MXSYM is a synonym for MXLAM.
LAMBDA is an array specifying the symmetries of the MXLAM different terms. It
is structured differently for the different values of ITYP; see the
documentation for an initialisation call to POTENL in Section 5.1 for
more details. Note that potential symmetries specified in LAMBDA may
be in any order and may be repeated, except for ITYP=2 and 7 for which
each symmetry (including vibration-rotation labels) must be unique.
LVRTP is a logical variable. If LVRTP = .FALSE. (the default), the potential
is specified in terms of its expansion in angular functions. If LVRTP
= .TRUE., subroutine VRTP is called as described in Section 4.3 to
evaluate the potential at specified collision coordinates. Note that
if MXLAM.LE.0 (which is the default value) LVRTP is forced to .TRUE..
LMAX When LVRTP=.TRUE. and MXLAM.LE.0, a nonnegative LMAX (default=-1)
MMAX may be input to generate the LAMBDA array internally (but only
L1MAX for non-IOS cases, i.e., ITYPE .lt. 100). LMAX specifies the
L2MAX highest Legendre term to include (ITYP=1,2). Similarly, for
ITYP=5,6 LMAX and MMAX indicate the highest spherical harmonic
to include; if MMAX is not specified (default = -1) MMAX is
set equal to LMAX. For ITYP=3, L1MAX and L2MAX are the highest
Legendre terms for the two molecules (L1MAX is a synonym for LMAX).
Note that IHOMO and ICNSYM (described below) may be used to exclude
values not allowed by symmetry.
IVMIN For ITYP=2, if LMAX is used with LVRTP=.TRUE. to generate the
IVMAX LAMBDA array internally, IVMIN must be given a nonnegative value
(default=-1) to indicate the lowest vibrational state. IVMAX, if
greater than IVMIN indicates the highest vibrational state;
otherwise IVMAX is set equal to IVMIN.
NPTS An array which contains the number of Gauss integration points
used in projecting the angular components of the potential.
NPTS(1) is the number of points for the theta integration common
to all supported ITYPE. NPTS(2) is for the integration over
(harmonic oscillator) vibrational functions for ITYP=2; for the
integration over phi for ITYP=5,6, and for the integration over
theta-2 for ITYP = 3. NPTS(3) is for the phi integration for
ITYP = 3.
NPT is a synonym for NPTS(1).
NPS is a synonym for NPTS(2).
NPT should be somewhat larger than the largest value of lambda
in the potential expansion. If IHOMO = 2 (see below)
only half the points will be used. For ITYP = 5 or 6, NPS should
be somewhat larger than the largest value of m in the potential
expansion; for ITYP = 2 NPS should be somewhat larger than the
largest value of v. If ICNSYM.gt.1 (see below), then NPS points
are used on the range 0 .lt. phi .lt. 2*pi/ICNSYM.
IHOMO If set to 2, indicates homonuclear symmetry (reflection about
theta=pi/2). Default value of 1 indicates heteronuclear symmetry.
For potentials which are not expanded in angular functions this
value should be set here or internally in VRTP; for potentials
which are expanded in angular functions, it is determined from
the LAMBDA array.
ICNSYM For nonlinear molecules denotes symmetry about z-axis; e.g., for
ammonia ICNSYM = 3 indicates three-fold symmetry. Default is 1.
For potentials which are not expanded in angular functions this
value should be set here or internally in VRTP; for potentials
which are expanded in angular functions, it is determined from
the LAMBDA array.
IHOMO2 and ICNSY2 have been added to support a 2nd molecule
in version 14. For historical reasons ICNSYM is used to
describe homonuclear symmetry of the second molecule in
ITYP = 3; however, input of IHOMO2 will produce the same effect.
Note: values of IHOMO and ICNSYM in &POTL may be overridden
by values set in the initialization call to a user supplied VRTP.
NTERM is an array of MXLAM integers, describing how to evaluate the
potential if LVRTP .EQ. .FALSE.. Each element of the NTERM array
corresponds to one element in the expansion of the potential.
If NTERM(I) is less than zero, this element of the potential array
will be evaluated using the VSTAR mechanism described below.
If NTERM(I) is positive, this element of the potential array will be
evaluated as a sum of NTERM exponential or inverse power terms,
specified by A, E and NPOWER.
4.1 Exponentials and powers; standard MOLSCAT POTENL routine
The following variables describe the input required when an element of the
potential array is to be represented as a sum of exponentials and inverse
powers.
NPOWER For each positive value in the NTERM array, there must be NTERM values
of NPOWER specified. If an element of NPOWER is 0, that potential term
is an exponential described by A and E. If an element of NPOWER is
negative, that potential term is an inverse power term described by A
and NPOWER.
E is an array of exponential factors. One value must be provided for
each zero in the NPOWER array. N.b. values should be negative.
A is an array of pre-exponential and pre-power factors. One value of A
must be supplied for each number in the NPOWER array. If NPOWER(I) is
0, the corresponding potential term is
A(I) * EXP(E(IE)*R)
while if NPOWER(I) is negative, it is
A(I) * R**NPOWER(I)
CFLAG For ITYP=5 and 6, there are two common definitions of the potential
expansion: the potential may be expanded either in terms of normalised
spherical harmonics Ylm or in terms of renormalised spherical
harmonics Clm. Internally, the program uses the former definition. If
the coefficients are input in terms of renormalised functions, CFLAG
should be set to 1 to instruct the program to multiply each input
coefficient by an appropriate factor.
4.2 The VSTAR mechanism
If any element of the NTERM array is negative, subroutines VINIT, VSTAR, VSTAR1
and VSTAR2 must be supplied by the user. However, VSTAR1 and VSTAR2 will be
called only if the VIVS or WKB propagator is being used in a mode requiring
explicit derivatives (see Section 2.12). The dummy routines supplied with
MOLSCAT, if called, will print an error message and terminate the program.
VINIT will be called once for each negative value in the NTERM array during the
initialisation call to POTENL. The syntax of the call is
CALL VINIT( I, RM, EPSIL )
On entry, I is the index of the potential element in the NTERM (and LAMBDA)
arrays, and RM and EPSIL have the values read in in &POTL. None of these
arguments should be altered by the routine, but it may be desirable to store
their values in local variables so that they are available on subsequent calls.
VINIT should also read any data required for the potential evaluation.
VSTAR is called many times during evaluation calls to POTENL. The syntax of the
call is
CALL VSTAR( I, RR, SUM )
On entry, I is the index of the potential element required, and RR is the
intermolecular distance (in units of RM). VSTAR must evaluate the I'th
potential element and return its value in SUM (in units of EPSIL).
The specifications for VSTAR1 and VSTAR2 are identical to that for VSTAR,
except that they must return the first and second radial derivatives of the
I'th potential element respectively. For most propagators, the radial
derivatives are not used: in these cases, it is adequate to supply VSTAR1 and
VSTAR2 routines that simply print an error message and exit.
4.3 The VRTP mechanism
It is often more convenient to specify V(R,angles) explicitly than to expand
the potential in angular functions. The general-purpose version of POTENL
allows this for IOS calculations for ITYPE = 101, 102, 103, 105, and 106 and
for other calculations for ITYP = 1, 2, 3, 5 and 6. The option is invoked by
specifying LVRTP = .TRUE. or MXSYM.LE.0; POTENL then expects
the potential to be supplied by a routine VRTP. (N.b. specifying MXSYM.gt.0
and LVRTP=.TRUE. for an IOS case will cause Legendre terms to be projected
via the VRTP mechanism, but this is generally a very inefficient way to
handlee IOS cases.)
The calling sequence is
CALL VRTP(IDERIV,R,V)
On an initialization call, IDERIV=-1, and RM(=R) and EPSIL(=V) can be set by
the routine or, if read in &POTL, used by the routine. The initialization call
may also set the variables IHOMO and ICNSYM in the common block ANGLES:
COMMON /ANGLES/ COSANG(7),FACTOR,IHOMO,ICNSYM,IHOMO2,ICNSY2
NOTE: dimensions in this common block were changed in version 14 to
anticipate more complex collision cases; EARLIER VERSIONS OF VRTP MUST
BE MODIFIED ACCORDINGLY when run with version 14 of MOLSCAT.
VRTP should set IHOMO=1 if the potential has heteronuclear symmetry or IHOMO=2
if it has homonuclear symmetry (with respect to theta), and ICNSYM to specify
the order of any rotational symmetry about the z axis. For ITYPE=103, ICNSYM is
used to specify heteronuclear or homonuclear symmetry for the second linear
molecule. Note that, for IOS calculations, IHOMO and ICNSYM would normally be
determined automatically by examining the LAMBDA array, but this can be done
only if the potential is expanded in angular functions.
Subsequent calls (IDERIV=0,1,2 for the potential and its first and second
derivatives respectively) require the routine to return in V the potential at
distance R and "angles" specified in the COSANG array.
Some propagators do not use POTENL/VRTP calls with IDERIV .gt. 0, and for these
it is adequate to supply a VRTP routine that traps calls with IDERIV .gt. 0,
prints an error message, and exits.
For ITYP = 1, 2, 3, 5, and 6, COSANG(1) is cos(theta).
For ITYP = 3 COSANG(2) is cos(theta-2) and COSANG(3) is phi.
For ITYP = 5 and 6, COSANG(2) is phi.
For non-IOS ITYP = 2 cases, COSANG(2) is q, the reduced harmonic oscillator
coordinate, such that the ground-state vibrational wavefunction of the diatom
is exp(-q**2/2); the code ASSUMES that the molecule is a harmonic oscillator.
For ITYPE = 102, VRTP does not use COSANG(2), and must return a vector of
matrix elements (v1|V(R,theta)|v2) in the order expected by POTENL.
Appropriate communication is provided by a negative MXLAM in &POTL input,
with v1,v2 values supplied as LAMBDA(3*i-1),LAMBDA(3*i),i=1,abs(MXLAM);
LAMBDA(3*i-2), which would normally have the order of the Legendre
polynomial, is not relevant and is ignored.
A sample VRTP routine is supplied with MOLSCAT which evaluates a simple
exponential/inverse power potential for a mock asymmetric top-atom case
and which is used in tests of the ITYP = 6 code. This routine must be
replaced in calculations using the VRTP mechanism.
5. User-supplied POTENL routines
___________________________________
As described above, many potentials can be specified either in terms of
exponentials and inverse powers or using a user-supplied VRTP or VSTAR routine.
It is best to use one of these mechanisms if it can be adapted to your needs.
However, in the rare cases where this is not possible, you will need to write a
complete POTENL routine. This section describes the specification that is
required of subroutine POTENL, in order to allow users to write their own
potential routines for complicated potentials.
In each run, POTENL is called once for initialisation purposes, and may on this
call read any data necessary to specify the potential. It returns information
about the terms present in the potential expansion for processing by MOLSCAT.
Subsequently, POTENL is called many times to evaluate the potential
coefficients for particular intermolecular distances.
The syntax of a call to POTENL is
CALL POTENL (IC,MXLAM,NPOTL,LAMBDA,RR,P,ITYP)
RR and P are REAL*8 arguments; the rest are INTEGER. LAMBDA and P are arrays
which should be dimensioned at LAMBDA(1), P(1) to switch off FORTRAN array
bound checking.
There are three basic types of call to POTENL;
IC = -1 Initialisation call. POTENL is called once with IC = -1, before any
other calls to it, to allow it to read any necessary data and set up
various parameters.
IC = 0 At subsequent calls to POTENL, IC is 0, and the routine must evaluate
the intermolecular potential.
IC = 1,2 The VIVS propagator is most efficient if radial derivatives of the
potential coefficients are available. If POTENL is called with IC = 1
or 2, it must evaluate the IC'th radial derivative of each potential
coefficient.
The WKB integrator also requires first derivatives to start the
integration (in locating turning points).
For most other propagators, IC = 1 and 2 are not used, and it is
adequate to supply a POTENL routine that prints an error message and
exits if it is called with IC = 1 or 2.
The specification of POTENL for initialisation and evaluation calls will be
described separately.
5.1 Initialisation
IC IC = -1 is the flag for an initialisation call.
MXLAM On entry, MXLAM specifies the size of the LAMBDA array which has been
provided externally. This value is used to check for array bound
errors.
On exit, MXLAM must specify the number of distinct terms in the
expansion of the potential (ie. the size of the P array that will
be returned by subsequent calls to POTENL.
NPOTL On entry, NPOTL is not set.
On exit, NPOTL must contain the appropriate value for the collision
partners involved, specified by ITYP = MOD(ITYPE,10). The required
values of NPOTL are
ITYP = 1, 3, 4, 5, 6: NPOTL=MXLAM
ITYP = 2 and 7: NPOTL=(KMAX+1), where KMAX is the order of the
largest Legendre component in the potential. See code in
COUPLE for precise use.
ITYP = 8: This is complicated, and depends on the particular version
of SURBAS in use. Basically, NPOTL = 2 for the usual version
of SURBAS, and NPOTL = 6 or 8 for special normal incidence
versions (depending on symmetry of lattice). Existing
versions of SURBAS place the necessary value in COMMON/NPOT/,
and POTENL picks it up from there.
LAMBDA Not set on entry.
On exit, the LAMBDA array must contain indices specifying the symmetry
of the potential terms that will be used. The LAMBDA array is used
differently for each collision type (ITYP). Although it is externally
a one-dimensional array, it is conceptually two-dimensional for some
collision types, and may be handled explicitly as a two-dimensional
array in POTENL if desired. Each element (or column) of LAMBDA
corresponds to an element of the P array returned by subsequent calls
to POTENL. MOLSCAT does not require that the symmetry terms be
supplied in any particular order, but just that the Ith column of the
LAMBDA array should correspond to the Ith element of the P array.
Detailed specifications of LAMBDA for the different ITYP values
follow. This should be read in conjunction with the documentation for
the P array for an evaluation call.
ITYP = 1:
Potential is expanded in Legendre polynomials (these are
unnormalized, such that P(2)=3*x*x-1).
LAMBDA is a one-dimensional array.
LAMBDA(I) contains the Legendre polynomial order of P(I)
ITYP = 2:
Potential is expanded in Legendre polynomials for each (v,v') pair.
LAMBDA is a two-dimensional array LAMBDA(3, ).
LAMBDA(1,I) contains the Legendre polynomial order.
LAMBDA(2,I) contains the vibrational quantum number v.
LAMBDA(3,I) contains the vibrational quantum number v'.
Note that the matrix elements are invariant to exchange of v and v',
and only one of these should be supplied.
ITYP = 3:
Potential is expanded as contracted spherical harmonics as described
by Green, JCP 62, 2271 (1975).
LAMBDA is a two-dimensional array LAMBDA(3, )
LAMBDA(1,I) is the index corresponding to l1
LAMBDA(2,I) is the index corresponding to l2
LAMBDA(3,I) is the index corresponding to the vector sum l = (l1+l2)
Note that, for identical molecules, the potential is invariant to
exchanging l1 and l2; nonetheless, for l1.ne.l2 BOTH terms MUST be
supplied.
ITYP = 4:
Potential is expanded in functions that are rotationally
invariant contractions of rotation matrices and spherical harmonics
as described by Phillips, Maluendes, McLean, and Green, JCP xxx,
nnnn (1994).
LAMBDA is a two-dimensional array LAMBDA(4, )
LAMBDA(1,I) is p1, the tensor order for the asymmetric rotor
LAMBDA(2,I) is q1, projection of p1 on the asymmetric rotor axis
LAMBDA(3,I) is p2, the tensor (harmonic) order for the linear rotor
LAMBDA(4,I) is p, the contraction of p1 and p2
ITYP = 5:
Potential is expanded as a sum over angular functions
[(Y(l,m) + (-)**m Y(l,-m))] / [(1 + delta(m,0))].
The Y(l,m) are normalized spherical harmonics; note that V(0,0) then
differs from the 'spherical average' by a factor. The angle theta is
measured from the z-axis, which is the symmetry axis of the top. The
angle phi is measured from the x-z plane, which MUST be a reflection
plane in the current implementation.
LAMBDA is a two-dimensional array LAMBDA(2, ).
LAMBDA(1,I) is l(I).
LAMBDA(2,I) is m(I).
ITYP = 6:
Specification of POTENL is the same as for ITYP = 5.
ITYP = 7:
LAMBDA(1,I) is the Legendre polynomial order;
LAMBDA(2,I), LAMBDA(3,I), and LAMBDA(4,I), LAMBDA(5,I)
are the diatom quantum numbers v,j and v',j' respectively.
Note that the matrix elements are invariant to exchange of vj and
v'j', and only one of these should be supplied.
ITYP = 8:
LAMBDA is a two-dimensional array LAMBDA(2, ).
LAMBDA(1,I) and LAMBDA(2,I) are the reciprocal lattice indices
corresponding to potential element P(I). They must be specified in the
unique form returned by subroutine ORDER. This last requirement makes
sure that the lattice symmetry is properly accounted for.
RR Not set on input.
On exit, RR must contain a length scaling factor which is applied to
(virtually) all the distances used internally in MOLSCAT. It is
specified in Angstroms. It is often convenient to set RR = 1.0D0 in
the initialisation call to POTENL, and to handle everything in
Angstroms thereafter.
P Not set on input.
On exit, P(1) must contain the energy scaling factor EPSIL to be used
internally in MOLSCAT, and subsequent calls to POTENL must return
energies in units of EPSIL. EPSIL is specified in 1/cm. It may
be convenient to set EPSIL = 1.0D0 in the initialisation call to
POTENL, and to handle everything in 1/cm thereafter. The value
given to EPSIL here does NOT affect the interpretation of energy
parameters input in &INPUT and &BASIS.
ITYP On entry, ITYP is the collision type.
If POTENL is coded specifically for a particular value of ITYP, it
should check that the correct value is passed, as a precaution against
the accidental use of the wrong executable version of MOLSCAT. The
value of this parameter should not be changed by POTENL.
5.2 Evaluation
For an evaluation call to POTENL, only the IC, MXLAM, NPOTL, RR and P arguments
are passed. LAMBDA and ITYP do NOT contain the values they were given in the
initialisation call, so these must be stored internally in POTENL if you need
them for an evaluation call.
IC IC = 0, 1 or 2 is the flag for an evaluation call:
IC = 0 evaluate the potential V
IC = 1 evaluate dV/dR
IC = 2 evaluate d2V/dR2
Note that calls to POTENL with IC = 1 or 2 only occur for the VIVS
propagator if IVP, IVPP, ISHIFT or IDIAG is set and for the WKB
integrator. Even then, there is an option (controlled by the &INPUT
logical variable NUMDER) which allows the derivatives to be evaluated
numerically without making an IC = 1 or 2 call to POTENL. It is thus
not altogether necessary for POTENL to cope with IC = 1 and 2 calls,
but it should at least trap an attempt to call it this way and print
an error message.
RR The intermolecular distance at which the potential is to be evaluated,
in units of RM (see the RR argument for an initialisation call to
POTENL)
P On exit, P(I),I=1,MXLAM must contain the potential array in the
order specified earlier by the LAMBDA array returned by the
initialisation call to POTENL. The P array is required in units of
EPSIL (see discussion of the initialisation call).
6. NAMELIST block &CONV
_______________________
If the variable ICONV in NAMELIST &INPUT is set to 1 on input, the program will
calculate S-matrices for successive values of the step size, RMIN, RMID or RMAX
and print out the root mean square (RMS) change in S-matrix elements and
transition probabilities each time. This provides an automated means of
choosing integration parameters capable of providing any required accuracy.
In this mode, the program performs NNRG calculations. Only a single energy
should be specified in the ENERGY array, and DNRG should be set to 0. Only a
single S-matrix is stored between calculations, so the program must not loop
over angular momenta or parity cases; thus JTOTU must be equal to JTOTL, and
MSET must be nonzero. Remember that, in some modes, MOLSCAT may determine RMIN,
RMID and RMAX internally, and it is safest to use the convergence testing with
these options switched off: IRMSET=0, RVFAC=0., IRXSET=0.
IVAR (default 0) selects which variable is to be altered. Note that the
first calculation is used as a "reference" calculation (unless ICONVU
.lt. 0, see below), so that the subsequent calculations must be less
accurate than the first, not more accurate. For successive
calculations, IVAR is used as follows: IVAR=0: step size is doubled
each time. IVAR=1: RMIN is increased by DR each time. IVAR=2: RMID is
increased by DR each time. IVAR=3: RMAX is decreased by DR each time.
DR is the amount by which RMIN or RMID is increased or RMAX is decreased
in successive calculations.
ICONVU (default 0) is used if convergence information is to be transferred
from one run of MOLSCAT to another. This is useful if the results
obtained using two different propagators are to be compared. If ICONVU
.gt. 0 on input, the first S-matrix is written out (unformatted) to
channel ICONVU; if ICONVU .lt. 0, a previously saved S-matrix is read
from channel |ICONVU|, and used as the reference S-matrix in
calculating RMS errors.
7. Array dimensioning
_____________________
The main program for MOLSCAT declares one large array, X, which is held in
COMMON/MEMORY/. DRIVER and other routines then partition up the X array
according to the size of the scattering problem being tackled. Thus, very
few of the arrays used internally by MOLSCAT are explicitly dimensioned, and it
is seldom necessary for users to concern themselves with array dimensioning. If
MOLSCAT finds that the X array is not big enough, it will (in most cases) print
an error message and exit tidily; it is then necessary to re-run the program
with a larger X array declared in the main program. Note that the error message
specifies the storage required for the current allocation; additional storage
may be needed for subsequent allocations.
There are a few exceptions to this, concerning arrays which are either input as
data or are in COMMON, and thus cannot be flexibly dimensioned. Restrictions
arising from the dimensions of these arrays have been described in this manual.
8. COMMON blocks
________________
MOLSCAT uses a number of COMMON blocks internally, and the names of these
should be avoided when writing subroutines to link with the MOLSCAT code. The
named COMMON blocks used in MOLSCAT are:
DRIVE MEMORY CMBASE PRBASE INTPAC ANGLES VLSAVE HIBRIN
LDVVCM WKBCOM EIGSUM LATSYM IOUTCM ASSVAR NPOT POPT
VLSAVE
9. Subroutines of the MOLSCAT system
____________________________________
Most subroutines in MOLSCAT are extensively commented; see the code for
details. The following descriptions list where the various subroutines are
called from, and in most cases give a brief description of their function.
AIRPRP Subroutine called by AXSCAT
Airy propagator routine
AIRYMP Subroutine called by SPROPN
Modulus and phase of Airy functions and derivatives
ASROT Subroutine called by SET6C
Calculates asymmetric rotor energies and wavefunctions
ASYME Entry point in subroutine SET6, called by COUPLE
ASYMF Entry point in subroutine SET6, called by COUPLE
ASYMG Entry point in subroutine SET6, called by MCGCPL
AXSCAT Subroutine called by STORAG
Controlling routine for the log-derivative/Airy hybrid propagator
BASE Subroutine called by DRIVER
Handles generation of basis set
BASE9, BAS9IN, SET9, CPL9, DEGEN9
Subroutine / entry points called from BASE etc.
User-supplied routines to specify a special coupling scheme
BASIN Entry point in subroutine BASE, called by DRIVER
Initialisation entry for BASE - reads &BASIS data
CHCK6I Subroutine called by SET6I
Check asymmetric rotor functions for orthogonality
CHECK6 Subroutine called by SET6
Check asymmetric rotor functions for orthogonality
CHKSTR Subroutine called by DRIVER, etc.
Checks that there is enough storage left in the X array
COLIM Subroutine called by RMTPRP
CONVRG Subroutine called by DRIVER
Performs convergence checking if &INPUT ICONV=1
CORR Subroutine called by AIRPRP
Correction terms for Airy propagator
COUPLE Subroutine called by BASE
Handles angular momentum coupling elements
CPLOUT Subroutine called by BASE
Output routine for angular momentum coupling elements
CPL3 Subroutine called by COUPLE
Coupling matrices for ITYPE = 3
CPL4 Subroutine called by COUPLE; entry CPL24 called from MCGCPL
Coupling matrices for ITYPE = 4, 24
CPL21 Subroutine called by MCGCPL
Coupling matrices for ITYPE = 21
CPL22 Subroutine called by MCGCPL
Coupling matrices for ITYPE = 22
CPL23 Subroutine called by MCGCPL
Coupling matrices for ITYPE = 23
CPL25 Subroutine called by MCGCPL
Coupling matrices for ITYPE = 25
CPL26 Subroutine called by MCGCPL
Coupling matrices for ITYPE = 26
DAPROP Subroutine called by AXSCAT (INTFLG = 8), DASCAT (INTFLG = 6)
Propagation routine for the diabatic modified log-derivative method
DASCAT Subroutine called by STORAG
Controls the diabatic modified log-derivative propagator
DASIZE, DAOPEN, DACLOS, DAWR1, DAWR2, DARD1, DARD2
Subroutines called by PRBR for in-core simulation of direct access
files
DEGENF Entry point in subroutine BASE, called by OUTPUT
Returns degeneracy factor for denominator of cross-section expressions
DELRD Subroutine called by VIVAS
Calculates length of next step
DERMAT Subroutine called by VIVAS WKB
Calculates first or second radial derivative of potential matrix
DGEMUL Subroutine called from AIRPRP QAPROP RMTPRP STABIL TRNSFM VIVAS
Matrix multiplication (ESSL routine)
DGESV Subroutine called from DVSCAT STABIL VIVAS YTOK
Solves linear equations (LAPACK routine)
DMSYM Subroutine called by ASROT
Symmetrizes spherical top wavefunctions
DRIVER Subroutine called by MOLSCAT main program
Driver reads &INPUT data and handles the calls to the major routines.
DSYFIL Subroutine called by AIRPRP LDVIVS RMTPRP TRNSFM WAVVEC YTOK
Fills in upper or lower triangle of symmetric matrix
DVFREE Subroutine called by DVSCAT
Asymptotic matching for DeVogelaere algorithm
DVSCAT Subroutine called by STORAG
Main routine for DeVogelaere algorithm
DVSTRT Subroutine called by DVSCAT
Starts DeVogelaere propagation
EAVG Subroutine called by DRIVER
Provides energy values suitable for Boltzmann averaging
ECNV Subroutine called by BASE DRIVER
Handles different possible energy units and returns conversion factor.
ECNVX Subroutine called by ECNV
Handles character string data values in EUNITS
EPSUM Subroutine called by OUTPUT
Calculates S-matrix eigenphase sum from K-matrix
ESYMTP Subroutine called by COUPLE SET6
F02AAF Subroutine called by WAVEIG
Eigenvalues of a real symmetric matrix
F02ABF Subroutine called by ASROT DMSYM EPSUM POTENT QAPROP
RMTPRP SHRINK VIVAS
Eigenvalues and eigenvectors of a real symmetric matrix
FIND Subroutine called by SURBAS
Finds potential term coupling two G-vectors in surface scattering
FINDRM Subroutine called by DRIVER
For IRMSET .gt. 0 option, finds suitable starting point for integration
FINDRX Subroutine called by DRIVER
Checks that RMAX is beyond centrifugal barrier and increments it if
necessary
FSYMTP Subroutine called by COUPLE SET6
GASLEG Subroutine called by GAUSSP
Returns points and weights for Gauss-Legendre quadrature.
GAUSHP Subroutine called by POTENL (VRTP case)
Returns points and weights for Gauss-Hermite quadrature.
GAUSSP Subroutine called by IOSBGP
Interface to GASLEG for Gauss-Legendre points and weights
GCLOCK Subroutine called by DRIVER IOSDRV IOSCLC
Timing routine
GDATE Subroutine called by DRIVER
Returns character string with current date
GET102 Subroutine called by IOSBIN
Special ITYPE=102 input; dummy version supplied with MOLSCAT
GSYMTP Subroutine called by MCGCPL SET6
GTIME Subroutine called by DRIVER
Returns character string with current time of day
HEADER Subroutine called by DRIVER
Writes or checks a header on the ISCRU scratch file to ensure that it
is compatible with the current run.
HERM Called by POTENL (VRTP case)
Generates Hermite polynomials
HRECUR Called by GAUSHP
IDPART Subroutine called by BASE
Handles identical particle symmetry for ITYP = 3
IOSBGP Entry point in subroutine IOSBIN, called by IOSDRV
IOSBIN Subroutine called by BASIN
Processes &BASIS data for IOS cases
IOS1 Entry point in IOSBIN, called by IOSDRV
Sets up points and weights for orientation quadrature
IOS2 Entry point in IOSBIN, called by IOSCLC
Initializes values for the propagator in IOS cases
IOSCLC Subroutine called by IOSDRV
Main control of IOS: loops over energies and orientations and
accumulates generalized IOS cross sections
IOSDRV Subroutine called by DRIVER
Sets up storage for IOS cases and then calls IOSCLC
IOSOUT Subroutine called by IOSCLC
Outputs IOS generalized and state-to-state cross
sections
IOSPB Subroutine called by IOSCLC
Calculates pressure broadening cross sections for IOS cases
IPASYM Subroutine called by SET6
Checks symmetry of asymmetric rotor coefficients (ITYPE = 6)
ISUTP Subroutine called by IOSCLC
IVCHK subroutine checks compatibility of symmetries with IV() indexing.
IXQLF Entry point in subroutine IOSBIN, called by IOSOUT IOSPB SIG6
J3J000 Subroutine called by COUPLE CPL21 SET6
Recursive routine for Wigner 3-j symbols with zero projections
J6TO4 Subroutine called by SET4
J6J Subroutine called by COUPLE J9J SET6 SIXJ
Recursive routine for Wigner 6-j symbol
J9J Subroutine called by XNINEJ
Recursive routine for Wigner 9-j symbol
KSYM Subroutine called by YTOK
Forces symmetry on K matrix
KTOS Subroutine called by DASCAT DVSCAT LDVIVS QASCAT RMTPRP
Converts the real symmetric K matrix into the S matrix
LDPROP Subroutine called by LDVIVS
Log-derivative propagator
LDVIVS Subroutine called by STORAG
Controls the hybrid log-derivative/VIVS propagator
MASK Called by DRIVER
Machine-dependent code to suppress floating-point underflows
MAXMGV Subroutine called by POTENT
Find largest element of vector
MAXMGV Called by POTENT
Maximum magnitude element of a vector
MCGCPL Subroutine called by BASE
Handles angular momentum coupling for McGuire-Kouri coupled states
approximation
MHAACK Subroutine called by DRIVER
Prints a citation request for the HIBRIDON propagator
NEXTE Subroutine called by DRIVER
Estimates next energy set in resonance search option
ODPROP Subroutine called by DASCAT
Single-channel implementation of diabatic modified log-derivative
propagator
ORDER Subroutine called by FIND
Takes account of lattice symmetry for Fourier components of
atom-surface potential
OUTERR Entry point in subroutine OUTPUT, called by DVSCAT
OUTINT Entry point in subroutine OUTPUT, called by DRIVER
Initialisation entry for OUTPUT
OUTMAT Subroutine called by AIRPRP
Read or write transformation matrix
OUTPCH Entry point in subroutine OUTPUT, called by DRIVER
OUTPUT Subroutine called by DRIVER
Outputs S-matrices etc. and calculates state-to-state cross sections
PARITY Function called by BASE COUPLE CSRTRT ESYMTP FCOEF
FSYMTP GSYMTP MCGCPL PRBR ROTROT SET6 THREEJ
Returns (-1)**argument
PERT1 Subroutine called by VIVAS
Calculates perturbation corrections for VIVAS propagator
PERT2 Subroutine called by VIVAS
Calculates perturbation corrections for VIVAS propagator
PLM Subroutine called by IOSBIN
Calculates associated Legendre polynomials
POTENL Subroutine called by DRIVER WAVMAT DERMAT ODPROP
Evaluates intermolecular potential at distance R
POTENT Subroutine called by AIRPRP
Calculates wavevector matrix and diagonalises average potential
POTIN9 Routine called by POTENL for ITP=9; dummy version must be
replaced with appropraite user-supplied routine.
PRBOUT Entry point in subroutine PRBR, called by DRIVER
PRBR Subroutine called by DRIVER
Handles calculation of pressure broadening cross sections
PRBR3 Subroutine called by PRBR
Pressure broadening code for ITYPE = 3
PRBR3R Entry point in subroutine PRBR3, called by PRBR
PRBRIN Entry point in subroutine PRBR, called by DRIVER
Initialisation entry for PRBR
QAPROP Subroutine called by QASCAT
Propagation routine for the quasiadiabatic modified log-derivative
method (INTFLG = 7)
QASCAT Subroutine called by STORAG
Controls the quasiadiabatic modified log-derivative propagator
QSYMTP Subroutine called by CPL4; matrix elements for ITYPE = 4
RBES Subroutine called by DASCAT DVFREE QASCAT RMTPRP
Generates Riccati-Bessel functions
RDPCH Subroutine called by OUTPUT
Outputs cross sections
RMSBF Subroutine called by YTOK
Ratio of derivative to function value for modified spherical bessel
functions of the third kind
RMTPRP Subroutine called by STORAG
Main routine for R-matrix propagator algorithm
RSYM Subroutine called by DVSCAT
Forces symmetry on K-matrix in DeVogelaere algorithm
RSYMTP Subroutine called by CPL4; matrix elements for ITYPE = 24
SCAIRY Subroutine called by AIRYMP
Scaled Airy function and derivatives
SET1 Entry point in subroutine SETBAS, called by BASE
Handles basis set and energy levels for ITYP = 1
SET2 Entry point in subroutine SETBAS, called by BASE
Handles basis set and energy levels for ITYP = 2, 7
SET3 Entry point in subroutine SETBAS, called by BASE
Handles basis set and energy levels for ITYP = 3
SET4 Subroutine for ITYP = 4 basis functions, called by BASE
SET5 Entry point in subroutine SETBAS, called by BASE
Handles basis set and energy levels for ITYP = 5
SET6 Subroutine called by BASE
Handles basis set and energy levels for ITYP = 6
SET6C Subroutine called by SET6
Generates asymmetric rotor basis sets from rotational constants
SET6I Subroutine called by IOSBIN
IOS version of SET6
SET8 Entry point in subroutine SURBAS, called by BASE
Handles basis set and energy levels for ITYP = 8
SETBAS Name of subroutine containing SET1 - SET5
SGNCHK Subroutine called by RMTPRP
Ensures consistency of eigenvector signs from one R step to the next
SHRINK Subroutine called by RMTPRP
Performs basis set contraction for R-matrix propagator algorithm
(ISHRINK .gt. 0 option)
SIG6 Subroutine called by IOSOUT
Handles ITYP = 6 cross sections for IOS case
SIXJ Function called by FCOEF FSYMTP PRBR PRBR3 ROTROT
Calculates Wigner 6-j symbol
SPROPN Subroutine called by AIRPRP
Diagonal propagator for Airy method
STABIL Subroutine called by DVSCAT
Stabilisation routine for DeVogelaere algorithm
STEFF Subroutine called by FINDRM
Steffensen iteration for accelerated convergence on RMIN
STORAG Subroutine called by DRIVER IOSCLC
Sets up storage for and calls a propagator; all calls to propagators
are currently done through STORAG
STRY Subroutine called by RMTPRP
Tests whether eigenvectors of potential are changing slowly enough to
end R-matrix propagation
STSRCH Called by ECNVX
SURBAS Subroutine called by BASE
Handles basis set for surface scattering (ITYPE = 8)
SWRITE Subroutine called by OUTPUT
Write unformatted S matrix
SYMINV Subroutine called by DAPROP KTOS LDVIVS LDPROP QAPROP RMTPRP
Calculates the inverse of a real symmetric matrix
THREEJ Function called by COUPLE CSRTRT FCOEF FSYMTP MCGCPL ROTROT
Calculates Wigner 3-j symbol with all projections zero
THRJ Function called by CSRTRT ESYMTP FSYMTP GSYMTP MCGCPL PRBR PRBR3
Calculates Wigner 3-j symbol with no restriction on projections
TRNSFM Subroutine called by AIRPRP POTENT QAPROP RMTPRP SHRINK VIVAS
Performs a similarity transform
TRNSP Subroutine called by AIRPRP POTENT QAPROP RMTPRP SHRINK
Transposes a matrix
VINIT Subroutine called by the general-purpose version of POTENL.
Initialises for calls to potential routines VSTAR, VSTAR1 and VSTAR2.
VIVAS Subroutine called by LDVIVS
VIVAS propagator
VRTP Subroutine called by general-purpose version of POTENL
User-supplied routine which may be called to evaluate an interaction
potential which is not expanded in angular functions.
VSTAR, VSTAR1 and VSTAR2
User-supplied routines which may be called by the general-purpose
version of POTENL to evaluate the interaction potential and its first
and second derivatives respectively.
WAVEIG Subroutine called by AIRPRP
Get eigenvalues of potential matrix
WAVMAT Subroutine called by DAPROP DERMAT DVSCAT DVSTRT HEADER
LDPROP POTENT QAPROP RMSET RMTPRP VIVAS WKB
Calls POTENL to obtain potential coefficients, and forms potential
matrix for scattering routines.
WAVVEC Subroutine called by WAVMAT DERMAT
Efficient routine to handle vector products required by WAVMAT
WKB Subroutine called by STORAG
WKB propagator for one channel problems
XNINEJ Function called by CSRTRT ROTROT
Calculates Wigner 9-j symbol
YRR Function called by IOSBIN
Bispherical harmonic for two linear rotors
YTOK Subroutine called by LDVIVS RMTPRP
Calculates the K matrix from the log-derivative matrix
ZBES Called by GASLEG
10. Machine-dependent features
______________________________
MOLSCAT is written in standard FORTRAN 77 as far as possible, but there are
inevitably a number of features which depend on the particular computer being
used. These will be described briefly here.
10.1 Main program
The MOLSCAT main program does not do any processing; it simply allocates
storage and calls DRIVER to do all the work.
If MOLSCAT terminates with the message
CHKSTR. CANNOT PROVIDE REQUESTED STORAGE.
then it is usually sufficient to modify the main program to increase the
parameter MXDIM and recompile.
10.2 NAMELIST
Most of MOLSCAT's input is in NAMELIST format. This is not standard FORTRAN 77,
but is implemented in most compilers and is too attractive a feature to forego.
It does, however, introduce some quirks; for example, older CRAY versions of
NAMELIST could not cope with CHARACTER variables, so the variable LABEL is
handled in a peculiarly complicated manner.
For compilers that do not provide NAMELIST, it is usually possible to simulate
it using other compiler extensions. Code for doing this is provided in
commented-out form in the routines that read data, and the extra routines
needed are available from JMH.
10.3 Integer length
MOLSCAT obtains working storage by partitioning an array of 8-byte real
elements. On most machines, integers occupy only 4 bytes, so it is possible to
pack 2 integers into each 8-byte element. The variable NIPR, set in subroutine
DRIVER, must be equal to the number of integers that may be packed into 8
bytes. NIPR should be 2 on most machines, but 1 on a CRAY.
10.4 Date and time routines
MOLSCAT obtains the date and time of a run for output in the header by calls to
routines GDATE and GTIME. These are not standard, and must be simulated. GDATE
must return the current date as a CHARACTER*11, and GTIME must return the time
of date as a CHARACTER*9. In the last resort, they can be replaced by routines
that just return spaces.
10.5 CPU time routines
MOLSCAT outputs information on the CPU time taken by various steps, which it
obtains by calls to subroutine GCLOCK. This must also be simulated in a
machine-dependent manner; it is required to return the CPU time taken so far,
in seconds. As a last resort, it may return a result of zero.
10.6 Floating point underflow
Various routines used by MOLSCAT require that the result of an arithmetic
operation causing floating point underflow should be zero. MOLSCAT calls
subroutine MASK in order to set this, and MASK should call an appropriate
machine-dependent routine to suppress underflow. Some machines do this
automatically, or use compiler options to achieve the effect; for example, a
VAX ignores floating-point underflow if the compiler is invoked with the
qualifier /CHECK=NOFLOATUNDER.
Take care that the routine used does not use excessive CPU time. For example,
the routine ERRSET can be very expensive indeed on some IBM machines, because a
complicated error-handling routine is called every time floating-point
underflow occurs; if available the IBM VS FORTRAN, CALL XUFLOW(0), is
generally preferable.
10.7 Linear algebra
LAPACK and BLAS
In version 12, MOLSCAT was modified to use LAPACK linear algebra routines
wherever possible. LAPACK is the successor to LINPACK and EISPACK, and is
designed to provide near-optimum performance for large problems on as wide a
range of architectures as possible. Suppliers such as NAG and CRAY have already
included substantial parts of LAPACK in their standard libraries, and will be
including more in the future.
If possible, you should run MOLSCAT using LAPACK routines that are optimised
for your particular computer. However, if this is not possible, you can obtain
FORTRAN versions of the LAPACK routines from a NETLIB server: simply send an
email message containing a line such as
send dsyevx from lapack
to a NETLIB server such as netlib@ornl.gov (Oak Ridge National Lab) to obtain
the source for the LAPACK routine DSYEVX.
The LAPACK routines use BLAS (basic linear algebra subroutines) as much as
possible. BLAS level 1, level 2 and level 3 routines exist. You should use BLAS
routines optimised for your particular computer if possible. However, if no
optimised routines are available, you can get FORTRAN versions of the BLAS from
a NETLIB server by sending an email message containing a line such as
send dblas2 from core
to obtain the double-precision level 2 BLAS routines.
It is important that any user who implements new options in MOLSCAT should
perform matrix operations by calls to the routines described below, both for
ease of maintenance and to simplify the creation of efficient versions for
other computers.
LINEAR ALGEBRA routines supplied with MOLSCAT
DGEMUL Matrix multiplication
DGESV Solve linear equations
SYMINV Invert symmetric matrix
F02AAF Diagonalise symmetric matrix without eigenvectors
F02ABF Diagonalise symmetric matrix with eigenvectors
MOLSCAT also calls BLAS (basic linear algebra subroutines) such as DAXPY, DDOT
etc (or single precision SAXPY, SDOT etc in the CRAY version) in many places.
Subroutine ODPROP (for efficient single-channel propagation) also requires
special treatment for vectorisation to be achieved, and there is a special
version of this routine for CRAYs and similar machines.
1) Matrix multiplication
MOLSCAT calls DGEMUL, which is a routine from the IBM ESSL library. In the
distributed version, DGEMUL calls the BLAS routine DGEMM. If the real ESSL
DGEMUL routine is available, use it; otherwise, use the best BLAS version of
DGEMM that you can find.
2) Symmetric matrix inversion
MOLSCAT calls SYMINV. Two versions of SYMINV are recommended. The first,
distributed with MOLSCAT, calls the LAPACK routines DSYTRF and DSYTRI to
carry out the inversion. The second is a pure FORTRAN routine, which on
most machines is faster than LAPACK equivalents for relatively small problems
(N .lt 70). For optimum performance, you need to test both routines on your
machine.
Note that MOLSCAT really does require matrix inversion, despite the usual rule
to use linear equation solvers instead. This is because the propagators
involved save information from one step to the next, and this advantage is lost
if the problem is formulated in terms of linear equation solvers.
3) Linear equation solver
The speed of this routine is not critical for most propagators. MOLSCAT calls
the LAPACK routine DGESV directly.
4) Eigenvalues and eigenvectors of symmetric matrices
These routines are important for INTFLG = 3, 4, 7 and 8. MOLSCAT calls
diagonalisers by the NAG names F02AAF and F02ABF, and the NAG routines give
acceptable performance on most machines. However, the distributed version of
MOLSCAT provides routines that simulate F02AAF and F02ABF by calls to the
LAPACK routine DSYEVX.
10.8 File handling
MOLSCAT adheres to the FORTRAN 77 standard in its use of READ and WRITE
statements (including direct access files).
The OPEN statements do not use FILE='fname' parameters and you will have to
provide files with the naming convention for your system. For example,
IBM OS/MVS and CMS use FTnnF001 as the filename for UNIT=nn; IBM AIX
uses filename fort.NN.
The following files may be used by MOLSCAT, depending on the value of
parameter in the &INPUT data set.
unit formatted use
------ --------- -------------------------------------
ISCRU no Propagator scratch unit used if ISCRU.ne.0
ISAVEU no S-matrix file written if ISAVEU .gt. 0
ISIGU yes,DA Updated cross section file if ISIGU .gt. 0
KSAVE yes Saved values for resonance search if KSAVE .gt. 0
In the VAX implementation, it is also convenient to OPEN the main data file and
a supplied ISCRU file with the keywords SHARED and READONLY, so that several
jobs can access them simultaneously. The resulting OPEN statements are
nonstandard, and are commented out in the distribution version of MOLSCAT.
MOLSCAT always writes state-to-state cross sections (card images) to unit 7
at the end of a run, and it is necessary to supply a (system dependent) dummy
data set if this output is not wanted.
11. References
______________
It would be difficult to give a complete set of literature references for all
the various capabilities of MOLSCAT. The following is a selected list for
various aspects; it tends to reflect publications of the authors as the program
often follows these works rather closely.
ITYPE = 1
_________
A.M. Arthurs and A. Dalgarno, Proc. Roy. Soc. A256 540 (1963).
ITYPE = 2
_________
S. Green, J. Chem. Phys. 70, 4686 (1979).
ITYPE = 3
_________
S. Green, J. Chem. Phys. 62, 2271 (1975); J. Chem. Phys. 67, 715 (1977);
T.G. Heil et al., J. Chem. Phys. 68, 2562 (1978).
ITYPE = 5,6
___________
S. Green, J. Chem. Phys. 64, 3463 (1976); J. Chem. Phys. 70, 816 (1979).
ITYPE = 8
----------
G. Wolken, J. Chem. Phys. 58, 3047 (1973).
J. M. Hutson and C. Schwartz, J. Chem. Phys. 79, 5179 (1983).
Effective Potential approximation
_________________________________
H. Rabitz, J. Chem. Phys. 57, 1718 (1972).
Coupled states (centrifugal sudden) approximation
_________________________________________________
P. McGuire and D.J. Kouri, J. Chem. Phys. 60, 2488 (1974).
Decoupled L-Dominant (DLD) approximation
________________________________________
A. E. DePristo and M. H. Alexander, J. Chem. Phys. 64, 3009 (1976).
IOS approximation
_________________
R. Goldflam et al., J. Chem. Phys. 67, 4149 (1977).
Pressure broadening cross sections
__________________________________
R. Shafer and R.G. Gordon, J. Chem. Phys. 58, 5422 (1973).
S. Green et al., J. Chem. Phys. 66, 1409 (1977).
S. Green, J. Boissoles, and C. Boulet, JQSRT 39, 33 (1988).
R-matrix propagator
___________________
J.C. Light et al., J. Chem. Phys. 69, 3518 (1978).
Modified log-derivative propagator
__________________________________
D. E. Manolopoulos, J. Chem. Phys. 85, 6425 (1986).
Hybrid modified log-derivative/Airy propagator
______________________________________________
D. E. Manolopoulos and M. H. Alexander, J. Chem. Phys. 86, 2044 (1987).
WKB semi-classical 'propagator'
_______________________________
R. T Pack, J. Chem. Phys. 60, 633 (1974).