Author: Chemistry Department
Last saved by: Theoretical Chemistry Group

VARIFLEX MANUAL
(Version 1.0; July 16, 1999)
Table of Contents
Overview	
More detailed description	
Input	
Rules for keyword input format	
Debugging options:	
Input Sections	
***** 1.  Title *****	
***** 2. Calculation Type *****	
A.  Ionic or neutral reactants	
B.  Reaction type.	
C.  Channels.	
D.  Level of Resolution in the Transition State Evaluations	
E.  Treatment of Collisional and Radiative Relaxation.	
F.  State-Counting.	
G.  Rotational Partition Functions.	
H.  Multiple Wells.	
I.  Collision Induced Dissociation Modeling.	13
J.  Various Restrictions.	13
***** 3.  Calculation Ranges  *****	
***** 4.  Experimental Data *****	18
***** 5.  Convolution Ranges *****	19
***** 6. CID *****	20
***** 7. Collision *****	21
***** 8. Complex *****	22
***** 9. Channel *****	23
***** 10. Fragments *****	24
***** 11. Transition State *****	27
A.  Input for keyword RRHOP	30
B.  Input for keyword PSTP	31
C.  Input for keyword VRCP	32
D.  Input for keyword LaplaceP	39
E.  Input for keyword ReadNP	40
***** 12. Channel 2, 3, etc. *****	41
***** End *****	41
***** Structure Data Module *****	42
***** Potential Input Module *****	45
1. CustomP.	45
2. VStretchP	46
3. VBendP	46
4. AtomAtomP	47
5. AnisotropicP	47
6. Ion-neutral interaction potential.	48
7. Inverse sixth power potential	48
***** Standard Hydrocarbon Radiative Association *****	49
Appendix 1	50
Definitions of coordinates for bending potentials in the built-in potential	50
Appendix 2	52
Output Contained in the Auxiliary Output Files	52
Appendix 3	54
Convergence Parameters in E/J and E Resolved Calculations	54

VARIFLEX MANUAL

Overview
	VariFlex is a program package for estimating rate constants of several types of gas phase reactions. Its focus is the application of transition state theory with an accurate treatment of the flexibility of the interfragment modes for variable definitions of the reaction coordinate.  Although aimed at calculating rates at this high level of sophistication for barrierless reactions, the package also allows convenient estimates of the rates using various alternative approaches such as phase space theory (PST), tight transition state theory, [i.e., with rigid rotor harmonic oscillator estimates (RRHO)], inverse Laplace transformations, and generic standard hydrocarbon estimates. 

The following reaction types may be considered:
I.  Bimolecular association reactions A + B => C.
II. Unimolecular dissociation reactions A => B + C.
III. Bimolecular reactions A + B => C + D.
IV. Competitive multichannel versions of each of the above, where more than one set of products is available, but involving only one complex or in the low pressure limit.

	For each of these reaction types one may consider either ionic or neutral species.

	The program focuses on the description of the temperature and/or pressure dependence of thermal reactions, where the reactants have a Boltzmann distribution of energies and angular momenta.  However, microcanonical and E/J resolved dissociation rates are also generally available as a subset of the thermal calculations (that is, rates for reactants of specified energy, or specified energy and total angular momentum).

Methods are described for calculating
I. High-pressure limiting thermal rates evaluated at the T, E, and/or E/J resolved levels.
II. Pressure-dependence of thermal rates using the RRKM or master equation or Troe's modified strong collision formalism.
III. Zero-pressure and low-pressure limit thermal rates.

More detailed description
This is a general program to calculate variable reaction coordinate/flexible transition state theory kinetic and dynamic estimates for barrierless reactions.  It implements the most recent algorithms for evaluating individual T, E, or E/J resolved rate constants.  These algorithms are coupled with master equation or Troe strong coupling analyses to provide kinetic estimates for a complete range of pressures. Procedures appropriate for both radical-radical and ion-molecule reactions are explicitly built in.
The basic variable reaction coordinate-transition state theory (VRC-TST) formalism is described in: 
[S. J. Klippenstein, Chem. Phys. Letts. 170, 71 (1990); J. Chem. Phys. 94, 6469 (1991); J. Chem. Phys. 96, 367 (1992); S. H. Robertson, A. F. Wagner, and D. M. Wardlaw, Faraday Disc. Chem. Soc. 102, 65 (1995).] 
This formalism is in turn based on the flexible transition state theory (FTST) formalism: 
[D. M. Wardlaw and R. A. Marcus, Chem. Phys. Letts. 110, 230 (1984); J. Chem. Phys. 83, 3462 (1985); J. Phys. Chem. 90, 5383 (1986); Adv. Chem. Phys. 70, 231 (1988).] 
The VRC-TST formalism has been implemented in this program at three different levels of resolution for evaluating the TS partition functions:
(i) The T resolved level [S. H. Robertson, A. F. Wagner, and D. M. Wardlaw, Faraday Disc. Chem. Soc. 102, 65 (1995); S. J. Klippenstein, J. Chem. Phys. 98, 6994 (1991)]. 
(ii) The E resolved level [S. J. Klippenstein, Chem. Phys.  Letts. 214, 418 (1993)] – in this first version this level is simply implemented as the integral over J of the E/J resolved results. 
(iii) The E/J resolved level [S. J. Klippenstein, J. Phys. Chem. 98, 11459 (1994)]. 
The following alternative procedures for estimating the TS partition functions are also available:
(a) Quantum phase space theory  (PST) [P. Pechukas and J. C. Light, J. Chem. Phys. 42, 3281 (1965); P. Pechukas, R. Rankin, and J. C. Light, J. Chem. Phys. 44, 794 (1966)].
(b) Tight transition state theory, i.e., rigid rotor harmonic oscillator (RRHO) based estimates.
(c) Inverse Laplace transformation [S. H. Robertson, M. J. Pilling, D. L. Baulch, and N. J. B. Green, J. Phys. Chem. 99, 13452  (1995)].
 (d) Classical phase space theory is not explicitly built in.  However, a close analogue may be obtained by simply incorporating the appropriate long range isotropic potential within the VRC-TST formalism and then considering the appropriate range of center-of-mass separation distances.  Alternatively, efficient PST estimates may be obtained by simply using large angular momentum step sizes in the quantum PST algorithm.
1-dimensional hindered and free rotor treatments of any particular mode are available for the separated fragments, the complex, and the transition state. Eckart and Wigner type tunneling corrections may be included for the transition state analyses.

The high pressure limiting thermal rate constants may be evaluated at any of the three levels of resolution and using any of the algorithms for evaluating the TS partition functions. When evaluated, the individual dissociation rate constants obtained at the E and E/J resolved levels are also printed out in secondary output files. 
Treatments for the pressure dependence require evaluations at the E and/or E/J resolved level and implement either: 
(i) The master equation formalism [R. G. Gilbert and S. C. Smith, Theory of Unimolecular Reactions, Blackwell, 1990; K. A. Holbrook, M. J. Pilling, and S. H. Robertson, Unimolecular Reactions, Wiley, 1996] or
(ii) Troe's modified strong collision format: [J. Troe, J. Chem. Phys. 66, 4745, 4758 (1977)]. 
The various implementations of the master equation incorporated here are based on work described in:
[M. A. Hanning-Lee, N. J. B. Green, M. J. Pilling, and S. H. Robertson, J. Phys. Chem. 97, 860 (1993); N. J. B. Green, S. H. Robertson, M. J. Pilling, J. Chem. Phys. 100, 5259 (1994); S. H. Robertson, M. J. Pilling, D. L. Baulch, N. J. B. Green, J. Phys. Chem. 99, 13452 (1995); N. J. B. Green, Z. A. Bhatti, Chem. Phys. Lett. 254, 358 (1996).] 

For ionic reactions we have also built in various components of the standard hydrocarbon analysis [R. C. Dunbar , Int. J. Mass Spec. Ion Proc. 160, 1 (1997); 100, 423 (1990)].  This procedure provides generic estimates for the zero-pressure radiative association rate constant.

Also available as subcomponents of the calculations are E or E/J resolved estimates for  
(i) Dissociation rate constants.
(ii) Rovibronic densities of states for the complex. 
(iii) Transition states numbers of  states. 
(iv) Radiative rate constants.
and T resolved estimates for the 
(v) Partition functions for the complex and for the primary  fragments. 
(vi) The equilibrium constant. 

Input
We try to provide a complete description of the data below, and this description should be at least perused prior to usage.  However, the best way to prepare the data, at least initially, will be to make modifications to that sample data file which most closely matches your own desired calculation.  We encourage you to follow this approach.  However, you must consider the convergence of your results with respect to E and/or J in any E or E/J resolved calculations.  A detailed description of such convergence procedures is provided in Appendix 3.  It is absolutely essential that one study this Appendix prior to running any E or E/J resolved calculations.

Rules for keyword input format
The input is keyword-based.  It is case insensitive, that is, uppercase and lowercase letters are not distinguished.  Different types of lines:  (i) Blank lines are permitted (except as noted below) and are ignored.  (ii) Comment lines beginning with ! are also ignored, but cannot immediately follow a keyword line with a data keyword.  (iii) Section header lines beginning with * denote the beginning of sections, and the program reads these to keep track of and ensure proper input flow.  (iv) Lines beginning with & can follow immediately after any section header line; these are read as title information, and are printed verbatim in the primary output.  (v) All lines not of these four types are read as keywords or data, depending on the context.  On a given keyword line (but not a data line), everything after a ! is comment, and is ignored.
Keyword lines contain first a single primary keyword followed on occasion by one or more qualifier keywords on the same line.  The primary keywords are of three types:
1. Program keywords (all end with P) give program instructions, and have no argument.  
2. Data keywords (all end with D) are followed by an argument consisting of a specific set of data read in free format (i.e., separated by spaces, or commas, or on separate lines).  This argument must be given on a separate line from the data keyword.
3. List keywords (all end with L) are followed by an argument consisting of multiple lines of data, where the number of items of data expected will usually be set by the keyword preceding the list keyword.  Unless noted otherwise, data items are again read in free format and may appear separated by commas or spaces on one line, or may be on separate lines. Again this argument must begin on a separate line from the data keyword.
The qualifier keywords are of the same three types and end in PQ, DQ, or LQ, with the P, D, and L components of these suffixes having the same meaning as for the primary keywords.  These qualifier keywords give instructions relevant to the primary program keyword, which appears on the same line.

Comments on input error messages
Error messages have been incorporated throughout suggesting which mandatory keyword was expected.  In some cases you may instead have been trying to input an optional keyword.  Thus, one should really read the error messages as an indication that something was wrong prior to the expected keyword point. Furthermore, these messages may in some instances be incorrect since it is difficult to test all possible circumstances for erroneous input information.

Debugging options:
Any of the debugging keywords may be inserted at any point in the input file, and the corresponding debugging output will be printed starting with that point in the program execution.  Currently there are four debugging keywords:
KeyDebugP is useful for debugging input files.  It causes the program to print out each keyword as it is read by the program.  This output is printed in the main output file variflex.out (if you are using the vrfx_run shell this file will instead be named ‘job’.out).
TSPFDebugP prints out information about the optimized Transition State Partition Function.  It is useful for debugging the actual VTST calculation.  It also provides additional information about the transition state that may be of interest to the user.  The output is in the file tspf.out (or ‘job’_tspf.out).  See Appendix 2 for more details.
TSPFRDebugP prints out information about the Transition State Partition Functions for each of the reaction coordinates (R) sampled in the VTST evaluations.  It is useful for debugging the actual VTST calculation.  It also provides additional information about the transition state that may be of interest to the user.  The output is in the file tspf.out (or ‘job’_tspf.out).  See Appendix 2 for more details.
EigVecDebugP prints out the steady state eigenvectors and the Boltzmann popultions when the eigenvector based solution to the master equation is employed.  These eigenvectors are normalized to the corresponding Boltzmann population for the lower energy considered in the master equation.  This output is printed in the main output file variflex.out (if you are using the vrfx_run shell this file will instead be named ‘job’.out).

