A difficulty of this approach is the fact that the colliding species have an infinite set of rotational states, and it is necessary to truncate the expansion of the total wavefunction to some finite number of these. In general, the wavefunction will converge by including higher energy basis functions. It is typically necessary to include all levels which are energetically accessible at the collision energy of interest (open channels) plus some energetically inaccessible levels (closed channels). Convergence is slower for more anisotropic interaction forces and for interactions with stronger attractive forces. Note that S-matrix elements are only defined between open channels.
This approach, which is exact except for the truncation of the basis set expansion, is called the close coupling method and is computationally feasible only for systems which have a rather small number of rotational and vibrational levels accessible at collision energies of interest. By making some approximations to the coupling terms it is possible to decouple the problem into smaller blocks and MOLSCAT is equipped to do this for several decoupling schemes. Of these, the coupled states approximation has been found to be reasonably accurate, especially for systems dominated by short-range forces and at higher collision energies. The infinite order sudden (IOS) approximation has been found to be useful in cases where the rotational energy spacings are small compared with the collision energy. For systems where coupled channel methods are not feasible other techniques such as classical trajectory calculations or semi-classical methods should be used, but MOLSCAT does not handle such calculations.
Besides expansion basis functions and an interaction potential, it is also necessary to provide input data which specify the collision energies, the method to use for integrating the coupled equations, approximate coupling scheme, optional processing, etc.
In MOLSCAT, the input data are divided into three sets of NAMELIST input. NAMELIST input is not standard FORTRAN but is implemented on most platforms and is too convenient to forego. In general it consists of data cards of the form, &<name> data1=value1, data2=value2, ... &END where <name> is the name associated with the input set; data1, data2, etc. are names of allowed variables in that set; and &END specifies the end of data for this set. &<name> generally must begin in column 2 and all other input must begin in column 2 or beyond, but, being nonstandard FORTRAN, the implementation specifics for a given platform may vary and local documentation should be consulted. One of the advantages of NAMELIST input is the ability to have default values for variables and many MOLSCAT input variables do, in fact, have sensible defaults and need not be included in the input data. The three sets of required data are &INPUT, &BASIS, and &POTL, and they are expected in this order. The &INPUT contains general control variables and &BASIS and &POTL describe the expansion basis set and interaction potential, respectively. Appropriate input for the latter two depends strongly on the kinds of collision species.
The NAMELIST data format is described more completely in Section 1.1 of the complete documentation.
V(R,theta) = sum-over-L V_L(R) * P_L(cos(theta)).(For more complicated collision systems there may be more than one sensible angular expansion available, but it is important to use the one expected by the MOLSCAT code; these are described in Section 5.1 of the full documention.)
Consider a simple "asymmetric Lennard-Jones" potential which was
used in early work on the
CO-He system (S. Green and P. Thaddues, Astrophys. J.
205, 766 (1976))
and which included only P_0, P_1, and P_2 Legendre polynomials.
The spherical term was
described by a Lennard-Jones function, and the anisotropic terms were
described by similar functional forms:
V_0(R)/EPS = (RM/R)**12 - 2 * (RM/R)**6 V_1(R)/EPS = -0.03 * (RM/R)**12 + 0.0073 * (RM/R)**7 V_2(R)/EPS = 0.2 * (RM/R)**12 - 0.34 * (RM/R)**6Here
RM
is the position of the minimum, approximately 3.5 Angstoms,
and EPS
is the well depth. approximately 21 1/cm.
The following &POTL parameters can be used to describe this potential.
First, MOLSCAT has to know how many angular terms there are; this is specified
by MXLAM
(MXSYM
is a synonym and may be used interchangeably):
MXLAM=3,It needs to know the indices for each of the terms, i.e. the degree of each Legendre polynomials. Note that the Legendre polynomials are each described by a single index; more complex systems may require several indices to describe each angular term. The input variable for these indices is
LAMBDA
and it should contain values for
the number of
terms specified by MXLAM
, in this case three symmetries times one term per
symmetry:LAMBDA=0,1,2,
Next, the radial coefficients, V_L(R) must be specified for each of these
angular terms. MOLSCAT has a built in mechanism for potentials which can be
described as sums of exponentials and inverse powers of the collision
distance, which is the case for this simple potential. The required
input variable is NTERM
which must be given a value for each of the angular symmetries
specified by MXLAM
.
The values specify how many exponential and/or inverse power terms are in each of
the angular functions. In this case there are two inverse powers in each:
NTERM=2,2,2,So far this has specified that there are six terms, two for each of the three angular functions. The inverse power for each of these terms is specified, sequentially, in the variable
NPOWER
:NPOWER=-12,-6,-12,-7,-12,-6,Note that the first two of these correspond to the P_0 term, the next two to the P_1 term, etc. (To specify an exponential, rather than an inverse power, set
NPOWER = 0
for that term.) The coefficient that multiplies each
inverse power must also be given, and these are input in the same order in
the variable A
. For the potential described above:A=1.,-2.,-0.03,0.0073,0.2,-0.34,A description of capabilities for potentials expressed as sums of exponential and inverse power terms is given in Section 4.1 of the full documentation.
Finally, note that there are scaling factors for both distance and energy
in the above potential and
these must be given in the variables RM
(in units of Angstroms) and EPSIL
(in units of 1/cm):
RM=3.5, EPSIL=21.,In cases where the potential terms do not have scaling factors like this, it is still necessary to describe the units of distance (in terms of Angstroms) and energy (in terms of 1/cm) in the variables
RM
and EPSIL
. Thus,
if distances are in Angstroms and energies are in 1/cm:RM=1.0, EPSIL=1.0,and, if distances are in Bohr radii (atomics units) and energies are in electron volts:
RM=0.529177, EPSIL=8065.7,Finally, for the above example, the complete required &POTL input is:
&POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2, NTERM=2,2,2, NPOWER=-12,-6,-12,-7,-12,-6, A=1., -2., -0.03, 0.0073, 0.2, -0.34, &ENDIn general NAMELIST variables can be given in any order and arrays can be filled with sequential values as shown. An equivalent input for LAMBDA would be
LAMBDA(1)=0, LAMBDA(2)=1, LAMBDA(3)=2,Because NAMELIST input is not part of standard FORTRAN, however, one should check local documentation for specific rules.
It is worth remarking that the angular symmetry terms may be listed in any order and (for most ITYP) they may be repeated. For example, suppose it were convenient for some reason to separate the repulsive and attractive terms in the above potential. One would have three repulsive terms (P_0, P_1, and P_2) and, similarly, three attractive terms, and each would consist of a single inverse power. One would then have the following input:
&POTL RM=3.5, EPSIL=21., MXLAM=6, LAMBDA=0,1,2,0,1,2, NTERM=1,1,1,1,1,1, NPOWER=-12,-12,-12, -6,-7,-6, A=1., -0.03, 0.2, -2., 0.0073, -0.34 &ENDwhich gives exactly the same interaction potential as the previous set. Note that the
NPOWER
and A
values have been
rearrange to correspond to the order in the LAMBDA
array.
It is generally less efficient to do the calculation with repeated
symmetries, and this discussion should be considered as an illustration
of the flexibility of the program and not as an indication of
recommended procedure. On the other hand, there is no penalty in
efficiency for specifying the symmetries in a different order.
ITYPE
,
is obtained
by adding to ITYP a number, IADD, which is a multiple of ten and which
specifies an approximate method.
In particular, for full close coupling IADD = 0, for coupled states,
IADD = 20,
and for the infinite order sudden approximation IADD = 100. For full close
coupling for the current case:ITYPE=1,The rotational basis set for this case requires specifying the rigid rotor quantum numbers, J, and this can be most easily done by setting
JMIN
, JMAX
,
and JSTEP
, which have rather obvious
meanings. Since the default values give JMIN = 0
and
JSTEP = 1
it is generally only necessary to specify
JMAX
.
The rotational energies can be specified by giving the usual
spectroscopic rotation constant (in 1/cm). For CO:BE=1.92265,Since rotational energies for a linear rigid rotor are given by
BE*J*(J+1)
it is readily verified that J = 5 is a closed channel at a collision
energy of 50 1/cm, which is the energy that will be chosen below for
this calculation. Therefore, to include all open channels and one
closed channel:JMAX=5,The following &BASIS data provide a complete description
&BASIS ITYPE=1, JMAX=5, BE=1.92265, &END
URED
. If the masses of the molecule and the atom are M1 and
M2, respectively, the reduced mass is M1*M2/(M1+M2). For CO-He,
with masses of 12, 16, and 4, respectively:URED = 3.5,It is necessary to specify the collision energies (these are total energy, kinetic plus internal rotor energy); generally
NNRG
specifies the
number of energies and the values (in 1/cm) are given in the array
ENERGY
.
To do calculations at a single, relatively low collision energy:NNRG=1, ENERGY=50.,Values for the starting and ending radial distances for the radial integration are generally required (in units of
RM
,
which is specified in &POTL), although
the default values are reasonable if RM
is approximately the distance of the
minimum, as it is here. Note that these values are normally used as guesses
and the program may adjust them, so it is not crucial to get precise
values. RMIN
should be well into the classically forbidden region,
however, and RMAX
should be far into the asymptotic region.
Reasonable values for the current system are:RMIN=0.7, RMAX=10.,It is necessary to choose a numerical method for solving the coupled equations. For many systems a good choice is the modified log-derivative method of Manolopoulos, J. Chem. Phys. 85, 6425 (1986), since it is a good compromise between speed and accuracy and is very easy to control. It is specified by:
INTFLG=6,Accuracy of this propagator is controlled by a single input variable,
STEPS
,
which is the number of steps per asymptotic wavelength of the radial
wavefunction. Good accuracy is achieved with values between 10. and 20.
unless the well depth is particularly large compared with the
collision energy.
Some other values which should generally be set because the defaults probably
do not provide what you desire include the following. PRNTLV
default is 0
which provides virtually no output; sensible values are probably from 3 to 5;
higher values give more detailed partial results. ISIGPR
default of 0 must be changed in order to print cross sections.
LABEL
can be set to a character string which will be printed
with the output to give some identification to the run. Thus:
LABEL='CO-He L-J potential of Green and Thaddeus', PRNTLV=3, ISIGPR=1,
It is generally necessary to do calculations for a range of partial waves; these
are the orbital angular momenta and correspond classically to an impact
parameter. They are labeled by the variable JTOT. JTOT = 0
corresponds to head-on collisions, and large JTOT
correspond to glancing collisions. For interactions which decrease sufficiently
rapidly with distance, contributions to the cross sections become negligible
for sufficiently large JTOT.
It is often desirable to specify the range of partial waves by specifying
JTOTL
, JTOTU
,
and JSTEP
as the lowest and highest JTOT and the increment,
respectively. However, the default values will cause the program to start
at zero and increment JTOT by steps of one
until state-to-state cross sections have achieved reasonable convergence.
Using automatic cross section convergence, the complete data set would be:
&INPUT URED = 3.5, NNRG=1, ENERGY=50., INTFLG=6, STEPS=10., RMIN=0.7, RMAX=10., LABEL='CO-He L-J potential of Green and Thaddeus', PRNTLV=3, ISIGPR=1, &END
&INPUT URED = 3.5, NNRG=1, ENERGY=50., INTFLG=6, STEPS=10., RMIN=0.7, RMAX=10., LABEL='CO-He L-J potential of Green and Thaddeus', PRNTLV=3, ISIGPR=1, &END &BASIS ITYPE=1, JMAX=5, BE=1.92265, &END &POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2, NTERM=2,2,2, NPOWER=-12,-6,-12,-7,-12,-6, A=1., -2., -0.03, 0.0073, 0.2, -0.34, &END
NNRG
and the input for
ENERGY
is expanded accordingly. For
example,NNRG=4, ENERGY = 120., 110., 100., 90.,It is generally best to put the highest energy first since the step size for some propagators is calculated from the first energy and an input variable,
STEPS
(as in the
INTFLG = 6
example above) and the highest energy will
give the most conservative estimate. When calculating multiple energies it
must be recalled that the expansion basis set should also be chosen to be
adequate for the highest energy as it will then likely be adequate for the
lower energies as well.
Most of the propagators can make use of a scratch data set, if supplied, to
reduce the work at energies subsequent to the first, and it is generally
advantageous to supply such a data set be using the variable ISCRU
.
For example, to request unit(2) as the scratch data set:
ISCRU=2,
MOLSCAT automatically accumulates state-to-state integral cross sections.
However, for many other purposes it is necessary to save the S-matrices on
a file for further processing. This is requested by setting the variable
ISAVEU
to the desired unit number. To request unit(13) as the save file:
ISAVEU=13,Both
ISCRU
and ISAVEU
are part of the &INPUT
data.
It is generally desirable to check that the parameters chosen for a propagator
give adequate accuracy. This can be done easily by running the program again
with a different value for STEPS
(in the case of INTFLG = 6
); larger values
should give higher accuracy. It also is easy to change to a different
propagator by simply changing the value of INTFLG
.
The propagators which appear to be most generally useful are
INTFLG = 6
, discussed above, INTFLG = 8
, and
INTFLG = 3
.
INTFLG = 8
is particularly useful for systems which have strong
long range interactions. It uses the same method as INTFLG = 6
at
short-range, and requires the STEPS
variable as above. It also
requires additional control variables, for switching to an Airy
propagator at long-range and for control
of the latter. RMID
is the distance at which it switches methods
(as usual, in units of RM
);
values somewhat beyond the potential minimum are recommended.
A value for TOLHI
, used to increase step size in the long-range
part, is also required, and values around TOLHI = 1.05
are useful.
INTFLG = 3
, the Walker-Light R-matrix propagator is another generally
useful method which also takes bigger steps at long-range. It is controlled
by DR
, which is the initial step size in units of
RM
. Values on the order
of 0.01 to 0.1 are generally useful. To control the increase of step size
at long-range a second parameter is required, RVFAC
.
Step size is increased
at distances greater than RVFAC
times the classical turning
distance; values of RVFAC
around 1.5 are generally useful.
If many calculations, or large calculations will be run, it is useful to try different propagators and vary the tolerance parameters to reach a compromise between accuracy and speed. See Section 2.12 of the full documentation for a more complete description of the different propagators and their required and allowed control parameters. See also Section 6 which describes MOLSCAT facilities for automating convergence testing.
Alternate mechanisms exist for specifying the rotational levels in the
&BASIS data. Rather than specifying a range with JMIN
,
JMAX
, and JSTEP
, one can specifically list the
levels using input parameters NLEVEL
to specify the number
of functions and JLEVEL
to list the rotational quantum numbers.
Thus, the following two data sets are equivalent:
&BASIS ITYPE=1, JMAX=5, BE=1.92265, &ENDand
&BASIS ITYPE=1, NLEVEL=6, JLEVEL=0,1,2,3,4,5, BE=1.92265, &ENDThis form can be used if one wanted, for some reason, to specify the basis in a different order. It also has the advantage of allowing the rotor energies to be specified in an input variable,
ELEVEL
, which is not allowed if the
basis set list is generated from limits like JMAX
.
For example, one could specify
&BASIS ITYPE=1, NLEVEL=6, JLEVEL=0,1,2,3,4,5, BE=1.92265, ELEVEL=0., 3.845, 11.5, 22.8, 38.4, 57.0, &ENDValues given by
ELEVEL
will override any value supplied by
BE
in this case.
This mechanism can be useful, for example, if the
rotational energies are perturbed for some reason from rigid rotor values.
For linear rotors excited by atoms (ITYPE=1
) use of
NLEVEL
, etc. is not generally helpful, but for more complex
collision systems, it often provides a desirable option. In general,
specifying a value for NLEVEL
will override values for
JMAX
, etc. and require input values for JLEVEL
.
To illustrate the use of this VSTAR mechanism, consider using it for
the V_1(R) and V_2(R) terms in the CO-He example. If NTERM
is
set to a negative value for a given angular term, the program will attempt to
obtain that coefficient from a routine called VSTAR (with entry
VINIT for initialization). The modified &POTL input would be
&POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2, NTERM=2,-1,-1, NPOWER=-12,-6, A=1., -2., &ENDNote that we still specify three symmetries, and must still give the Legendre function indices. Also, we specify that V_0(R) is still obtained from two inverse power terms. However, the 2nd and 3rd symmetry terms will be obtained by the VSTAR mechanism and the
NPOWER
and A
values
for these have been eliminated from the &POTL input.
The VSTAR routine supplied with MOLSCAT is a dummy. If called it will print an error message and terminate the program. Therefore it must be replaced with a user written routine. It is necessary to relink or reload the program so that the user supplied routines replace those supplied with MOLSCAT.
The following simple FORTRAN routine will provide the same potential as specified for the second and third angular terms of the CO-He potential described above.
SUBROUTINE VSTAR(I,R,V)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IF (I.EQ.2) THEN
V= -0.03*(1/R)**12 + 0.0073*(1/R)**7
ELSEIF (I.EQ.3) THEN
V= 0.2*(1/R)**12 - 0.34*(1/R)**6
ELSE
WRITE(6,*) ' VSTAR CALLED FOR ILLEGAL SYMMETRY',I
ENDIF
RETURN
ENTRY VINIT(I,RM,EPSIL)
IF (I.EQ.2.OR.I.EQ.3) THEN
WRITE(6,*) ' VSTAR INITILIZED FOR SYMMETRY',I
RETURN
ELSE
WRITE(6,*) ' VINIT CALLED FOR ILLEGAL SYMMETRY',I
STOP
ENDIF
ENTRY VSTAR1
ENTRY VSTAR2
WRITE(6,*) ' VSTAR1 OR VSTAR2 CALLED; NOT SUPPORTED'
STOP
END
Several things should be noted.
Calls to VINIT and VSTAR will specify the symmetry number in the first
argument and it should be checked to see if calls for this symmetry
are expected.
There will always be a VINIT initialization call for each symmetry for
which the VSTAR mechanism is being used and it will be called before any
calls to VSTAR. The
present routine simply prints an informative initialization message.
One could do other
initialization processing, read input data, etc.; the values of RM
and EPSIL
from &POTL are made available during the initialization call.
The value of R passed to VSTAR will be in units of RM
and the value of V returned will be multiplied by EPSIL
;
note that these scaling factors are not applied
by the VSTAR routine, itself. For some propagators first and/or second
derivatives are required and these are provided by entries VSTAR1 and VSTAR2;
it is generally not necessary to deal with these cases, but such calls
should be trapped, as is done here.
For coupled channel calculations MOLSCAT actually requires that the potential be expanded in terms of angular functions, but can do this expansion automatically via an appropriate Gauss quadrature using results from the VRTP routine. For some IOS calculations, the angle dependent potential may be used directly.
Details for specifying the VRTP routine are given in Section 4.3 of the full documentation. As an example, the following simple FORTRAN code will generate the CO-He potential discussed above.
SUBROUTINE VRTP(IDERIV,R,P)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION P(1)
C
C THE FOLLOWING NAMED COMMON BLOCK MUST BE PRESENT
C IT PROVIDES COMMUNICATION WITH POTENL
C THE FOLLOWING HAS DIMENSIONS FOR VERSION 14
COMMON /ANGLES/COSANG(7),FACTOR,IHOMO,ICNSYM,IHOMO2,ICNSY2
C
C LEGENDRE POLYNOMIAL FUNCTIONS
P0(CT)=1.D0
P1(CT)=CT
P2(CT)=0.5D0*(3.D0*CT*CT-1.D0)
C
IF (IDERIV.LT.0) THEN
C THIS IS THE 'INITIALIZATION' CALL
WRITE(6,*) ' VRTP INITIALIZED FOR CO-HE TUTORIAL POTENTIAL'
RETURN
ELSEIF (IDERIV.EQ.0) THEN
C THIS CALL REQUESTS CALCULATION OF POTL AT R, COS(THETA)
C GET INVERSE POWERS
RM1=1.D0/R
RM6=RM1**6
RM12=RM6*RM6
C ASSEMBLE V0, V1, V2
V0=RM12-2.D0*RM6
V1=-0.03*RM12+0.0073*RM6*RM1
V2=0.2*RM12-0.34*RM6
C COMBINE WITH LEGENDRE POLYNOMIALS
C FOR ITYP=1 COSANG(1) CONTAINS COS(THETA)
CT=COSANG(1)
P(1)=(V0*P0(CT)+V1*P1(CT)+V2*P2(CT))*FACTOR
C FACTOR IS AN ITYP DEPENDENT FACTOR SUPPLIED BY MOLSCAT IN /ANGLES/
ELSE
WRITE(6,*) ' VRTP NOT SUPPORTED FOR IDERIV',IDERIV
STOP
ENDIF
RETURN
END
Several things should be noted. During the initialization call to VRTP
(signalled by IDERIV < 0) R may contain RM
and P(1) may contain EPSIL
from the &POTL input;
alternatively, values for RM
and EPSIL
may be set here and returned during this
initialization call. Note that the latter will override the former.
The initialization call may do other processing, for example,
reading input data or scaling parameters to &POTL
input values for RM
or EPSIL
.
In subsequent calls to VRTP, IDERIV = 0 requests the potential and
IDERIV = 1 or 2 requests the first or second derivative; the latter two cases
generally do not have to be supported but such calls should be trapped as
indicated. For an evaluation call (IDERIV = 0) R contains the distance, in
units of RM
,
and the values for the appropriate angles are in COSANG; use of COSANG by
different ITYP is detailed in Section 4.3
of the full documentation. Finally, the calculated
value of the potential should be multiplied by FACTOR, an ITYP-dependent
geometrical factor which is passed in COMMON/ANGLES/, and
the resulting value should be returned in P(1).
To invoke the VRTP mechanism, appropriate values for
&POTL data must be specified. In particular LVRTP
must have the
value .TRUE.; it will be set to this value if MXLAM
is less than
or equal to 0 (default is 0). However, to specify projection of appropriate
angular terms (Legendre polynomials for ITYP = 1), MXLAM
and
LAMBDA
generally will be specified as above, so LVRTP
must be specifically included in the &POTL data:
LVRTP=.TRUE.,Also, the number of Gauss points for the numerical integration to project the Legendre components must be specified. A conservative value for this case is:
NPTS=7,Thus, the revised &POTL data to invoke the VRTP mechanism for this case is:
&POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2, LVRTP=.TRUE., NPTS=7, &ENDNote that
NTERM
, NPOWER
, and A
are no longer relevant.
For cases like this an alternative form for specifying the angular terms is
also available. One generally
desires all the allowed symmetries up to some maximum and it is possible
to specify this using LMAX
instead of MXLAM
and
LAMBDA
; then the default value of MXLAM = 0
automatically sets LVRTP
to .TRUE. so that input equivalent
to the above is given by
&POTL RM=3.5, EPSIL=21., LMAX=2, NPTS=7, &END
In general, the only change required to request a decoupling approximation is to modify ITYPE in &BASIS. For the CS approximation one adds 20 and for IOS one adds 100. However, these options do allow certain additional parameters which can be useful in a real calculation and some of these are discussed here.
The input parameter JZCSMX
can be used to
limit calculations to m blocks less than or equal to JZCSMX
.
As a simple example, in the CO-He calculations at 50 1/cm collision energy
and a basis designated by JMAX = 5, the j=5 level is closed; it is
sensible to then specify
&BASIS ITYPE=21, JMAX=5, JZCSMX=4, BE=1.92265, &ENDSometimes only cross sections out of the lowest rotational levels are required and
JZCSMX
can be used to avoid unnecessary
calculations.
The number of generalized cross sections to compute
is determined by input parameters in the &INPUT set. In particular,
for a linear rotor excited by an atom, LMAX
determines the
highest Q(L) value which will be computed. (If LMAX
is
not specified, the program will try to chose a value consistent with
the rotational levels specified in the &BASIS input, but this is
a less desirable method for controlling this parameter.)
The IOS approximation does not actually expand a system wavefunction
in terms of rotational basis functions; it is, in fact, equivalent
to expanding in the complete, infinite set of basis functions. Rather,
it calculates S-matrices as a function of rotor orientations and
it is the orientation dependence of the IOS S-matrices which determines
the generalized IOS cross sections. This orientation dependence is
handled with Gauss quadratures and, rather than specify a basis set,
it is necessary to specify the number of points to use in the Gauss
quadratures (in fact, the number of points for each dimension of
the orientation dependence must be specified for more complex
collision systems). This value is specified by the input paramter
IOSNGP
in the &BASIS set. It should generally be
somewhat larger than the value of LMAX
.
For the CO-He example discussed initially, the following input would be appropriate for an IOS calculation:
&INPUT URED = 3.5, NNRG=1, ENERGY=50., INTFLG=6, STEPS=10., RMIN=0.7, RMAX=10., LABEL='CO-He L-J potential of Green and Thaddeus', PRNTLV=3, ISIGPR=1, LMAX=10, &END &BASIS ITYPE=101, JMAX=5, IOSNGP=16, &END &POTL RM=3.5, EPSIL=21., MXLAM=3, LAMBDA=0, 1, 2, NTERM=2,2,2, NPOWER=-12,-6,-12,-7,-12,-6, A=1., -2., -0.03, 0.0073, 0.2, -0.34, &ENDPoints to note:
LMAX
was chosen to provide all the Q(L)
that will be required to calculate state-to-state cross sections
among the levels specified by JMAX
, and IOSNGP
was chosen to give reasonable accuracy for this LMAX
. The
rotation constant has been removed from the &BASIS data; this is
not necessary, but this value will be ignored in an IOS calculation.
Finally, it should be noted that the WKB propagator INTFLG=-1
is often sufficiently accurate for IOS calculations and is generally
very fast. It essentially approximates a radial integral from the
classical turning point to infinity using Gauss quadrature. The number
of quadrature points is specified as a triplet of values in the
input parameter NGMP
in the &INPUT data set. The
first value is an initial number of points, the second value is an
increment, and the final value the maximum number of points. Calculations
are done first with NGMP(1) points and then with NGMP(1)+NGMP(2) points,
and continue by incrementing with NGMP(2) until two successive S-matrices
are converged to a specified tolerance, STEST, or until NGMP(3) is reached.
Typical effective values are:
NGMP=20, 3, 29, STEST=0.002,It should be recalled that the WBK propagator is only applicable to problems with a single classical turning point and should therefore be used only when the collision energy is large compared with the potential well depth.