1476
J . Phys. Chem. 1988, 92, 1476-1488
AlgorRhms and Accuracy Requirements for Computing Reaction Paths by the Method of Steepest Descent Bruce C. Garrett, Michael J. Redmon, Chemical Dynamics Corporation, 9560 Pennsylvania Avenue, Upper Marlboro, Maryland 20772
Rozeanne Steckler, Donald G. Truhlar,* Department of Chemistry and Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 55455
Kim K. Baldridge, David Bartol, Michael W. Schmidt, and Mark S. Gordont Department of Chemistry, North Dakota State University, Fargo, North Dakota 58105-5516 (Received: July 6. 1987; In Final Form: October 6, 1987)
We present several algorithms for computing minimum energy reaction paths by the method of steepest descents of the potential energy and for interpolating reaction-path potentials. The algorithms are suited for use with global analytic potential surfaces or ab initio calculations of energies and gradients. The efficiencies of the various algorithms are demonstrated and compared by applications to several test cases.
1. Introduction
Recent advances in electronic structure theory allow the direct computation of gradients, and sometimes higher derivatives, of potential energy surfaces with only a small increase in computational effort over the calculation of the energy This development has led to increased interest in minimum energy paths computed by the method of steepest descent5-8 and in dynamics approximations based on potentials expanded about reaction path^.^,^'^ In the present article we address the question of the accuracy with which the reaction path must be calculated, for a given level of electronic structure theory, in order to converge the calculation of thermal rate constants based on reaction paths and reaction-path potentials. The simplest way to calculate a reaction path is to start at a saddle point and take successive small steps in the direction of the negative gradient. As the step size is decreased, this leads to a converged minimum energy path (MEP). If the coordinate system is mass-weighted or mass-scaled,s this is called an intrinsic reaction coordinate (IRC). We have formulated variational transition-state theory and multidimensional semiclassical transmission coefficients in terms of the IRC and the expansion of the potential in its vicinity.I3J4 We have found, however, that the gradient-following step size must be very small for the error in the rate constant due to lack of convergence of the IRC to be negligible.13J4J9 In particular, in calculating reaction rates using analytic potential energy surfaces, we have usually tried to converge all calculated rate constants with respect to numerical parameters to l% or better, and we have used mass-weighted step ul/*a. or less (often much less), where u is sizes of about the atomic mass unit and a. the abbreviation for bohr. Such small step sizes are not too expensive when analytic potential energy surfaces are used, but they are generally prohibitive for ab initio calculations. Therefore, we decided to make a systematic study for some prototype systems of the maximum acceptable step size for a given required level of accuracy in the rate constant. In addition we decided to investigate whether more sophisticated methods of calculating the IRC might allow for larger step sizes to be used, and we have tried to devise an efficient general procedure for future applications. The results of these studies are presented here. 'Visiting Fellow, University of Minnesota Supercomputer Institute, 1985-1986.
0022-3654/88/2092-1476$01.50/0
Most previous work on calculating IRC's has used the simplest method, Le., stepping in the direction of the negative gradient, ( 1 ) Pulay, P. Applications of EIectronic Structure Theory; Schaefer. H. F., Ed.; Plenum: New York, 1977; p 153. (2) Pople, J. A,; Krishnan, R. A,; Schlegel, H. B. Int. J. Quantum Chem. Symp. 1979, 13, 225. (3) Morokuw, K.; Kato, S. Potential Energy Surfaces and Dynamics Calculations; Truhlar, D. G., Ed.; Plenum: New York, 1981; p 243. Morokuma, K.; Kato, S.; Kitaura, K.; Obara, S.;Ohta, K.; Hanamura, M. New Horizons of Quantum Chemistry; Reidel: Dordrecht, 1983; p 221. (4) Caw, J. F.; Yamaguchi, Y.; Schaefer, H. F. J . Chem. Phys. 1984, 81, 6395. (5) Shavitt, I. University of Wisconsin Theoretical Chemistry Laboratory Technical Report WIS-AEC-23, Madison, WI, 1959 (unpublished). Shavitt, I.; Stevens, R. M.; Minn, F. L.; Karplus, M. J. Chem. Phys. 1968,48, 2700. (6) Marcus, R. A. J . Chem. Phys. 1968, 49, 2610. (7) Truhlar, D. G. J . Chem. Phys. 1970,53, 2041. Truhlar, D. G.; Kuppermann, A. J. Am. Chem. SOC.1971, 93, 1840. Garrett, B. C.; Truhlar, D. G. J . Phys. Chem. 1979, 83, 1052; J . Phys. Chem. 1983,87, 4553. (8) Fukui, K. J. Phys. Chem. 1970, 74, 4161. Fukui, K. The World of Quantum Chemistry; Daudel, R.; Pullman, B., Eds.; Reidel: Dordrecht, 1974; p 113. Fukui, K.; Kato, S.; Fujimoto, H. J . Am. Chem. SOC.1975, 97, I . (9) Hofacker, G. L.2.Naturforsch., A : Astrophys., Phys., Phys. Chem. 1963, 18A, 607. Hofacker, G . L. Int. J. Quantum Chem. Symp. 1969,3, 33. Hofacker, G. L.; Levine, R. D. Chem. Phys. Lett. 1971, 9, 617; Chem. Phys. Lett. 1972.15, 165. Hofacker, G. L.; RBsch, N . Ber. Bunsenges. Phys. Chem. 1973, 77,661. Hofacker, G. L.; Michel, K. W. Ber. Bunsenges. Phys. Chem. 1974, 78, 174. (10) Marcus, R. A. J . Chem. Phys. 1968,45,4500;Discuss. Faraday SOC. 1968,44,7; J . Chem. Phys. 1968,49,2610,2617. Wu, S.-F.;Marcus, R. A. J. Chem. Phys. 1970, 53, 4026. (11) Light, J. C. Adu. Chem. Phys. 1971, 19, I . (12) Kato, S.; Kaot, H.; Fukui, K. J . Am. Chem. SOC.1977, 99, 684. Yamashita, K.; Yomabe, T.; Fukui, K. Chem. Phys. Lett. 1981, 84, 123. (13) Garrett, B. C.;Truhlar, D. C. J. Phys. Chem. 1979,83, 1079; J. Phys. Chem. 1980,84,682;J . Phys. Chem. 1983,87,4553. Garrett, B. C.;Truhlar, D. G.; Grev, R. S.;Magnuson, A. W. J . Phys. Chem. 1980,84, 1730; J . Phys. Chem. 1983,87,4554;Skcdje, R. T.; Truhlar, D. G.; Garrett, B. C. J. Chem. Phys. 1982, 77, 5955. (14) Isaacson, A. D.; Truhlar, D. G . J . Chem. Phys. 1982, 76, 1380. Truhlar, D. G.; Isaacson, A. D.; Garrett, B. C. Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; Vol. 4, p 65. (15) Garrett, B. C.; Truhlar, D. G. J. Phys. Chem. 1982, 86, 1136. Truhlar, D. G.; Kilpatrick, N . J.; Garrett, B. C. J. Chem. Phys. 1983, 78, 2438. Steckler, R.; Truhlar, D. G.; Garrett, B. C.; Blais, N . C.; Walker, B. J. Chem. Phys. 1984, 81, 5700. (16) Miller, W. H.; Handy, N. C.; Adams, J. E. J . Chem. Phys. 1980, 72, 99. Gray, S. K.; Miller, W. H.; Yamaguchi, Y.; Schaefer, H. F., I11 J . Chem. Phys. 1980, 73, 2733. Miller, W. H. Potential Energy Surfaces and Dynamics Calculations; Truhlar, D. G., Ed.; Plenum: New York, 1981; p 265. Cerjan, C . ;Shi, S.; Miller, W. H. J. Phys. Chem. 1982, 86, 2244. Miller, W. H. J . Phys. Chem. 1983, 87, 3811.
0 1988 American Chemical Society
The Journal of Physical Chemistry, Vol. 92, No. 6, 1988 1477
Reaction Paths by the Method of Steepest Descent which is equivalent to Euler integrationZoof the defining differential equation, but there have been some attempts to use higher order methods. For example, Hase and Duchovic used a Runge-Butta method with a step size of 9.4 X lo4 u1IZa0.l8 Our goal is to make a more systematic survey to see if any higher order algorithms allow appreciably larger step sizes. 2. Theory
We consider a reacting system composed of N atoms, A, B, etc., and we denote its geometry by a mass-scaled N-dimensional vector where p is arbitrarily chosen as is the reduced mass of the reactants and ZA, 3B, ... are the coordinates of A, B, .... Then the IRC is the solution of the differential equation d2/ds = v'
(2)
where s is the mass-scaled reaction coordinate and v' is a unit vector defined as the negative of the normalized gradient of the potential: a = -bV/lbvl (3) where 9 is the N-dimensional gradient operator in mass-scaled coordinates and V is the Born-Oppenheimer potential energy (which equals the total electronic energy plus nuclear-nuclear repulsion). Note: to convert distances along s to mass-weighted values (as quoted in the Introduction), one multiplies by p ' / 2 . For the studies presented here of bimolecular reactions with analytic or potential energy surfaces, we set I.C equal to mOHmH2/mOH3 mcH3mH2/mCHJ, both of which equal 1.8 u. For the studies based on a b initio gradients and Hessians we set p = 1 u. In practice one uses the initial conditions that 3(s=O) is the saddle-point geometry. One then integrates from s = 0 to s = s-, where s- is negative (Le., on the reactant side), and from s = 0 to s = s+, where s+ is positive (Le., on the product side). The range (s-, s+) must be wide enough to converge all dynamics calculations of interest. As mentioned in the Introduction, computing the IRC by stepping in the direction of the negative gradient is equivalent to the Euler methodZo(also called the Euler one-step method in some books) for integrating ( 2 ) . In addition to this simple method we consider several more sophisticated methods, including the fixed-step-size Adams-Moulton method, denoted FAM, the variable-step-size Adams-Moulton method, denoted VAM, and the fixed-stepsize Adams predictor method, denoted FAP. (The FAP method of order n is denoted FAPn.) The Adams-Moulton algorithm used is the same for the FAM and VAM methods; it is a fourth-order predictor-corrector integrator and is described In general an integrator of order p uses information from the current and p previous steps; Le., it is a (p 1)-point integrator. Specification of a method as order p does not imply that the error varies with the step size 6s as (8s)'; in particular, ~ . FAP method of order 0 the error for p = 4 varies as ( 6 ~ ) The is equivalent to the Euler one-step method. We will also be especially interested in the two-point (first-order) FAP method which is discussed below in terms of the quantities available from ab initio gradient calculations. This method utilizes the gradient from the previous integration point as well as the geometry and gradient from the current integration point. One might consider using the geometry at the previous integration point plus the geometry and gradient a t the current point. This yields the quadratic version of the point-slope formula, which is known to be unstable.22 It is also unstable to use both the geometry and
+
(17) Kato, S.;Morokuma, K. J . Chem. Phys. 1980, 73, 3900. (18) Hase, W. L.; Duchovic, R. J. J . Chem. Phys. 1985,83,3448. Mondro, s. L.; Linde, s. V.; Hase, W. L. J . Chem. Phys. 1986, 84, 3783. (19) Truhlar, D. G.; Brown, F. B.; Steckler, R.; Isaacson, A. D. The Theory of Chemical Reaction Dynamics; Clary, D. C., Ed.; Reidel: Dordrecht, 1986; p 285. (20) See, e.g., Carnaham, B.; Luther, H. A.; Wilkes, J. 0.Applied Numerical Methods; Wiley: New York, 1969; pp 344ff. (21) Miller, W. H.; George, T. F. J . Chem. Phys. 1972, 56, 5668.
gradient at both points. We also consider a method employing a specialized corrector to improve the convergence of the Euler algorithm. This method, which is called the Euler stabilization (ES) method, is a refined version of the algorithm of Schmidt, Gordon, and D ~ p u i s ?which ~ is itself a refined version of an earlier algorithm proposed and used by Ishida et al.24 Having obtained the IRC, we calculated transmission coefficients and rate constants by methods described elsewhereI4 and in section 3.2. We consider two approximations to the transmission ~ corresponding ~ ~ , to the minimum-energy-path coefficient: K semiclassical adiabatic ground-state a p p r o x i m a t i ~ n , and ' ~ ~ ~K ~ ~ ~ corresponding to the small-curvature adiabatic ground-state approximation?6 We note that involves an s-dependent effective reduced mass that depends on the reaction-path curvature, and K~~~ is obtained by replacing this with the constant p. We consider three approximations to the rate constant: kCVT, corresponding to canonical variational t h e ~ r y ; ' ~ . ' kCVTIMEP, ~,'~ defined byI3 kCVT/MEP = ,MEPkCVT (4) and kCVTfSC, defined byz6 kCVT/SC = ,sCkCVT
(5)
(Note: the transmission coefficients are sometimes called K ~ and K ~ ~ in other ~ articles ~ because / ~ they ~ / must be consistent with the CVT version of variational transition state theoryI3 and because they are semiclassical (S) calculations based on the adiabatic ground-state (AG) potential curve, but here we omit the CVT, S , and AG since no confusion should result in the present article.) To calculate these transmission coefficients, we perform generalized transition-state normal-mode analyses at a grid of points along the reaction path, and all subsequent calculations are performed by interpolation off this grid.I4 These analysis points are evenly spaced with separation As, where As is larger than the step size 6s used to converge the IRC calculation. To distinguish As from 6s, we call 6s the reaction-path step size and As the Hessian grid size. A converged calculation is obtained by decreasing both the Hessian grid size and the reaction-path step size.
3. Algorithms 3.1. Calculating the Reaction Path. 3.1 .I. Euler and Adams Predictor Methods. In the Euler method,z0 discussed above, for calculating the reaction path only the geometry and gradient at the current point on the reaction path are used to find the next point. The higher order FAP methods and the FAM and VAM methods also use the knowledge of the gradient at other previous points on the reaction path to predict the next point. All three of these methods use an Adams-like predictor algorithm to integrate eq 2. The Adams-like predictor method of order p uses the current geometry and gradient and the previous p gradients to calculate the next geometry byz1
where 3, and 6,are the geometry and negative of the normalized gradient at point n on the MEP, the coefficients C, are defined by C,,, = fs*'ds (s - s,)(s J
-)'s,
,
...(s - s,+~-,)
(7)
S"
(22) Davis, P. J.; Polonsky, I. Handbook of Mathematical Functions: Abramaowitz, M.; Stegun, I. A., Eds.; National Bureau of Standards: Washington, DC, 1965; p 896. (23) Schmidt, M. W.; Gordon, M. S.; Dupuis, M. J. Am. Chem. Soc. 1985, 107,2585. (24) Ishida, K.; Morokuma, K.: Komornicki, A. J. Chem. Phys. 1977.66, 2153. (25) Truhlar, D. G.; Kuppermann, A. J . Chem. Phys. 1972, 56, 2232. (26) Skodje,R. T.; Truhlar, D. G.; Garrett, B. C. J . Phys. Chem. 1981, 85, 3019. Truhlar, D. G.; Isaacson, A. D.; Skodje, R. T.; Garrett, B. C. J . Phys. Chem. 1982, 86, 2252.
~~
~
Garrett et al.
1478 The Journal of Physical Chemistry, Vol. 92, No. 6, 1988 TABLE I: Effect of Interpolation Method for p&) Reaction R1 at 200 K" no. of points used
in Lagrangian interpoln 2 5 7 742 742 742 741 74 1 741 740 740 740 735 736 736 726 729 729 713 717 717 698 703 706 677 685 687 657 664 659 622 636 635 599 60 1 598 456 450
As, a. 0.005 0.010 0.020 0.030 0.040 0.050 0.060 0.070 0.080 0.090 0.100 0.300
on
K~
TABLE II: Effect of Interpolation Method for p&) Reaction R2 at 200 K" no. of points used
for
in Lagrangian interpoln
spline
interpoln
As, a0
2
742 74 1 740 736 729 718 705 687 666 637 597 496
0.01 0.02 0.03 0.04 0.05 0.06 0.08 0.10 0.15 0.20
53.0 52.8 52.7 52.0 51.4 50.5 48.3 44.7 33.5 27.5
and the quantity Jim) is the divided difference function of 5,,,, , C, ..., zj,,-* defined by = 5(s,)
(8)
and
a i m ) = (Jim-,) - JSI))/(s,
- s,-m)
+
- s,)5(s,)
53.0 52.8 52.7 51.9 51.3 50.5 48.1 44.3 33.4 28.2
given by Schmidt, Gordon, and Dupuk2j The methods are similar to a predictor-corrector method in that they involve two steps. First a step of length 6 , is specified by the Euler method as
2iy1 = 2,
(10)
The higher order FAP methods also have very simple geometrical interpretations; e.g., the FAP method of order 1 is equivalent to fitting a quadratic to the current geometry and gradient and the previous gradient.
+ 6,5,
(13)
where 6 , 6s. Then a corrector (stabilization) step is taken to a point at the minimum of a quadratic fit to the ptential along a line through parallel to a "bisector" vector d,+, defined by
(9)
Notice that the coefficients in the predictor are independent of the order p . The FAP method consists of only using the Adams-like predictor just described. As stated above, the FAP method of order 0 is equivalent to the Euler method, which steps along a straight line fitted to the current geometry and gradient:
Ziy, = 2,
spline interpoln
7 53.0 52.8 52.7 51.9 51.3 50.5 49.1 44.3 33.5 27.0
" Euler integrator, 6s = 0.001 a,.
"Euler integrator, 6s = 0.001 a,,.
620)
5 53.0 52.8 52.7 51.9 51.3 50.5 48.0 44.4 33.3 27.1
on K~ for
;in+l= (5, - i @ I ) / I C ,
-
$y,l
(14)
(Note: in ref 23 and 24 this vector is stated to bisect the two successive gradients v, and viyl. While strictly speaking this vector is perpendicular to the bisector, the statement is intended to convey bisection in geometrical construction in which the gradients are considered to be placed head t,o tail, as in Figure 1 of ref 24. It is in this sense that we refer to d,, as the "bisector".) The explicit algorithm for the corrector step is
Z,+, = Ziy, + &+, where the step length along the unit vector the parabolic fit to be
f = f(62/2)[1
(15)
&+,
is predicted by
* @vn+l/~2~~+l)l-1
(16)
where 3.1.2. A d a w M o u l t o n Methods. The algorithms for the other two classes of methods, FAM and VAM, also include a corrector of the Moulton type. The Moulton corrector of order p is defined by Z,+, = X'lp+', C,(u7p,'I- 5, - h,@,') - h,(h, + h,,)di2)- ... h,(h, h,-,)...(h, + h,, + ... h,+l,)d,@))[h,(h, + h p l )...(h, h,, ... + h,+,-,)]-l (12)
+
+
+
+
+
where C, is the coefficient defined by eq 7, h, is the distance is the gradient evaluated a t Zi!,. between s,+~ and s, and In the FAP and FAM methods the distance h, is kept constant for all n (i.e,, the step size is fixed). In the VAM method the step size is varied according to an algorithm2I that limits the truncation error to a specified value e. Whenp 2 1, the program must start the IRC calculation with lower orders. In particular the first step uses a p = 0 algorithm, and the order is increased by 1 at each step until it reaches the desired value. This appears to be more efficient than using a Runge-Kutta starter. In the VAM method it is possible that a reaction-path step takes one past the next scheduled Hessian grid point. In this case we still take the large reaction-path step and we perform a normal-mode analysis at the end of this step. The distance As to the next Hessian grid point is now measured from this new irregular grid point. The order p is fixed at 4 for all FAM and VAM calculations. 3.1.3. The Euler Stabilization Method. We use two versions of the Euler stabilization method: one (ES1) is a modification of and the other (ES2) is the original of the stabilization method
AV,+, = V(Ziyl)- V(Zi'$,
+ fi2;,,+,)
(17)
and
G,+,= ;i"+,.VV(2iyl) (18) In this equation, G,+,9 the component of the analytic derivative of the potential along d,, at 2iyl, AVn+l/62 is a finite-difference approximation to this derivative, 62 is the step size parameter used in the finite-difference approximation, and the signs are both chosen opposite to the sign of G,+l so that the step in (15) is in the direction of decreasing potential. Normally both signs are negative, giving a positive [. Negative ['s are very rare and are almost always near zero. become parallel, the "bisectorn In the event that 5, and is indeterminate. In this case, the stabilization step is unnecessary anyway (the IRC is linear) and can simply be skipped, taking as Z,+,. In practice, the stabilization may be skipped whenever exceeds some angle, say 175". the "elbow" between 5, and Skipping the stabilization step is frequently necessary when small step sizes 61 are being used. In the calculations presented here € the criterion for skipping the stabilization step was 15, lo-* au. Instead of being nearly linear, an IRC may be rapidly curving. This situation is associated with fairly sizeable values of [ such that the minimum energy is located rather far from x$!, along the "bisector". If f lies outside the range < 262,the stabilization step predicted by eq 16 is deemed unsatisfactory. It should be mentioned that if reasonably chosen values of 6, are being used, the failure of eq 16 is rare. The handling of such a failure depends on the economics of analytic versus ab initio surfaces, which dictate two different strategies.
Reaction Paths by the Method of Steepest Descent
The Journal of Physical Chemistry, Vol. 92, No. 6. 1988 1419
TABLE 111: Effect of Changing Reaction-Path Step Size for Euler and ESl Methods for Reaction R1
T = 200 K as, a0
N.
SW,
a0
kCVT cm3 molecule-l
T = 300 K
KMEP
s-I
KSC
kCVT, cm3 molecule-'
s-I
KMEP
KSC
1.18 (-15) 1.18 (-15) 1.18 (-15) 1.18 (-15) 1.17 (-15) 1.18 (-15) 1.35 (-15) 1.49 (-15) 1.85 (-15)
4.5 4.5 4.5 4.5 4.5 4.5 12 15 11
20 20 20 20 20 17 b b b
1.20 (-15) 1.20 (-15) 1.20 (-15) 1.20 (-15) 1.19 (-15) 1.17 (-15) 1.17 (-15) 1.31 (-15) 1.86 (-15)
4.5 4.5 4.5 4.5 4.5 4.6 4.6 6.1 126
20 20 21 379 93 21 24 b b
Euler (FAPO) 0.0010 0.0025 0.0050 0.0100 0.0125 0.0250 0.0500 0.0600 0.1000
1500 600 300 150 120 60 30 25 15
0.0010 0.0025 0.0050 0.0100 0.0250 0.0500 0.0600 0.1000 0.1500
4494 1800 912 504 254 116 112 82 84
a
a a a
a
-1 .o -0.25 -0.18 -0.20
9.98 (-18) 9.97 (-18) 9.95 (-18) 9.91 (-18) 9.89 (-18) 9.88 (-18) 1.23 (-17) 1.42 (-17) 2.05 (-17)
50 50 50 50 50 51 325 355 160
717 715 71 1 704 695 616 b b b
Euler Stabilization (ESl)c a
a a a a a
-0.89 -0.73 -0.53
1.02 (-17) 1.02 (-17) 1.02 (-17) 1.02 (-17) 1.01 (-17) 9.87 (-18) 9.79 (-18) 1.16 (-17) 2.08 (-17)
49 49 49 49 50 53 58 128 126
713 713 736 40800 10500 868 1220 c c
'In these cases the error in is less than 1 kcal/mol for the whole (-1.2 a. I s 5 +0.3 ao) interval studied. these cases the SC calculation was not completed because of too small an overlap of successive eigenvectors on the Hessian grid. eFor these runs we set 6, = 6s and 6, = & / l o . TABLE IV: Results Obtained by Fixed-Step Adams Predictor Method (FAP) for Reaction R1 as a Function of Order p
T = 200 K
bs, a,
NB
0.0010 0.0025 0.0050 0.00625 0.0100 0.0 125 0.0250 0.0500
1500 600 300 240 150 120 60 30
0.0010 0.0025 0.0050 0.00625 0.0100 0.0125 0.0250 0.0500 0.0600
1500 600 300 240 150 120 60 30 25
0.0010 0.0025 0.0050
1500 600 300
0.0005 0.00 10
3000 1500
sM, a. a a
-1.20 -1.15 -1.05 -0.90 -0.25 -0.15 a
-1.20 -1.05 -0.95 -0.60 -0.15 -0.15 -0.15 -0.12 a
-1.05 -0.05 a
-0.05
kCVT cm3 molecule-I
s-I
T = 300 K
KMEP
KSC
kCVT, cm3 molecule-' s-l
KMEP
Ksc
9.99 (-18) 9.99 (-18) 9.99 (-18) 9.99 (-18) 9.98 (-18) 1.00 (-17) 1.30 (-17) 1.29 (-17)
p=l 50 50 50 50 51 55 123 5 26
719 719 b b 802 826 5270 b
1.18 (-15) 1.18 (-15) 1.18 (-15) 1.18 (-15) 1.18 (-15) 1.19 (-15) 1.39 (-15) 1.38 (-15)
4.5 4.5 4.5 4.5 4.5 4.5 8.8 26
20 20 b b 22 20 97 b
9.99 (-18) 9.99 (-18) 1.03 (-17) 1.18 (-17) 1.12 (-17) 1.10 (-17) 1.22 (-17) 9.78 (-18) 1.20 (-17)
p = 2 50 50 49 43 162 470 208 216 125
719 b b 1050 17300 b 2730 b b
1.18 (-15) 1.18 (-15) 1.20 (-15) 1.18 (-15) 1.25 (-15) 1.25 (-15) 1.33 (-15) 1.12 (-15) 1.32 (-15)
4.5 4.5 4.4 4.0 7.9 16 17 19 5.5
20 b b 25 b 178 101 b b
9.99 (-18) 6.27 (-18) 1.81 (-17)
p = 3 49 80 110
716 1270 4058
1.18 (-15) 8.79 (-16) 1.69 (-15)
4.5 6.5 6.2
20 30 54
9.99 (-18) 6.21 (-18)
p = 4 44 64
682 1476
1.18 (-15) 8.69 (-16)
4.2 5.7
20 34
'In these cases the error in Lf is less than 1 kcal/mol for the whole (-1.2 a. I s 5 C0.3 ao) interval studied. was not completed because of too small an overlap of successive eigenvectors on the Hessian grid.
If the energy and gradient may be obtained at little expense, as in the analytic surface case, the best procedure is to restart the step. The Euler step size 6, is halved, and the step is repeated. This procedure is iterated until a step is accepted, or after j step halvings, 6 , / 2 j I6> In the latter case an Euler step is taken with the reduced step size 6 , / 2 j and is not stabilized. At the beginning of each new step 61 is reset to the original 6s. This method will be referred to as ES1. If the energy and gradient are obtained by ab initio calculations, the abandonment of the calculated point a t x$!, may be undesirable since the calculation at that point may have required a large amount of computer time. However, the quadratic fit to
these cases the SC calculation
the energy along the "bisector" probably is unreliable (since the fitted data are the energy and gradient projection at x'ip!l and an energy at step along the "bisector"). In this case, computation of the energy at 2S2,4s2,Sa2, ... is carried out, until three successive energies bracket the minimum. One, two, or three additional energy evaluations are more economical than recomputation of the step with a smaller a1, which costs two energies and a gradient. The minimum of a quadratic fit to the last three energies is taken as the next IRC point, This method is referred to as ES2. 3.1.4. Initiating an IRC. None of the procedures described above can be used for the first step away from a transition state, for these structures have a zero gradient of the energy. Instead
Garrett et al.
1480 The Journal of Physical Chemistry, Vol. 92, No. 6, 1988 TABLE V Results Obtained by FAM and VAM Methods for Reaction R1 at 200 K no. of e 6s, a0 integrn steps N B sMt FAM 3000 b a 0.0010 1500 b 600 1200 a 0.0025 300 600 b a 0.0050 b a 0.00625 240 480 150 300 -1.20 a 0.0100 120 240 -1.10 a 0.0125 60 120 -0.65 a 0.0250 30 60 -0.25 a 0.0500
10-6 0.00001 0.0001 0.0010 0.0025 0.0050 0.0 100 0.015 0.020 0.050
24 1 221 208 187 163 139 117 96 94 70
0.0005-0.028 0.001-0.032 0.001-0.040 0.001-0.040 0.001-0.038 0.001-0.035 0.001-0.036 0.001-0.042 0.001-0.031 0.001-0.037
520 479 453 44 1 355 294 240 213 192 142
VAM~ b b b b b b b b -1.14 -0.71
kcvT,
cm3 molecule-’ s-l
KMEP
9.99 (-18) 9.99 (-18) 9.99 (-18) 9.99 (-18) 9.98 (-18) 1.00 (-17) 1.01 (-17) 1.91 (-17)
50 50 50 50 50 51 45 27
719 719 723 721
9.99 (-18) 9.98 (-18) 9.99 (-18) 9.99 (-18) 9.99 (-18) 1.05 (-17) 9.51 (-18) 9.44 (-18) 1.25 (-17) 1.53 (-17)
50 50 50 50 50 48 46 38 23 150
720 722 787 2090 6630 33050
Ksc
C C C
c
C
2.8 (5) 4.9 (5) C
is less than 1 kcal/mol for the whole interval (-1.2 a, 5 s 5 0.3 ao). CInthese cases the SC these cases the error in “Does not apply. calculation was not completed because of too small an overlap of successive eigenvectors on the Hessian grid. “In these runs the initial stepsize bo was lo-’ a,, and the 6s column gives the range of 6s that resulted from a given e. TABLE VI: Number of Gradient Evaluations Necessary in Integration Procedure To Obtain 15%Accuracy for Reaction R1 T = 300 K T = 200 K integrn
method Euler (FAPO) ES 1 FAP 1 FAP2 FAP3 FAP4 FAM VAM
kcw 5 6 8 17 68 142 14 24
K~~~
59 112 120 240 1250 1990 112 228
K=
59 912 600 1500 1320 2730 480 453
kcvT 4 6 9 4 58 127 10 24
K~~~
57 101 120 267 1200 2030 62 184
Sc 60 912 600 1500 1230 2680 480 451
the IRC is initiated by following the unique normal mode whose frequency is imaginary. In the harmonic approximation, the energy lowering AE for a given step 6s is approximately AE = (1/2)k6s2
approximation on the Hessian grid. In order to compute this derivative, we must order the eigenvectors of each symmetry in order of decreasing eigenvalue. This is done by taking overlaps of eigenvectors at successive points on the Hessian grid to maintain the correspondence of modes at these successive points. If the overlaps of successive eigenvectors are too small, the run terminates. This may result when mode-mode couplings are large and one or both step sizes is not small enough. We include tunneling in the total energy range from tRG,which is the ground-state energy of reactants, up to VAG,which is the maximum of the vibrationally adiabatic ground-state potential curve The small-curvatureapproximation semiclassical tunneling probability at energy E is computed as P c ( E ) = {l
+ exp[2~9(E)])-~
(20)
where the imaginary-action integral is
(19)
where the negative force constant k is obtainable from the magnitude of the imaginary frequency. This formula can be used in two ways, to take a known step 6s no matter how much the energy is lowered or to take whatever 6s produces the desired energy lowering. We use the former strategy for our calculations with analytic surfaces, and the latter choice is used in the ab initio calculations. 3.2. Calculating the Reaction-Path Curvature and the Effective Reduced Mass. All the methods employed to calculate rate constants and transmission coefficients are the same as those described in detail elsewherei4 except for the calculation of the curvature component^'^.'^ BmF(s)and the effective reduced massz6 kefds). These are not needed for kCVTand K ~ but ~if they ~ are , not computed carefully they became the quantities that most limit the numerical precision of #c. The f ~ r m u l for a ~the ~ curvature ~~~ , ’ derivative ~ components involves dZ/ds. In previous ~ o r k ’ ~the of the gradient was,computed by a two-point backward difference formula. In the present work we replaced this by the analytic derivative a t the current point of a three-point quadratic fit to 0’ at current, previous, and next points on the reaction-path grid. For fixed-step-size integrations this is equivalent to using the central difference formula with step size 6s. Having obtained the curvature components, we calculated the effective reduced mass for tunneling at every point on the Hessian grid. The effective reduced mass depends on Em&), tm(s),and dt,/ds, where t , is the distance along normal mode m to the ground-state vibrational turning point on the concave side of the MEP. The derivative dt,/ds is computed by a finite difference
J
s