Input Sections
The sections must appear in the order given below, but (except as noted for Mandatory sections) sections may be omitted.  The end of each section is detected by the appearance of a new * line, and the input file must end with the line *END (and should include the carriage return or end of line command).

***** 1.  Title *****
(Mandatory section).
This section begins with the line
*Title
which is the first line expected in the input file.  The only input expected in this section (other than comment lines and debugging keywords) are title lines beginning with &, which will be read and echoed in the output files.

***** 2. Calculation Type *****
(Mandatory section)
Begins with the line
*CalculationType
The basic choices here fall into the following categories.  Choose one keyword from each group.  But some groups may be omitted, depending on other choices.  See the general discussion above for guidance regarding which combinations of choices are currently implemented.  The order of keyword entry in this section is unimportant, aside from the qualifier keywords which must be on the same line as their corresponding keyword.

A.  Ionic or neutral reactants
	IonP
	NeutralP

B.  Reaction type.
	AssociationP  Treat only the association direction.
	DissociationP  Treat only the dissociation direction.
	ReversibleP  Treat both association and dissociation simultaneously.

C.  Channels.  
If more than one reaction channel is to be considered, use the keyword
	ChannelsD  Allows more than one product channel.  Requires the number of 				product channels as argument.  (Default is one channel).

D.  Level of Resolution in the Transition State Evaluations
	TResP Evaluate the TS partition function at the T resolved level
	EResP Evaluate the TS partition function at the E resolved level
	EJResP Evaluate the TS partition function at the E/J resolved level
StHydroP generic standard hydrocarbon estimate -- only appropriate for ionic association reactions.  If this keyword is given, no other Section 2 keywords will be used.  The input should include the temperature and binding energy parts of Section 3, and then should proceed to the special format specified in the Standard Hydrocarbon Radiative Association input section given below.

E.  Treatment of Collisional and Radiative Relaxation. 
HighPresP or equivalently HighEffP Consider only the high pressure/high efficiency limit for the reaction.
LowEffP or equivalently LowPresP  Consider only the low pressure/low efficiency limit for the reaction.
LimitJPQ Qualifier for LowEffP and LowPresP calculations which limits the range of J for certain integrals. (Default is no limit which can also be specified with NoLimitJPQ.)
	StrongCollisionP  Perform a strong collision (collision efficiency = 1) RRKM analysis.
TroeP  Perform a Troe style analysis incorporating both the Strong Collision and Weak Collision broadening factors as described in [J. Troe, Ber. Bunsenges. Phys. Chem. 87, 161 (1983); R. G. Gilbert, K. Luther, and J. Troe, Ber. Bunsenges. Phys. Chem. 87, 169 (1983)] 
ModifiedStrongP  Perform a modified strong collision RRKM analysis evaluating the collisional efficiency in terms of the Troe formalism described in [J. Troe, Ber. Bunsenges. Phys. Chem. 87, 161 (1983)]; but with no weak collision correction.
1DMasterP  Employ 1-D Master equation in treating the pressure dependence.  In all cases a simple exponential down model is employed for the energy transfer coefficients.  The radiative relaxation rates are evaluated from the individual populations in each of the E,J levels and the Einstein coefficients for absorption and emission as obtained from the IR intensities.  
PCompPQ, PReacPQ: Qualifiers for 1DMasterP which specify the particular form of J averaging employed in deriving the 1D master equation from E/J resolved rate constants. (Default is PCompPQ.)
PCompPQ implies the use of the Boltzman probability for the J distribution in the master equation 
k(E) = ∫ dJ k(E,J) P_comp (E,J;T) / ∫ dJ P_comp (E,J; T);
P_comp(E,J) = rho_comp(E,J)*exp(-E/k_B T) / Q_comp(T)
PReacPQ implies the use of the product of the Boltzman probability for the complex with the stabilization probability for the J distribution in the master equation as described in [S. C. Smith and R. G. Gilbert, Int. J. Chem. Kinet. 20, 307,979 (1988)].  
k(E) = ∫ dJ k(E,J) P_comp(E,J;T) P_stab(E,J) / ∫ dJ P_comp P_stab;
P_stab(E,J) = ω / [ω + k(E,J)]
The present implementation differs however in the choice of total E and total J as the separable variables instead of E_vib and E_rot.

EigenvaluePQ, EigenvectorPQ, InversionPQ, DiffusionPQ: Qualifiers for 1DMasterP which specify the algorithm employed in solving the master equation.  For single channel calculations any may be used and the default is EigenvectorPQ.  For multichannel reactions the dissociation part is always treated with the eigenvector based routine and the association part can be treated with either the eigenvector or the inversion routine.  Note that the inclusion of the f_ne factor [cf. Eq. 6.8.26 of Gilbert and Smith’s book] for the association rate requires the use of the eigenvector routine.  This factor differs from unity for high temperatures.
EigenvaluePQ implies the use of an eigenvalue solver based approach.
EigenvectorPQ implies the use of an eigenvector based approach.
	InversionPQ implies the use of a matrix inversion based approach.
	DiffusionPQ implies the use of a diffusion equation based approach (this approach is not currently functioning)
It is generally wise to try both an eigen based routine and the inversion based routine to ensure that the results you are obtaining are numerically stable.
The steady state population distributions can be printed out by adding the keyword EigVecDebugP
	2DMasterP (Not currently functioning)

F.  State-Counting. 
In the evaluation of vibrational state densities, the default approach is direct state counting by the Beyer-Swinehart algorithm.  The steepest-descents approximation may be used instead for either the density of states for the complex (by entering the keyword RhoCSteeP) and/or for the conserved mode contribution to the transition state number of states (by entering the keyword RhoTSSteeP).  The primary advantage of the steepest descent approach is that it provides a smooth state density even for cases where the direct count vibrational density fluctuates due to a sparsity of the vibrational levels.  Sparse vibrational state spacing may be a problem in calculations involving low binding energy, small molecules, or vibrational frequency sets with many degeneracies and few distinct frequencies. Using the steepest descents option will probably never give significant errors for the complex density of states, and if its use gives large changes in results compared with exact counting, one should consider the physical meaning of the sparsity.  Quite often the smoothing provided by the steepest descent algorithm will be appropriate.  However, some of the algorithms for the conserved mode contribution to the transition state number of states rely on the proper fluctuations in the vibrational density and so RhoTSSteeP should only be used with great caution.  It may also be important to use the direct count algorithm when performing master equation analyses of the pressure dependence.  A secondary advantage of the steepest descent approach is that there are some instances where it will provide a significantly faster calculation, such as when many (e.g., 1000) vibrational modes are considered in a radiative rate calculation. 
Note: The steepest descent algorithm will not function properly when some of the modes are treated as free-rotors or as hindered rotors.
G.  Rotational Partition Functions.  
The rotational contributions to the canonical partition functions for the complex and for the fragments are generally evaluated classically.  Direct quantum counts may be employed for these partition functions by entering the keywords QuanPFCompP (for the complex), QuanPFFrag1P (for the first fragment) and/or QuanPFFrag2P (for the second fragment).   Note, however, that no nuclear spin (e.g., ortho vs. para) will be considered unless specified in the fragment input.
H.  Multiple Wells.  
Low pressure RRKM analyses for bimolecular reactions passing through multiple complexes may be performed by indicating the connectivity of all the different complex wells to each other and to the reactants and products. This requires first the keyword 
MultiWellP
Followed by the keyword
NWellD whose argument is the number of complexes to be considered (NWELL) and the number of product channels to consider (NPROD).  
Then give the keyword
ReacWellConnectD whose argument is a list of NWELL integers.  For each of the integers, a value of 1 implies that the reactants are directly connected to that well and a value of 0 implies that they are not.   Explicit entrance channel transition state evaluations will be performed for each of those wells that are connected to the reactants.  Typically, all but the first value will be 0. 
The next keyword must be 
WellProdConnectD whose argument is a list of NWELL x NPROD integers describing the connectivity of each well to each of the products.  In this input, the outer loop is over the wells and the inner loop is over the products.  I.e., one first inputs the connectivity for the first well to each of the products.  Then the connectivity for the second well is given, and so on, until the connectivity for each of the wells has been given.  Explicit exit channel transition state evaluations will be performed for each of the wells connected to products.
The next keyword must be
WellWellConnectD whose argument is a list of NWELL x (NWELL – 1) / 2 integers describing the connectivity of each well to each of the other wells.  In this input, the lower triangular portion of the matrix is cycled over.  I.e., the connectivities are read in the order W_2W_1, W_3W_1, W_3W_2 W_4W_1, …Transition state evaluations will performed for each of the pair of wells which are connected.  
The final keyword given is 
ReacProdConnectD whose argument is a list of NPROD integers describing the connectivity of the reactants to each of the products. Such connectivities would be appropriate for abstraction reactions for example.
The input of the data for the transition states should follow the order of the abovedescribed input for the connectivities.

I.  Collision Induced Dissociation Modeling
The keyword CIDP results in a collision induced dissociation modeling calculation along the lines of Eq. (2) in Rodgers, Ervin, and Armentrout and Ervin (J. Chem. Phys. 106, 4499 (1997)], but with the current estimate of the rate constants.  Furthermore, the rate constants may be evaluated at the angular momentum J resolved level followed by an integration within Eq. (2) over the full range of J values. 

J.  Various restrictions
Standard Hydrocarbon calculations (StHydroP) are available only for Ionic Associations (IonicP and AssociationP).  Only one channel may be considered in T Resolved calculations (TResP) and in High Pressure calculations (HighPresP or HighEffP).   T Resolved (TResP) calculations can only consider the High Pressure limit (HighPresP or HighEffP). 2D Master (2DMasterP) equations must employ an E/J Resolved (EJResP) determination of the number of states.

***** 3.  Calculation Ranges  *****
(Mandatory section)
Begins with the line
*CalculationRanges
This section specifies the temperatures, pressures, radiative rates, binding energies, and energy transfer parameters that will be ranged over in the calculations. The parameters may be specified in any order.  Also, any number of these specifications may be skipped in which case the defaults will be employed.  The basic format is generally the same for each of the parameters and so we describe in detail only that for the temperature.  The radiative rate involves additional possibilities beyond the basic format.  The other exception is that the ranges of temperature and pressure may also be input as a list of pairs 

Temperatures:
Default is one temperature of 300 K
Keyword TRangeP
Followed by one of the following three qualifiers (on the same line):
StepDQ
Specifies a range of equally spaced temperatures.  The argument line consists of 3 numbers, which are the minimum temperature, the step interval between temperatures, and the total number of temperatures. 
ScaleDQ
Specifies a range of temperatures increasing by a constant factor.  The argument line consists of 3 numbers, which are the minimum temperature, the factor by which each temperature is larger than the preceding one, and the total number of temperatures.
ListDQ
Specifies that a series of temperatures will be read in.  Its argument (one number) is the number of temperature values to be read.  It must be followed by the keyword
ValuesL
after which is given the series of temperatures to be used.

Pressures:
Default is one pressure of 760 torr
Keyword PRangeP
Followed by one of the following three qualifiers (on the same line):
StepDQ     or
ScaleDQ    or
ListDQ     followed by
ValuesL

Temperature and Pressure pairs:
If this keyword option is used, the keywords TRangeP and PRangeP should not be used.  Default is one T,P pair of 300 K, 760 torr
Keyword TPRangeP
Give the keyword NTPPairsD whose argument is the number of sets of T,P pairs of data.
 Then give the keyword TPPairsL followed by the given number of pairs of T,P values

