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, MMAX

(default 0). For IOS calculations, these are the highest L and M 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 completed (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 (total energy equals collision kinetic energy plus asymptotic channel energy) at which scattering calculations are performed are controlled by the array ENERGY and the variables NNRG, NTEMP, TEMP, DNRG, EUNITS, and EUNITC. 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) or by the character*4 variable EUNITC (default is blanks) as follows:

 
      EUNITS = 1 or  EUNITC = 'CM'    1/cm
      EUNITS = 2 or  EUNITC = 'K'     Kelvin
      EUNITS = 3 or  EUNITC = 'MHZ'   MHz
      EUNITS = 4 or  EUNITC = 'GHZ'   GHz
      EUNITS = 5 or  EUNITC = 'EV'    eV
      EUNITS = 6 or  EUNITC = 'ERG'   erg
      EUNITS = 7 or  EUNITC = 'AU'    Hartrees
      EUNITS = 8 or  EUNITC = 'KJ'    kJ/mol
      EUNITS = 9 or  EUNITC = 'KCAL'  kcal/mol

If a nonblank character string (length 4 or less) is input via EUNITC, it will override any value input via EUNITS, and an attempt will be made to match one of the allowed units and set the conversion factor. This matching is insensitive to case and extraneous characters are generally ignored, so that EUNITC='1/CM' has the same effect as EUNITC='cm-1' or EUNITS=1. Likewise, both EUNITC='e.v.' and EUNITC='eV' cause conversion from electron volts.

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

These are described separately in Section 2.9 and Section 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,
 
      ((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 because there are no open channels for that energy, because it was skipped owing to IFEGEN.gt.1 in a pressure broadening calculation, 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 for pressure broadening (line shape) cross sections are controlled by the following variables in the &INPUT data: NLPRBR (synonym is IFLS), IFEGEN, LINE, and LTYPE.

Calculations 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:

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.

Parameters common to all propagators

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.

Parameters for specific propagators

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 a 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 willnot 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 than 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 synonym 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.

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.

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:
  • 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.
  • 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.

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.

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

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 STEST

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

STEST

(default = 0.0001) is the convergence tolerance for the WKB phase shift; values on the order of 0.001 are probably more suitable than the default value.

ONWARDOn to Section 3

TOPBack to the Table of Contents


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