Electronic structure computations and theoretical chemical kinetics: developments at the interface

N•H Computer Physics Communications ELSEVIER

Computer Physics Communications 84(1994) 115—130

Electronic structure computations and theoretical chemical kinetics: developments at the interface Michael Page Department of Chemirtiy, North Dakota State University, Fargo, ND 58105, USA Received 1 March 1994

Abstract

The past quarter century has seen considerable advances both in the accurate ab initio calculation of features of the Born—Oppenheimer potential energy surface and in the statistical dynamical calculation of chemical reaction rates. Although each of these areas has in the past played a significant role in fueling advances in the other, the relationship between them is now more immediate. These areas are now, and will be for the forseeable future, inextricably linked.

1. Introduction About the time the first issue of Computer Physics Communications was published, 25 years ago, Pulay [1] first demonstrated the practical implementation of a class of computation that has become a cornerstone of the application of ab initio electronic structure methods to problems in chemical kinetics. Pulay computed the derivative of the ab initio Hartree-Fock Born—Oppenheimer potential energy with respect to nuclear coordinates. Because the negative of this potential energy gradient represents the vector of forces on the atoms in a molecule, Pulay’s new technique was called the force method. Although the equivalent technique implemented within a semiempirical molecular orbital theory framework [2] rapidly blossomed into a viable tool for studying organic reaction mechanisms [31,the ab initio approach to reaction potential energy surfaces did not truly blossom for about a decade. At about the same time as Pulay’s early work, research in theoretical chemical kinetics focused largely on the investigation of the dynamical foundations of statistical reaction rate theories [4—111.This focus represented a renewed interest in the foundations of transition state theory, a theory that over the prior

three or four decades had emerged as an enormously popular and useful tool [12] for explaining and predicting trends among activated rate processes. Although the substantial progress made by these researchers paved the way for later developments, the ultimate goal of ab initio prediction of elementary reaction rate constants of useful accuracy was unrealized. The primary difficulty was that the potential energysurface information necessary for even a simple application of conventional transition state theory was largely unavailable. Not only was there no generally applicable computational approach to obtain such potential energy surface information, but the severe dependence of the predicted reaction rate on O010-4655/94/$07.OO 0 1994 Elsevier Science B.V. All rights reserved SSDI 0010-4655(94)00031-V

M. Page / Computer Physics Communications 84 (1994)115—130

116

the calculated barrier height factors of 6 and 31 in the rate accompanying errors of one and two kcal/mol, respectively, in the barrier height made the development of these approaches a daunting task. The onus through much of the 1970s was thus on the electronic structure community to (1) develop methods capable of describing electronic wavefunctions and energies in transition state regions of the potential energy surface regions often characterized by partially formed and partially broken chemical bonds, and (2) develop and refine techniques for searching and characterizing important regions of these many dimensional surfaces. Two achievements that had a profound influence on this endeavor are the efficient evaluation of derivatives of two-electron integrals [13] and the efficient implementation of optimization techniques for multiconfiguration self-consistent field wavefunctions [14]. With the enthusiasm surrounding these and other developments in electronic structure calculations as well as the anticipation of future capabilities, Truhiar, Garrett and others in the late 1970s and early 1980s developed practical theoretical methods to compute chemical reaction rates [15,16].The severe dependence of the reaction rate on the barrier height notwithstanding, these new methods reflected an awareness based on the earlier careful tests that the rate and its temperature dependence were also profoundly influenced by a much broader region of the potential energy surface. The new practical Generalized Transition State Theory or Variational Transition State Theory methods centered on the idea of separating out one degree of freedom (the path of steepest descent in mass-scaled Cartesian coordinates that passes throught the saddle point) as a reaction coordinate and approximating the potential energy surface as the energy profile along this reaction path plus a harmonic approximation to the potential energy for displacements perpendicular to this reaction path. Miller, Handy and Adams [17] derived the classical Hamiltonian for this reaction path potential and in doing so showed how to determine parameters of use in the variational transition state approaches. Throughout the 1980s, variational transition state theory methods and semiclassical tunneling corrections based on the reaction path potential were refined and improved though detailed tests [181. Although by the mid 1980s algorithms sufficient to perform variational transition state theory calculations directly using results of ab initio electronic structure calculations as input were in place, most actual applications used a more cumbersome two-step process: fit results of ab initio electronic structure calculations to an analytic functional form, and then determine variational transition state theory input from the analytic form of the potential energy surface. Because there is no general, widely applicable approach to performing this intermediate fitting step, substantial recent emphasis has been placed on the direct determination of variational transition state theory potential energy surface information from ab initio electronic structure calculations [19—22]. We discuss here some of the aforementioned advances in somewhat greater detail, with emphasis on the recent and current interplay between the capabilities of electronic structure computations and theoretical approaches to the calculation of chemical reaction rates. The discussion of this active and rapidly developing field is by no means exhaustive, and no doubt there are many important contributions that go unmentioned. We end with a brief prognosis for the future. —

—

—

2. Potential energy surfaces Central to the theoretical description of chemical phenomena and implicit in all of the preceding discussion is the Born—Oppenheimer approximation [23],which asserts that the motion of the electrons in a molecule is uncoupled from the motion of the nuclei. The Born—Oppenheimer approximation leads to the concept of a potential energy surface (PES). For each possible arrangement of the nuclei in a molecule, a potential energy can be determined by solving the quantum mechanical equations of motion for the electrons. The PES is the function that describes how this potential energy changes as the nuclei

117

M Page /Computer Physics Communications 84 (1994) 115—130

move relative to one another. Important regions of the PES include local minima, which represent equilibrium molecular configurations, simple saddle points, which represent to a first approximation regions of dynamical instability separating one local minimum or valley from another, and the broad vicinity of the low-energy valley connecting these minima and passing through the saddle point region [24]. Thermodynamic parameters, i.e., molecular heats of formation, depend on energies and other properties of molecules at or near their equilibrium configurations on the PES, while kinetic parameters are most sensitive to the PES in the broad vicinity of saddle-point regions, which generally are characterized by partially formed and partially broken chemical bonds. We are not concerned with efforts to fit ab initio, semiempirical and experimental data to global functional forms to represent the PES; rather, our focus is the direct determination of important features of the potential energy surface from ab initio electronic structure calculations. In choosing the coordinate system for describing these features, we note that one of the most useful and enduring mental constructs in all of theoretical chemistry is the representation of chemical dynamics the complex simultaneous motion of several atoms in a molecule by the rolling or sliding of a single particle on a many dimensional surface. This concept is brought to life through our everyday experience with gravitational potential energy. Words and phrases like hills, valleys, mountain passes and plateaus are pervasive in discussions involving PESs, and, in fact, these phrases can all be found in Eyring’s paper at the 1937 meeting of the Faraday Society on reaction kinetics [25]. For this gravitational potential energy construct to hold, the coordinate system must be isoinertial; that is, each coordinate direction must have the same reduced mass. This can be accomplished by mass-scaling the 3N Cartesian coordinates of the N atoms, —

—

x1

=

(1)

X3N

where, for example, ma and Za are the mass and Z coordinate of atom a, respectively. Collecting the mass-weighted Cartesian coordinates into a vector, x, much of what we discuss in the next several sections is based on the following local quadratic approximation to the energy, tFo(x_x E=E0+got(x_xo)+(x_xo) 0)+H.O.T. (2) Here the vector g0 and the matrix F0 represent the first (gradient) and second (Hessian or force constants) derivatives of the energywith respect to mass-weighted Cartesian coordinates evaluated at the point x0. The quantities E0, g0, and F0 represent the basic tools of PES exploration from the point of view of ab initio electronic structure calculations. Although curvilinear coordinates are often more useful, ab initio electronic structure calculations are wedded to rectilinear displacements. At critical (stationary) points, the vector g0 vanishes, and the eigenvectors and eigenvalues of F0 represent normal-mode displacements and squares of the harmonic vibrational frequencies, respectively. At local minima, the matrix F0 is positive definite, while at simple saddle points it has one negative eigenvalue. The path of steepest descent in mass-weighted Cartesian coordinates that passes through the saddle point is known as the minimum energy path (MEP) [26]or intrinsic reaction coordinate (IRC) [27,28]and is a convenient definition of a reaction path for use in variational transition state theory calculations. This path corresponds physically to the infinitely damped trajectory initiated (in both directions) from the saddle point with zero kinetic energy. The MEP satisfies the differential equation dx(s) —g(s) ds = c (3) ‘

118

M Page / Computer Physics Communications 84 (1994) 115—130

where c, the norm of the gradient vector, insures that the right hand side represents a vector of unit length pointing in the steepest downhill direction. The unit vector, ri(s), is thus directed along the tangent to the MEP. The solution to Eq. (3), x(s), is a vector of the mass-weighted Cartesian coordinates as a function of the distance (arc length) parameter, s. Because the energy is not known as an analytical function of the Cartesian coordinates, the solution to Eq. (3) must be determined numerically. 3. Born—Oppenheimer potential energy derivatives In 1986, Yaniaguchi and Schaefer [29]wrote a review article titled “A New Dimension to Quantum Chemistry,” in which they state “Analyticderivative methods have, over the last fifteen years, revolutionized the way in which quantum chemistry is done.” The new dimension these authors referred to is the calculation of analytic derivatives of the Born—Oppenheimer potential energy with respect to nuclear coordinates. The revolution concerns the resolution of the dimensionality dilemma for exploring certain aspects of many dimensional potential energy surfaces. To illustrate the dilemma and the resolution, consider methanol (CH3OH), with six atoms and 18 Cartesian degrees of freedom. For some arbitrary orientation of the molecule, we can perform an electronic structure calculation and obtain a point on the PES, but we learn nothing about the dynamics of the molecule unless we consider how the PES changes as the atoms move. The gradient vector in Eq. (2) can be obtained by a finite-difference procedure: compute the energy for very small displacements in each of the coordinate directions and see how much it changes. Using a two-point difference method for acceptable accuracy, the Cartesian gradient takes, in this case, 36 energy evaluations (reduced to 24 by exploiting the invariance of the energy to overall translations and rotations [301).The matrix of force constants can be obtained by finite differencing the gradients. For methanol, however, if an energy evaluation takes one unit of effort, then a finite-difference gradient evaluation takes 24 and a force constant calculation about 600 units. In contrast, determining these quantities analytically, i.e., by the direct evaluation of the algebraically differentiated electronic energy expressions, takes roughly one, three and ten units of effort for evaluating the same quantities [31] *• Further, the ratio of these efforts force constants to gradient or gradient to energy is only weakly dependent on the number of atoms for analytic evaluation, in contrast to being directly proportional to the number of atoms for the finite-difference procedure. In analytic derivative evaluation, the energy expression for a particular electronic structure method, the expectation value of the Hamiltonian, is first differentiated and then the resulting expression is evaluated. Although straightforward in principle, the algebra tends to be cumbersome and there can be data handling problems. There are fortunately some significant simplifications for wavefunctions that have been determined by the variational principle. If we separate the parameters determining the wavefunction into those that are determined by the variational principle, {q}, and those that are not variational, {p}, then the derivative of the Born—Oppenheimer potential energy with respect to a nuclear coordinate, x (or in general any parameter in the Hamiltonian), can be written as d aEa~ ~—{E(x, p(x), q(x))} = + E-~-~. (4) —

—

(~_)+ ~

The difficult parts of the calculation are evaluating the derivatives of the parameters {p} and (q} with

These ratios depend substantially on the electronic structure method and its specific implementation. The numbers given are for a 1760-configuration complete active space (CAS) self-consistent field calculation on CH3COOH with a correlation consistent polarized valence double zeta basis set using the MESA L31] program on an IBM RS6000/580. The force constants take about one hour of CPU time.

*

M. Page / Computer Physics Communications 84 (1994) 115—130

119

respect to x. Fortunately, because the parameters {q} are determined through the variational principle by the requirement that they minimize the energy, dE/dq = 0, the last term in the above expression vanishes. For example, the derivatives of the molecular orbitals with respect to nuclear coordinates do not have to be determined to calculate gradients for Hartree—Fock or multiconfiguration SCF wavefunctions, since the orbitals are determined variationally ~ Likewise, the derivatives of the linear configuration expansion coefficients for MCSCF or configuration interaction (CL) wavefunctions do not have to be determined. For CL wavefunctions however, for which the MOs have not been determined variationally, the simplification involving the MOs is no longer available. The derivatives of the molecular orbital parameters are determined though the coupled-perturbed Hartree—Fock (CPHF) or coupled perturbed MCSCF equations [32,33].These equations require the solution of a large linear system of equations for each nuclear degree of freedom. Pople and coworkers [34] introduced an efficient iterative method for the solution of these equations. Handy and Schaefer [35] showed that by rearranging the order of summations in the derivative expression, the full coupled perturbed equations need not be solved for the determination of gradients in methods for which the orbitals are not determined variationally. Instead of 3N linear systems of equations, only one system must be solved. Rice and Amos [36] also present the general correlated gradient expression in a form that significantly reduces the complexity of the calculation. A major part of the derivative evaluation for all the wavefunction formulations considered here is the evaluation of the derivatives of the two-electron integrals with respect to nuclear displacements. These are included in the above expression as part of dp/dx. These integrals are the basic building block of the electronic energy expressions. The Hartree—Fock energy expression for a calculation with order of 102 basis functions, for example, has 0(108) terms, each involving orbital expansion coefficients multiplied by an integral over up to four basis functions. Since the basis functions are centered on the atoms, these integrals depend on the positions of the atoms and must be differentiated. The derivative of a two-electron integral is itself a small linear combination of new two-electron integrals involving slightly more complicated functions. Second derivatives of the energy involve second derivatives of the integrals which leads to more integrals that need to be evaluated. Recent advances in the efficient evaluation of these derivative integrals [13]are a major catalyst for the rapid development of derivative methods. A lucid and detailed review of analytic derivative techniques through 1987 has been presented by Pulay [32],and many of these techniques are elaborated on in Ref. [33].Some more recent advances include gradients for Moller—Plesset perturbation theory to fourth order [37];gradients for second order Moller—Plesset perturbation theory based on a spin-restricted Hartree—Fock reference function [38]; gradients for an implementation of the coupled cluster singles and doubles method [39,40], and its perturbative triples correction [41]; and gradient for the direct (i.e., no storage of the two-electron integrals) implementation of second-order Moller—Plesset perturbation theory [42]. 4. The reaction path Ilamiltonian Shortly following the first useful implementation of analytic second derivatives for Hartree—Fock wavefunctions [35],Miller, Handy and Adams [17]further entrenched the reaction path idea with the derivation of the classical Hamiltonian for a simple potential based on the MEP. A number of dynamical

*

Actually optimum molecular otbitals do change to first order, but only to preserve orthonormality as the metric of the space

spanned by the atom-centered basis functions changes. The parameters (q) in Eq. (4) mix the molecular orbitals while preserving orthonormality.

120

M. Page/Computer Physics Communications 84 (1994) 115—130

models based on this Hamiltonian have since been proposed. The reaction path Hamiltonian (RPH) is an outgrowth of a model developed in rotation-vibration spectroscopy theory [43]to treat systems where one degree of freedom undergoes large-amplitude motion. The idea, conceptually, is to consider the potential as a trough or a streambed plus 3N 7 harmonic walls that are free to close in or widen out as one proceeds along the trough. The PES is approximated as the potential energyof the MEP (V0(s)) plus —

a quadratic approximation to the energy in directions perpendicular to the MEP, V(s, Q1, Q2, Q3N—~)= V0(s) +

(5)

...,

k

Here the {QJ are referred to as generalized normal coordinates and the (w} are the associated harmonic frequencies. They are obtained at each point on the path by diagonalizing the force constant matrix for which the reaction path direction as well as directions corresponding to rotations and translations have been projected out [17].Thus the projected force constant matrix is obtained as K(s) = (I-- P(s))F(s)(I P(s)). (6) Here F is the matrix of force constants expressed in mass-weighted Cartesian coordinates, and P is a projection matrix expressed in terms of the unit reaction path tangent, v0, and six vectors corresponding to overall translation and rotation, v~as —

3N—7

~

v1v~.

(7)

Diagonalization of K gives the generalized normal-mode eigenvectors, L., and the transverse harmonic vibrational frequencies, o~. K(s)L1(s) =o~L~(s). (8) The generalized normal modes thus change with arc length along the reaction path. If (~}represent Cartesian displacements from the reaction path at a particular value of the arc-length parameter, s, then the kth generalized normal coordinate, Q,~,can be represented as 3N

Q~(s)= E~IL~k(s).

(9)

To obtain the Hamiltonian function for this reaction path potential, it is necessary to express the kinetic energy in terms of the momenta congugate to the reaction path coordinates. The necessary coordinate transformation, as well as a second transformation to the harmonic action-angle variables, are presented in the original paper [17].The resulting classical Hamiltonian for the zero-angular-momentum case in terms of the action-angle variables is H(p~,s, n, q)

=

F-i ~

(n,~+

~-)Wg(S)

+

V0(s)

k

2nk. + +

—

l)ük~/Wk~

cos

~ ~/(2n~+ 1)( k,k’

[1 + ~ sl(2nk + l)/Wk

2

.

(10)

5~fl~kh3k,F(5)]

k

This second transformation is particularly instructive because the Hamiltonian reduces to a simple Hamiltonian for just the reaction path degree of freedom with an effective one-dimensional potential if all of the (B) coupling terms are set to zero. 3N—7

H(s,p~)= ~ k~1

(nk+~)wk(s)+Vo(s)+fp~.

(11)

M. Page / Computer Physics Communications 84 (1994) 115—130

121

Because this approximate form of the reaction path Hamiltonian is independent of the angle variables {q}, the conjugate momenta, {n}, are constants of the motion. The classical action variable ~k can be identified with a quantum number associated with the kth generalized vibrational mode perpendicular to the path, and neglecting the (B) coupling elements can be identified with a vibrationally adiabatic approximation. Note that although the energy in these vibrations is not conserved, the vibrational action, which is the ratio of the energy to the frequency, is conserved. As the transverse frequencies change along the MEP, energy leaks in or out of these degrees of freedom giving rise to an effective one-dimensional potential. The effective potential obtained with all these vibrational modes in their ground state, all nk = 0, is the vibrationally adiabatic ground-state potential, often referred to as Va°(S) or VAG [44].A family of vibrationally adiabatic potentials can be obtained by considering the different sets of values for the quantum numbers ~ These considerations provide a powerful model for predicting rate enhancement associated with reactant vibrational excitation [18,44]. The coupling terms that appear in the reaction path Hamiltonian can be given a geometrical interpretation. The most important coupling terms are those that couple the transverse vibrational modes directly to the reaction path, BkF, where the index k runs over the 3N 7 transverse modes and F denotes the reaction path degree of freedom. These are called the curvature coupling elements because they are a measure of the degree to which the reaction path curves into a particular transverse mode as the reaction coordinate is traversed. The curvature coupling elements are also important ingredients in the description of quantum mechanical tunneling effects for the reaction path degree of freedom in variational transition state theory calculations [16,45]. The curvature of the reaction path can be obtained by differentiating Eq. (3) with respect to s. The result is [191 —

dy(s) ds __[Fv_(vtFv)v]/c. (12) The scalar curvature is just the length of this vector, and the curvature coupling elements in the reaction path Hamiltonian are just elements of the curvature vector expressed in the basis of generalized normal modes. Bl,F =

—

( )

dy(s) ~ ds L1

=

—v~F0L~/c.

(13)

Note that the curvature can be calculated from only first and second derivatives at a point on the path. The classical notion that a trajectory will overshoot the path and climb the wall if the path curves on the way down the hill is a reflection of this curvature coupling. Climbing the wall in a transverse direction is tantamount to exchanging energy between the reaction path and the transverse vibration. The curvature coupling elements are also very important for the treatment of quantum effects in reaction coordinate motion. If the path is curved in the vicinity of the saddle point, then a quantum particle might be able to cut the corner and tunnel through a narrower section of the potential. It should be cautioned however that Eq. (12) becomes indeterminate at the saddle point and the curvature must be obtained from a limiting equation based on l’Hospital’s rule [191.Determining the curvature at the saddle point as well as the remaining coupling elements in the reaction path Hamiltonian require knowing limited information about the third derivatives on the PES. Specifically, they require the derivative of the force constant matrix with respect to the reaction coordinate direction, which can be obtained by finite differences, ~

dF

F(s0+8)—F(s0—8) 28

(14)

These remaining coupling elements can be identified with the direct coupling of transverse vibrational

M. Page/Computer Physics Communications 84(1994) 115—130

122

modes (the mode—mode coupling elements) and with the derivatives of the transverse frequencies with respect to s, dK(s) B~,1=J4 ds L~/(w?—cofl, iv’j, and d(w~)

(15)

dK(s)

B~,1= ds =L~( ds )L1. (16) Here dK/ds can be obtained from dF/ds (Eq. (14)) and the derivative of the projection matrix [19].All of the terms in the reaction path Hantiltonian can thus be determined from ab initio electronic structure calculations given x(s), g(s), F(s) and F’(s). 5. VarIational transition state theory The most fruitful use of reaction path potentials has been as input to variational transition state theory (VTST) and semiclassical tunneling calculations. Truhlar, Garrett, and coworkers have over the past 15 years developed VTST into a practical computational method primarily for computing rates of gas-phase chemical reactions and also more recently for the study of solvation effects [461and heterogeneous reactions on surfaces [47].Many of the wide variety of approaches are implemented and available in a computer program called POLYRATE [22].A comprehensive treatment of much of the theory is given in Ref. [431. We present here a brief outline that illustrates the basic ideas and shows how the necessary input is related to the information available from ab initio electronic structure calculations. The fundamental assumption of conventional transition state theory (TST) (in classical mechanics) is that a plane in configuration space passing through the saddle point on the potential energy surface and perpendicular to the non-bound normal coordinate provides a perfect dynamical bottleneck for the reaction: no trajectory that passes through this surface from reactants to products ever returns. If trajectories do recross, then the TST rate is too high. The variational approach recognizes that whenever this fundamental assumption is not satisfied, any alternative dividing surface (in phase space) that has fewer one-way trajectories passing through it provides a better approximation to the actual one-way rate, and the best dividing surface (transition state) is the one with the least rate of passage. In microcanonical VTST, the variationally best transition state is found at each energy, and then a thermal rate constant can be obtained by appropriate averaging with a Boltzmann probability term. Canonical VTST (CVT) is easier to apply; one simply finds the best compromise transition state at each temperature. In the usual hybrid approach to quantizing transition state theory [48],classical partition functions are replaced with quantum partition functions for all bound modes, while the reaction coordinate is treated classically. The application of CVT then requires a search at each temperature over a range of possible dividing surfaces to find the dividing surface that gives the lowest rate. This search in practice is limited to planar configuration-space dividing surfaces perpendicular to the MEP and identified by the arc length, s, k~(T)

=

min k(s, T).

(17)

Here the generalized TST rate constant at location is given by k(s, T)

=~r~—

~~‘~s)

e’~~~’IcBT.

(18)

M. Page / ComputerPhysics Communications 84 (1994) 115—130

123

Vibrational Frequencies ..U.

3000 ..UU

2500

.......u.....

2000

E

1500 saelsISsIsIIaes.m..$$*:$::::$$. ~

1000 500 0 -0.8 -0.6 -0.4 -0.2 -0.0 0.2 0.4 0.6 0.8

arc length [sqrt(amu)bohr] Fig. 1. Harmonic vibrational frequencies LComplete Active Space (11 electrons in 9 orbitals) correlation consistent Valence Double Zeta basis set] as a function of arc length along the minimum energy path for the reaction H + HNO —~ H 2NO (see Ref. [51]).

In Eq. (18), o~is a symmetry factor accounting for ‘two or more symmetry-related reaction paths [491 and is used in place of symmetry numbers in all rotational partition functions, kB is the Boltzmann constant, h is Planck’s constant, VMEP(s) is the potential along the MEP relative to reactants, Q~(T,s) and ‘~P’~(T) are the partition function of the generalized transition state at s and the reactant partition function per unit volume, respectively. The partition functions in Eq. (18) typically are approximated as products of translational, rotational, vibrational, and electronic partition functions. Furthermore, a separable mode formalism is generally applied to compute the vibrational partition function [161,and the individual modes are often treated as harmonic oscillators or sometimes hindered rotors [50]. Aside from a possible barrier to hindered rotation, all of the information necessary to apply Eqs. (17) and (18) using ab initio potential energy surface information is available from x(s), g(s), and F(s). From the perspective of reaction path potentials, variational effects the deviation between the CVT and TST rate coefficients may arise because the harmonic walls close in as one moves from the saddle point toward reactants or products. This narrowing of the passageway has greater influence on the reaction at higher temperatures. The CVT approach, in fact, can be formulated as locating a maximum in the Gibbs free energy of activation along the path. The origin of the variational effect, then, corresponds to a change in the entropic contribution to the free energy of activation along the path. Figs. 1 and 2 help to illustrate these ideas. Fig. 1 shows the harmonic vibrational frequencies along the MEP for the reaction H + HNO H2NO. These vibrational frequencies were obtained at several points along the MEP by diagonalizing the projected force constant matrix computed at the MCSCF level [51]. Large negative arc lengths correspond to the H + HNO reactants. Here, as can be seen from the leftmost portion of Fig. 1, there are three non-zero vibrational frequencies, corresponding to the NH stretch, the NO stretch, and the HNO bend of HNO. The other two modes are called transitional modes because although they are bonafide vibrations in the vicinity of the transition state (and in this case at the product too), they correlate with relative translations freesmall rotations These two modes t overand a fairly rangeofof the the reactants. reaction path, corresponding increase in frequency by nearly 1000 cm to the approaching hydrogen atom moving about 0.5 angstroms. This tightening of these incipient —

—

—‘

124

M Page/Computer Physics Communications 84 (1994) 115—130

Location of Transition State 3 2

Ii

600K//.

s~ox

/ OK’

~

•.

-0.8 -0.6 -0.4 -0.2 0.0 arc~ngth (sqt(anii)bofrj

0.2

0.4

Fig. 2. Temperature dependence of the location of the variationally best transition state for the reaction H+HNO —~H

2NO (see

Ref. [51]).

bending degrees of freedom implies increasingly severe geometrical restrictions on the low-energy structures. Fig. 2 shows the ground state vibrationally adiabatic energy profile for the H + HNO H2NO reaction. Note that the maximum in this curve is not at s = 0.0. This is for two reasons: (1) the energy profile here is calculated at a higher level of theory than that used to determine the MEP and the vibrational properties [51], and (2) the zero-point vibrational energy contributions are included in the profile, so even if this profile was not determined at a higher level of theory, the top of the hill would not correspond to s = 0.0. Fig. 2 identifies the location of the variationally determined transition state at several temperatures. Note that at low temperatures, the maximum in the vibrationally adiabatic profile corresponds to the variational transition state, while as the temperature is raised the transition state moves to the right, i.e., toward tighter structures where the frequencies of the transitional modes are greater. —~

6. Semiclassical tunneling corrections The previous discussion indicates that variational effects do not lead to large deviations of the variational state location from the saddle point at low temperatures. It is tempting to conclude that detailed information about the shape of the energy profile and properties along the MEP are unimportant at low temperature. On the contrary, often reaction path Hamiltonian parameters are needed over an even wider portion of the MEP to obtain accurate low-temperature rate constants. This is because at low temperatures the approximation that the passage through the transition state dividing surface can be treated as a classical translation breaks down. Whereas a classical treatment of passage over the barrier requires no potential energy surface information, a reasonably accurate quantum treatment requires a considerable amount of potential energy surface information. A quantum treatment of the reaction coordinate degree of freedom must consider quantum mechanical tunneling for energies below the barrier height and quantum mechanical reflection for energies above

M Page/Computer Physics Communications 84 (1994) 115—130

125

the barrier height. When these energy dependent transmission probabilities are integrated with a temperature-dependent Boltzmann probability term, a temperature-dependent multiplicative factor is obtained. Thus a quantum mechanical treatment of reaction coordinate motion simply introduces a transmission factor, K(T), and the rate constant is obtained by first computing k~’T(T) and then correcting it with a transmission factor, k(T)==K(T) k~(T).

(19)

Computing K(T) is not so easy. The simplest approach is to consider transmission through an inverted parabolic barrier determined from the one negative eigenvalue of the force constant matrix at the saddle point. But this is bad for two reasons: first, this approximation to the potential energy profile is only good close to the barrier top, whereas the dominant tunneling contribution is often from energies well below the top of the barrier; and second, to the extent that a vibrationally adiabatic approximation is valid, this is not even a good local approximation to the effective one-dimensional potential. A better approximation is to consider tunneling through the actual shape of the ab initio vibrationally adiabatic potential energy barrier. This requires that x(s) and {w(s)) be determined accurately over a fairly wide range on either side of the barrier. This is called the vibrationally adiabatic zero-curvature approximation. The barrier penetration calculation can be performed numerically through the actual potential energy curve or it can be performed analytically by first fitting the vibrationally adiabatic potential energy curve to a simple functional form such as an Eckart potential [52] for which the barrier penetration problem is solved. Because the MEP is a curved path, it is possible to have corner cutting effects, and thus the vibrationally adiabatic zero-curvature approximation underestimates the tunneling contribution. If the curvature is small, then a limiting equation that still uses only reaction path Hamiltonian information can be used. In particular the small-curvature semiclassical adiabatic ground-state (SCSAG) method [53] requires the reaction path curvature over a wide range of the MEP. A newer version, called the centrifugal dominant (CD) [22]small-curvature approximation requires the same information. When the reaction path curvature is large, such as in the exchange of a light atom between two heavy atoms, a broader swath of the potential energy surface than a quadratic expansion about the MEP is required. Other methods, such as the least-action ground-state (LAG) method [54] and the large curvature ground-state version 3 (LCG3) method [45]have been developed for these problems. The agreement between these semiclassical tunneling methods and accurate or best available quantum rate constants has been excellent [55]. 7. Determining the reaction path The most computationally intensive step in VTST calculations or in dynamical studies based on the reaction path Hamiltonian is the determination of the reaction path, x(s), by numerical integration of Eq. (3) and the evaluation of potential energy derivatives along the path. Many reaction paths are obtained using methods based on simple Euler integration of Eq. (3). Beginning at the saddle point, one first steps along the path tangent which, at the saddle point, is an eigenvector of the force constant matrix. Subsequent small straight line steps are taken along the local steepest descent direction, which requires only the gradient. x(s + 6s) —x(s)

—

(g/c)06s.

(20)

Because the MEP is actually a curved path, the Euler step of finite length starting on the path will lead to a new point that deviates from the true path on the convex side. The negative of the gradient vector at the new point then has a component that leads back toward the MEP as well as a component

126

M. Page / Computer Physics Communications 84 (1994) 115—130

that leads downhill along the MEP. A subsequent step may cross the MEP and end on the concave side. Thus oscillations about the MEP are not uncommon. A number of devices based on constrained energy minimization have been proposed to return to the path after too large an Euler step has been taken. Ishida et al. [56]considered a one-dimensional search along the bisector of the old and new gradient. Schmidt et al. [57]refmed this procedure by computing one additional energy point along the bisector and taking the correction step to be the minimum along a parabolic fit to this information. The efficiency of this algorithm depends on the specification of two parameters, the distance along the bisector used to perform the extra energy evaluation, and the criterion for skipping the stabilization step for locally straight MEPs [58].Melissas et a]. [21]present best compromise values for these parameters and call the resulting method ES1* (Euler stabilization method, version 1). The ES1* method does remarkably well [211in comparison to schemes that require more potential energy surface information in tests where the reaction path information is used to obtain VTST and semiclassical tunneling results. Other methods [59—611include N 1 dimensional constrained energy minimizations as correction steps. The method of Gonzalez and Schlegel [60]is particularly attractive. The new point on the path is chosen such ‘that the old point and the new point lie on the arc of a circle, and the old and new gradients are tangent to this arc. In the limit of small step lengths, both the path tangent and the path curvature are predicted correctly. Also the constrained optimization step requires only a small number of additional gradient evaluations because the well developed minimization technology [62],such as useful guesses to the Hessian matrix and efficient Hessian updating schemes can be used profitably. Gonzalez and Schiegel [611have shown that their method is one of a general class of implicit methods for the numerical solution of ordinary differential equations. Implicit here means that the value of the function at the endpoint is required. The contribution of Gonzalez and Schlegel is the use of the constrained optimization as the manner of obtaining the value of the function at the endpoint. A host of methods in this general class have been implemented and tested on model two-dimensional surfaces. These authors, as well as Garrett, et al. [581and Baidridge et al. [20]have concluded that higher-order algorithms are not necessarily more efficient than lower-order algorithms. It is important to note that there exists no local criteria for determining whether a given point is on the MEP. Every point on the PES is on some solution to Eq. (3), but the MEP is the one solution that passes through the saddle point. Therefore the constraints in the constrained-optimization approach depend on the validity of reaction path characteristics for previous points along the path. Another class of methods includes no stabilization steps and can be discussed in terms of a Taylor series representation of the path in the arc length parameter, s, about a point on the path at arc length ~v~1)(s ~ (21) x(s) =x0(s) + ~~°~(s s0) + fv~(s so)2 + —

...

—

—

—

The coefficients, which are all vectors, can be interpreted [19,63,64].The first-order coefficient is the unit tangent to the MEP; it is the negative of the normalized gradient vector. The second-order coefficient is the curvature vector evaluated at s 0. It can be computed (Eq. (12)) using g0 and F0. The third-order coefficient describes the first-order change in the curvature vector and can be computed with g0, F0, and P0. (Eq. (14)). A general expression for all of the coefficients in Eq. (21) in terms of energy derivatives have been given in Ref. [191.The third- and higher-order terms in Eq. (21) require third- and higher-energy derivatives evaluated at the point of expansion, but even without these higher energy derivatives, the contributions to these terms from first-, second-, and partial third-derivative information can be obtained. Thus the LQA method [191takes a step along the exact solution on a locally quadratic surface and therefore includes the second-order contribution to all the terms in Eq. (21). Use of a local quadratic approximation to the energy in the solution of Eq. (3) was introduced by Pechukas as a

M Page / Computer Physics Communications 84 (1994) 115—130

127

pedagogical tool to study the conservation of symmetry along paths of steepest descent [651.It has since been used computationally as an aid to function minimizations, both to optimize wave function parameters [66]and to optimize geometries [67].The LQA method has been used to calculate the MEP for several systems [51,68—75]. If F~is available, then the third-order term can be computed exactly, and the LQA contribution to the remaining terms can be retained. This is the corrected LQA (CLQA) method [641.The CLQA method is a significant improvement on the LQA method and can be applied approximately with no additional electronic structure calculations by estimating F~from successive force constant calculations along the path [64].The significance of these methods is that they are based on the same potential energy surface information as that necessary to compute the reaction path Hamiltonian, so the two problems are solved in concert. Sun and Ruedenberg [76] recently have shown that the integration of Eq. (3) using local quadratic approximations to the energy can be significantly improved by relaxing the requirement that the expansion point lie on the path. These authors also introduce methods based on trust regions for assessing the accuracy of the integration, a very important issue because, as discussed earlier, there exist no local criteria for determining this accuracy. An advantage of this method is that larger step sizes can be taken for the reaction path determination. A disadvantage is that force constants on the path for determining vibrational frequencies along the path are not a natural byproduct of the method. Therefore for use in VTST calculations, either additional force constants will need to be determined, or force constants from points close to the path can be used to approximate vibrational frequencies on the path. Recently [21]a number of techniques for determining the MEP have been distinguished according to the effort necessary to determine VTST rate constants with semiclassical tunneling corrections to a given accuracy. The methods were tested for two different reactions calculating rate constants at two different temperatures. The computational effort was evaluated in each case for different hypothetical values of the ratios of efforts to calculate F0, g0, and E. Over the twelve test cases, the most efficient of the tested algorithms for cases with relatively low force-constant cost was the LQA method and with relatively high force-constant cost was the ES1 method. The CLQA method [63],the higher-order implicit methods of Gonzalez and Schlegel [61]and the new method of Sun and Reudenberg [76]have all not yet been tested in actual dynamical calculations. *

8. Prognosis The accuracy of ab iitio potential energy surface information is of paramount importance for the prediction of reaction rates. There is currently intense activity in developing new methods for computing dynamical electron correlation effects as well as developing ways to extend existing methods to larger systems. Two recent developments that will play a significant and increasing role in developing and testing new electronic structure methods are the full-CL benchmark calculations of Bauschlicher, Langhoff, Taylor and coworkers [77]and the basis set benchmark calculations of Woon, Peterson, Kendall and Dunning [78,79]. These complementary sets of calculations are both designed to assess existing methods and provide a framework to assess new methods as they are developed. Two sources of error in ab initio electronic structure calculations are truncation of the one-electron basis set and truncation of the n-electron basis set. These errors, unfortunately, are not uncoupled. The full-CL benchmark calculations [77] attempt to choose a one-electron basis set that is large enough to reduce the most signfficant coupling between the two sources of error. Then by performing CL calculations that include all possible configurations of a given spin and space symmetry, the n-electron basis set error is eliminated. These calculations, which have now been performed over a wide variety of

128

M. Page / Computer Physics Communications 84 (1994) 115—130

electronic structure environments, have been used and will continue to be used to provide general criteria by which to judge approximate treatments of the n-electron space and to provide guidance regarding which methods are appropriate for which problems. The basis set benchmark calculations [78,79],by contrast, attempt to choose a n-electron basis set that is large enough to reduce the most significant coupling between the two sources of error. Then by performing a wide variety of calculations with systematically improved one-electron basis sets, the convergence of the calculations with respect to the one-electron space can be established. Because there is no reliable internal measure of the accuracy of electronic structure calculations, these benchmark calculations and extensions of these that are sure to be available in the future will play a crucial role in developing and assessing new electronic structure techniques. A similar situation exists for assessing and refining statistical reaction rate theories: comparison of computed rates to experimental measurements suffers from the inability to separate errors in the rate determination from errors in the potential energy surface. Therefore, as mentioned previously, comparison of variational transition state theory and semiclassical tunneling models to converged (exact) calculations of the quantum dynamics [18,55,80] is of paramount importance. There is currently vigorous activity in the computation of quantum dynamics for many-degree-of freedom systems [81]. One area that is sure to receive attention in the future is the development of a more detailed understanding of the sensitivity of the computed VTST/ST results to specific regions of the potential energy surface. Such an understanding would complement current work in interpolation methods [82] and in semiempirical or hybrid methods [83] for obtaining the potential energy surface information. In interpolation methods, the path dependence of the basic reaction path ingredients (energy, frequencies, moments of inertia, curvature coupling elements) are fit to simple functional forms using limited input from ab initio electronic structure calculations. In the semiempirical approach [83], approximate molecular orbital theory methods are parameterized for each problem to give desired results for certain features of a specific potential energy surface. These features can be designed to agree with experiment or can be chosen as more accurate ab initio results of limited portions of the potential energy surface. The fields of ab initio electronic structure theory and reaction rate theory have over the last quarter century become inextricably merged. Now and in the future, new developments in each of these areas will be spawned by and will spawn new developments in the other.

Acknowledgements This work was supported by the Mechanics Division of the Office of Naval Research (grant No. N000149310045) and by the NSF EPSCoR program (grant No. R11861075).

References [1] P. Pulay, MoL Phys. 17 (1969) 197; 18 (1970) 473; 36 (1971) 2470. P. Pulay and W. Meyer, J. Mo!. Spectrosc. 41(1971) 59; J. Chem. phys. 57 (1972) 3337. [2] J.W. Mclver Jr. and A. Komornicki, Chem. Phys. LetI. 10 (1971) 303. [3] J.W. Mclver Jr. and A. Komornicki, J. Amer. Chem. Soc. 94 (1972) 2625. A. Komornicki and J.W. Mclver Jr., J. Amer. Chem. Soc. 95 (1973) 4512; 96 (1974) 5799. [4] J.C. Keck, J. Chem. Phys. 32 (1960) 1035; Adv. Chem. Phys. 13 (1967) 85. [5]D.L. Bunker, J. Chem. Phys. 37 (1962) 393; 40 (1964) 1946. [6] R.A. Marcus, J. Chem. Phys. 41 (1964) 610; 45 (1966) 3489. [7] E.M. Mortensen, J. Chem. Phys. 48 (1968) 4029. [8] R.E. Wyatt, J. Chem. Phys. 51(1969) 3481.

M. Page/ Computer Physics Communications 84 (1994) 115—130

129

[9] D.G. Truhlar and A. Kuppermann, Chem. Phys. Lett. 9 (1971) 269; J. Chem. Phys. 59 (1973) 395. [10] J.M. Bowman, A. Kuppermann, J.T. Adams and DO. Truhiar, Chem. Phys. Lett. 20 (1973) 229. [11] P. Pechukas and F.!. McLafferty, J. Chem. Phys. 58 (1973) 1622. [12] S.W. Benson, Thermochemical Kinetics, 2nd ed. (Wiley, New York, 1976). [13] M. Dupuis and H.F. King, Tnt. J. Quantum Chem. 11 (1977) 613. P. Saxe, Y. Yamaguchi and H.F. Schaefer, J. Chem. Phys. 77 (1982) 5647. H.B. Schiegel, J. Chem. Phys. 77 (1982) 3676. [141E. Da!gaard and P. Jorgensen, J. Chem. Phys. 68 (1978) 3833. B.H. Lengsfield III, J. Chem. Phys. 73 (1980) 382. K. Ruedenberg, L.M. Cheung and S.T. Elbert, Tat. J. Quantum Chem. 16 (1979) 1069. L.M. Cheung, K.R. Sundberg and K. Ruedenberg, tnt. J. Quantum Chem. 16 (1979) 1103. B.O. Roos, P.R. Taylor and P.E.M. Siegbahn, Chem. Phys. 48 (1980) 152. [15] B.C. Garrett and D.G. Truhlar, J. Phys. Chem. 83 (1979) 1052, 1079, 2921; 84 (1980) 805; J. Amer. Chem. Soc. 101 (1979) 4534; J. Amer. Chem. Soc. 101 (1979) 5207. D.G. Truhlar and B.C. Garrett, Accts. Chem. Res. 13 (1980) 440. [16] A.D. Isaacson and D.G. Truhiar, J. Chem. Phys. 76 (1982) 1380. D.G. Truhlar, A.D. Isaacson, R.T. Skodje and B.C. Garrett, J. Phys. Chem. 86 (1982) 2252. (17] W.H. Miller, N.C. Handy and I.E. Adams I. Chem. Phys. 72(1980) 99. [18] D.G. Truhlar and A.D. Isaacson, J. Chem. Phys. 77 (1982) 3516. D.C. Clay)’, B.C. Garrett and D.G. Truhiar, J. Chem. Phys. 78 (1983) 777. D.G, Truhlar, N.J. Kilpatrick and B.C. Garrett, J. Chem. Phys. 78 (1983) 2438. D.L. Bondi, J.N.L. Conner, B.C. Garrett and D.G. Truhlar, J. Chem. Phys. 78 (1983) 5981. D.G. Truhlar, B.C. Garrett, P.G. Hipes and A. Kuppermann, J. Chem. Phys. 81 (1984) 3542. A.D. Isaacson, M.T. Sund, S.N. Rai and D.G. Truhiar, J. Chem. Phys. 82 (1985) 1338. B.C. Garrett, N. Abusalbi, D.J. Kouri and D.G. Truhiar, I. Chem. Phys. 83 (1985) 2252. B.C. Garrett, D.G. Truhlar and G.C. Schatz, J. Amer. Chem. Soc. 108 (1986) 2876. R. Steckler, D.G. Truhlar and B.C. Garrett, J. Chem. Phys. 84 (1986) 6712. [19] M. Page and J.W. Mclver Jr., J. Chem. Phys. 88 (1988) 922. (20] K.K. Baldndge, M.S. Gordon, R. Steckler and D.G. Truhiar J. Phys. Chem 93 (1989) 5107. [21) V.S. Melissas, D.G. Truhiar and B.C. Garrett, J. Chem. Phys. 96 (1992)5758. [22] POLYRATE 4: A new Version of a Computer Program for the Calculation of Chemical Reaction Rates for Polyatomics, D.H. Lu, T.N. Truong, V.S. Melissas, G.C. Lynch, Y.P. Liu, B.C. Garrett, R. Steckler, A.D. Isaacson, S.N. Rai, G.C. Hancock, J.G. Lauderdale, T. Joseph and D.G. Truh!ar, Comput., Phys. Comm. 71(1992) 235. [23] H. Eyring, J. Walter and G.E. Kimball, Quantum Chemistry (Wiley, New York, 1944) pp. 190—192. [24] D.G. Truhiar and M.S. Gordon, Science 249 (1990) 491. (25] H. Eyring, Trans. Faraday Soc. 33 (1937) 3. [26] 1. Shavitt, J. Chem. Phys. 49 (1968) 4048. R.E. Weston Jr., J. Chem. Phys. 31(1959) 892. R.A. Marcus, J. Chem. Phys. 45 (1966) 4493; 49 (1969) 2610, 2617. H.F. Schaefer III, Chem. Brit. 11(1975) 227. [27] K. Fukui, in: The World of Quantum Chemistry, R. Daude) and P. Pullman, eds. (Reidel, Dordrecht, 1974) pp. 113—141. [28] K. Fukui, Acc. Chem. Res. 14 (1981) 368. [29] H.F. Schaefer and Y. Yamaguchi, J. Molec. Struct. (Theochem) 135 (1986) 391. [30] M. Page, P. Saxe, G.F. Adams and B,H. Lengsfleld III, Chem. Phys. Lett. 107 (1984) 587. [31] MESA (Molecular Electronic Structure Applications) P. Saxe, B.H. Lengsfield Ill, R, Martin, and M. Page, copyright 1990, University of California. [32] P. Pulay, Adv. Chem. Phys. 69 (1987) 241. [33] P. Jorgensen and J. Simons, eds., Geometrical Derivatives of Energy Surfaces and Molecular Properties, (Reidel, Dordrecht, 1986). (34] l.A. Pople, K. Raghavachari, H.B. Schlegel and J.S. Binkley, tnt. J. Quantum. Chem. 13 (1979) 225. [35] N.C. Handy and H.F. Schaefer 111, J. Chem. Phys. 81 (1984) 5031. [36] J.E. Rice and R.D. Amos, Chem. Phys. Lett. 122 (1985) 585. [37] J. Gauss, D. Cremer, Chem. Phys. Lett. 153 (1988) 303. E.A. Salter, G.W. Trucks and R.J. Bartlett, J. Chem. Phys. 90 (1989) 1752. [38] D.T. Tozer, J.S. Andrews, R.D. Amos and N.C. Handy, Chem. Phys. Lett. 199 (1992) 229. [39] N.C. Handy, J.A. Pople, M. Head-Gordon, K. Raghavachari and G.W. Trucks, Chem. Phys. Lett. 164 (1989) 185. [40] R. Kobayashi, N.C. Handy, R.D. Amos, G.W. Trucks, MJ. Frisch and J.A. Pople, J. Chem. Phys. 95 (1991) 6723.

130

M Page/ComputerPhysics Communications 84 (1994)

115—130

[41] R. Kobayashi, R.D. Amos and N.C. Handy, Chem. Phys. Lett. 184 (1991)195. [42] Mj. Frisch, M. Head-Gordon and l.A. Pople, Chem. Phys. Lett. 166 (1990) 275, 281. [43] J.T. Hougen, P.R. Bunker and J.W.C. Johns, J. Mo!. Spectrosc. 34 (1970) 136. [44] T.H. Dunning Jr., E. Kraka and R.A. Eades, Faraday Discuss. Chem. Soc. 84(1987)427. [45] D.G. Truhiar, A.D. Isaacson and B.C. Garrett, in: Theory of Chemical Reaction Dynamics, Vol. IV, M. Baer, ed. (CRC Press, Boca Raton, FL, 1985) pp. 65—137. [46] S.C. Tucker and D.G. Truhlar, J. Amer. Chem. Soc. 112 (1990) 3347. D.G. Truhlar, G.K. Schenter and B.C. Garrett, J. Chem. Phys. 98 (1993) 575. [47] J.G. Lauderdale and D.G. Truhlar, Surf. Sci. 164 (1985) 558. T.N. Truong and D.G. Truhiar, J. Phys. Chem. 91(1987) 6229. (48] See, for example, 1.M.W. Smith, Kinetics and Dynamics of Elementary Gas Reactions (Butterworth, London, 1980). [49] V.S. Melissas and D.G. Truhiar, J. Phys. Chem. 98 (1994) 875. [50] D.G. Truhiar, I. Comp. Chem. 12 (1990) 266. [51] M. Page and M.R. Soto, J. Chem. Phys. 99 (1993) 7709. [52) H.S. Johnston, Gas Phase Reaction Rate Theory (Ronald, New York, 1966). [53] R.T. Skodje, D.G. Truhiar and B.C. Garrett, J. Phys. Chem. 85 (1981) 3019. [54] B.C. Garrett and D.G. Truhlar, J. Chem. Phys. 79 (1983) 4931. [55] B.C. Garrett and D.G. Truhiar, J. Phys. Chem. 95 (1991) 10374. [56] K. Ishida, K. Morokuma and A. Komornicki, I. Chem. Phys. 66(1977) 2153. [57) M.W. Schmidt, M.S. Gordon and M. Dupuis, J. Amer. Chem. Soc. 107 (1985) 2585. [58] B.C. Garrett, M.J. Redmon, R. Steckler, D.G. Truhlar, K.K. Baidridge, D. Bartol. M.W. Schmidt and M.S. Gordon, J. Phys. Chem. 92 (1988) 1476. [59] K. Muller and L.D. Brown, Theor. Chim. Acta. 53 (1979) 75. [60] C. Gonzalez and H.B. Schiegel, I. Chem. Phys. 90 (1989) 2154; 94 (1990) 2553. [61] C. Gonzalez and H.B. Schiegel, J. Chem. Phys. 95 (1991) 5853. [62] H.B. Schlegel, Adv. Chem. Phys. 67 (1987) 249. [63] M. Page, C. Doubleday and J.W. Mclver Jr., J. Chem. Phys. 93 (1990) 5634. [64] E. Kraka and T.H. Dunning, in: Advances in Molecular Electronic Structure Theory, T.H. Dunning, ed. (JAI Press, Greenwich, CT, 1990) pp. 129—173. [65] P. Pechukas, J. Chem. Phys. 64 (1976) 1516. [66] R.N. Camp and H.F. King, J. Chein. Phys. 77 (1982) 356. [671J.M. McKelvey and J.F. Hamilton, J. Chem. Phys. 80 (1984) 579. [68] 1. Ischtwan and M.A. Collins, J. Chem. Phys. 89 (1988) 2881. [69] S. Koseki and M.S. Gordon, I. Phys. Chem. 93 (1989) 118. (70] C. Doubleday, J.W. Mclver Jr. and M. Page, J. Phys. Chem. 92 (1988) 4367. [71] N.J. Caidwell, J.K. Rice, H,H. Nelson, G.F. Adams and M. Page, J. Chem. Phys. 93 (1990) 479. [72] B.C. Garrett, M.L. Koszykowski, C.F. Melius and M. Page, J. Phys. Chem. 94 (1990) 7096. [73] M.R. Soto, M. Page and M.L. McKee, Chem. Phys. 153 (1991) 415. [74] M. Page, M.C. Lin, Y. He and T.K. Choudhury, J. Phys. Chem 93 (1989) 44414. [75] M.R. Soto and M. Page, J. Chem. Phys. 97 (1992) 7334. [76] J.Q. Sun and K. Ruedenberg, J. Chem. Phys. 99 (1993) 5257; 5269. [77] C.W. Bauchlicher, S.R. Langhoff and P.R. Taylor, Adv. Chem. Phys. 77 (1990) 103, and refererences therein. [78] D.E. Woon and T.H. Dunning Jr., J. Chem. Phys. 99 (1993)1914. [79] K.A. Peterson, R.A. Kendall and T.H. Dunning, J. Chem. Phys. 99(1993)1930, 9790. [80] D.C. Clary, J. Chem. Phys. 83 (1985) 1685. G.C. Lynch, D.G. Truhlar and B.C. Garrett, I. Chem. Phys. 90 (1989) 3102. K. Haug, D.W. Schwenke, D.G. Truhiar, Y. Zhang, J.Z.H. Zhang and D.J. Kouri, I. Chem. Phys. 87(1987)1892. Y. Zhang, J.Z.H. Zhang, D.J. Kouri, K. Haug, D.W. Schwenke and D.G. Truhiar, Faraday Discuss. Chem. Soc. 84(1987) 381. G.C. Lynch, P. Halvick, D.G. Truhlar, B.C. Garrett and D. Schwenke, Naturforsch. 44a (1989) 427. B.C. Garrett and D.G. Truhiar, I. Phys. Chem 95 (1991) 10374. [81] See, for examples, D.C. Chatfield, R.S. Friedman, D.W. Schwenke and D.G. Truhiar, J. Phys. Chem. 96 (1992) 2414. U. Manthe, T. Seideman and W.H. Miller, J. Chem. Phys. 99 (1993) 10078, and references therein. [82] A. Gonzalez-Lafont, T.N. Truong and D.G. Truhiar, J. Chem. Phys. 95 (1991) 8875. [83] A. Gonzalez-Lafont, T.N. Truong and D.G. Truhlar, I. Phys. Chem. 95 (1991) 4618. A.A. Viggiano, J. Paschkewitz, R.A. Morris, J.F. Paulson, A. Gonzalez-Lafont and D.G. Truhiar, J. Amer. Chem. Soc. 113 (1991) 9404. Y.P. Liu, D. Lu, A. Gonzalez-Lafont, D.G. Truhlar, B.C. Garrett, J. Amer. Chem. Soc. 115 (1993) 7806.