Radiative rates:
Default is one k_r of zero value
Keyword KRRangeP
Followed by one of the following five keywords (on the same line):
CalcPQ
Specifies that the needed radiative rates will be calculated subsequently from the infrared mode intensities.
StepDQ
ScaleDQ
ListDQ          followed by
ValuesL
StHydPQ indicates that the radiative emission rate is to be estimated from generic standard hydrocarbon rules
The keyword MINPROBDQ may also be added as a third keyword on the radiative keyword lines.  This keyword specifies the minimum mode excitation probability considered in the evaluation of the net radiative emission rate (default = 1.0e-05). 

Zero-point corrected binding energies D_0.
Energy to go from the vibrational/rotational ground state of the desired electronic state of the complex to the vibrational/rotational ground states of the fragments for channel 1.  Default is one D_0  value of 20000 cm^-1
Keyword DZeroRangeP
Followed by one of the following three keywords (on the same line):
StepDQ
ScaleDQ
ListDQ          followed by
ValuesL

E transfer parameters: 
Defines the model and parameters for the collision induced transition between initial and final  levels.  The default is an exponential down model for the energy transfer with an α  value of 100 cm^-1. 
Keyword ETransferRangeP
Followed by one of the following four keywords which specify the model used for the downwards energy transfer (ΔE =E_final-E_initial< 0; α>0):
ExpDownP implies the standard exponential down model; [P(ΔE) = exp(ΔE/α)]
GaussianDownP implies a Gaussian down model; [P(ΔE) = exp(-ΔE^2/α^2)]
PoissonDownP implies a Poisson down model; [P(ΔE) = (-ΔE/α)exp(ΔE/α)]
TwoExpDownP implies a sum of two exponentials down; [P(ΔE) = A_1exp(ΔE/α_1) + A_2exp(ΔE/α_2)]
VarExpDownP implies a sum of two exponentials down; [P(ΔE) = exp{-(-ΔE/α)^n}]

and the keyword NETransferDQ on the same line.  
One of the keywords EDependPQ and TDependPQ may also be provided on the same line to yield an energy transfer coefficient α which varies with E; α = α_0 + α_1*E;  or T; α = α_0 + α_1*T.  The reference energy for the energy dependent form is the threshold energy of the primary fragments.  The same reference energy is also employed in the printing of the calculated <ΔE_down>,  <_,ΔE_tot>, and  ΔE_rms  = [<ΔE^2>]^1/2 .

The next line should be 
ValuesL
followed by a list of values appropriate for the given keyword combination

α_0 (cm^-1)		for	(ExpDownP) 
α_0(cm^-1),  α_1 		for	(ExpDownP,EdependPQ)
α_0(cm^-1),  α_1 (cm^-1/K^-1)		for	(ExpDownP,TDependPQ)

α_0 (cm^-1)		for	(GaussianDownP)  
α_0(cm^-1),  α_1 		for	(GaussianDownP,EdependPQ)
α_0(cm^-1),  α_1 (cm^-1/K^-1)		for	(GaussianDownP,TDependPQ)

α_0 (cm^-1)		for	(PoissonDownP)  
α_0(cm^-1),  α_1 		for	(PoissonDownP,EdependPQ)
α_0(cm^-1),  α_1 (cm^-1/K^-1)		for	(PoissonDownP,TDependPQ)

α_0 (cm^-1), n		for	(VarExpDownP) 
α_0(cm^-1),  α_1, n 		for	(VarExpDownP,EdependPQ)
α_0(cm^-1),  α_1 (cm^-1/K^-1) , n		for	(VarExpDownP,TDependPQ)

A_1, α_1,0(cm^-1), A_2, α_2,0(cm^-1) 		for	 (TwoExpDownP) 
A_1, α_1,0(cm^-1), A_2, α_2,0(cm^-1), α_1,1, α_2,1   for	 (TwoExpDownP,EdependPQ) 
A_1, α_1,0 (cm^-1), A_2, α_2,0 (cm^-1), α_1,1 (cm^-1/K^-1), α_2,1 (cm^-1/K^-1)			
for   	(TwoExpDownP,TdependPQ) 

***** 4.  Experimental Data *****
(Optional Section)
Begins with line
*ExperimentalData
This section allows for the input of experimental data for the temperatures and/or pressures at which the rate constants are to be calculated.  These data are used in evaluating root-mean-square deviations between theory and experiment for each of the parameter ranges.   

For comparison with high pressure experimental data give one of the optional keywords RevHighPExptD or BiHighPExptD or UniHighPExptD followed by the corresponding set of experimental data
RevHighPExptD is followed by pairs of data corresponding to the high pressure association and dissociation rate constants for each of the temperatures to be considered.
BiHighPExptD is followed by the bimolecular high pressure association rate constants for each of the temperatures to be considered.
UniHighPExptD is followed by the high pressure unimolecular dissociation rate constants for each of the temperatures to be considered.

For comparison with pressure and temperature dependent experimental data give one of the optional keywords RevExptD or BiExptD or UniExptD followed by the corresponding set of experimental data
RevExptD is followed by sets of data corresponding to the association and dissociation rate constants for each of the channels and each of the temperature/pressure pairs to be considered.  For each temperature and pressure combination, pairs of association and dissociation rate constants are given for first channel one and then channel 2, etc. 
BiExptD is followed by the bimolecular rate constants for each of the channels and each of the temperature/pressure pairs to be considered. Again the inner loop of this input corresponds to the different channels and the outer loop corresponds to the different temperature and pressure pairs. 
UniExptD is followed by the unimolecular dissociation rate constants for each of the temperature/pressure pairs to be considered. Again the inner loop of this input corresponds to the different channels and the outer loop corresponds to the different temperature and pressure pairs. 

***** 5.  Convolution Ranges *****

(Mandatory for all E Resolved and E/J Resolved calculations.)
Begins with line
*ConvolutionRanges
This section specifies the ranges of energy and angular momentum to be used in the convolutions over the thermal distribution of reactants. No defaults are currently included.Thus an ERangeP keyword is mandatory for E Resolved calculations and both an ERangeP and a JRangeP keyword are mandatory for E/J Resolved calculations.  The basic format is the same as in the calculation ranges section.The ordering of the E and J data is immaterial.

Energies:
Keyword ERangeP
Followed by one of the following three keywords:
StepDQ
ScaleDQ
ListDQ       followed by
ValuesL

Total Angular momenta: must be input as integer values
Keyword JRangeP
Followed by one of the following three keywords:
StepDQ
	ScaleDQ
ListDQ      followed by
ValuesL

One must generally determine the appropriate ranges of E and J by trial and error.  I.e., the calculation must be run for different ranges and then the output should be examined for convergence.  The output file intgl.out, which illustrates the contributions to the E and J integral for high pressure or strong collision or modified strong collision calculations, is often of some use in examining this convergence.  See Appendix 3 for further discussion.

***** 6. CID *****
[Only used in Collision Induced Dissociation Modeling] 
Begins with line
*CID
This modeling follows that described in [M. T. Rodgers, K. M. Ervin, and P. B. Armentrout, J. Chem. Phys. 106, 4499 (1997); K. M. Ervin and P. B. Armentrout, J. Chem. Phys. 83, 166 (1985)].
Give the line TauGamCIDD
Followed by the pair of values for τ and γ describing the experimental observation timescale and the width of the kinetic energy distribution of the ion P(ΔE) ∝ exp(-ΔE^2/γ^2).  A related neutral parameter is taken as M_reactants/(M_reactants + M_collider)
Then give the keyword NCIDSetsD followed the number of sets of CID parameters.
Then give the keyword CIDPairsL followed by the list of sets of n and σ, where σ implicitly has units of Å^2/(eV^n-1).  

***** 7. Collision *****
[Mandatory except in the High Pressure limit (HighPresP or HighEffP)]
Begins with line 
*Collision
Give the mass of the collider (bath gas) with keyword 
	CollMassD
	whose argument is the mass in g mol^-1.
Neutrals and ions are different:
	A. For ions (keyword IonP given in Section 2) give the following:
		Collider polarizability, using keyword 
			PolzD
whose argument is the polarizability in Å^3, 
		and the collider dipole moment (if non-zero), using keyword 
			DipoleD
whose argument is the dipole moment in Debye. 
	B. For neutrals, (keyword NeutralP given in Section 2) give two sets of Lennard-			Jones parameters for the self interactions of the complex and the collider.
			SigmaD, 
		whose argument is the sigma parameters for the complex and the collider in Å.
EpsilonD, 
whose argument is the epsilon parameters for the complex and the collider in units of cm-1.
For master equation analysis a pressure for performing the low pressure limit analysis may be specified with the keyword PlowPD followed by its argument on the next line.  The default pressure for the low pressure limit analysis is 1.0e-10 torr for master equation analyses.  For other pressure treatments the low pressure limit is obtained analytically.

***** 8. Complex *****
Mandatory except in the High Pressure limit of an Association calculation ((HighPresP or HighEffP) + AssociationP)
Begin with line
*Complex
Then
1.	 Input the structure data for the complex following the standard form in the Structure Data Module instructions (given below).
2. 	The vibrational step size for the complex density of states calculation may be input here.  This is done by giving the keyword 
	EVCStepD
followed by the value. (Default is 1.0 cm^-1)
3. 	A correction factor for the complex density of states and partition function may be supplied here if desired.  Use the keyword 
	RhoCorrD
followed by the value.  (Default is 1.0)
4. 	Either spherical top or symmetric top expressions may be used in evaluating the density of states for the complex by giving either the keyword 
RhoSpherTopP  
or the keyword
RhoSymmTopP 
The default is to use the symmetric top expressions. 
5.	Similarly, either spherical top or symmetric top expressions may be used in evaluating the radiative emission rates by giving either the keyword  
KRadSpherTopP  
or the keyword
KRadSymmTopP
The default is to use the spherical top expression. 
In some instances, use of the spherical top expressions for the energies may save a substantial fraction of the CPU, particularly in the radiative emission rate evaluations. However, when the energy is low and/or the molecule is very nonspherical then the results can differ substantially. Geometric means of the rotational constants are used in the symmetrizations.
The spherical top and symmetric top keywords are irrelevant for linear tops.

***** 9. Channel *****
(Mandatory Section.)
This section begins the input which is specific for each reaction channel. Begin with the keyword for the appropriate section, that is, input one of the section keywords
*Channel1  or
*Channel2  or
*Channel3  etc.

This is followed by the input of the fragment data (as necessary), the transition state data, and the potential data (as necessary) for the given channel, as specified in the following sections. 

***** 10. Fragments *****
This section inputs geometries and other properties of the two fragments for the given channel. 
For the first channel this section is mandatory with the exception of high pressure dissociation calculations which do not involve either PST or VRC. That is, the requirements for omitting Section 10 for channel 1 are (1) HighPresP or HighEffP, and (2) DissociationP, and (3) no PSTP, VRCP, or VRCEffP transition state treatments (cf. Sec. 11, below).  For the other channels this section is mandatory only for PST (PSTP) or VRC (VRCP or VRCEffP) transition state treatments.  This section may always be included and any irrelevant information will simply be ignored.
 Begin with the section header

