MOLSCAT User's Manual -


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:

Approximate methods are specified as follows:

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, EUNITC - EUNITS is an integer specifying the energy units for ELEVEL and for other molecular energy constants input in &BASIS (for example, BE, WE, WEXE, etc., discussed below). EUNITC is a character*4 variable which may be used for the same purpose. Both are interpreted in the same way as the EUNITS and EUNITC parameters of the &INPUT NAMELIST block (see Section 2.4) including the fact that a nonblank value for EUNITC will override any value input for EUNITS. Although these &BASIS input variables have the same names as values input via the &INPUT data set, they are distinct; a value input in &INPUT EUNITC, for example, has no effect on ELEVEL values.

IVLU - (default 0) if nonzero specifies a scratch data set unit used to hold the coupling array. This can provide a great savings in required memory at the cost of a great deal of I/O. This option is currently implemented only for ITYP = 6 and ITYP = 9.

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 functions 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) may be used to limit the k values included: if KSET is negative, only levels with k=|J2MAX| are included; if KSET is zero or positive, only levels with k less 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); this might not correspond with 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 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 .le. 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 choose specific levels from 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):

PRTYeven/odd kepsilon

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.

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 - 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 to L(I) = IABS(JTOT + JZCSFL*J(I)), with allowed values of JZCSFL = -1, 0, 1).

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,/code> 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.

Forward to Section 4

Back to Section 2

Back to the Table of Contents