CCP6 BOUND User manual Version 5 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 BOUND is a general-purpose package for performing calculations of bound state energies in weakly bound molecular systems. It is an outgrowth of the MOLSCAT package, which performs molecular scattering calculations. It is written in near-standard FORTRAN 77, and versions of it are available for IBM, VAX, SUN and CRAY computers. This manual describes the capabilities and modes of use of the program, and explains the input data required. JMH would be grateful to be informed of any bugs you encounter, either in BOUND or in this documentation. Version of 18 March 1994 Table of Contents 1. Introduction 1.0 Acknowledging use of BOUND 1.1 Monomer types 1.2 Propagators 1.3 Data format 2. NAMELIST block &INPUT 2.1 General purpose variables 2.2 Print level control 2.3 Angular momentum 2.4 Specifying energies 2.5 Calculating expectation values 2.6 Propagator scratch file 2.7 Surface band structure calculations 2.8 Control of propagators - Units of length - Range of integration - Step size 2.9 Automatic convergence checking 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 potential routines 6.1 Initialisation 6.2 Evaluation 6. Array dimensioning 7. Subroutines of the BOUND program 8. Machine-dependent features 9. References 1. Introduction _______________ The quantum-mechanical bound state problem can be formulated as a set of coupled differential equations similar to those encountered in molecular scattering theory. The only difference between the two cases is in the boundary conditions that must be applied. The BOUND program uses the technology developed for the scattering problem, and applies it to bound states. The discussion below is thus often expressed in scattering terminology. Bound states are located using the method originally developed by Johnson, J. Chem. Phys. 69, 4678 (1978), Appendix A, and further developed by Manolopoulos (Ph. D. Thesis, Cambridge University, 1988) and Hutson (Comp. Phys. Commun., to be published (1994)). A matching point RMID is chosen in the classically allowed region, and the log-derivative matrix is propagated to RMID, both outwards from the short-range classically forbidden region and inwards from the long-range classically forbidden region. At a bound state energy, one of the eigenvalues of the difference between the two log-derivative matrices is zero. The BOUND program locates bound states by performing a one-dimensional search for such zeroes of the matching eigenvalue. The BOUND program does not yet calculate wavefunctions explicitly, although this is planned for the next version. However, some information on the makeup of the bound state is available from the eigenvector corresponding to the smallest eigenvalue of the matching matrix, which may be shown to be the (unnormalised) wavefunction at the matching distance RMID. This clearly depends on the matching distance chosen, but can nevertheless give useful information when attempting to assign quantum numbers to the levels found by BOUND. The BOUND program also calculates a "node count" as described by Johnson. For a calculation at a specified energy, the node count is the number of eigenvalues of the coupled equations which lie below that energy. BOUND sets up the coupled equations for one total angular momentum JTOT and parity case M at a time (see below), and the combination of the node count and the matching eigenvector usually provides enough information to assign quantum numbers to the levels found. Unfortunately, the algorithms used are sometimes prone to showing spurious nodes in the deep classically forbidden regions; if this happens with your problem, I know of no solution other than trying another propagator! 1.0 Acknowledging use of BOUND Any publication that uses BOUND should reference it as: J. M. Hutson, BOUND computer code, version 5 (1993), 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. 1.1 Monomer types BOUND 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 systems ITYPE = 2 Atom - vibrating diatom systems ITYPE = 3 Linear rigid rotor - linear rigid rotor systems ITYPE = 5 Atom - symmetric top systems (also handles near symmetric tops) ITYPE = 6 Atom - asymmetric top systems (also handles spherical tops) ITYPE = 7 A generalised form of ITYPE = 2, allowing the intermolecular potential to depend on the diatom j quantum number ITYPE = 8 Surface band structure calculations. At present, the code is restricted to centrosymmetric lattices, for which the potential matrices are real. ITYPE = 9 Allows the user to specify the quantum numbers and coupling matrix elements to define a new collision type, by supplying a set of routines BASE9, BAS9IN, CPL9, etc. This feature is not yet documented in detail, as the interface to these routines is liable to change slightly in future versions. However, examples of these routines and a template for creating your own are available from JMH on request. So far, the feature has been used for atom-diatom close-coupling calculations in a body-fixed representation, for diatom-symmetric top close-coupling calculations, for atom-open shell diatom bound states and for triatomic bound states in hyperspherical harmonics. The computer time required to solve a set of N coupled equations is proportional to N**3, and in practice the limit on N is from 50 to 200, depending on the speed of the computer. The basis sets necessary for converged close-coupling calculations quickly exceed this limit as rotational constants decrease, particularly for monomer types 2 to 7. However, BOUND 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 Helicity decoupling (centrifugal decoupling) calculations, neglecting off-diagonal Coriolis terms. IADD = 30 Decoupled l-dominant 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 helicity decoupling approximation for atom - vibrating diatom systems. Not all these approximations are supported for all monomer types, though the most common ones are. The program will print a warning message and exit if you request an option that is not supported. 1.2 Propagators The coupled equations may be solved using any one of several methods, controlled by the input parameter &INPUT INTFLG (see section 2.1). 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, and is generally recommended except for problems with very strong coupling or where many energies are required. 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, especially at the first energy. 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 modified log-derivative propagators, some of the remaining work is also avoided at subsequent energies, so that subsequent energies may cost only 30% as much as the first. For the basic log-derivative propagator, the saving in time is much less spectacular unless the potential is very expensive to evaluate. 1.3 Data format BOUND's main data file is read on channel 5, and consists of 3 blocks of NAMELIST data: &INPUT, read in subroutine DRIVER &BASIS, read in entry BASIN of subroutine BASE &POTL, read in initialisation call of subroutine POTENL if the general-purpose version of POTENL is used (see below). Each NAMELIST block starts with &, where is INPUT, BASIS, or POTL as appropriate, and is terminated by &END. Between these delimiters, the NAMELIST input data consist of "=value" entries in the data file, where is the variable name. Note that all lines of data must start in column 2 or later; characters in column 1 will cause errors. On CRAYs, 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. The following sections will describe the variables which can be input in &INPUT, &BASIS, &POTL data, classified according to function. Note that there are also a number of variables in the NAMELIST blocks which are associated with MOLSCAT rather than BOUND. These should not be altered from their default values, and are not described here. 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.2 above) LASTIN (default 1). If LASTIN is set to 1, the run will terminate after this calculation has been completed. If LASTIN is 0, another complete set of data will be read once this set has been processed. 2.2 Print level control The amount of printed output produced by BOUND is controlled by the integer variable PRNTLV; sensible values of PRNTLV vary from 1 (when only details of energy convergence are required) to 35 (when debugging). PRNTLV = 3 prints out the matching determinant for each set of coupled equations solved, and some information on the makeup of the bound state is printed if PRNTLV is 5 or more. Voluminous debugging output starts appearing at PRNTLV = 12. The number of page throws generated is controlled by the parameter ITHROW (default 0). If ITHROW = 1, each new energy calculation starts on a new page; otherwise, page throws occur only between different parity cases. 2.3 Angular momentum BOUND has built-in loops over total angular momentum JTOT and symmetry types ("parity cases") 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. It should be noted that bound state energies can shift quite substantially between different angular momenta (and parity cases). By default, the program loops over all allowed parity cases M for each JTOT. 2.4 Specifying energies The total energies at which calculations are performed are controlled by the variables EMIN and EMAX and EUNITS. The program begins by performing calculations at EMIN and EMAX (for each JTOT and parity case) to find the number of eigenvalues that lie between them. It then tries to locate the eigenvalues. If NODMIN or NODMAX is set, the program only tries to find eigenvalues whose node counts lie between these limits. The program converges on each eigenvalue until the step size taken is less than DTOL. Since the convergence is usually quadratic, the resulting energies are usually much more precise than +/-DTOL. For each JTOT and parity case, the program carries out a maximum of MXCALC calculations (default 100) in searching for the eigenvalues. It may be necessary to increase MXCALC if you are searching for many eigenvalues in a single run. The units of EMIN, EMAX and DTOL 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 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 accordingly. The input energies are converted immediately into 1/cm using the appropriate conversion factor, and all energies subsequently printed by BOUND are in 1/cm. 2.5 Calculating expectation values In addition to calculating bound state energies, the bound program implements the calculation of expectation values using the finite-difference approach (Chem. Phys. Lett. 151, 565 (1988)). The present version can calculate expectation values of any function which is made up of a product of one of the angular functions in the potential function and a power of R; see the documentation on potential expansions for each collision type below. More complicated functions could easily be handled by modifying subroutine PERTRB. The expectation values are specified by a variable NPERT and arrays IPPERT, NPOW, DELTA and FACTOR in NAMELIST block &INPUT. For each bound state located, the program attempts to calculate NPERT (up to 20) expectation values; if the angular function associated with potential term I is A(I), the Jth value calculated is the expectation value of FACTOR(J) * A(IPPERT(J)) / R**NPOW(J) The program uses a finite difference calculation with a perturbation of DELTA(J) * A(IPPERT(J)). Small values of DELTA give smaller higher-order contributions to the result but poorer numerical stability; DELTA=0.001 cm-1 usually gives good results. The program uses the secant method, without a bisection step, in searching for the eigenvalues of the finite difference problem. It is thus possible for it to fail to converge, or to converge on the wrong eigenvalue. If this is a problem, try modifying RMID to a value that gives better convergence, or decrease DELTA. 2.6 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 modified log-derivative propagators. If ISCRU is greater than 0 (default 30), this save file is created on channel ISCRU and used for any subsequent energies. It can be very large. Note that this option saves CPU time at the expense of disc i/o. It is always advantageous for the modified log-derivative propagators, but for the basic log-derivative propagator 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. 2.7 Surface band structure calculations Atom - surface band structure calculations are not plagued by complications arising from total angular momentum and parity. However, the parallel component of the momentum must be specified. BOUND calculates parallel momenta from a polar angle THETA (measured from the surface normal) and an azimuthal angle PHI (measured relative to the surface reciprocal lattice vector g1), assuming an incoming beam with a wave vector of 100 inverse Angstroms. 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: DO 840 JTOT = JTOTL,JTOTU,JSTEP THETA = THETLW + THETST*JTOT DO 830 M = 1,MXPHI PHI = PHILW + PHIST*(M-1) .. .. Calculations for parallel momentum corresponding to angles THETA, PHI and wave vector 100 A-1. .. .. 830 CONTINUE 840 CONTINUE This rather peculiar method of specifying parallel momenta is used because the basis set generating routines were originally designed for MOLSCAT, which handles atom-surface diffractive scattering. 2.8 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 BOUND. For the log-derivative propagators used in BOUND version 5, the error in the eigenvalues due to the use of a finite step size is proportional to DR**4 in the limit of small step size. It is thus possible to obtain an improved estimate of the bound state energy by performing calculations with two different step sizes and using Richardson extrapolation. See the section on convergence testing below. Units of length _______________ BOUND operates with units of lengths specified by RM, which is set in subroutine POTENL (see section 6) 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. Range of integration ____________________ BOUND propagates the coupled equations outwards from RMIN to RMID, and inwards from RMAX to RMID. RMIN (default 0.8) is used exactly as input. RMAX (default 10.0) is used exactly as input. RMID must be specified, and to give acceptably fast convergence must be within the classically allowed region at the energy concerned. Although the value of RMID does not affect the calculated bound state energies, it does affect the rate at which the program converges to them. A value of RMID slightly outside the inner classical turning point usually gives rapid convergence from a reasonable distance away. Step size _________ In the BOUND program, all the propagators use a constant step size given by the input parameter DR. This may be slightly modified internally if necessary, to give an integer number of steps between RMIN and RMAX. RMID may also be modified slightly for the same reason. 2.9 Automatic convergence checking If the &INPUT variable NCONV is set greater than zero on input, BOUND performs NCONV extra calculations of each eigenvalue, with different values of DR, RMAX or RMIN. This is useful in checking the convergence with respect to these parameters, or in estimating a correction for the error due to the use of a finite step size. The aspect of convergence to be checked is selected with the variable ICON as follows: ICON selects the parameter to be varied; for each eigenvalue calculation: ICON = 1 doubles the step size ICON = 2 decreases RMAX by DRCON ICON = 3 increases RMIN by DRCON DRCON is the amount by which RMAX is to be decreased or RMIN is to be increased. BOUND prints the difference between the first calculation done (which should be the most accurate provided DRCON is positive) and subsequent calculations. For the case of DR convergence (ICON=1), BOUND uses Richardson extrapolation to extrapolate the results to zero step size: since the limiting error for the log-derivative propagators is known to be proportional to DR**4, BOUND can use the results of two calculations with different step sizes to estimate the actual error, and correct for it. If you need very high precision (step size errors less than 10**-3 cm-1), it is usually far better to use Richardson extrapolation based on two relatively large step sizes (with NCONV=1, ICON=1) rather than decreasing DR until the desired precision is achieved directly. 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 monomer 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 = 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. 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 NLEVEL If NLEVEL > 0, it is 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 as described for the individual monomer 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 400, although it is conceptually two-dimensional for ITYP > 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 200. 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. 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. 3.2 ITYP = 2 For atom - vibrating diatom systems, 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 is not used for ITYP = 2. Note that energies calculated from WE, etc. are set with the zero of energy at v = j = 0. 3.3 ITYP = 3 For linear rigid rotor - linear rigid rotor systems, 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. It is also possible to specify the statistical weights to be applied to the different symmetry combinations. 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. This information is actually intended for the MOLSCAT program, rather than BOUND; its only effect for BOUND is that the program will not perform a 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. 3.4 ITYP = 4 Not used in this version of BOUND 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=|KSET| are included; if KSET is zero or poisitive, 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); 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. 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 NLEVEL is 0 and all three of 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. Alternatively, the program (subroutine SET6) may obtain the wavefunctions by reading data from unit IASYMU (specified in &BASIS, default = 5, the standard input unit). If IASYMU is the standard input unit, data cards should follow the &BASIS ... &END set and precede the &POTL set. If NLEVEL is set in &BASIS, this number of rotor descriptions will be expected. (If IASYMU is an auxiliary unit, it may be convenient to read levels until an end-of-file condition occurs.) Each level is described 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 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. 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 version 11, 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 >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 helicity decoupling 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. 3.7 ITYP = 7 The basis set for ITYP = 7 is specified in the same way as that for ITYP = 2. 3.8 ITYP = 8 For atom - corrugated surface systems, 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 BOUND's methods of generating basis sets for surface systems, it is important to realise that the operation is carried out in two steps. At the time that the &BASIS NAMELIST block is read, BOUND sets up a master list of basis functions which MAY be included in subsequent 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, BOUND 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 > 0 and the JLEVEL array is read in explicitly. Subsequently, for particular values of THETA, PHI and ENERGY, BOUND 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 BOUND does not allow the master list to contain more than 200 functions. For the common case of surface systems with parallel momentum along a symmetry direction, BOUND 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 cases with zero parallel momentum involving 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 on request. 3.9 ITYP = 9 Selecting this option makes the program call the user's own basis-set generating routines; the programmer decides how to use the various arrays. Further information is available from JMH on request. 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 JZCSMX is set > -1 on input, m will be limited by JZCSMX instead of JMAX. 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 desirable to include the diagonal part of the Coriolis coupling exactly, and to restrict the basis set to functions with a body-fixed projection number not greater than the total angular momentum JTOT. This is achieved if IBOUND = 1 on input. Note that this is NOT the default. This option is supported only for ITYPE = 21, 22, 23 and 27; for other values of ITYPE, setting IBOUND not equal to zero is likely to produce incorrect results. Another effect of IBOUND .EQ. 1 is to suppress calculations for projection quantum numbers greater than JTOT. 4. NAMELIST block &POTL _______________________ Subroutine POTENL is used by BOUND to evaluate the intermolecular potential function at a given distance R as an expansion in appropriate angular functions. The general-purpose version of POTENL supplied with BOUND is described in this section; it is adequate for simple potentials whose coefficients can be expressed as a sum of exponentials and inverse powers, but for more complicated potential forms users must supply their own potential routine(s). This may be done either using the VSTAR mechanism described in section 4.2, or by replacing the POTENL routine completely. The latter method is often more efficient, since it is often necessary to repeat a lot of computation if the different potential terms must be evaluated separately. The specification required from POTENL is described in detail in section 5. The parameters that may be input in &POTL NAMELIST are as follows: RM is the length scaling factor that will be used by BOUND. 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. 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 BOUND'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 6.1 for more details. Note that potential symmetries specified in LAMBDA may be in any order and may be repeated, except for ITYP=7. LMAX IF LMAX .GT.0 a LAMBDA array is generated internally instead of using the input array. LMAX is interpreted as the highest-order Legendre polynomial or spherical harmonic to include. POTENL also reads IHOMO, and uses it as the step size for this list. Note that IHOMO may be overwritten in VRTP. MMAX For ITYP = 5 or 6, MMAX is the maximum projection of the spherical harmonics to be included in the expansion. POTENL also reads ICNSYM and uses it for the step size in M. Note that ICNSYM may be overwritten in VRTP. IVMAX For ITYP = 2, IVMAX is the largest diatom vibrational quantum number for which potential matrix elements are required. 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 is specified at particular angles. NPT is the number of Gauss-Legendre points (in theta) used in projecting out the angular components of the potential if LVRTP = .TRUE.. NPT should be somewhat greater than the largest value of LAMBDA required in the potential expansion. Note that, if IHOMO .EQ. 2 (see below) then only the first (NPT+1)/2 quadrature points are actually used. NPS is the number of secondary quadrature points (in r or phi) used in projecting out the components of the potential if LVRTP = .TRUE.. For ITYPE = 5 or 6, NPS should be somewhat larger than the largest value of m in the potential expansion; for ITYPE=2, NPS should be somewhat larger than the largest value of v (but note that the version of the Gauss-Hermite routine supplied with the program at present supports only NPS = 2, 3, 4, 5 and 10; this could be generalised simply by substituting another routine for supplying Gauss-Hermite points and weights). If ICNSYM .GT. 1 (see below), then NPS points are used on the range 0 < phi < 2*pi/ICNSYM. NTERM is an array of MXLAM integers, describing how to evaluate the potential if LVRTP is .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 BOUND 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 ITYPE=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 propagator is being used in a mode requiring explicit derivatives (see section 2.12). 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 ITYP = 1, 2, 5 and 6. The option is invoked by specifying LVRTP = .TRUE.; POTENL then expects the potential to be supplied by a routine VRTP. 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 should also set the variables IHOMO and ICNSYM in the common block ANGLES, if they have not already been set by data read in subroutine POTENL: COMMON /ANGLES/ COSANG(3),FACTOR,IHOMO,ICNSYM 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. 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. COSANG(1) is cos(theta) and, for ITYP = 5 and 6, COSANG(2) is phi. For 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). POTENL projects out the components of the potential using Gaussian quadratures; note that the ITYP=2 code ASSUMES that the diatomic molecule is a harmonic oscillator. Some propagators do not use POTENL/VRTP calls with IDERIV .GT. 1, and for these it is adequate to supply a VRTP routine that traps calls with IDERIV .GT. 0, print an error message, and exits. The VRTP mechanism is not supported for ITYP=3 cases in the present version of BOUND. A sample VRTP routine is supplied with BOUND. 5. User-supplied potential 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 BOUND. 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 is 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 to 6 NPOTL=MXLAM ITYP = 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. BOUND 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. 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 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) 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 BOUND. 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 BOUND, 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 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 BOUND. 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 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. Array dimensioning _____________________ The main program for BOUND declares one large array, X, which it passes to the principal routine DRIVER. DRIVER and other routines then partition up the X array according to the size of the problem being tackled. Thus, very few of the arrays used internally by BOUND are explicitly dimensioned, and it is seldom necessary for the user to concern himself with array dimensioning. If BOUND 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. 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. 7. Subroutines of the BOUND program ___________________________________ Most subroutines in BOUND 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. However, the list is not entirely up-to-date. ASROT Subroutine called by SET6C Calculates asymmetric top energy levels 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 BASE Subroutine called by DRIVER Handles generation of basis set BASIN Entry point in subroutine BASE, called by DRIVER Initialisation entry for BASE - reads &BASIS data BASE9, BASIN9, SET9, CPL9 User-supplied subroutines called by BASE etc. to specify the details of a special coupling case BDPROP Subroutine called by DRIVER Calls propagators and evaluates matching determinant BINOM Function called by THREEJ Provides binomial coefficients CHKSTR Subroutine called by DRIVER Checks that there is enough storage left in the X array COUPLE Subroutine called by BASE Handles angular momentum coupling elements CPLOUT Subroutine called by BASE Output routine for angular momentum coupling elements CPL21 Subroutine called by MCGCPL Coupling matrices for ITYPE = 21 CPL25 Subroutine called by MCGCPL Coupling matrices for ITYPE = 25 CPL26 Subroutine called by MCGCPL Coupling matrices for ITYPE = 26 CPL3 Subroutine called by COUPLE Coupling matrices for ITYPE = 3 CSRTRT Subroutine called by MCGCPL Handles angular momentum for coupled states rotor-rotor collisions DAPROP Subroutine called by BDPROP Propagation routine for the diabatic modified log-derivative method (INTFLG = 6). DEGENF Entry point in subroutine BASE, called by Returns degeneracy factor for denominator of cross-section expressions DGEMUL Subroutine called by QAPROP TRNSFM Matrix multiplication DMSYM Subroutine called by ASROT Symmetrizes spherical top wavefunctions DRIVER Subroutine called by BOUND main program Driver reads &INPUT data and handles the calls to the major routines. ECNV Subroutine called by BASE DRIVER Handles different possible energy units and returns conversion factor. ESYMTP Subroutine called by COUPLE SET6 F02AAF Subroutine called by BDPROP Eigenvalues of a real symmetric matrix F02ABF Subroutine called by ASROT DMSYM QAPROP ZERVEC Eigenvalues and eigenvectors of a real symmetric matrix FCOEF Subroutine called by COUPLE Calculates Percival-Seaton angular momentum coupling coefficient FIND Subroutine called by SURBAS Finds potential term coupling two G-vectors in atom - surface systems FSYMTP Subroutine called by COUPLE SET6 GCLOCK Subroutine called by DRIVER Timing routine GSYMTP Subroutine called by MCGCPL SET6 IDPART Subroutine called by BASE Handles identical particle symmetry for ITYP = 3 IPASYM Subroutine called by SET6 Checks symmetry of asymmetric rotor coefficients (ITYPE = 6) J3J000 Subroutine called by COUPLE CPL21 SET6 Recursive routine for Wigner 3-j symbols with zero projections 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 LDPROP Subroutine called by BDPROP Log-derivative propagator MATPRN Subroutine called by BDPROP DRIVER Matrix printing routine MCGCPL Subroutine called by BASE Handles angular momentum coupling for McGuire-Kouri coupled states approximation and helicity decoupling ORDER Subroutine called by FIND Takes account of lattice symmetry for Fourier components of atom-surface potential PARITY Function called by BASE COUPLE CSRTRT ESYMTP FCOEF FSYMTP GSYMTP MCGCPL ROTROT SET6 THREEJ Returns (-1)**argument POTENL Subroutine called by DRIVER WAVMAT Subroutine to evaluate intermolecular potential at distance R QAPROP Subroutine called by BDPROP Propagation routine for the quasiadiabatic modified log-derivative method (INTFLG = 7) ROTROT Function called by COUPLE Calculates angular momentum matrix element for rotor-rotor collisions 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 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 top basis sets from rotational constants 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 SIXJ Function called by FCOEF FSYMTP ROTROT Calculates Wigner 6-j symbol SURBAS Subroutine called by BASE Handles basis set for atom - surface systems (ITYPE = 8) SYMINV Subroutine called by DAPROP LDPROP QAPROP 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 Calculates Wigner 3-j symbol with no restriction on projections TRNSFM Subroutine called by QAPROP Performs a similarity transformation. TRNSP Subroutine called by QAPROP Transposes a matrix VINIT Subroutine called by the general-purpose version of POTENL. Initialises for calls to potential routines VSTAR, VSTAR1 and VSTAR2. 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. WAVMAT Subroutine called by DAPROP HEADER LDPROP QAPROP RMSET Calls POTENL to obtain potential coefficients, and forms potential matrix for propagator routines. WAVVEC Subroutine called by WAVMAT Efficient routine to handle vector products required by WAVMAT XNINEJ Function called by CSRTRT ROTROT Calculates Wigner 9-j symbol ZERVEC Subroutine called by DRIVER Evaluates and prints the wavefunction at the matching point 8. Machine-dependent features _____________________________ BOUND 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. 8.1 Main program The BOUND main program does not do any processing; it simply allocates storage and calls DRIVER to do all the work. If BOUND 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. 8.2 NAMELIST Most of BOUND'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. 8.3 Integer length BOUND 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. 8.4 Date and time routines BOUND 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. 8.5 CPU time routines BOUND 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 used so far, in seconds. In the last resort, it may return a result of zero. 8.6 Floating point underflow Various routines used by BOUND require that the result of an arithmetic operation causing floating point underflow should be zero. BOUND 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. 8.7 Linear algebra BOUND is usually distributed with FORTRAN versions of all the linear algebra routines. However, on many machines it is preferable to use specially optimized versions. Wherever possible, remove the FORTRAN versions and use calls to locally available library routines. This is particularly important on vector processors, where it is crucial to use vectorised versions of the BLAS routines. It is important that any user who implements new options in BOUND 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. The important routines called by BOUND and their functions are: DGEMUL Matrix multiplication SYMINV Invert symmetric matrix F02AAF Diagonalise symmetric matrix without eigenvectors F02ABF Diagonalise symmetric matrix with eigenvectors BOUND also calls BLAS (basic linear algebra subroutines) such as DAXPY, DDOT etc (or SAXPY, SDOT etc in the CRAY version) in many places. LAPACK and BLAS In version 5, BOUND 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 BOUND 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. The linear algebra routines supplied with BOUND 1) Matrix multiplication BOUND 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 BOUND calls SYMINV. Two versions of SYMINV are recommended. The first is a pure FORTRAN routine, which on most machines is faster than LAPACK equivalents for relatively small problems (N<70). The second calls the LAPACK routines DSYTRF and DSYTRI to carry out the inversion. For optimum performance, you need to test both routines on your machine. Note that BOUND 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) Eigenvalues and eigenvectors of symmetric matrices These routines are important for INTFLG = 7. BOUND calls diagonalisers by the NAG names F02AAF and F02ABF, and the NAG routines give acceptable performance on most machines. However, the distributed version of BOUND provides routines that simulate F02AAF and F02ABF by calls to the LAPACK routine DSYEVX. 8.8 File handling BOUND adheres to the FORTRAN 77 standard in its use of READ and WRITE statements (including direct access files). 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 BOUND. 9. References ______________ It would be difficult to give a complete set of literature references for all the various capabilities of BOUND. The following is a selected list for various aspects. General theory ______________ B.R. Johnson, J. Chem. Phys. 69, 4678 (1978). G. Danby, J. Phys. B 16, 3393 (1983). For a review of the coupled-channel method in general, see J. M. Hutson, Computer Physics Communications, to be published (1994). Preprints are available from JMH. For an introduction to the coupled equations as they arise for Van der Waals complexes, see J. M. Hutson, Advances in Molecular Vibrations and Collision Dynamics 1A, 1 (1991), published by JAI Press. ITYPE = 2 _________ S. Green, J. Chem. Phys. 70, 4686 (1979). ITYPE = 3 _________ S. Green, J. Chem. Phys. 66, 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). Coupled states (centrifugal sudden) approximation _________________________________________________ P. McGuire and D.J. Kouri, J. Chem. Phys. 60, 2488 (1974). Bound state applications ________________________ J.M. Hutson, J. Chem. Phys. 81, 2357 (1984). Expectation values __________________ J.M. Hutson, Chem. Phys. Lett. 151, 565 (1988).