*Fragments
Then:
1. Input structure data for fragment 1 using the Structure Data module input form. See instructions for keywords and format.  (If one fragment is an atom it must be fragment 2.) (See notes at the end of this section for special instructions pertaining to fragment structure inputs).
2. Next, structure data for fragment 2 is input in exactly the same way, using the Structure Data module again. (See notes at the end of this section for special instructions pertaining to fragment structure inputs).
3. Now the energy of this channel must be input.  Use keyword ChanEnergyD, followed by the energy relative to the primary channel.  (This keyword must be omitted for channel 1, which is defined to be the reference state with an energy of 0.0.)
4. If considering a neutral reaction then the coefficient for an r^-6 interaction potential between neutral fragments may be input using the following keyword.  This value may be used here in any PST calculations.  If this property is to be used in a VRC calculation it will need to be input again in the Potential Input Module.
PotSixD followed by the coefficient (cm^-1 Å^6).
5. If considering an ionic reaction, the default is for fragment 1 to be the neutral fragment.  If it is desired to make fragment 2 the neutral, give the keyword NeutFragD here followed by the value 2.
6. Again, if considering an ionic reaction, then the polarizability (and dipole if any) of the neutral fragment should be input using the following keywords.  These values are used to obtain reference values for the orbiting collision rate and possibly in any PST calculations.  If these properties are to be used in a VRC calculation employing, for example, an ion-induced dipole and/or ion-dipole potential, then they will need to be input again in the Potential Input Module.  
PolarizP followed by one of the following three qualifier keywords, which specify the choice for the neutral polarizability assignment:
EstimatePQ  The polarizability is crudely estimated as Mass/10.
CalculatePQ  Approximate calculation based on atom additivity.
ReadDQ  Input the neutral polarizability value (Å^3).
Note that both EstimatePQ and CalculatePQ assume that the first fragment is the neutral fragment unless you specified otherwise in part 4 above. 
	The dipole moment for the neutral fragment may be input with the keyword
DipoleD
These quantities are used in obtaining reference Langevin and dipole corrected Langevin rate constants.  If they are not input then the polarizability for the neutral fragment will be taken to be equal to that for the collider as input in Sec. 5, and the dipole will be set to zero.
7. For the first channel one may input a correction factor for the partition function for the primary fragments using the keyword
	FragCorrD

Note 1 on coordinate system orientation:  It is required that both fragment structures be input in Cartesian coordinates, using a particular principal axis coordinate frame. If the user does not know the orientation of the "correct" principal axis system (or even any principal axis system), an arbitrary coordinate system should be used to input the atom positions, and the program should be run.  The only requirement for this arbitrary coordinate system is that the first atom corresponds to the "bonding" atom.  The program will calculate and print out the correct principal axis coordinates of the atoms, and will then stop so that the user can reenter the fragment’s atom coordinates in the correct orientation before rerunning the program.  Note that this will often require repositioning the Pivot Points defined in Section 10C for a VRC calculation, since these will have been input initially using the user's original (arbitrary) coordinate frame.
It may be useful to know the rules the program uses for selecting the axes of its principal coordinate frame, which are as follows:
(a) the first atom read in for each fragment is the "bonding" atom and this atom is always in the first octant (x, y, z >=0);
(b) for linear tops, the z-axis coincides with the molecular axis;
(c) for symmetric tops, the z-axis coincides with axis of highest symmetry, i.e., the axis whose moment is non-degenerate;
(d) for asymmetric tops, the z-axis has the most unique moment (relative difference of other moment pair is smallest).

Note 2 for input of linear components of nonlinear fragments:  If a nonlinear fragment is input and the general potential is to be used then the first three atoms input must not be linear.
Note 3 for use of anisotropic ion-neutral potentials:  Three anisotropic ion-neutral potential functions are given as options in the Potential Input Module below, namely anisotropic ion-polarizability, ion-dipole, and ion-quadrupole.  The orientation of these terms in the neutral fragment is determined by dummy atoms which must be input in the present Fragments section. If any of these three potentials are used, the user must add the corresponding dummy atom to the end of the atom coordinate list for the neutral fragment.  The anisotropic polarizability dummy atom (if any) should be first, followed by the dipole dummy atom (if any), followed by the quadrupole dummy atom (if any).  There are thus anywhere from zero to three dummy atoms at the end of the neutral fragment atom list, depending on the potential functions chosen, and (for this ion potential) to it is important to include only those dummy atoms which are required for the potential functions specified.  The dummy atoms must have mass 0.0, and may be located anywhere except at the center of mass. Dummy atoms must be included in the total atom count for the fragment.
Note 4 for use of diplaced charge:  For ion-neutral calculations, the charge of the ion is located at its center of mass by default.  As described in the Potential Input Module, the charge may be displaced to any desired location in the ionic fragment.  The position of the displaced charge is designated by a dummy atom of mass 0.0 placed at the end of the atom coordinate list for the ionic fragment.  This dummy atom should be added only if the DisplacedChargeP keyword is used in the Potential Input Module. Dummy atoms must be included in the total atom count for the fragment.

***** 11. Transition State *****
[This section is omitted only for a low-efficiency case designated by the LowEffP (or equivalently LowPresP) keyword above.]
Begin with the header
*TransState
In some cases one may wish to consider the effect of more than one well distinguished TS location for a given channel.  In this case, use the following procedure: (1) designate the form of the effective transition state treatment with one of the following five keywords; (2) complete the input for transition state 1 by going all the way through Section 10; (3) return to the point designated *SecondTS or *ThirdTS below, and give the input for transition state 2 or 3.
The default, if none of the following set of 5 keywords are provided, is to use a single transition state for the channel in question, which will be located at the position of minimum flux over the sampled reaction coordinates.  This will be the usual approach. In this case, skip the following set of 5 keywords, and go through the rest of Section 10 to the end.

Set of special keywords for multiple transition state input:
MINTSP denotes the minimum of the partition function evaluated for two separate TS locations
EFFTSP denotes the inverse of the sum of the inverses of the partition functions for two separate TS locations
USTTSP denotes the unified statistical theory expression of Miller (W. H. Miller, J. Chem. Phys. 65, 2216 (1976)) which requires the consideration of three separate TS locations (1/N = 1/N_1 + 1/N_2 – 1/N_3).
SUMTSP denotes the sum of the partition function evaluated for two separate TS locations
RATIOTSP denotes the product N_1*N_2 / N_3 and thus requires three separate TS evaluations.  This ratio may be useful, for example, in correcting PST calculations of ortho/para effects for the loss in phase space due to the conversion from free rotors to hindered rotors.  I.e., a PST calculation for N_1 could be multiplied by N_2 for a VTST calculation employing classical state counts for the full potential divided by N_3 for the potential which is implicit in PST.  For such a correction to be valid the PST potential should be more attractive than the full potential for each orientation.   

Proceed to the next paragraph to continue input. If multiple transition states are used, then you will return to the present paragraph to begin input for the second and third TS locations.  In this case you must at this point give the mandatory section header *SecondTS or *ThirdTS, specifying the transition state for which the data is to be entered.  Then give the TS input for the corresponding TS location, .by proceeding through the rest of Section 10 to the end.

