A General Chemical Dynamics Computer Program
William L. Hase, Ronald J. Duchovic,a Xiche Hu,b Andrew Komornicki,c
Kieran F. Lim,d Da-hong Lu,e Gilles H. Peslherbe,f Kandadai N. Swamy,g
Scott R. Vande Linde, Antonio Varandas,h Haobin Wang, and Ralph J. Wolf,i
Department of Chemistry,
Wayne State University
Detroit, MI 48202
a. Current address: Department of Chemistry, Indiana University-Purdue University at Fort Wayne, Fort Wayne, IN 46805-1499
b. Current address: Department of Chemistry, University of Illinois, Urbana, IL 61801
c. Current address: Polyatomics Research Institute, 1101 San Antonio Road, Suit 420, Mountain View, CA 94043
d. Current address: Biological and Chemical Sciences, Deakin University, Geelong, VIC 3217, AUSTRALIA
e. Current address: Department of Chemistry, Fitchburg State College, Fitchburg, MA 01420
f. Current address: Department of Chemistry and Biochemistry, University of Colorado at Boulder, Boulder, CO 80309-0215
g. Current address: IBM Corporation, Kingston, New York 12401
h. Current address: Departamento de Quimica, Universidade de Coimbra, P-3049 Coimbra Cedex, Portugal
i. Current address: Building 773-42A, Scientific Computations Section, Westinghouse Savannah River Col, Aiken, SC 29802
TABLE OF CONTENTS
I. INTRODUCTION ..................................................................... 3
II. STATIONARY POINT SEARCH (NSELT=-3) .................................. 5
III. THE REACTION PATH (NSELT=-2) ............................................. 6
A. Properties of the Reaction Path ................................................. 6
B. Free Energy Maximum .......................................................... 7
IV. NORMAL MODE ANALYSIS (NSELT=-1) ..................................... 11
V. FIND THE EQUILIBRIUM GEOMETRY (NSELT=1) ........................ 13
VI. SELECTION OF INITIAL CONDITIONS (NSELT>2) ........................ 14
A. Polyatomic Rotational Energy .................................................. 14
B. Orthant Sampling ................................................................ 15
C. Microcanonical Normal Mode Sampling ...................................... 16
D. Fixed Normal Mode Sampling ................................................. 17
E. Local Mode Sampling ........................................................... 17
F. Thermal Sampling ............................................................... 18
G. Barrier Excitation ................................................................ 18
H. Diatomic Molecule ............................................................... 19
I. Orientation of Reactants ......................................................... 21
VII. TRAJECTORY RESULTS .......................................................... 23
VIII. POTENTIAL ENERGY FUNCTIONS ............................................ 25
IX. ENTERING THE DATA FILE ...................................................... 27
A. Potential Energy Surface Parameters .......................................... 27
B. Calculating the Trajectories and Selecting the Initial Conditions ........... 34
C. Parameters for Classifying Reaction Events .................................. 39
D. Printing Trajectory Results ...................................................... 40
E. Initial Coordinates and Momenta for NSELT≠ 2 ............................. 42
F. Parameters for Calculating the Reaction Path and Locating the
Canonical Variational Transition State ......................................... 42
X. FILES FOR WRITING OUTPUT .................................................. 44
XI. INSTALLING AND RUNNING THE COMPUTER PROGRAM ............. 45
A. Array Dimensions and Parameters ............................................. 45
B. List and Description of the Subroutines and Functions ...................... 46
C. Sample Input Datafiles ........................................................... 49
This computer program will perform seven different types of calculations, which are specified by the parameter NSELT:
NSELT = -3 Search for stationary points, given an initial guess.*
NSELT = -2 Determine the reaction path and canonical variational transition state properties. NSELT = -1 Perform a normal mode analysis at a geometry which is read in.
NSELT = 0 Calculate a trajectory from coordinates and momenta which are read in.
NSELT = 1 Determine the minimum energy geometry for a specified potential energy surface.
NSELT = 2 Calculate a trajectory for one or two reactants.
NSELT = 3 Calculate a trajectory starting at a saddlepoint.
For NSELT > 2 polyatomic initial conditions may be chosen for different types of excitation, which are specified by the parameter NACT:
NACT = 0 The reactants are atoms and/or diatoms.
NACT = 1 Choose polyatomic reactant initial conditions with orthant sampling.
(for NSELT=2 only).
NACT = 2 Choose polyatomic reactant initial conditions with microcanonical normal mode sampling.
NACT = 3 Choose polyatomic reactant initial conditions with fixed normal mode energies.
NACT = 4 Choose polyatomic reactant initial conditions with local mode sampling.
NACT = 5 Choose polyatomic reactant initial conditions with a Boltzmann distribution of normal mode energies.
The methodology of classical trajectory calculations for triatomic systems has been reviewed on several occasions.1-4 Descriptions of trajectory methods for polyatomics are limited,5 so a rather detailed summary of the techniques used in this program is given below. Integration of the classical equations of motion is standard in that combined fourth-order Runge-Kutta and sixth-order Adams-Moulton algorithms are used.1,6
* Not available in this version. Code under development.
1. D. L. Bunker, Meth. Comput. Phys. 10, 287 (1971).
2. R. N. Porter and L. M. Raff, in Dynamics of Molecular Collisions, Part B, edited by W. H. Miller, (Plenum, New York, 1976), p. 1.
3. M. D. Pattengill, in Atom-Molecule Collision Theory, edited by R. B. Bernstein (Plenum, New York, 1979), p. 359.
4. D. G. Truhlar and J. T. Muckerman, in Atom-Molecule Collision Theory, edited by R. B. Bernstein (Plenum, New York, 1979), p. 505.
5. W. L. Hase, in Aspects of the Kinetics and Dynamics of Surface Reactions, AIP Conference Proceedings No. 61, edited by U. Landman, AIP, New York, 1980), p. 109.
6. C. A. Parr, Ph.D. Dissertation Thesis, California Institute of Technology.
II. STATIONARY POINT SEARCH. NSELT = -3
This option is not available yet, and still under development. Calculations with this option would determine properties of stationary points (minima, transition states, or higher-order saddlepoints) on a potential energy surface by zero gradient search. Two algorithms that can be interfaced to Venus for root finding are the Newton-Raphson and Broyden methods.1 Users can contact the authors for more technical information.
1. W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes: the At of Scientific Computing, Second Edition (Cambridge University Press, New York, 1992).
III. THE REACTION PATH. NSELT = -2
A. Properties of the Reaction Path
The symbol R is used in this presentation to denote the reaction path. The position of the transition state along the reaction path is identified by R‡. The reaction path, independent of the coordinate system, is defined by a trajectory that moves from reactants towards products in infinitesimal steps with its kinetic energy removed after each step.1 For an association reaction without a saddlepoint the calculation of the reaction path is initialized with a large separation between the reactants. If the initial separation is sufficiently large the resulting path is independent of the starting separation.
In mass-weighted Cartesian coordinates, the reaction path trajectory traces out the steepest descent path and can be found by following the negative gradient vector of the 3N Cartesian coordinates from the initializing separation.2-4 The negative gradient vector is given by
Here i = 1,2,...N labels the atoms and g = x,y,z labels the Cartesian components, are the mass-weighted Cartesian coordinates, mi is the mass of atom i, V the potential energy surface, and c is a constant to normalize the 3N-dimensional vector to unity;
The partial derivatives are obtained analytically from the potential energy function. The 3N first order differential equations given by Eq. (II.1) which define the reaction path, are solved numerically; e.g. by a Runge-Kutta integration algorithm. A normal reaction path step-size DR for the numerical integrations is 0.00002 - 0.001 amu1/2-Å.
Properties found along the reaction path by solving Eqs. (II.1) and (II.2) are the classical potential energy VMEP(R) and Cartesian coordinates for each of the atoms. Internal coordinates such as bond lengths, valence angles, and torsion angles can be determined as a function of R from the atomic Cartesian coordinates. The Cartesian coordinates are also used to calculate the principal moments of inertia Ig(R), with g = x,y,z and g = x,y for nonlinear and linear systems, respectively.
In addition to the above properties, generalized normal coordinates may be defined for the vibrational motions orthogonal to the reaction path.2-4 Harmonic frequencies for these normal coordinates are found by diagonalizing the projected force constant matrix FP,
where F is a 3N x 3N matrix of the mass-weighted Cartesian force constants ,
and P the projector. The role of P is to project out of F the direction along the reaction path and the directions corresponding to infinitesimal rotations and translations. As a result of this projection, FP will have seven zero eigenvalues corresponding to infinitesimal rotations, translations, and motion along the reaction path. There will also be 3N-7 non-zero eigenvalues which give the harmonic frequencies ni(R) for the vibrations orthogonal to the reaction path.
In an actual calculation for an analytic potential energy surface, the ni(R) are determined during the numerical integration of Eq. (II.1).5,6 The Cartesian force constants and the projector are evaluated at intervals along the reaction path. The Cartesian force constants are calculated numerically from the analytic first-derivatives . Analytic expressions are used for determining components of the projector.2 The F and P matrices are then inserted into Eq. (II.3) to form FP, which is then diagonalized to obtain the 3N-7 harmonic frequencies. The vibrationally/rotationally adiabatic ground-state potential is found from the harmonic frequencies and the classical potential; i.e.
where s equals 3N-6 if every degree of freedom is treated quantum mechanically. Zero-point energy is not included in Eq. (II.5) for modes treated classically (see the following section).
B. Free Energy Maximum
The variational transition state theory rate constant is given by7-11
where kb is Boltzmann's constant, T is temperature, h is Planck's constant, s is the symmetry factor Q‡(T,R‡) is the partition function for the transition state, QR(T) is the partition function for the reactant(s), and is the difference in the reference energies between the transition state and the reactants. Conventional transition state theory locates the transition state at the saddlepoint. Canonical variational transition state theory, which is applicable to reactants with Boltzmann energy distributions, chooses the transition state at the point along the reaction coordinate R where the thermal rate constant is minimized,
In terms of the free energy of activation DG‡, the thermal rate constant k(T,R‡) is expressed by
where DG‡ is the difference in free energy between the transition state and the reactants. Thus, the maximum value of DG‡ will result in the minimum value of k(T,R).
Partition functions are assumed to be separable into contributions from translation, rotation, vibration, electronic and internal rotation. Rotational partition functions are determined by the classical approximation and vibrational partition functions are determined by the harmonic independent normal mode approximation. Partition functions for two-dimensional internal rotations are discussed below.
The free energy along the reaction path is calculated from the equation
The first three terms account for contributions from vibration, external rotation, and internal rotation. is the harmonic adiabatic ground state potential given by Eq. (II.5).
For reactions like Li+ + DME ® Li+ (DME) and Li+ + H2O ® Li+(H2O), rocking motions along the reaction path change from a free rotation of a neutral molecule about its center of mass at infinite separation to low frequency vibrations at the equilibrium separation [DME = (CH3)2O]. These rocking motions are treated as either quantum harmonic independent normal modes, or 2-dimensional classical hindered internal rotors.12
The 2-dimensional hindered internal rotor treatment allows these rocking motions to gradually change from a free rotation to a hindered rotation, and finally to a vibration at short separations. The Hamiltonian utilized for this motion is
where q has the limits 0,π and j has the limits 0,2π. The classical partition function13 for this Hamiltonian is
where dt = dpjdpqdjdq, and s is the symmetry number. If V(q,j) is zero, the partition function is that for free rotation and equals .
The reduced moments of inertia Iq and Ij for the internal rotation are calculated as described by Pitzer and Gwinn.14 As example, consider H + CH3 ® CH4, where
is either the Iq or Ij moment of inertia for CH3, and is the same for the H---CH3 reactive system. Both and are calculated from the reaction path Cartesian coordinates. The integral in Eq. (II.11) is evaluated by treating pq analytically and using ten point Gaussian quadrature to numerically integrate over pj, j and q.
1. K. Fukui, S. Kato, and H. Fujimoto, J. Am. Chem. Soc. 97, 1 (1975); S. Kato, H. Kato, and K. Fukui, ibid. 99, 684 (1977); H. F. Schaefer, Chem. Brit. 11, 227 (1974); K. Ishida, K. Morokuma, and A. Komarnicki, J. Chem. Phys. 66, 2153 (1979).
2. W. H. Miller, N. C. Handy, and J. E. Adams, J. Chem. Phys. 72, 99 (1980).
3. S. Kato and K. Morokuma, J. Chem. Phys. 73, 3900 (1980); K. Morokuma and S. Kato, in Potential Energy Surfaces and Dynamics Calculations, edited by D. G. Truhlar (Plenum, New York, 1981), p. 243.
4. A. D. Isaacson and D. G. Truhlar, J. Chem. Phys. 76, 1380 (1985).
5. W. L. Hase and R. J. Duchovic, J. Chem. Phys. 83, 3448 (1985).
6. S. R. Vande Linde, S. L. Mondro, and W. L. Hase, J. Chem. Phys. 86, 1348 (1987).
7. W. L. Hase, J. Chem. Phys. 64, 2442 (1976); W. L. Hase, Acc. Chem. Res. 16, 258 (1983).
8. M. Quack and J. Troe, Ber. Bunsenges. Phys. Chem. 81, 329 (1977).
9. D. G. Truhlar and B. C. Garrett Acc. Chem. Res. 13, 440 (1980).
10. D. G. Truhlar , W. L. Hase, and J. T. Hynes, J. Phys. Chem. 87, 2644 (1883).
11. D. G. Truhlar, A. D. Isaacson, and B. C. Garrett in Theory of Chemical Reaction Dynamics, edited by M. Baer (Chemical Rubber, Boca Raton, FL, 1985).
12. W. L. Hase, S. L. Mondro, R. J. Duchovic, and D. M. Hirst, J. Am. Chem. Soc. 109, 2926 (1987).
13. D. A. McQuarrie, Statistical Thermodynamics (Harper and Row, New York, 1973).
14. K. S. Pitzer and W. D. Gwinn, J. Chem. Phys. 10, 428 (1942).
IV. NORMAL MODE ANALYSIS. NSELT = -1
Normal mode coordinates and frequencies are computed for the specified potential energy surface by diagonalizing the mass-weighted Cartesian force constant matrix.1 The derivatives of the potential with respect to the Cartesian coordinates, ∂V/∂xi, are calculated analytically. Second derivatives and force constants are determined numerically. The Cartesian force constants are then mass weighted.
In mass-weighted Cartesian coordinates q, the kinetic and potential energy in matrix notation are
where ƒ is a 3N x 3N matrix of the mass-weighted Cartesian force constants for the N atoms. The matrix [L] diagonalizes [ƒ] to give the diagonal eigenvalue matrix [L],
The components of [L] are the vibrational frequency parameters .
The transformation between normal coordinates [Q] and Cartesian coordinates [q] is given by
[q] = [L] [Q] (III.4)
and the inverse
Thus, a column of [L], an eigenvector, is a representation of normal mode coordinate Qi in mass-weighted Cartesian coordinates. A more convenient representation is in Cartesian displacement coordinates x;
where [M] is a diagonal matrix whose elements are the atomic masses. From (III.5) and (III.6)
The computer program determines the frequencies ni in cm-1 and the Cartesian coordinate eigenvectors given by the columns of the matrix [Lx].
1. S. Califano, Vibrational States (Wiley, New York, 1976), p.23.
V. FIND THE EQUILIBRIUM GEOMETRY. NSELT = 1
To determine the geometry and energy at a minimum on a potential energy surface initial Cartesian coordinate and momentum vectors, x and px, are entered. A trajectory is then integrated for NIP time steps from the phase space point given by x and px, and the resulting momenta are then set to zero. This process is repeated until the trajectory is integrated for a total of NC time steps. Coordinates and energy are printed at each NIP time step interval, and convergence on the minimum energy geometry is determined by inspection. If the potential energy surface has multiple minima, each can be found by judicious choices for the initial x and px. For most cases the best choice for px is zero.
VI. SELECTION OF INITIAL CONDITIONS. NSELT > 2
Initial conditions for a polyatomic reactant are chosen by orthant sampling,1-3 microcanonical normal mode sampling,4 fixed normal mode energies,5-7 local mode sampling,14 or sampling a Boltzmann distribution. A diatomic reactant can be either a rotating harmonic or Morse oscillator.
For a trajectory started at a saddlepoint (NSELT=3), the initial conditions can be chosen by microcanonical normal mode sampling,4 fixed normal mode energies,5-7 or thermal sampling for all the modes but the reaction coordinate. The reaction coordinate energy is either fixed or chosen from a thermal distribution.
These sampling methods are described below. [Ri is a fresh random number uniformly distributed in the interval (0-1)].
A. Polyatomic Rotational Energy
Rotational energy and angular momentum for the polyatomic molecule is selected by assuming separability of vibrational and rotational motion. The options are specified by the parameters NROT and NLIN:
NROT = 0 Choose rotational energy from a thermal distribution by assuming a symmetric top,
Ix < Iy = Iz or Ix = Iy < Iz.8
NROT = 1 Rotational energy about each axis is RT/2.
NLIN = 0 Molecule is nonlinear.
NLIN = 1 Molecule is linear.
The program automatically rotates each reactant molecule to an axis frame defined by its principal moments of inertia. Either the x- or z-axis defines the symmetry axis for a symmetric top molecule. The x-axis defines the molecular axis for a linear molecule. For the NROT = 1 option the angular momentum about each axis is ±(IkT)1/2 and the rotational energy = 3kbT/2. If the molecule is linear Jx = 0 and = kbT.
For NROT = 0 the total angular momentum and its x or z component are sampled from the probability distributions,8
Jx(z) is sampled from P(Jx(z)) by the rejection method. J is sampled by the cumulative distribution function formula.
For linear molecules Jx = 0 and J is sampled from eq. (V.2). The components Jy and Jz are
The rotational energy is
B. Orthant Sampling. NACT = 1
Orthant sampling was developed for cases where the number of dimensions in the phase space is large and where the number of samplings is much less than the orthants in the phase space (orthants are in n-dimensional space what quadrants and octants are in two and three dimensions).1 It produces an approximate microcanonical ensemble, and is exact for a collection of harmonic oscillators.3,4 The steps are:
(1) The dimensions of the energy hypersurface at fixed H(x,px) - V0 = + are approximated. This is done by varying the coordinates and momenta to find their maximum and minimum values, , and (V0 is the minimum potential energy).
(2) To obtain the correct relationship between the average potential and kinetic energy, the are uniformly scaled by a parameter pscale (e.g., = for harmonic oscillators).2,3 For anharmonic and coupled systems does not equal .13
(3) The initial momenta and coordinates are then chosen from
where is the equilibrium position and the vi are components of a 6n-dimensional unit random vector (vi ranges between 0 and 1).1 Any center of mass velocity is then subtracted from the molecule.
(4) A spurious angular momentum arises from step 3. This angular momentum is computed and saved.
(5) The desired initial angular momentum Ji is added to the molecule by forming the vector
and adding the rotational velocity to each of the atoms, where
and I-1 is the inverse of the inertia tensor.
(6) The actual internal energy E is calculated by using the correct Hamiltonian and compared with the intended energy . If they do not agree to within 0.1 percent the coordinates and momenta are scaled according to
Any center of mass translational energy is then subtracted from the molecule and the procedure loops back to step 4.
C. Microcanonical Normal Mode Sampling. NACT = 2
This is an exact microcanonical sampling technique for harmonic oscillators.4 It is a complementary approximate method to orthant sampling for anharmonic oscillators.
The energies for n classical harmonic oscillators in a microcanonical ensemble are chosen from
where the sampling begins at i = 1 and . Cartesian coordinates and momenta are then chosen according to the following steps.
(1) The normal mode amplitudes are found from where .
(2) n different random numbers Ri are chosen to give each normal mode an initial random phase 2πRi. The initial normal mode coordinates and their time derivatives are then
(3) A transformation to Cartesian coordinates x and momenta px is then made via
where is a matrix of the equilibrium coordinates. Since normal modes are approximate for finite displacements a spurious angular momentum arises following this transformation.
Steps (4) - (6) are the same as those in section B above.
D. Fixed Normal Mode Energy Sampling. NACT = 3
For this option the normal mode quanta ni are read in and these values are used for all the initial conditions. Initial conditions are chosen according to steps (1) - (6) in section C above.
E. Local Mode Sampling. NACT=4
For this option14 the first step is to add an internal vibration-rotation energy to the reactant as described in the previous Section D. The next step is to add energy to a particular bond (i.e. local mode) so that the total energy of the bond given by
is equal to the energy of the quantum level n of a Morse oscillator; i.e.
In choosing the analytic potential energy function for the reactant, a Morse potential must be chosen for the local mode. The atoms for this local mode must be numbered 1 and 2, and one of these two atoms must be a terminal atom.
Cartesian coordinates and momenta for the reactant are chosen according to the following steps.
(1) Using E(n), Eq. (V.13), values for the local mode bond length r0 and momentum are chosen following the prescription of Porter, Raff, and Miller.10
(2) Calculate the current energy Er, Eq. (V.12), in the local mode. (Initially, adding E0 puts some energy in the local mode). Compare Er with E(n). If they do not agree to within 0.5 percent, continue with the following steps.
(3) The current bond length and momentum are r and pr. Add (r0-r) and ( - pr), respectively, to the bond length and momentum.
(4) Subtract any center of mass velocity on the reactant after step 3. Any spurious angular momentum which may have been added to the reactant is subtracted as described in step (5) of Section B above. The procedure loops back to step (2).
F. Thermal Sampling. NACT=5
The normal mode vibrational energies are sampled from a Boltzmann distribution at temperature Tvib, in which the quantum number ni in the ith mode is sampled using the probability distribution
where ni is the vibrational frequency of the ith mode.
G. Barrier Excitation. NSELT=3
With the reaction coordinate held fixed, the energy of the reaction coordinate is either a constant value or selected from a thermal distribution. In the latter case, the reaction coordinate momentum is sampled using the cumulative distribution function11
where R is a random number and Trc is the temperature chosen for the reaction coordinate. NACT of 2, 3, and 5 may be used to excite the vibrational modes at the barrier. NROT may be either 0 or 1.
H. Diatomic Molecule
The diatomic molecule is treated as a rotating oscillator, for which the Hamiltonian is
H = p2/2µ + V(r) + J2/2µr2 (V.16)
The oscillator internal energy Enj is calculated semiclassically (see Section VII for details) as a function of vibrational and rotational quantum numbers n and J, with an iterative (fixed-point) procedure. The linear momentum p for the rotating oscillator is chosen from
The following rejection method is used to choose the interatomic separation r. It is based upon the probability for an internuclear distance r
The steps are as follows:
(1) Determine the classical turning points r- and r+ (done in the semiclassical quantization procedure.
(2) Choose a value of r from r = r- + (r+-r- )R1; where R1 is a random number.
(3) Compute |p| and accept if P(r)/P0(r) is greater than a second random number R2. If not, repeat step 2. P0(r) is the relative probability for the most probable r.
The diatomic molecule lies along the x-axis and the angular momentum components Jy and Jz are determined in accordance with Eq. (V.3). They are added to the molecule by vectorially adding the angular velocity to each of the atoms.
I. Orientation of Reactants
If there is only one reactant, the above sampling is sufficient. However, for a collision between two reactants the following initial conditions are also needed.7
The impact parameter b for the collision is either kept constant or chosen randomly by
b = bmaxR1/2 (V.19)
where R is a random number. bmax must be chosen by tests.11 The Cartesian coordinates x and velocities for each reactant are randomly oriented by rotation through Euler's angles within the reactant's space-fixed center-of-mass Cartesian coordinate frame:12
The q, j, and c are chosen randomly according to
where R1, R2 and R3 are three different random numbers. A different set of random numbers is used for each reactant.
J. Relative Translational Energy (Erel)
The relative translational energy Erel can either be fixed or chosen from the Boltzmann distribution15 at temperature Trel
If using this distribution, the rate constant versus temperature Trel should be calculated via
where s(Erel) is the reactive cross section at Erel, m is the reduced mass of the two reactants, bmax is the maximum impact parameter, Nt is the total trajectory number, and Nr is the reactive trajectory number.
Finally, the centers of mass for the two reactants are separated by the distance s and a relative translational energy Erel added with the constraint that the total system center-of-mass remain at rest.
1. D. L. Bunker and W. L. Hase J. Chem. Phys. 59, 4621 (1973).
2. W. L. Hase, R. J. Wolf, and C. S. Sloane, J. Chem. Phys. 71, 2911 (1979).
3. R. J. Wolf and W. L. Hase, J. Chem. Phys. 72, 316 (1980).
4. W. L. Hase and D. G. Buckowski, Chem. Phys. Lett. 74, 284 (1980).
5. S. Chapman and D. L. Bunker, J. Chem. Phys. 62, 2890 (1975).
6. C. S. Sloane and W. L. Hase, J. Chem. Phys. 66, 1523 (1977).
7. W. L. Hase, D. M. Ludlow, R. J. Wolf and T. Schlick, J. Phys. Chem. 85, 958 (1981).
8. D. L. Bunker and E. A. Goring-Simpson, Faraday Discuss. Chem. Soc. 55, 93 (1973).
9. I. N. Levine, Molecular Spectroscopy (Wiley, New York, 1975), p.155.
10. R. N. Porter, L. M. Raff and W. H. Miller, J. Chem. Phys. 63, 2214 (1975).
11 D. L. Bunker, Methods. Comput. Phys. 10, 287 (1971).
12. E. B. Wilson, Jr., J. C. Decius, and P. C. Cross, Molecular Vibrations (McGraw Hill, New York, 1955), p. 285.
13. G. Nyman, K. Rynefors, and Leif Holmlid, J. Chem. Phys. 88, 3571 (1988).
14. D.-h. Lu and W. L. Hase, J. Phys. Chem. 92, 3217 (1988).
15. R. N. Porter and L M. Raff, in Dynamics of Molecular Collision, Part B, edited by W. H. Miller (Plenum, New York, 1976), p.1.
VII. TRAJECTORY RESULTS
All cartesian coordinates and momenta, and specified internal coordinates may be listed after each integration step. For a reaction A + B ® AB ® C + D the following results are listed at the conclusion of the trajectory.
1. The internal energy E, the vibrational and rotational energy Ev and Er (E = Ev + Er), the standard deviation in the rotational energy, and the total rotational angular momentum j are given for each product C and D. The rotational energy is computed by averaging over the longest vibrational period tvib of the product
is found from the trajectory by the summation
where n equals tvib divided by the trajectory integration stepsize. The standard deviation in the Er(i) is also computed. The rotational energy is found from
Vibrational and rotational quantum numbers are determined for a diatomic product molecules. The rotational quantum number J is found from the diatomic’s rotational angular momentum, i.e.
The vibrational quantum number n is found by Einstein-Brillouin-Keller (EBK)1 semiclassical quantization of the action integral
2. Relative properties between C and D are given. The relative translational energy is
Also listed is the relative translational energy plus the centrifugal potential
is the orbital angular momentum. At large separations L = µbvr .
3. Angles between the following pairs of vectors are computed: (Lf,jc),(Lf,jD), (jC,jD), (vi,vf), (Li,Lf), (Li,jC), (Li,jD), (jA,Li), (jA,Lf), (jA,jC), (jA,jD), (jA,jB), (jB,Li), (jB,jC) and (jB,jD).
1. M. C. Gutzwiller, Chaos in Classical and Quantum Mechanics (Springer, New York, 1990).
VIII. POTENTIAL ENERGY FUNCTIONS
The computer program allows the user to "build" an arbitrary potential energy surface from the following types of potential energy functions.
Type 1: Harmonic stretch, V = 1/2 ƒr (r-ro)2
Type 2: Morse stretch, V = D[1 - exp(-bDr)]2, b = c1 +c2 Dr + c3Dr2 + c4Dr3, Dr = r-ro
Type 3: Harmonic bend, V = 1/2 ƒq(q-qo)2, ƒq = S(rj) S(rk) where rj and rk define the bend. The S(r) attenuation terms are unity if Dr < 0. Otherwise, S(r) = exp[-C Dr2].
Type 4: Harmonic alpha bend (wag), V = 1/2 ƒa (a-π)2
Type 5: Generalized Lennard-Jones, V = a/rn + b/rm + c/r
Type 6: 2-fold torsion (double bond), V = Vo (1-cos2t)/2 = Vo sin2t
3-fold torsion (methyl group), V = Vo cos2(3t/2)
Type 7: Generalized exponential repulsion and attraction, V = a exp(-br) + c/rn
Type 8: Ghost pair interaction, V = q1/r12 + q2(1/r1 + 1/rln + 1/r2i + 1/r2j)
G. C. Lie and E. Clementi, J. Chem. Phys. 62, 2195 (1975)
Type 9: Tetrahedral center
R. J. Duchovic, W. L. Hase, and H. B. Schlegel, J. Phys. Chem. 88, 1339 (1984). Revised according to W. L. Hase, S. L. Mondro, R. J. Duchovic, D. M. Hirst, J. Am. Chem. Soc. 109, 2916 (1987)] and X. Hu and W. L. Hase, J. Chem. Phys. 95, 8073 (1991).
Type 10: Non-diagonal stretch-stretch interaction, , . The S(r) terms have the same functional form as in Type 3 above.
Type 11: Non-diagonal stretch-bend interaction,
V = ƒrq (rij - r) (qk - q), ƒrq = S(rij) S(rkm) S(rlm). The S(r) terms have the same functional form as in Type 3 above. S(rkm) and S(rlm) have the same exponential constant C as in Type 3. For S(rij) the constant C may be different. If ij = km or lm then S(rij) is set to unity.
Type 12: Non-diagonal bend-bend interaction,
, . The S(r) terms are defined in Type 3 above and have the same exponential constant C as in Type 3.
Type 13: Torsion for dihedral angle,
Type 14: Axilrod-Teller three-body potential, B.M. Axilrod and E. Teller, J. Chem. Phys. 11, 299 (1943).
Type 15: SN2 potential energy functions.
S. R. Vande Linde and W. L. Hase, J. Phys. Chem. 94, 2778 (1990) and H. Wang, L. Zhu and W. L. Hase, J. Phys. Chem. 98, 1608 (1994).
Type 16: Rydberg potential, where .
Type 17: Hartree-Fock Dispersion potential, R. Ahlrichs, R. Penco and G. Scoles, Chem. Phys. 19, 119 (1977).
and the Hartree-Fock repulsion is given by .
Type 18: Extended LEPS Potential for triatomic system, P. J. Kuntz, E. M. Nemeth, J. C. Polanyi, S. D. Rosner and C. E. Young, J. Chem. Phys. 44, 1168 (1966).
si is the Sato parameter. The 12 adjustable parameters in this potential are: 1Di, 1bi, 1ri0, and si (i = 1,2,3).
Type 19: Generalized LEPS Potential for triatomic system, P. J. Kuntz, in Dynamics of Molecular Collisions, Part B, edited by W. H. Miller (Plenum, New York, 1976) p.53.
This is similar to the LEPS function of Type 18
There are 18 adjustable parameters 1Di, 3Di, 1bi, 3bi, 1ri0, 3ri0 (i = 1,2,3). If one sets 1bi = 3bi, 1ri0 = 3ri0, and 3Di = 0.5[(1-si)/(1+si)] 1Di (si is the Sato parameter in Type 18), one gets an expression equivalent to Type 18.
Type 20: Double Many-Body Expansion (DMBE) IV potential energy surface, parameterized for
HO2 (so far); see M.R. Pastrana, L.A.M. Quintales, J. Brandao and A.J.C. Varandas, J. Phys. Chem. 94, 8073 (1990). A generalized version of this potential energy function and the EHFACE2 models [A.J.C. Varandas and J.D. Silva, J. Chem. Soc. Faraday Trans. 88, 941 (1992)] used for describing some two-body interactions in DBME IV will be available as standard options in future versions of Venus .
IX. ENTERING THE DATA FILE
Free format is used in the computer program. Time and date will be printed in the output file at the beginning of the execution time (this feature is implemented in the computer program for systems which support the Unix™ operating system).
Enter two lines of title.
A. Input Group One. Potential Energy Surface Parameters.
0. NATOMS (number of atoms); the following give the number of each type of potential energy function, NST (harmonic stretches), NM (Morse stretches), NB (harmonic bends), NA (harmonic wags), NLJ (generalized Lennard-Jones potentials), NTAU (torsions), NEXP (generalized exponential repulsion and attraction), NGHOST (ghost pairs) NTET (tetrahedral centers), NVRR (non-diagonal stretch-stretch interactions), NVRT (non-diagonal stretch-bend interactions), NVTT (non-diagonal bend-bend interactions), NANG (dihedral angles), NAXT (Axilrod-Teller), NSN2 (index number for the SN2 potential), NRYD (Rydberg potential), NHFD (Hartree-Fock Dispersion), NLEPSA (LEPS potential, Type A), NLEPSB (LEPS potential, Type B), NDMBE (DMBE potential)
For items (1-20) data is not entered if the number of a surface type (NST, NM, etc.) is zero.
1. N1J(I),N1K(I),RSZ(I), and FS(I) for I = 1, NST. This is information for the harmonic stretches. N1J is the J-atom, N1K is the K-atom, RSZ is the equilibrium bond length in Å and FS is the force constant in mdyn/Å. N1K must be greater than N1J.
2. N2J(I), N2K(I), RMZ(I), B(I), and D(I) for I = 1, NM. This is information for the Morse stretches. N2J is the J-atom, N2K is the K-atom, RMZ is the equilibrium bond length in Å, B is the exponential constant in Å-1, and D is the bond strength in kcal/mol. N2K must be greater than N2J. If B(I) is negative, the Morse beta term is represented by a cubic polynomial in r. Coefficients CM1(I), CM2(I), CM3(I), and CM 4(I) must then be entered.
3. N3J(I), N3K(I), N3M(I), THETAZ(I), FBZ(I), RJZ, CJ, RKZ, and CK for I = 1, NB. This is information for the harmonic bends. N3M is the middle atom. N3K must be greater than N3J. THETAZ is the equilibrium angle in degrees and FBZ is the force constant in mdyn-Å/rad2. The RZ (Å) and C(Å-2) terms are the equilibrium bond lengths and scale parameters for the J and K bond-stretching attenuation. If C<0 the attenuation term is not included.
4. N4J(I), N4K(I), N4M(I), N4N(I) and FA(I) for I = 1, NA. This is information for the harmonic wags. See the diagram below for the definition of the J, K, M and N atoms. FA is in the force constant in mdyn-Å/rad2.
5. N5J(I), N5K(I), NREP(I), MREP(I), LREP(I), ALJ(I), BLJ(I) and CLJ(I) for I = 1, NLJ. This is information for the generalized Lennard-Jones potentials. N5J and N5K are the atoms with N5K greater than N5J. NREP, MREP, and LREP are the powers in the potential energy expression. If a power equals zero the appropriate term is skipped in the calculations. ALJ, BLJ, and CLJ are in kcal/mol-Åpower.
6. N6I(I), N6J(I), N6K(I), N6L(I), N6M(I), N6N(I) and VZTAU(I) for I = 1, NTAU. This is information for the torsions. See the diagram below for the definition of the I, J, K, L, M and N atoms. VZTAU is the torsional barrier in units of kcal/mol. The bond vectors I®L and M®N define one plane. The bond vectors I®L and J®K define another plane:
The potential is defined in terms of the dihedral or torsional angle, t, which is the angle subtended by the two planes.
If VZTAU is set less than zero, then a symmetric 3-fold torsion is used, with the torsional barrier equal to the absolute value of VZTAU.
Note that two bond vectors may share an atom. Thus an ethyl H1-C2-C3-H4 grouping can be expressed by
Hydrogen (1) atom J
Carbon (2) atoms I and K
Carbon (3) atoms L and N
Hydrogen (4) atom M
Note also that if vectors J®K and M®N are parallel, t = 0 degrees. If J®K and M®N are anti-parallel, t = 180 degrees.
Thus, the above example (for a 3-fold torsion, i.e. VZTAU < 0) gives the staggered configuration as the equilibrium geometry for ethane. The following example gives the eclipsed configuration as the equilibrium geometry:
Hydrogen (1) atom K
Carbon (2) atoms I and J
Carbon (3) atoms L and N
Hydrogen (4) atom M
7. N7J(I), N7K(I), NPOW(I), AEX(I), BEX(I) and CEX(I) for I = 1, NEXP. This is the information for the generalized exponential repulsion and attraction potential functions. N7J is the J-atom and N7K is the K-atom, where N7K is greater than N7J. NPOW is the power for the bond length and BEX is the exponential constant in Å-1. AEX is in kcal/mole and CEX is in kcal/mol-Ånpow. If npow = 0 the c/rn term is skipped. If AEX < 0 the exponential term is skipped.
8. N8I(I), N8J(I), N8K(I), N8L(I), N8M(I), N8N(I), GCl(I), GQ1(I), and GQ2(I) for I = 1, NGHOST. This is the information for the ghost pair interaction. The atoms are defined by
1 and 2 are locations for fictitious charges, and are located on the angle bisector GC1 angstroms from the central atom. GQ1 and GQ2 are in kcal-Å/mol. N8J > N8I and N8N > MBL.
9. Information for the tetrahedral centers. Read in the following stacked input for I = 1,NTET.
N9I(I), N9J(I), N9K(I), N9L(I), and N9M(I). These are the atoms defining the tetrahedral center. The M-atom is at the tetrahedral center. The remaining four atoms, I, J, K, and L form bonds with M.
FT0 (I,J), J = 1,6. These are the quadratic force constants for the six valence angles at the tetrahedral center in units of mdyn-Å/rad2. Force constants for the angles are entered in the following order: IMJ, IMK, IML, JMK, JML, KML, angle.
FT2(I,J), J=1,6. Values for the valence angle quadratic force constants when one of the bonds not defining the angle dissociates.
GT0(I,J), J=1,6. These are the cubic force constants in units of mdyn-Å/rad3 for the six valence angles at the tetrahedral center.
GT2(I,J), J=1,6. Values for the valence angle cubic force constants when one of the bonds not defining the angle dissociates.
HT0(I,J), J=1,6. These are the quartic force constants in units of mdyn-Å/rad4 for six valence angles at the tetrahedral center.
HT2(I,J), J=1,6. Values for the valence angle quartic force constants when one of the bonds not defining the angle dissociates.
THT(I,J), J=1,6. Equilibrium values for the angles, in degrees, at the tetrahedral center.
THT1(I,J), J=1,6. Asymptotic equilibrium values for the valence angles, when one of the bonds defining the angle dissociates.
THT2(I,J), J=1,6. Asymptotic equilibrium values for the valence angles, when one of the bonds not defining the angle dissociates.
FD1(I,J), J=1,4. Out of plane quadratic force constants for the four-atom polyatomic product. J=1,2,3,and 4 are for dissociation of the I,J,K, and L atom, respectively.
HD1(I,J), J=1,4. Out of plane quartic force constants for the four-atom polyatomic products. J=1,2,3 and 4 are for dissociation of the I,J,K, and L atom, respectively.
GNO(I,J), J=1,5. Values for the five non-diagonal cubic force constants.
RO(I,J), J=1,4. Equilibrium values for the MI, MJ, MK, and ML bonds, respectively.
10. N10I(I), N10J(I), N10K(I), N10L(I), FKRR(I), RIJZ(I), CIJ(I), RKLZ(I), and CKL(I) for I=1,NVRR. This is information for the non-diagonal stretch-stretch interaction. Atoms N10I and N10J define one bond, and N10K and N10L the other; N10I>N10J and N10K>N10L. FKRR is the quadratic force constant in mdyn/Å. RIJZ and RKLZ are the equilibrium bond lengths in Å. The CIJ and CKL are scale parameters for the attenuation switching functions. If C < 0 the attenuation term is not included.
11. N11I(I), N11J(I), N11B(I), FKRT(I), R11Z(I), and CRT(I) for I=1,NVRT. This is information for the non-diagonal stretch-bend interaction. Atoms N11I and N11J define the stretch, with N11I>N11J. N11B is the index for the bend, from Input Type 3 above. The switching functions, which depend on the two bonds that define the angle, have the same C scale parameters as in Input Type 3. The switching function that depends on rij has the scale parameter CRT. If CRT<0 the switching function is fixed at unity. FKRTZ is the quadratic force constant in mdyn/rad. R11Z is the equilibrium bond length for the rij.
12. N12B(I), N12BB(I), NTT(I), FKTTZ(I) for I=1, NVTT. This is information for the bend-bend interaction. The index for the two valence bends are N12B and N12BB. If NTT = -1 the bend-bend interaction is not attenuated as any of the bonds defining this interaction are stretched. If NTT = 0 the bend-bend interaction is attenuated in the same way as are the harmonic bends in the Type 3 potential above. FKTTZ is the quadratic force constant in mdyn-Å/rad2.
13. Information for the dihedral angles. Read in the following stacked input for I=1, NANG. N13I(I), N13J(I), N13K(I) and N13L(I) for I=1, NANG. These are the indices for atoms defining the dihedral angle I. The atoms are defined in the figures below.
NDH(I). This is the order for the dihedral angle.
FDH(I,J) and GDH(I,J) for J=1, NDH(I). FDH (kcal/mol) and GDH (degrees) are the potential energy and phase angle for the Jth term in the potential energy function for the dihedral angle.
14. N14I(I), N14J(I), N14K(I), ZAXT(I) for I = 1, NAXT. This is information for the Axilrod-Teller interaction. N14I, N14J and N14K are the atom indices. ZAXT is the coefficient in the Axilrod-Teller potential function, in kcal/mol-Å9.
15. SN2 potential functions for the Cl- + CH3Cl and Cl- + CH3Br reactions. No information needs to be added at this point. The choice of the SN2 potential function depends on the NSN2 index, as follows:
16. N16J(I), N16K(I), RYDZ(I), DRYD(I), ARYDZ(I) for I = 1, NRYD. This is information for the Rydberg potential. N16J and N16K are the atom indices. RYDZ is the re equilibrium distance in Å, DRYD is the D depth of the potential in kcal/mol and ARYD is the dimensionless a exponential parameter.
17. N17J(I), N17K(I), RHFD(I), AHFD(I), BHFD(I), C6HFD(I), C8HFD(I), C10HFD(I) for I = 1, NHFD. This is information for the Hartree-Fock dispersion potential. N17J and N17K are the atom indices. RHFD is the rm equilibrium distance in Å, AHFD and BHFD are the Hartree-Fock repulsion energy A and a coefficients in kcal/mol and Å, respectively. C6HFD, C8HFD and C10HFD are the parameters of the dispersion expansion energy in kcal/mol-Å6, kcal/mol-Å8 and kcal/mol-Å10, respectively.
If NHFD is -1 then read in N17J, N17K, RHFD, AHFD, BHFD, C6HFD, C8HFD, C10HFD. All atoms with labels from N17J to N17K will have two-body interactions with the same listed parameters (convenient for clusters).
If NHFD is -2 then read in J, K, RHFD, AHFD, BHFD, C6HFD, C8HFD, C10HFD, N17. Atom N17 will have two-body interactions with atoms J to K, with the listed parameters.
If NHFD is -3 then first read one line corresponding to NHFD=-1 and then one line corresponding to NHFD=-2.
18. N18J1(I), N18K1(I), RLZ1(I), BL1(I), DL1(I), DELTA1(I); N18J2(I), N18K2(I), RLZ2(I), BL2(I), DL2(I), DELTA2(I); N18J3(I), N18K3(I), RLZ3(I), BL3(I), DL3(I), DELTA3(I); for I = 1, NLEPSA. This is information for the LEPS potential, type A. Similar to Morse function, N18J is the J-atom, N18K is the K-atom, RLZ is the equilibrium bond length in Å, BL is the exponential constant in Å-1, DL is the bond strength in kcal/mol, and DELTA is the Sato parameter. N18K must be greater than N18J. Since LEPS potential is three-body potential, there are 3 pairs: (J1, K1), (J2, K2), AND (J3, K3).
19. N19J1(I), N19K1(I), RLZS1(I), BLS1(I), DLS1(I), RLZT1(I), BLT1(I), DLT1(I); N19J2(I), N19K2(I), RLZS2(I), BLS2(I), DLS2(I), RLZT2(I), BLT2(I), DLT2(I); N19J3(I), N19K3(I), RLZS3(I), BLS3(I), DLS3(I), RLZT3(I), BLT3(I), DLT3(I); for I = 1, NLEPSB. This is information for the LEPS potential, type B. Similar to Morse function, N19J is the J-atom, N19K is the K-atom. RLZS, BLS, and DLS are parameters for groud-state Morse functions: RLZS is the equilibrium bond length in Å, BLS is the exponential constant in Å-1, and DLS is the bond strength in kcal/mol. RLZT, BLT, and DLT are parameters for excited-state “anti-Morse” functions: RLZT in Å, BLT in Å-1, and DLT in kcal/mol. N19K must be greater than N19J. Since LEPS potential is three-body potential, there are 3 pairs: (J1, K1), (J2, K2), AND (J3, K3).
20 No additional data needs to be entered for the DMBE IV potential energy surface of HO2. If NDMBE is a nonzero integer in input line 1, the H and O atoms will be identified by their masses.
21. VZERO. A parameter, in kcal/mol, for arbitrarily setting the potential energy zero.
B. Input Group Two. Calculating the Trajectories and Selecting the Initial Conditions.
1. NSELT, NCHKP. NSELT is the option for the type of calculation and is described in Sections I-V. NCHKP is an option for (re)starting calculations using checkpoint file information. If NCHKP≠0, calculations are restarted from a checkpoint file (any type of Venus calculation produces a checkpoint file in Unit 50):
NCHKP=1 Calculations are restarted where left, in case of a computer crash, or under-estimation of the cycle number in reaction path following, number of trajectories, etc. (the code follows input updates);
NCHKP=-1 New calculations are started with cartesian coordinates and momenta read from the checkpoint file. For example, one can perform a sequence of calculations with this option: minimum energy search, (NSELT=1) refinement by reaction path following (NSELT=2) normal mode analysis, (NSELT=-1) where coordinates and momenta can be passed through the checkpoint file. This option is only available for NSELT<2;
NCHKP=-2 The random number generator information from the previous calculation is read in from the checkpoint file (option available for NSELT=2 only).
2. W(I) for I = 1, NATOMS. Masses of the atoms in amu's.
In the following, the parameters to be entered are defined according to the value of NSELT.
Parameters for NSELT = -3, -2 and -1.
3. HINC and NPTS. HINC is the cartesian coordinate displacement interval for calculating force constants numerically. HINC must be chosen by trial, but for many cases the best value of HINC is approximately 0.001Å. NPTS is the number of displacements about the minimum. Choices are 1 and 2; with 2 the optimum choice.
Parameters for NSELT = -3 (Not available in this version)
4. NT, NS, NIP, SCALE. NT is the number of calculations to be performed (and the number of initial guess structures to be read in). NS is the maximum number of iterations during the energy gradient root search. NIP is the number of cycles between intermediate printouts. SCALE is a scaling factor for the convergence criteria. If SCALE is set to zero, the convergence criterion will be the computer limiting precision.
If NSELT = -3 or -1 continue reading data at D.
Parameters for NSELT = -2.
5. NATOMA(1), NLINA. NATOMA(1) is the number of atoms for reactant A. NLINA=0 if the molecule is nonlinear and =1 if the molecule is linear.
NATOMA(1) ≠ 0.
NATOMB(1), NLINB. Same as above, but for reactant B.
The reactant A atoms are labeled 1, 2, ..., NATOMA, and the reactant B atoms are labeled NATOMA+1, NATOMA+2, ... NATOMA+NATOMB. The reactants define NPATH = 1.
Continue reading data at D for NSELT = -2
Parameters for NSELT = 0, 1, 2 and 3.
6. NT (the number of trajectories to be calculated), NS (the total number of integration steps per trajectory), NIP (the number of steps between intermediate printout), and NCROT (the number of cycles over which the rotational energy in the products is averaged for analyzing final results). NIP must be greater than six. NCROT cannot be greater than 500.
For NSELT = 1, where a minimum in the potential is found, the kinetic energy is (re)set to zero after every NIP cycles.
7. TIME. TIME is the numerical integration stepsize in units of 10-14 sec.
If NSELT=0 continue reading data at C.
If NSELT=1 continue reading data at D.
If NSELT = 2 continue reading data at 9.
8. NBAR, EBAR or NBAR, TBAR. NBAR is an option for the barrier excitation. If NBAR=1, the reaction coordinate energy is fixed and must be entered as EBAR. If NBAR=2, the reaction coordinate energy is chosen from a thermal distribution, at fixed temperature, which has to be entered as TBAR. NATOMB(1) must equal zero for barrier excitation.
9. NACTA, NACTB, ISEED. NACTA and NACTB are options for the initial conditions of the trajectories (NSELT=2) for reactants A and B, respectively. NACT is described in Sections I and VI.
ISEED is a random number integer which initializes the selection of a chain of random numbers that are used in generating the initial coordinates and momenta. ISEED must be an integer between 0 and 231-1 and the last integer in ISEED must be odd. Larger values containing odd numbers are better.
If both NACTA and NACTB are 0 or 1, continue reading data at 11.
10. HINC and NPTS. HINC is the Cartesian coordinate displacement interval for calculating force constants numerically. HINC must be chosen by trial, but for many cases the best value of HINC is approximately 0.001Å. NPTS is the number of displacements about the minimum. Choices are 1 and 2; with 2 the optimum choice.
11. NATOMA(1), NLINA. NATOMA(1) is the number of atoms for reactant A. NLINA=0 if the molecule is nonlinear and =1 if the molecule is linear.
NATOMA(1) ≠ 0.
12. QZA(1,J), J=1,K. Equilibrium Cartesian coordinates for reactant A, with K=3* NATOMA(1).
If NATOMA(1) = 1 continue entering data at 20.
If NATOMA(1)>2 continue entering data at 14.
13. N, J. N is the diatomic number of vibrational quanta and J is the diatomic number of rotational quanta.
Continue entering data at 20.
14. Parameters for NACTA = 1 or 2.
ENMTA, PSCALA. ENMTA is the vibrational energy in kcal/mol. PSCALA is the momenta scale factor for orthant sampling, and is not used for NACT=2.
15. Parameters for NACTA = 3 or 4.
ANQA (I), I=1,J. Quantum numbers for reactant A's normal modes of vibration. J is the number of normal modes, sorted according to increasing frequency.
16. Parameter for NACTA = 5.
TVIB. Vibrational temperature for Boltzmann distribution.
17. NROTA, TROTA. IF NROTA = 0, A will have a thermal rotational energy distribution. IF NROTA = 1 the rotational energy about each axis of A will be RT/2. TROT is the temperature in Kelvin. The thermal sampling assumes a symmetric top, with Ixx @ Iyy ≠ Izz.
If NACT ≠ 4 continue entering data at 20. NATOMB(1) must equal zero if NACT = 4.
18. NEXM, NLEV. NEXM is the index for the Morse oscillator local mode to be excited (N2J for the Morse oscillator is the light atom). NLEV is the local mode vibrational level.
19. MPLOT, NPLOT, NLM, MNTR. These four quantities are needed to determine mode energies and mode populations. MPLOT = 1 if average mode energies every NPLOT cycles are to be printed on unit 15. If not, MPLOT = 0. Enter zero for NLM. This is a parameter for an option, to be added, in which energies in local modes different than the one initially excited are monitored. MNTR is the index for the normal mode, whose population is monitored. Average energies are determined every NPLOT cycles for each normal mode. Three or more trajectories must be run for averages.
20. NATOMB(1), NLINB. Parameters for reactant B as for reactant A in (11) above.
If NATOMB(1) = 0 go to C.
If not, enter parameters for B as in steps (12) - (17) above.
21. NABJ(1), NABK(1), RMAX(1), RBAR(1), DELH(1). These are reactant A and B internuclear parameters. They are defined in sections C.2 and C.3 below, and are also parameters for the first reaction path. If there is only one reactant NABJ(1) is set equal to NABK(1).
22. NREL, SEREL, S. IF NREL = 0, the initial A + B relative translational energy will be selected from a Boltzmann distribution at temperature SEREL (K). IF NREL = 1 the relative translational energy is fixed at SEREL (kcal/mol). S is the initial A + B separation in Å.
23. NOB, BMAX. If NOB = 0 the impact parameter is sampled randomly between 0 and BMAX. If NOB = 1 the impact parameter equals BMAX. BMAX is in Å.
C. Input Group Three. Parameters for Classifying Reaction Events.
Two fragments, A and B, are assumed to be formed by the reaction. NPATHS is the number of different reaction paths allowed in addition to that for the reactants.
If NPATHS = 0, go to D.
Items 2 through 7 are entered consecutively for I = 2, NPATHS + 1.
2. RMAX(I), RBAR(I), NATOMA(I), NATOMB(I), DELH(I). RMAX is the distance between A and B fragments in Å at which the trajectory is halted. RBAR is an intermediate A-B distance at which trajectory results are printed. It could correspond to a saddlepoint (barrier) distance. NATOMA and NATOMB are the number of atoms in the A and B fragments. DELH (kcal/mol) is the difference in the potential energy at the minimum and for A and B at infinite separation.
3. NABJ(I), NABK(I). Indices for two atoms J and K. The distance between the atoms defines RBAR and RMAX. K > J.
4. LA(I,J) for J = 1, NATOMA(I). An array of indices for the atoms of fragment A.
5. QZA(I,J) for J = 1, 3* NATOMA(I). Equilibrium coordinates for the atoms of fragment A.
6. LB(I,J) for J = 1, NATOMB(I). An array of indices for the atoms of fragment B.
7. QZB(I,J) for J = 1, 3*NATOMB(I). Equilibrium coordinates for the atoms of fragment B.
D. Input Group Four. Printing Trajectory Results.
1. NFQP and NFCOOR. NFQP is an option for listing coordinates and momenta. NFQP = 1 if they are to be listed. Otherwise NFQP = 0. NFCOOR is an option for saving coordinates in Unit 8 to prepare animation. NFCOOR = 0 or 1 as above for NFQP.
2. NFR and NUMR. NFR is the option for listing interparticle distances and NUMR is the number of distances to be listed. NFR = 1 if distances are to be listed, otherwise NFR = 0
If NFR = 0 do not include item 3.
3. JR(I) and KR(I) for I = 1, NUMR. JR and KR are the J- and K- atoms with KR > JR.
4. NFB and NUMB. NFB is the option for listing theta angles and NUMB is the number of theta angles to be listed. NFB = 1 if theta angles are to be listed, otherwise, NFB = 0.
If NFB = 0 do not include item 5.
5. IB(I) for I = 1, 3*NUMB, which are the atoms which define the theta angles. The central atom is listed third.
6. NFA and NUMA. NFA is the option for listing the alpha angles and NUMA is the number of alpha angles to be listed. NFA = 1 if alpha angles are to be listed, otherwise, NFA = 0.
If NFA = 0 do not include item 7.
7. IA(I) for I = 1, NUMA, which are the indices for the alpha angles to be listed.
8. NFTAU and NUMTAU. NFTAU is the option for listing the torsion angles and NUMTAU is the number of torsion angles to be listed. NFTAU = 1 if torsion angles are to be listed, otherwise, NFTAU = 0.
If NFTAU = 0 do not include item 9.
9. ITAU(I) for I = 1, NUMTAU, which are the indices for the torsion angles to be listed.
10. NFTET and NUMTET. NFTET is the option for listing the valence angles at the tetrahedral center. NFTET = 1 if the angles are to be listed, otherwise, NFTET = 0.
If NFTET = 0 do not include item 11.
11. ITET(I) for I = 1, NUMTET, which are the indices for the tetrahedral valence angles to be listed.
12. NFDH and NUMDH. NFDH is the option for listing the dihedral angles and NFDH is the number of dihedral angles to be listed. NFDH = 1 if torsion angles are to be listed, otherwise, NFDH = 0.
If NFDH = 0 do not include item 13.
13. IDH(I) for I = 1, NUMDH, which are the indices for the dihedral angles to be listed.
If NSELT = 2, stop entering data at this point.
E. Input Group Five. Initial Coordinates and Momenta for NSELT≠ 2.
1. Enter Q(I) for I = 1,3*NATOMS followed by P(I) for I = 1,3*NATOMS. The Q(I) and P(I) arrays are for the coordinates and momenta, respectively. Q(I) and P(I) can be stacked for NSELT ≠ -2. Only one Q(I) and P(I) array are entered for NSELT = -2. In using the NSELT = - 2 option to calculate the reaction path and free energy maximum, the axis system for the coordinates must be that defined by the principal moments of inertia, if there are 2-dimensional hindered rotors, i.e. N2DR ≠ 0 (see below). In addition if N2DR ≠ 0, the reaction path must lie along the y-axis.
If NSELT ≠ -2, stop entering data at this point.
F. Input Group Six. Parameters for Calculating the Reaction Path and Locating the Canonical Variational Transition State; NSELT = -2
1. NS, NIP, DS. NS is the number of integration steps along the reaction path. Free energy is calculated every NIP steps. DS is the integration stepsize in units of amu1/2‑Å.
2. WM(I), I=1, 3*NATOMS-6. These are the product vibrational frequencies in cm-1. The product is the polyatomic formed by association.
3. NTEMP. Number of temperatures for which free energy will be calculated along the reaction path.
4. TEMP(I) for I = 1, NTEMP. Values for the temperatures in K.
5. N1DR, N2DR. N1DR is the number of 1-dimensional free rotations along the reaction path. N2DR is the number of 2-dimensional hindered rotations along the reaction path. N2DR may not exceed 2. If N2DR is 1, fragment A is the polyatomic and fragment B is the atom.
If N1DR=0, continue entering data at 8.
6. AI1D(I) for I = 1, NIDR. Reduced moments of inertia in amu-Å2 for the 1-dimensional free rotations.
7. SYMM(I) for I = 1,N1DR. Symmetry numbers for the 1-dimensional free rotations.
If N2DR = 0, stop entering data
8. SYMA. Symmetry number for the 2-dimensional hindered rotation of the A reactant.
9. SYMB. Same as above, but for the 2-dimensional hindered rotation of the B reactant.
X. FILES FOR WRITING OUTPUT
File 6: general listing.
File 7: used in normal mode analysis.
File 8: save coordinates for making computer animations.
File 9: relative properties for products A and B.
File 10: properties for product A.
File 11: properties for product B.
File 12, 13, 14: populations of the n, n-1 and n+1 levels of an excited local mode versus time.
File 15: average energies of normal modes versus time.
File 16: population of an unexcited local mode versus.
File 50: checkpoint file.
XI. INSTALLING AND RUNNING THE COMPUTER PROGRAM
The standard Venus96 package - for Unix™ systems - comes as a compressed archive file called venus.tar.Z. The archive has been created using the standard Unix™ tar utility and it has been compressed using the standard Unix™ compress utility. It includes the Fortran 77 source code of the Venus96 computer program, venus96.f, a collection of input datafiles in a subdirectory called datafiles, a job script file which compiles the program and runs the different test inputs, and a PostScript version of this Venus96 manual. One should be able to compile and run on any workstation, but this code is primarily developed for Unix™ supported systems. The job script file has to be modified according to the user's system requirements.
It is safest to compile VENUS96 with the static memory allocation option. On the IBM RISC/6000 and SGI platforms it is “-qsave” and “-static”, respectively. Less memory will be required if the static option is not used. However, users should test the accuracy of the code, compiled in this manner, by comparing its output with the output obtained from the code compiled with the static option.
A vectorized version of Venus96 - called Vecf - is also available from the authors for use on vector capable workstations and Supercomputers. The latter code uses the similar input files as the scalar Venus96. For reference, see X. Hu, W.L. Hase and T. Pirraglia, J. Comp. Chem. 12, 1014 (1991).
The trajectories can be displayed, using the Sybyl® molecular modeling software. A C-language computer program utility is available from the authors to convert the output file in Unit 8 to the Sybyl® standards. Sybyl® is a registered trademark of Tripos Associates Inc.
The - adjustable - array dimensions and parameters are reviewed in the following Section. A list and description of the subroutines and functions is then given for the programmer's use. A list and description of the input datafiles is provided next and a hard copy of some of them is given in the last Section.
A. Array Dimensions and Parameters
Two parameters governing the array dimensions can be adjusted in VENUS:
ND1 maximum number of atoms for the calculations (default = 100).
NDP maximum number of reaction paths to be considered. See Input Group Three (default = 10).
The arrays for the potential energy functions can be adjusted by appropriate global substitutions. However, they are set to the following default dimensions:
1. Harmonic stretches 100
2. Morse stretches 100
3. Harmonic bends 200
4. Harmonic wags 20
5. Generalized Lennard-Jones 200
6. Torsions 20
7. Exponential repulsions and attractions 100
8. Ghost pair interactions. 20
9. Tetrahedral centers 20
10. Non-diagonal stretch-stretch interactions 100
11. Non diagonal stretch-bend interactions 100
12. Bend-bend interactions 100
13. Dihedral angles 20
14. Axilrod-Teller interactions 300
15. SN2 potential energy functions ----
16. Rydberg functions 100
17. Hartree-Fock dispersion functions 100
18. LEPS potential functions, type A 100
19. LEPS potential functions, type B 100
20. DMBE IV potential energy functions ----
The next version of Venus96 will have all array dimensions and parameters stored and adjustable in an include file.
B. List and Description of the Subroutines and Functions
The main Venus96 program makes use of the following subroutines and functions:
RUNGEK perform one cycle of Runge-Kutta-Gill integration of the equations of motion
ADAMSM perform one cycle of sixth-order Adams-Moulton integration of the equations of motion.
WLBOND evaluate and write out local mode population arrays
ENMODE calculate normal mode energies
WENMOD evaluate and write out normal mode population arrays
EBOND calculate local mode (Morse oscillator) energies
WEBOND evaluate and write out excited local mode population arrays
GWRITE write relevant information in output file
SELECT select initial conditions for coordinates and momenta
RANDST set up data for random number generator.
TEST check for intermediate and final events
ENERGY calculate potential, kinetic and total energy of the molecular system
PARTI driver routine for calculating energy partial derivatives
DVDQ calculate potential energy derivatives with respect to coordinates
STRET calculate harmonic stretch potential energy derivatives
MORSE calculate Morse potential energy derivatives
HBEND calculate harmonic bend potential energy derivatives
HALPHA calculate harmonic alpha bend potential energy derivatives
LENJ calculate Lennard-Jones potential energy derivatives
HTAU calculate 2 and 3-fold torsion potential energy derivatives
HEXP calculate general exponential repulsion or attraction potential energy derivatives
GHOST calculate ghost pair interaction potential energy derivatives
TETRA calculate tetrahedral center potential energy derivatives
CUBEND part of the tetrahedral center potential energy derivatives evaluation. called by subroutine TETRA
VRR calculate non diagonal stretch-stretch potential energy derivatives
VRT calculate non diagonal stretch-bend potential energy derivatives
VTT calculate non diagonal bend-bend potential energy derivatives
DANGLE calculate dihedral angle torsion potential energy derivatives
AXT calculate Axilrod-Teller potential energy and derivatives
VSN2 calculate SN2 potential energy and derivatives
PARSN2 set parameters for the selected SN2 potential energy surface
RYDBG calculate Rydberg potential energy and derivatives
HFD calculate Hartree-Fock Dispersion potential energy and derivatives
LEPS1 calculate type A LEPS potential energy and derivatives
LEPS2 calculate type B LEPS potential energy and derivatives
DMBE Driver for the DMBE IV potential function energy and derivatives. Calls a series of routines and functions (ho2sur, ho2der, dexdis, delect, threbq, voh, ehfoh, disoh, dvoh, voo, ehfoo, disoo, dvoo, cef, exdis, elect, disp, ddisp) and uses data block HO2DAT
FINAL calculate the product final properties
CENMAS calculate the center of mass momenta and coordinates
ROTN calculate angular momentum, moment of inertia tensor, angular velocity, and rotational energy
ANGVEL subtract off a given angular velocity from a molecular system
GFINAL write out trajectory final analysis results in the output file
INITQP initialize coordinates and momenta from normal mode parameters
LMODE select the desired bond length for a diatom.
LMEXCT calculate energy in (Morse oscillator) local mode.
INITEBK initialize parameters for an oscillator with given quantum numbers n and j by semiclassical EBK quantization
FINLNJ calculate vibrational and rotational quantum numbers for a product diatom
GLPAR calculate parameters for Gauss-Legendre quadrature
HOMOQP initialize coordinates and momenta for a diatom
ROTEN select angular momentum and rotational energy
THRMAN calculate vibrational quantum numbers from a thermal (Boltzmann) distribution
ORTHAN initialize coordinates and momenta with orthant sampling
ROTATE randomly rotate a molecule about its center of mass by Euler's angles.
NMODE driver for normal mode analysis
FMTRX evaluate the force constant matrix numerically
FGMTRX construct an effective Wilson -FG- matrix. the normal modes and spectroscopic frequencies are then evaluated.
EIGOUT write out eigenvalues and eigenvectors
EIGN diagonalize a matrix using the Givens-Householder algorithm
MPATH perform reaction path following
GRCONV calculate mass-weighted and normalized gradient for reaction path following
GPATH calculate free energy along the reaction path
GINROT determine classical hindered rotor partition function by Gaussian quadrature
POTEN rotate molecule on its center of mass and calculate resulting potential energy
BAREXC choose momentum for the reaction coordinate (barrier excitation)
SYBMOL write out information for the Sybyl® .mol file for molecular graphics animation
FDATE get current time and date (version for Unix™ systems)
GDATE get current time and date (version for Cray-Unicos™ systems)
RAND0 generate a random number
RAND1 generate a random number chain using the multiplicative congruential method. called by RAND0.
GAMA generate a gamma-distributed random number.
PL evaluates Legendre polynomial and derivative.
C. Sample Input Datafiles
The sample input datafiles on the datafiles subdirectory are:
al13.dt5 Al13 icosahedral cluster dissociation. Trajectory calculations. Microcanonical sampling. Lennard-Jones + Axilrod-Teller potential energy surface from L.G.M. Pettersson, C.W. Bauschlicher and T. Halicioglu, J. Chem. Phys. 87, 2205 (1987).
al3.dt5 Al3 triangular cluster dissociation. Trajectory calculations. Orthant sampling. Rydberg potential functions.
ar13nm.dt5 Ar13 clusters. Norrmal mode analysis (NSELT=-1). Hartree-Fock Dispersion potential functions (cluster-type input).
arch4m.dt5 Ar + CH4 collision reaction. Trajectory calculations. Microcanonical sampling.
arh2o.dt5 Ar + H2O collision reaction. Trajectory calculations. Fixed normal mode sampling.
benzene.dt5 Benzene. Trajectory calculations. Local mode excitation. Normal mode energy monitoring.
cf3cn.dt5 CF3 + CN minimum energy path. Reaction path following.
cf3h.dt5 Trifluoromethane. Trajectory calculations. Local mode excitation.
cluster.dt5 H surrounded by Ar12 + CH3 collision. Trajectory calculations. Fixed normal mode sampling.
diamond.dt5 Two-layer diamond model. Trajectory calculations. Fixed normal mode sampling.
h2o.dt5 H2O unimolecular dissociation. Trajectory calculations. Fixed normal mode sampling.
h2odim.dt5 Water dimer equilibrium geometry. Minimum energy search by annealing procedure (NSELT=1).
h2oh2o.dt5 H2O + H2O collision dynamics. Trajectory calculations. Fixed normal mode sampling.
h2onm.dt5 H2O molecule. Normal mode analysis. Harmonic force field potential. Use of the checkpoint file is made (need to run the h2orp.dt5 datafile first).
h2orp.dt5 H2O molecule. Geometry optimization by reaction path following. Harmonic force field potential.
ho2rp.dt5 HO2 DMBE IV potential. Part of the isomerization reaction path, from the transition state to one of the HO2 isomers.
hch3.dt5 H + CH3 reaction. Reaction path following. 2-D hindered internal rotor treatment of the free rotations along the path. Potential from J. Am. Chem. Soc. (1987) 109, 2916.
heh2.dt5 He + H2 collision dynamics. Trajectory calculations. Diatomic (harmonic oscillator) sampling.
ho+hbr.dt5 Trajectory calculations for HO + HBr diatomics collisions, with an hypothetical potential model
lidme.dt5 Li+ + (CH3)2O minimum energy path. Reaction path following. 1-D and 2-D hindered internal rotor treatments of the free rotations along the path.
lmch4.dt5 Methane local mode excitation. Trajectory calculations. Normal mode energy monitoring.
mimo.dt5 (Li+)H2O + H2O minimum energy path. Reaction path following. 1-D and 2-D hindered internal rotor treatments of the free rotations along the path.
nacrown.dt5 Na+ + D3d 18-Crown-6 complexation reaction. Reaction path following.
nmch4.dt5 Methane vibrational frequencies. Normal mode analysis.
nono.dt5 Van der Waals NO-NO dimer. Trajectory calculations. Fixed normal mode sampling.
propyn.dt5 Propyne vibrational frequencies. Normal mode analysis.
sn2bar.dt5 Cl- + CH3Br SN2 reaction complex. Trajectory calculations, started at the saddlepoint. PES1(Br). Use of the checkpoint file is made for the geometry (need to run the sn2ts.dt5 file first).
sn2nm.dt5 Cl- + CH3Cl SN2 reaction complex. Normal mode analysis. PES1.
sn2rp.dt5 Cl- + CH3Cl SN2 reaction. Reaction path following. PES2.
sn2sel.dt5 Cl- + CH3Br SN2 bimolecular collision. Trajectory calculations. PES1(Br). Erel selected from Boltzmann distribution at 300 K.
sn2trj.dt5 Cl- + CH3Br SN2 bimolecular collision. Trajectory calculations. PES1(Br). Erel fixed at 0.5 kcal/mol.
xe2xe.dt5 Xe2 + Xe minimum energy path. Reaction path following.