You may give one of the optional keywords specifying the form of the conserved mode vibrational convolution: NoVibConvP, or VibAdiabaticP, or VibStatisticalP (Default is VibStatisticalP, which will usually be correct for canonical TST calculations. VibAdiabaticP will normally be specified for microcanonical TST calculations.  NoVibConvP may be used to obtain solely the transitional mode contribution).
Next you may specify two energy grids to be used in (i) the convolution over the conserved modes and the transitional modes and (ii) in the transitional mode number of states evaluations.  There are two choices for this specification.  The generally preferred specification employs the keyword StepScaleD follows by three numbers.  The first is a grid size (in cm^-1) for the evaluation of the conserved mode (fragment vibrations) state.  The second and third numbers describe the starting step size (in cm^-1) and number of steps for a variable step size grid of energies at which the transitional mode number of states will be evaluated [E(j) = ETSTP*{(j-1) + [(EMAX-N*ETSTP)/ETSTP]**{(j-1)/(N-1)}].  Alternatively, the keyword StepSizeD may be given in which case the transitional mode number of states is evaluated on a grid of evenly spaced energies.  This keyword is followed by two numerical values with the first again defining the step size (in cm^-1) for the conserved mode state density and the second now defining a constant step size (in cm^-1) for the transitional mode number of states.  The default is to employ the StepSizeD set of grids with default parameter values of (1.0 cm^-1, 50.0 cm^-1).  However, it is worth noting that StepScaleD will generally provide sufficient convergence with a greatly reduced grid size and correspondingly reduced CPU times.
In the interest of efficiency, one should generally use a transitional mode step size which is ≥ the vibrational step size.  The convolution integral is carried out using the vibrational step size, with interpolated values of the transitional mode number of states.  The program may revise both the default and input values to larger values if this would result in certain matrix dimensions being exceeded (cf. output in flow.out).  If desired you can bypass this by increasing the corresponding matrix dimensions in param.f and then recompiling the program.  These stepsizes are irrelevant for T resolved calculations (TResP).

Next is an optional keyword for the electronic state properties.
Give either keyword ElecAdiabaticP or ElecStatisticalP to choose the type of folding in of the electronic state contributions (Default is ElecStatisticalP).
Give the keyword NElecStatesD followed by the number of electronic states of the transition state.
Now the cycle over the electronic states begins
For the first electronic state give the keyword ElecStatesL, followed by the degeneracy and the energy for the first electronic state. 
For subsequent electronic states, precede this line with the keyword NextElecState
Now specify the algorithm to be employed in evaluating the TS partition function.
RRHOP ==> Rigid Rotor Harmonic Oscillator/tight transition state treatment; employing various simple assumptions for tunneling (i.e., no tunneling, Wigner tunneling or Eckart tunneling).
PSTP ==> phase space theory evaluation of  N_EJ^‡ or N_E^‡ or Q_T^‡ with a range of possibilities for the conserved-transitional mode convolutions.
VRCP ==> Variable Reaction Coordinate evaluation of  N_EJ^‡ or N_E^dagger or Q_T^‡ employing various conserved-transitional mode convolution procedures.
LaplaceP ==> read in modified Arrhenius form for k(T) and use inverse laplace transform techniques to generate N(E)
ReadNP  ==> read in N_EJ^‡ or N_E^‡ from nej.dat 
The input stream is different for each of these choices as described in the following subsections.

When the input for the current electronic state is complete, return to the point above designated "cycle over the electronic states," and repeat for the next state, if more than one electronic state is available to the fragments. 
Special input format for multiple transition states:
After finishing the specification of the first transition state with each of the fragment electronic states, you should return to the paragraph above specifying the section headers “*SecondTS or *ThirdTS” if input is needed for a second or third transition state.  
Now proceed to Section 11.

***** Section 11, subsection A *****
A.  Input for keyword RRHOP
Start by inputing the structure data for the transition state following the standard form in the Structure Data Module instructions.
The treatment of tunneling may be specified by giving one of the following keywords (default is NoTunnelP) 
NoTunnelP  No tunneling correction.
WigTunD  Wigner tunneling correction. The argument is first the frequency of the imaginary mode at the barrier top followed by a positive cutoff energy defining how far below the barrier top the consideration of tunneling will extend. 
EckTunD  Eckart tunneling correction.  The argument is first the frequency of the imaginary mode followed by the zero-point corrected energies of the barrier top relative to the reactant and product energies for the tunneling process. 
The tunneling corrections are currently only implemented for E/J and E resolved calculations
Next, for the first electronic state of a given channel, the energy of the transition state relative to the separated fragments for channel 1 is given by
TSEnergyD,  whose argument is the energy (cm^-1).
Finally, if desired give the correction factor for the transition state partition function with the optional keyword TSCorrD (Default is 1.0).

***** Section 11, subsection B *****
B.  Input for keyword PSTP

You may give the optional keyword JStepsD, followed by five numbers specifying the integration step sizes for the angular momentum variables of the two fragments, j_1, k_1, j_2, and k_2; and l, the orbital angular momentum (default is 1,1,1,1,1).
Now give the rotational symmetry factor for the transition state with keyword SigRotD.
A correction factor for the transition state  partition function may be supplied here if desired using the optional keyword TSCorrD  (Default is 1.0).
Finally, either a polarizability (for ion/neutral reactions) or a coefficient for the inverse sixth power potential (for neutral/neutral reactions) may be given.  
For ionic reactions the polarizability may be supplied with the PolarizP keyword followed by one of the qualifier keywords EstimatePQ, CalculatePQ, or ReadDQ, as specified in Sec. 8 for the fragments.  The default value for this parameter is to use that specified in Sec. 8 for the fragments.
For a neutral reaction the coefficient for the r^-6 potential (C6) may be input here with the keyword:  
PotsixD 
and the C6 potential parameter entered on the following line in cm^-1 Å^6 (The default is to use the value input for the fragments in Section 9.)

***** Section 11, subsection C *****
C.  Input for keyword VRCP
The reaction coordinate for the Variable Reaction Coordinate approach is defined as the distance between a Pivot Point in fragment 1 and a Pivot Point in fragment 2.  The Pivot Point locations are defined in each molecule relative to a reference point, which may be either the center of mass or a chosen ("bonding") atom.  With respect to the reference point in each fragment, the Pivot Point is located using spherical polar coordinates.  The reference frame for these spherical polar coordinates is the principal axis system used for the fragment structure in Section 10 above.  The θ = 0 direction is the z direction, and the φ = 0 direction is the x direction.
Give keyword BondRefP or CMRefP to choose whether the reference points in the two fragments are the bonding atoms or the centers of mass.  If BondRefP is used, the bonding atom in each fragment is defined as the first atom in the list of atom coordinates given in Section 10.
Next specify the ranges of reaction coordinates to consider.  In the most general case this will correspond to specifying R (the separation between the Pivot Points); the position of Pivot Point 1 (specified by Dist1, Theta1 and Phi1, which are the spherical polar coordinates for the vector from the reference point in fragment 1 to the pivot point in fragment 1); and the position of Pivot Point 2, specified by the corresponding quantities Dist2, Theta2, and Phi2 for fragment 2. Note that there are restrictions on the combinations of R, Dist1, and Dist2.  In particular, R must be greater than the sum of the magnitudes of the vectors from the center-of-mass to the Pivot Point for each of the fragments.  If the reference points are the center-of-mass then this restriction corresponds to R > Dist1 + Dist2.  If instead the reference points are the bonding atom the restriction may be written as R > Dist1’ + Dist2’; where Dist1’ is the distance from the center-of-mass of fragment 1 to the pivot point of fragment 1 which depends on Dist1, Theta1, and Phi1, as well as the vector from the center-of-mass of fragment 1 to its bonding atom. 

There are five choices (listed below as A-E) of input systems for specifying the ranges of these quantities.  In each instance distances should be given in Å, and angles in degrees.

A.  RRangeP with the qualifier NoPPSetsPQ.  The program assumes that the Pivot Points are fixed at the corresponding reference points and only the separation R between the reference points is to be varied.  The only input required is a set of R values.  The scheme for this case is identical to that for the Ranges input, employing the usual choice of one of the following three qualifiers:
StepDQ
ScaleDQ or
ListDQ (with the number of elements in the list as its argument) and followed by 
ValuesL
with the list of R values.

B.  RRangeP with the qualifier TotalPPSetsPQ.  A range of values for R is given.  Separately a series of pairs of Pivot Point positions (Dist1, Theta1, Phi1, Dist2, Theta2, Phi2) is given.  A calculation is done for each R value with each given pair of pivot points.
Input the range of R values in the standard format (using StepDQ, ScaleDQ or ListDQ as with option A above).  Then give the keyword 
NPPSetsD followed by the number of pivot point sets to be read in and then
PPSetsL followed by the specified number of sets of (Dist1, Theta1, Phi1, Dist2, Theta2, Phi2), either six numbers to a line, or in free format.
C.  RRangeP with the qualifier TwoPPSetsPQ.  A range of values for R is given.  Separately a series of pivot point 1 positions (Dist1, Theta1, Phi1) is given.  Finally a series of pivot point 2 positions (Dist2, Theta2, Phi2) is given.  A calculation is done for each R value with each possible combination of pivot points.
Input the range of R values in the standard format (using StepDQ, ScaleDQ or ListDQ as with option A above).  Then give the keyword 
NPP1SetsD followed by the number of pivot point 1 sets to be read in, and then
PP1SetsL followed by the specified number of sets of (Dist1, Theta1, Phi1), either three numbers to a line, or in free format.  Then give
NPP2SetsD followed by the number of pivot point 2 sets to be read in, and then
PP2SetsL followed by the specified number of sets of (Dist2, Theta2, Phi2), either three numbers to a line, or in free format.
D.  TotalSetsP.  With this option, a series of seven-number sets (R, Dist1, Theta1, Phi1, Dist2, Theta2, Phi2) is read in. Give the keyword
NSetsD whose argument is the number of sets of values to read, followed by
SetsL which is followed by the given number of sets of values (R, Dist1, Theta1, Phi1, Dist2, Theta2, Phi2), either seven numbers to a line, or free format.
E.  IndividualRangesP.  A range of values for each of the seven parameters is specified separately. A calculation is done for every possible combination of the given values of the seven parameters. The input scheme is again the standard Ranges input, using StepDQ, ScaleDQ or ListDQ for each of the parameters.  Cycle through the parameters in the following order using the keywords
RRangeP
Dist1RangeP
Theta1RangeP
Phi1RangeP
Dist2RangeP
Theta2RangeP
Phi2RangeP

Next is the optional keyword for input of R-dependent vibrational frequencies of the conserved modes. (Default is NoRDepVibP).
Give RDepVibP to select R dependence of the vibrational frequencies, or NoRDepVibP to use the same frequencies (those input for the fragments) at all R values.
If RDepVibP was specified then give the keyword
	NRDepModesD followed by the number of reaction coordinate R values at which vibrational frequencies will be input, and a multiplicative scaling factor for those frequencies.  The R values must be input in monotonically increasing or decreasing order.
Follow this with the keyword
	RDepModesL.  Cycle over each of the R points (with R either monontically increasing or decreasing), specifying R and the vibrational frequencies for each of the conserved modes for each of the R's.  The definition of R used here to label the sets of mode frequencies should correlate with the definition of the Rmep variable in the potential.  The default definition of R as specified in the standard potential is the "bonding" separation defined as the distance between the first atom of fragment 1 and the first atom of fragment 2.
Next is the keyword for input of R-dependent geometries of the fragments.
Give RDepGeoP to select R dependence of the fragment geometries, or NoRDepGeoP to use the fragment geometries (those input for the fragments) at all R values.  The only change in fragment geometries allowed is one where there is a change in the geometry from the asymptotic value but that change does not vary along the reaction path.  For example, CH3 is a planar radical in isolation but will splay away from the approach of a second reactant.  Using the RdepGeoP option, it is possible to specific a planar geometry for the CH3 asymptotic reactant and a splayed geometry with a fixed splay angle for the CH3 as a fragment along the reaction path.
If RDepGeoP was specified, then give the keyword PositionL followed by the masses and positions of the atoms for fragment 1 (unless fragment 1 is an atom).  On each line, give the mass (g mol^-1) and the x, y and z coordinates (Å) of each atom.  This describes the fixed fragment geometry that the fragment will assume in all variations along the reaction path.  Repeat, starting from the PositionL keyword, for fragment 2.  If fragment 1 or fragment 2 is an atom, then no input (including the PositionL keyword) should be supplied, as there is no internal geometry specification for an atom.  If neither fragment is an atom but only one fragment’s geometry changes from its asymptotic value, input is still required for both fragments, even though one fragment’s geometry is unchanged from what has already been entered elsewhere.
The keyword MinEPathD is given next followed by the data for the determination of the minimum energy path.  The minimum energy path is determined via a random sampling over relative orientations of the two fragments, for each definition of the reaction coordinate.  This requires an integer seed and a number of sampling points.  These two data values are given on the line following this keyword.

Now give keyword IntegrationParP followed by one of the following three keywords:
MonteCarloEJD. [only to be used for E or E/J resolved calculations (EResP or EJResP)]
It has a one-line set of arguments, consisting of three integers:
1.	The starting seed for the Monte Carlo integration 
2.	Total number of configurational sampling points for the integration 
3.	Number of momentum samplings per configuration 

Or

GridIntD  [only to be used for T resolved calculations (TResP)]
It has a two to six-line set of arguments.  Many of the test examples show values for the arguments listed below that are appropriate for about 1% relative accuracy in the canonical rate constant.  For each line, the arguments are:
Line 1 consists of three integration parameters specifying features of a customized Romberg quadrature evaluation of the transitional mode partition function:
1. eps: the targeted value of overall relative accuracy of the integration.  Romberg integration calculates the integral at successively finer (by a factor of two) integration grids.  It then takes the last Kextrp grids (where Kextrp is input to be described shortly) and extrapolates the integral to an infinitely fine grid.  The integration is considered converged when the inclusion of the integral for the most recent grid does not make a relative change in the extrapolated value of the integral by more than eps.  If a transitional mode partition function requires a multiple integration, the integrals are nested and each one dimensional integration will be converged according to eps.
2. ecut: cutoff value of ratio of the minimum to maximum value of integrand in a subsection of the integration range.  Because of the exponential character of the integrand, over much of its range, an integrand can be essentially approximated as zero.  For any grid, a contiguous region of grid points is identified which encloses the maximum integrand value and whose edge values are less than this maximum value by a relative amount of ecut.  If all points outside this region are ecut less in a relative sense than the maximum, the integrand outside this region is set to zero.  All evaluations of integrands on finer grids is only carried out within this region.  If a transitional mode partition function requires a multiple integration, the integrals are nested and each one dimensional integration will have its region of non-zero integrands determined according to ecut.

3. Kextrp :number of most recent grids whose integral value will be used in the extrapolation to the value on an infinitiely fine grid.  Each successive grid is a factor of two finer than the previous grid.  If a transitional mode partition function requires a multiple integration, the integrals are nested and each one dimensional integration will be converged according to Kextrp.
Line 2 is specific to integration over the theta1 (T1) orientaton angle of the first fragment. It consists of three parameters:
1. nsym: the number of symmetrically equivalent integration ranges for the given angle.  Each orientation angle has a full range of 0° - MaxValue.  The actual range of the integration is 0° - MaxValue/nsym.  Note that pivot points of low symmetry degrade the symmetry of the overall integration.  Therefore take care when making nsym not unity with low symmetry pivot points.
2. eterm: the ratio by which the current estimate of the transition mode partition function must exceed the lowest previous value in order to terminate the integration.  In canonical calculations, the lowest rate constant at a given temperature is the best rate constant.  If the program is searching over an input grid of pivot point locations and separation distances, it keeps track of the lowest rate constant for each temperature as it progress through the grid. It thus begins the current transitional mode partition function integration with knowledge of the maximum value the final integral can be in order for the current pivot point selection and separation distance to produce the lowest rate constant.  By design, Romberg integration produces for each level of grid refinement an estimation of the integral.  If that estimate exceeds the maximum value by the inverse of eterm, then the integration is terminated.  The coarsest grid for Romberg integration is just two points at the extremes of the range.  The integration will be terminated at that grid using min{eterm,0.1) in place of eterm.
3.	nspline: number of contiguous points beyond which a spline fit replaces explicit evaluation of the integrand.  In essentially all applications, the determination of the integrand is the most computer intensive aspect of the integration.  Since Romberg integration involves systematic refinements by a factor of the grid, it is possible to approximate the integrand on all finer grids by a spline fit to the integrand on the current grid.  Because of the dynamic range of the integrand, the spline is not fit to the integrand directly but to its natural log.  nspline controls the number of continuous points in the current grid at which this approximation is applied.  If nspline = 0, this approximation is never applied.  For technical reasons, non-zero values of nspline should be 1 plus a positive integer power of 2.  If a transitional mode partition function requires a multiple integration, the integrals are nested such that the theta1 and, if necessary, the theta2 integration are the highest level integrations.  Thus activating this approximation eliminates evaluating an integrand that itself may require numerical integrations. 

Line 3 thru 6 may or may not be required depending on the number of orientation angles. The following table indicates what orientation angle is required for each line:
| Case | Line 3 | Line 4 | Line 5 | Line 6
| Linear Top + Atom | absent | absent | absent | absent
| Nonlinear Top + Atom | phi1(P1) | absent | absent | absent
| Linear Top + Linear Top | theta2 (T2) | chi (CHI) | absent | absent
| Nonlinear Top + Linear Top | phi1(P1) | theta2(T2) | chi(CHI) | absent
| Nonlinear Top + Nonlinear Top | phi1(P1) | theta2(T2) | phi2(P2) | chi(CHI)
If required, each line has input similar to Line 2 only for a different orientation angle:
1. nsym: the number of symmetrically equivalent integration ranges for the given angle.  Each orientation angle has a full range of 0° - MaxValue.  The actual range of the integration is 0° - MaxValue/nsym.
2. eterm: the fraction that controls whether the current estimate of the one dimensional integral in this orientation angle needs to be improved upon.  The orientation angle integration covered here is an interior integral in a multiple integral and therefore must be repeated many times in the evaluation of the integrand of the outer integrals.  Only the largest values of the integrand of these outer integrals are important.  So if the current integration over this orientation angle is smaller than MinVal, a value determined from past evaluations during the determination of the same multiple integral, the integration is terminated at its current value.  If and only if the integration is not terminated by the MinVal test, MinVal for this orientation angle is set to the minimum of past values of MinVal or eterm times the final current value of the integrand.
3.	nspline: number of contiguous points beyond which a spline fit replaces explicit evaluation of the integrand.  ***Only present if the orientation angle is theta2.  For all other orientation angles, there is no nspline entry.***  In essentially all applications, the determination of the integrand is the most computer intensive aspect of the integration.  Since Romberg integration involves systematic refinements by a factor of the grid, it is possible to approximate the integrand on all finer grids by a spline fit to the integrand on the current grid.  Because of the dynamic range of the integrand, the spline is not fit to the integrand directly but to its natural log.  nspline controls the number of continuous points in the current grid at which this approximation is applied.  If nspline = 0, this approximation is never applied.  For technical reasons, non-zero values of nspline should be 1 plus a positive integer power of 2.  Theta2 integration occurs only when the transitional mode partition function requires a multiple integration.  The integrals are nested such that the theta2 integration has an integrand that always contains an interior integral. Thus activating this approximation eliminates evaluating an integrand that itself requires numerical integrations.

Or
MonteCarloTD  [only to be used for T resolved calculations (TResP)]
It has a six-line set of arguments:
Line 1 consists of four integer parameters specifying features of a Monte Carlo integral evaluation of the transitional mode partition function.
1.The type of Monte Carlo integration: 1 for crude; 2 for stratified (selected angles).
2. Random number seed: must be a negative integer.
3. The number of integration subunits. Cumulative results are printed at end of each subunit.
4. The number of (multidimensional) Monte Carlo points per subunit. The total number of points is the product of this number and number of subunits.
Lines 2-6: Same as for GridIntD except that each line consists of only one integer specifying number of symmetrically equivalent integration ranges for the given angle. The second and third entries included under GridIntD are not needed here but, if included, are ignored and will cause no harm.

Now give the rotational symmetry factor for the transition state with keyword SigRotD.
Then, if desired, give the correction factor for the transition state with the optional keyword TSCorrD (Default is 1.0).
Finally, the potential function must be specified.  Begin with the header *Potential.  Then input the potential parameters using the Potential Input Module given below.

***** Section 11, subsection D *****
D.	Input for keyword LaplaceP

The inverse Laplace transform procedure requires as input the parameters A_∞ (in cm^3 s^-1) n_∞, E_∞, and T_∞ for the modified Arrhenius form k = A_∞(T/T_∞)^n_∞ exp[-E_∞/(k_BT)].  These parameters should be specified on the line following the LaplaceP keyword.

Note:  Quantum evaluations of the fragment partition functions (QuanFrag1P) cannot be used in combination with an inverse Laplace transform evaluation for the first channel.   Also, the inverse Laplace transform procedure requires a constant grid spacing (i.e., do not use StepScaleD for the corresponding transition state). 

***** Section 11, subsection E *****
E.	Input for keyword ReadNP

The parameter ReadNP implies that the number of states will be read in directly.  

For the E resolved case this involves reading first the number of energies.  This is followed by a loop over the energies reading pairs of energies (in cm^-1) and transition state number of states (E, N_E^‡).  The program requires that the energies match those in the initial grid of total energies from the convolution section.   

For the E/J resolved case this involves reading first the number of energies E and the number of total angular momenta J.  This is followed by a nested loop over J and E (with E as the inner loop) reading E, J, N_EJ^‡.  Again, the program requires that the E and J read in match those in the initial grid of E and J from the convolution section.   

***** 12. Channel 2, 3, etc. *****
If there is more than one channel, you should return to Section 9, input the next channel section keyword, and repeat all the input of sections 9, 10 and 11 for the next channel.

***** End *****
At the end of all input (which will be after the completion of either the Potential Input Module or the TS data for the last channel), give the keyword *End.  

***** Structure Data Module *****
The following input is given in order each time the structure input module is executed.

Molecule type: Give one of the six molecule type keywords first:
AtomP
LinearP
OblateP  (A=B>C, where A,B,C are the rotational constants)
ProlateP  (A>B=C)
NearOblateP  (A~B>C, where A,B,C are the rotational constants)
NearProlateP  (A>B~C)

Electronic states: Next give information for the electronic states (For RRHO TS’s only one electronic state can be given and so the following two lines are skipped.  If you desire more than one electronic state for a RRHO TS simply include them as additional channels.):
NElecStatesD followed by the number of electronic states considered. After this must come the following keyword with the list of state degeneracies and energies:
ElecStatesL, with succeeding lines giving the degeneracy and energy (cm^-1) of each state relative to that for the ground electronic state of the particular channel or complex.  Use 0.0 for each ground electronic state. 
If the species is an atom, give the keyword AtomMassD, followed by the atomic mass, but omit this keyword if the species is a molecule.  
If the species is an atom, this is the end of the Structure Input module.

Vibrational modes:
First give the keyword NModesD, followed by two parameters, which are the number of vibrational modes of the molecule, and the correction factor by which the given mode frequencies are to be multiplied.
Some or all of the mode frequencies may be estimated automatically by the program according to the Standard Hydrocarbon pattern.  To do this, give the optional keyword StdModesD here, followed by the number of modes you wish to have the program estimate for you.  This number should not be less than 15.
Give the keyword ModesL, followed by the list of mode frequencies (cm^-1).  If some frequencies were estimated according to the preceding paragraph, input only the remaining mode frequencies here.  If all were estimated, omit this keyword.

Radiative Rates:
If the data is for the complex and the radiative rate keyword in Section 3 is CalcPQ, the radiative intensities must be assigned here using the keyword ModesIntensL.  There are two choices:
1.  The radiative intensities of the modes (km/mol) may be read in from an input list.
2.  The qualifier StdIntensPQ may be given after ModesIntensL, in which case the intensities will be estimated automatically from the standard hydrocarbon model.
If desired, a multiplicative correction factor for the radiative rates may be given with keyword KRCorrD.

Free Rotors:
Some of the modes may be treated as classical 1-d free rotors by giving the keyword FreeRotorD at this point in the input stream, followed by the number of 1-d free rotors on the next line. Then give the keyword FreeRotorL followed by the list of rotational constants (cm^-1) for each of the free rotors.  The corresponding modes should be removed from the list of vibrational modes (that is, the number given for NModesD above should be the number of harmonic oscillators, not including free or hindered rotors).

Hindered Rotors:
Some of the modes may be treated as classical 1-dimensional hindered rotors by giving the keyword HindRotorD at this point followed by the number of 1-d hindered rotors on the next line. The corresponding modes should be removed from the list of vibrational modes (that is, the number given for NModesD above should be the number of harmonic oscillators, not including free or hindered rotors).  Then cycle through each of the hindered rotors providing the following set of data:
Give the keyword HindTypeD followd by an integer which specifies the format of the potential for that hindered rotor.  Then give the keyword HindParL followed by the list of parameters for that format of potential (all in cm^-1 or cm^-1/rad^2). In each case the first parameter specifies the effective rotational constant (cm^-1) for the given torsional motion.
HindTypeD = 1.  V = V_0 + V_1 cos (nΔφ).  There are three parameters. The second parameter specifies the value of the potential at Δφ=π/n, and the third parameter specifies the value of n.
HindTypeD = 2.  V = V_0 + V_1 cos (Δφ) + V_2 cos (2Δφ). There are three parameters.  The second parameter specifies V(Δφ=180) and the third parameter specifies the force constant at Δφ=0.
HindTypeD = 3.  V = V_0 + V_1 cos (Δφ) + V_2 cos (2Δφ). There are three parameters.  The second parameter specifies V(Δφ=180) and the third parameter specifies V(Δφ=90).

HindTypeD = 4.  V = V_0 + V_1 cos (Δφ) + V_2 cos (2Δφ) + V_3 cos (3Δφ).  There are four parameters. The second parameter specifies V(Δφ=180), the third and fourth parameters specify the force constants at Δφ=0 and Δφ=180, respectively.
HindTypeD = 5.  V = V_0 + V_1 cos (Δφ) + V_2 cos (2Δφ) + V_3 cos (3Δφ).  There are four parameters. The second parameter specifies V(Δφ=180), the third parameter specifies the force constants at Δφ=0 and the fourth parameter specifies V(Δφ=90).respectively.
HindTypeD = 6.  V = V_0 + V_1 cos (Δφ) + V_2 cos (2Δφ) + V_3 cos (3Δφ) + V_4 cos (4Δφ).  There are five parameters.  The second parameter specifies V(Δφ=180), the third and fourth parameters specify the force constants at Δφ=0 and Δφ=180, respectively, and the fifth parameter specifies V(Δφ=90).
In each instance V(Δφ= 0) is taken to be the equilibrium point with 0 potential value.
Note that the relation between frequency and rotational constant is B = ω^2/(2 * force const.)  as long as all are in the same energy units (e.g., all in cm-1).
Rotational Symmetry Number: 
Next give the rotational symmetry number using the keyword SigRotD, followed by the value given as real number.
Ortho/para option (designed for systems involving H_2): One may force the fragment partition function, and also the corresponding PST TS evaluation in a phase-space-theory calculation, to include a sum over only the even or odd fragment rotational states.  At this point, give the keyword OrthoParaP followed by either EvenPQ or OddPQ.  For this to work properly one must include the appropriate QuanPFFrag1P or QuanPFFrag2P keyword in the Calculation Type section.  You will also probably need to revise your fragment symmetry numbers, which will generally now be unity.
Geometry:
Now the geometry is input.  Give one of the two keywords RotConstP or CartesianP to specify input of either the rotational constants or of the masses and positions of the individual atoms of the molecule.
A. If the keyword is RotConstP, the total mass is first input as the argument to the MassD keyword.  Then the rotational constant(s) are given as arguments to the RotConstD keyword. One value for linear (LinearP), three values for nonlinear (OblateP or ProlateP).
B. If the keyword is CartesianP, then the input is as follows:
Use keyword NAtomsD, whose argument is the number of atoms in the molecule.  If dummy atoms are used in the atom list, they should be included in this count.
Then use keyword PositionL followed by the masses and positions of the atoms.  On each line, give the mass (g mol^-1) and the x, y and z coordinates (Å) of the atom.
Linear species must lie along the z-axis.
Nonlinear species should not have the first three atoms as linear if the general potential is to be used.

***** Potential Input Module *****
	Each channel for which a VRC (VRCP) transition state calculation is to be performed needs potential  input, which is done as follows. 
	Certain standard potential forms are built in, and their parameters are read from the input file.  If desired, following input of any standard potential parameters, the user can call a custom potential input, which is read in from a separate input data file "pot.dat".  The user must put the necessary potential-reading routine computer code in a subroutine called potdatcust to read his/her potential parameters.  The user must also modify the potential-creation routine mcpot to process these custom potential parameters.
After the *Potential header is read in the appropriate "TransState" section, the following standard potential forms are available by specifying the keywords below, any combination of which may be used to construct the total potential.  There are some preliminary inputs which may be required:
1.  If you are using an ion-neutral potential, the center of charge is located by default at the center of mass of the ion.  It may be displaced, if desired, to an arbitrary position in the ionic fragment.  Designate this by giving the optional keyword DisplacedChargeP here.  The position of the center of charge will have been located in the Fragments input section using a dummy atom in the atom coordinate list of the ionic fragment.
2. If you are using one of the anisotropic potentials, and if one (or both) of the fragments is a non-linear polyatomic, first you must inform the program whether the "bonding" atom is terminal or central for fragment 1 (keyword TerminalAtom1P or CentralAtom1P) and/or fragment 2 (keyword TerminalAtom2P or CentralAtom2P). 

The following standard potential choices may be employed.  Note that for choices 2-7 the contributions are added together to form the final potential.

1. CustomP.  
	This keyword is set up to execute a separate user-written parameter-input subroutine called potdatcust, located in file potread.f.  The parameters for the custom potential terms should be read from a separate input file called pot.dat, unit 15. A sample form for the subroutine potdatcust has been included in the distribution files.  The user may use this form with the appropriate data in the file "pot.dat", or substitute his/her own subroutine using the name potdatcust. 
		The calculation of the potential from these parameters must be carried out in the subroutine mcpot, located in file potcalc.f.  The program actually uses the total potential given by variable vtot located in subroutine mcpot.  The user must add code to mcpot to process his custom potential parameters and add the additional potential energy (in atomic units) to the variable vtot.  The arrays nipotpar(nchmx), ipotpar(nchmx,nippmx), npotpar(nchmx) and potpar(nchmx,nppmx) are available (included in the common blocks in commonpot.f) for passing data from potdatcust to mcpot.
		A variety of examples of customized potentials are found in the examples section of the manual.

2. VStretchP
The stretching potential is defined with respect to the "bonding" separation defined as the distance between the first atom of fragment 1 and the first atom of fragment 2.  Specify one of the following three stretching potential forms:
MorseD, whose argument is the three Morse parameters D_e (cm^-1), beta (Å^-1) and R_e (Å).  (Positive D_e).
VarshniD, whose argument is the three Varshni potential parameters D_e (cm^-1), beta (Å^-2) and R_e (Å).  The Varshni potential is a potential form which may give a more faithful representation of the long-range part of the bond stretching curve than the Morse curve.  See S. J. Klippenstein et al., J. Chem. Phys. 89, 4761 (1988); Y. P. Varshni, Rev. Mod. Phys. 29, 664 (1957).
InterpD.  This form interpolates between a set of given potential points.  The argument is the number of points.  This is followed by PointsL giving a list of the potential points, given in the form of pairs of numbers consisting of distance (Å) and potential (cm^-1) (in ascending order of distances).

3. VBendP
This keyword is followed by one of two valence bending potential parameter formats (extrapolation or interpolation).  The bending potential is calculated for any given bond length based on the equilibrium bending angles and bending force constants input in this section using either a sinusoidal hindering (SinHindPQ) or a quadratic hindering (QuadraticPQ). (Angles in degrees, bending force constants in hartree/rad^2, bond lengths in Å.) See Appendix 1 for definitions.
ExtrapP,  extrapolation of bending force constants from the values at the equilibrium bonding distance r_e. Give on the same line either the qualifier keyword SinHindPQ or QuadraticPQ.  Then on the next lines provide the following information:
NAnglesD, the number of angular variables (call this n_a).
PeriodD, the period of the potential in fractional multiples of 2π.  E.g., 0.5 for a potential which has a period of π.  A default value of 1.0 is assumed.
AnglesL,  the set of equilibrium values of the angular variables (list of n_a values).
AlphasL,  the exponential extrapolation decay constants for each of the angular variables (list of n_a values).
ForcesL, the matrix of bending force constants at r_e.  (Square matrix of size n_a x n_a).
InterpP,  interpolation of bending force constants between the values given for a series of r values. Give on the same line either the qualifier keyword SinHindPQ or QuadraticPQ.  Then on the next lines provide the following information:
NAnglesD, the number of angular variables (n_a).
PeriodD, the period of the potential in fractional multiples of 2π.  E.g., 0.5 for a potential which has a period of π. A default value of 1.0 is assumed.
NPointsD, the number of r values for which the force constant information is going to be given (call this n_r).
PointsL, the list of bond length values for which information is to be given (n_r values).
AnglesL, the equilibrium angle values at each of the given r values (n_r lines having n_a values on each line).
ForcesL, the force constant matrix at each given r value (a series of n_r lower triangular matrices, each of which has dimension n_a).

4. AtomAtomP  
	Pairwise atom-atom interactions.  Choice of hard-sphere or Lennard-Jones pairwise potentials.  Many default values are built in.
HardSphereP
Diameters1L.  Diameters of fragment 1 atoms (Å).  Enter -1.0 for the default value for a given atom. (Default values are present for H, D, C, N, O, F, Si, P, S, Cl, Ge, As, Se, Br, Sb, I).  For dummy atoms, values must be entered: use a diameter of -10.0 for these.
Diameters2L.  Diameters of fragment 2 atoms (Å).  Enter -1.0 for default value.
LennardJonesP  
Parameters1L.  Epsilon (cm^-1) and sigma (Å) values for fragment 1 atoms, entered in pairs.  Enter a negative sigma value for each atom you wish to use default parameters for. (Default values are present for H, C, N, O, F, S, Cl, Br). For dummy atoms, values must be entered: for these enter 0.0, 0.0.
Parameters2L.  Epsilon (cm^-1) and sigma (Å) values for fragment 2 atoms, entered in pairs.  Enter a negative sigma value for each atom you wish to use default parameters for. For dummy atoms, values must be entered: for these enter 0.0, 0.0.

5. AnisotropicP
	A potential anisotropy form assuming a bonding potential which is cylindrically symmetric with respect to each fragment.  It is assumed that the stretching potential form, number 2 above, is used in conjunction with this potential form.  See Appendix 1 for definitions.  Specify the following keywords:
	PowerD, the exponent of the angular dependence function (normally 2).
It must be indicated whether the bonding site in each fragment is one-sided or two-sided.  Input one of the two keywords OneSided1P or TwoSided1P, and one of the two keywords OneSided2P or TwoSided2P (unless fragment 2 is an atom).
Give the two equilibrium bonding angles:
Theta1D, the equilibrium value of the fragment 1 bonding angle (degrees)
Theta2D, the equilibrium value of the fragment 2 angle

6. Ion-neutral interaction potential.  
	The potential may include any or all of three terms, namely ion-induced-dipole, ion-dipole, and ion-quadrupole potentials.  The polarizability may be isotropic or anisotropic.  The reference axes for the anisotropic potentials are determined by the direction from the center of mass to the appropriate dummy atoms in the neutral fragment input.  That is, the "parallel" direction for the anisotropic polarizability (assumed to be cylindrical with parallel and perpendicular components) is in the direction from the center of mass to the polarizability dummy atom;  the dipole (positive end) points in the direction from the center of mass to the dipole dummy atom (for positive ions; for negative ions the dummy atom corresponds to the negative end of the dipole); and the quadrupole (assumed to be cylindrical) is oriented in the direction from the center of mass to the quadrupole dummy atom.  If the polarizability is isotropic, no dummy atom is used for it.  To include any one or more of these terms, give the appropriate keyword group(s) from the following:
i)  IonIndDipP  Ion-induced-dipole potential term.  There are four choices for supplying the polarizability of the neutral fragment (give exactly one of these keywords):
EstimateP.  The polarizability is crudely estimated as Mass/10.
CalculateP,  Approximate calculation based on atom additivity.
PolzD.  Input the isotropic neutral polarizability value (Å^3).
AniPolzD.  Input two numbers which are the parallel and perpendicular components of the anisotropic polarizability (Å^3).  The parallel component is oriented along the reference axis described by the polarizability dummy atom.  The "parallel" direction will normally be assigned along the cylindrical axis, or the highest symmetry rotation axis, of a molecule having symmetry.
ii) IonDipoleD.  Ion-permanent-dipole potential. The argument is the permanent dipole of the neutral (Debye). The dipole is oriented (positive end) along the reference axis described by the dipole dummy atom.
iii) IonQuadD.  Ion-permanent-quadrupole potential.  The argument is the permanent quadrupole of the neutral (Debye- Å). The quadrupole is oriented (positive field) along the reference axis described by the quadrupole dummy atom.

7. Inverse sixth power potential
	For phase space calculation with neutrals. Use the keyword Inverse6D, whose argument is the coefficient of the r^-6 potential term (in units of  cm^-1A^6).

***** Standard Hydrocarbon Radiative Association *****

Giving the keyword StHydroP in section 2 invokes a self-contained module which provides a semiquantitative estimate of the radiative association rate constant.  No other calculations will be performed when the input file contains this keyword.
The ranges of temperature and binding energy to be considered are input using Section 3.  Then the rest of the input takes the special form given below:
After the Section 3 input, give the section header *StHydro
Then use the following keywords (any order):
DegreesD, the number of vibrational degrees of freedom of the association complex.
MassesD,  the mass of the ion and the mass of the neutral (g/mole).
PolzD,  the polarizability of the neutral (A^3).  Optional keyword: if it is omitted, an estimated value will be used.
If one of the reactants is an atom, a correction will be applied to the association rate estimation; this is specified by giving the keyword AtomicP, which is required if one reactant is monatomic.
After completing the required input, type *End.

Appendix 1
Definitions of coordinates for bending potentials in the built-in potential

Read this Appendix if you are using the built-in potential forms designated as Bending (Number 3, VBendP) or Anisotropic (Number 5, AnisotropicP) in the Potential Input Module.
The term General Potential used here refers to the set of built-in potential functions supplied with Variflex.  Parameters which specify the General Potential are read in by subroutine potdat located in file Potread.f.
The potential function generated by the General Potential is calculated as a function of relative orientation of the two fragments by subroutine mcpot, located in file Potcalc_gen.f.
Note (as described in the Potential Input Module instructions above) that Variflex may be run using solely the potential functions provided by the General Potential; or a user-defined potential function may be used; or the General Potential may be used and supplemented by additional user-supplied potential terms.
The angularly anisotropic potential functions (Potential forms 3 and 5) are calculated based on the deviation of the relative orientation of the two fragments from an equilibrium (or reference) orientation.  The input supplied by the user must specify the orientation of the fragments in this reference orientation, using the correct intermolecular coordinates.  This appendix describes these coordinates and their relation to the Variflex input.
The angular potentials are defined with reference to a bonding atom in each fragment.  (These atoms need not actually participate in bonding, but serve as reference positions.)  The bonding atom in fragment 1 is designated A_1, and it is the first atom in the atom list for fragment 1.  Similarly bonding atom B_1 is the first atom in the atom list of fragment 2.  The line R joining A_1 and B_1 is a convenient reference axis for keeping track of the relative orientations.  The next two atoms in the atom list for fragment 1 are A_2 and A_3, which are used to define angles.  Similarly B_2 and B_3 are the second and third atoms in the atom list for fragment 2.  The torsions are defined in the form Torsion(ABCD), using the usual Z-matrix convention, viewing from the first atom given, back along the bond to the second atom given, and so forth back along the remaining two atoms.
First define the bonding angles Theta1 and Theta2, which are the angles between R and reference axes in the molecules.  These are defined by 
Theta1 = Angle(B_1A_1A_2)
Theta2 = Angle(A_1B_1B_2)
These are the angles used in the AnisotropicP potential function input (Potential form 5), and their equilibrium values are to be specified in the input.
The more general potential of Potential form 3 requires the specification of a set of angular coordinates which depends on the specific case.  These are described individually below.
1. Fragment 1 linear, Fragment 2 atomic.  One coordinate.  NAngles = 1.  The angle is Theta1 = Angle(B_1A_1A_2).
2. Fragment 1 non-linear, Fragment 2 atomic.  Two coordinates. NAngles = 2.  First angle is Theta1 = Angle(B_1A_1A_2).  Second angle is the torsion Tau1 = Torsion(B_1A_1A_2A_3).
3. Fragment 1 linear, Fragment 2 linear.  Three coordinates. NAngles = 3.  First angle is Theta1 = Angle(B_1A_1A_2).  Second angle is Theta2 = Angle(A_1B_1B_2).  Third angle is the torsion Tau1 = Torsion(B_2B_1A_1A_2).
4. Fragment 1 non-linear, Fragment 2 linear.  Four coordinates.  Nangles = 4.  First angle is Theta1 = Angle(B_1A_1A_2).  Second angle is the torsion Tau2 = Torsion(B_1A_1A_2A_3).  Third angle is Theta2 = Angle(A_1B_1B_2).  Fourth angle is the torsion Tau1 = Torsion(B_2B_1A_1A_2).
5. Fragment 1 linear, Fragment 2 non-linear.  Four coordinates.  Nangles = 4.  First angle is Theta1 = Angle(B_1A_1A_2).  Second angle is Theta2 = Angle(A_1B_1B_2)..  Third angle is the torsion Tau1 = Torsion(B_2B_1A_1A_2).  Fourth angle is the torsion Tau3 = Torsion(B_3B_1B_2A_1).
6. Fragment 1 non-linear, Fragment 2 non-linear.  Five coordinates.  Nangles = 5.  First angle is Theta1 = Angle(B_1A_1A_2). Second angle is the torsion Tau2 = Torsion(B_1A_1A_2A_3).  Third angle is angle is Theta2 = Angle(A_1B_1B_2).  Fourth angle is the torsion Tau1 = Torsion(B_2B_1A_1A_2).  Fifth angle is the torsion Tau3 = Torsion(B_3B_1B_2A_1).

Appendix 2
Output Contained in the Auxiliary Output Files

In addition to the principal output file variflex.out, several other files contain useful output information.  flow.out, pf.out, and trate.out are always printed, erate.out and intgl.out, are always printed when an E or E/J resolved calculation is performed.  The files tspf.out and tspfr.out are only printed when the corresponding TspfDebugP and TspfRDebugP keywrods are given in the Title section.  (The names given here for these files are the generic names.  If the vrfx_run script file is used then these files will be renamed with ‘job_’ preceding them, except for the Variflex.out file which will be renamed ‘job’.out.) 

1. flow.out.
This file prints messages at many points during program execution to help track the progress of the calculation.  No timing information is provided in the standard distribution version, but the file Numeric.f contains a subroutine timer that includes several timing routines that have been used on various machines. The routine Unixtime works on some UNIX systems to give machine time information.  The routines date_and_time and cpu_time work with Fortran 90 and 95, respectively to give real elapsed time and CPU time, respectively.  The routine MCLOCK works on IBM machines.  To activate any one of these, remove the initial c character from all lines in the appropriate routine, and recompile.  The timing information is printed in flow.out as well as the main output.

2.  pf.out
	This file contains the canonical partition functions (at each temperature used) for the complex and for each of the fragments.  These are the same values as given in Variflex.out, which are collected together here for convenience.

3.  trate.out.
	This file gives a summary of the final rate constant values.  These are the same values as given in Variflex.out, which are collected together here for convenience.

4.  erate.out.  
	The first part of this output gives the density-of-states ρ(E,J) and the microcanonical unimolecular dissociation rate (for each channel) for the complex.  These values are listed for each value of E and J for which they are calculated.
The second part of the output gives the vibrational photon emission rate from the complex for each value of E and J.

5.  intgl.out.
	This file displays the accumulating values of the numerator of the thermal rate constant expression, which are accumulated during the numerical integration of the Boltzmann weighted contributions over all values of J and E.

6.  tspf.out.
	The first part of this output gives the number of states function N^‡(E,J,transitional) for the transitional degrees of freedom in the transition state.  Sometimes this information will be given for each electronic state and each channel.  Other times it will be given for the sum over the electronic states.  Also given in each case is the reaction coordinate (R) value found for the transition state, and the positions of the two pivot points (either defined, or optimized, depending on whether a variable reaction coordinate functional form is being used.)
	The second part of the output gives the total number-of-states function N^‡(E,J) obtained by convoluting N^‡(E,J,transitional) with the density-of-states function of the conserved modes.

7.  tspfr.out.
	This output gives the number of states function N^‡(E,J,transitional) as in tspf.out but for each different transition state dividing surface considered (i.e. each R, d1, d2 combination). 

Appendix 3
Convergence Parameters in E/J and E Resolved Calculations 

Obtaining converged results for E/J or E resolved calculations requires the consideration of a number of different convergence parameters.  The following convergence parameters must be considered: 
(I)	For the convergence of the final thermal integrals over E and/or J within the RRKM style analyses (cf Section 5 of this manual, and Eq. (12) of J.Chem.Phys. 104, 4502, 1996).  These parameters also play a related role in the master equation analyses.  The ERangeP parameters determine the minimum E, the number of E’s and the step size in E for the thermal convolution.  Similarly, the JrangeP parameters determine the minimum J, the number of J’s, and the step size in J for the thermal convolution.  Both of these sets of parameters are key parameters and one should generally check the convergence of the final results with respect to their variation.  For room temperature evaluations an energy step size of about 40 cm^-1 and a maximum E of 3000 cm^-1 will be in about the right range.  For RRKM analyses the minimum energy should be near 0, while for master equation analyses –1000 cm^-1 would be a typical minimum energy.  The range and step sizes for J are system dependent (i.e., due to variations in the rotational constants) but steps of 5 to 20, and maximum values of 100 to 400 are fairly typical.  The minimum J value should generally be near zero.  For convenience the convergence of the thermal integral is printed out in intgl.out for the high-pressure limit of the first and last temperatures considered.  For RRKM analyses the calculation time will generally be linearly proportional to the total range of energy and to the total number of J values considered.  For master equation analyses the matrix dimension will be linearly related to the total number of E values and the total number of J values.  Even beginning users must adjust these parameters to obtain convergence
(II)	For the estimate of the transitional mode contribution to the transition state number of states for a given energy (cf. MonteCarloEJD in Section 10C).  This contribution is evaluated via Monte Carlo integration and there are three integer parameters for this integration. The first integer simply determines the initial seed for the random numbers.  Any value for this number is appropriate.  The second integer is a key parameter which defines the number of configurations sampled in the random sampling. Numbers for this parameter will typically be in the range from 100 to 10000; a value of 1000 will generally provide convergence to about 10-15% which is often satisfactory.  A value of 100 is useful for obtaining first estimates of transition state locations and final results, but is not usually satisfactory for final results.  Convergence with only 100 points would be limited to about 40%.  A value of 10000 is useful for detailed examinations of the transition state, such as the dependence of its location on E or J. Final estimates can be converged to about 5% or better with a value of 10000.  The third parameter defines the number of momentum samplings for each configuration.   A value of unity is most common.  However, convergence rates for reaction coordinates that are different from the center-of-mass separation can sometimes be improved by choosing values between 2 and 5. For beginning users we recommend simply choosing these three parameters as 99, 2000, 2.

(III)	For the convolution of the transitional and conserved modes (cf. StepSizeD or StepScaleD in section 10 of this manual, and Eq. (13) of J.Chem.Phys. 104, 4502, 1996).  The vibrational step size determines the step size actually used in the evaluation of the convolution integral.  The speed of the calculation seldom depends on the value for this parameter so a small number such as 1.0 cm^-1 is typically used for this parameter.  The code may increase this number above the input value in order to keep certain matrices within the limitations imposed in the program.  Such increases are reported in flow.out.  If necessary, such increases can be shortcircuited via the increase of mxdme in param.f, and then recompiling the program.  The transitional mode contribution is first evaluated over a grid of energies whose spacing is determined by the second parameter, which is again a key parameter.  The total time for the calculation is generally linearly proportional to the number of energies in this grid.  Thus, it is generally wise to employ as large as spacing as is feasible.  We have found that 50 cm^-1 often provides an adequate spacing.  However, some situations, such as low temperatures, may require much smaller spacings. For beginning users we recommend simply using the defaults of 1.0 cm^-1, 50.0 cm^-1 for these two parameters.  For StepScaleD we recommend 1.0 cm^-1, 20.0 cm^-1, 20 as a reasonable set of starting parameters.

(IV)	For the estimate of the partition functions for the complex (cf. EVCStepD in Section 7).  This parameter determines the vibrational step size employed in the evaluations of the state density (and when necessary radiative rates) for the complex.  Such evaluations are generally not the slow part of the calculation and so a small value of 1.0 cm^-1 is generally appropriate.  Again, the code may internally reset this variable to avoid exceeding certain matrix dimensions and messages to this effect are printed in flow.out.  To bypass this resetting you will need to increase mxdmec in param.f and recompile the program. 

(V)	For the transition state location (cf. Section 10C). The location of the variational transition state must lie within the range of reaction coordinate values specified, and the spacing of these values must be close enough to give adequate accuracy for the transition state location.  If there is doubt about these questions, the transition state results given in the tspf.out output file should be surveyed.

Sample Procedure for Testing Convergence
First analyze the range of E and J that are required to obtain converged results.  Use defaults for parts (III) and (IV) above.  Using 1000 Monte Carlo configurational points set up a guess for the E and J ranges and step sizes.  Run the calculation.  Then increase the E range by about 50% and run again.  Check for any variation in the final result.  Repeat until the results are sufficiently converged.  (Note there may be some minor variation from one calculation to the next due to random variations in the Monte Carlo integrations.  One way to avoid worrying about these minor variations is to consider the ingtl.out file.  The variation with E of the terms in this file may be used to determine whether the convergence is complete.  Care should be taken to look at the integral for the highest temperature of interest.  Note that the range of E required for convergence to be complete in the intgl.out file may be larger than necessary if one is only interested in low pressures.)   Next do the same thing for the J range.  Now increase the number of Monte Carlo points by factors of 2 to 10 until you are satisfied with the variation in this parameter.  Then decrease the E step size by about 50% and again check for variation.  Similarly, decrease the J step size by about 50% and check for variation.  At this point one might test convergence with respect to the transitional mode step size (paragraph III above) by decreasing it by a factor of two; this would be especially important if one is interested in temperatures below 300 K or in low energies.  Finally, one should examine the tspf.out file to make sure that the optimized TS locations are never at the limits of the grid of reaction coordinate values considered.  If they are then the range should be increased. The variation of the optimized TS results with reaction coordinate step size(s) should also be tested.