J. Phys. Chem. 1992,96, 531-537 and Stan Miklavic for making data available to us prior to publication. We also thank Frans Leermakers, Gregor Cevc, Christiane Helm, Wolfgang Knoll, Clay Radke, and Joe Zasadzinski for helpful comments and the NSF for supporting J.N.I. under Grant CTS-9015537.
Glossary force law
decay length
variation of force For pressure P with distance D between two surfaces. Since many force laws have a near expoor P = P$-D/A, nential distance dependence, F = these are usually plotted on a log linear (semi-log) graph where the data falls on a straight line whose inverse slope equals the decay length, A. for any exponential function of the type F = Foe-D/Athe decay length X is a measure of the rate of decrease of F with distance D F falls by a factor of 2.7 as D changes by A. The smaller the value of X the more rapid is the
a X
x, P e,
6 Yi
Q
531
decay and thus the shorter is the range of the force. area (m-’), Hamaker constant (J) distance (m) force (N), pressure (N bilayer compressibility modulus (J m-’), bilayer bending modulus (J) thickness of brush layer or headgroup layer thickness (m) density of molecular groups per unit area (m-’) mean distance between surface binding sites or between headgroups (m) protrusion energy per unit length (J m-I) decay length of the density of molecular groups protruding out from a surface decay length of the force number density (m3) electronic charge (C), surface charge density (C m-’) molecular diameter or width (m) interfacial energy (J m-’)
ARTICLES Global Potential Energy Surfaces from Limited ab Initio Data Alan D. Isaacson Department of Chemistry, Miami University, Oxford, Ohio 45056 (Received: May 1 , 1991)
-
Two new potential energy surfaces for modeling the dynamics of the reaction OH + H2 H 2 0 + H are presented. While the form of these surfaces is relatively simple, varying the surface parameters along the reaction path allows these surfaces to reproduce the ab initio information of Walch and Dunning as well as either the bamer shape computed by Dunning, Harding, and Kraka or the experimental rate constants of Ravishankara et al. By comparing the results obtained from these surfaces and an earlier one by Schatz and Elgersma, we are able to study the sensitivity of the predicted dynamics on both the form of the surface and the data to which the surface is fit. In particular, the results indicate that both the barrier shape and the degree of reaction path curvature can strongly influence the predicted rate constants at lower temperatures.
1. Introduction Any detailed theoretical study of reaction dynamics relies on a knowledge of the potential energy surface (PES) for the reacting system. However, since the PES for an N-atom system is a function of 3 N - 6 variables, the accurate calculation and representation of a PES poses a major challenge in dynamics studies even for simple systems. As a result, a common method for obtaining a PES of a reacting system involves the inversion of a limited amount of theoretical and/or experimental data on reactant, product, and transition-state properties as well as experimental rate constants and equilibrium constants. Such an inversion is always based on an assumed functional form for the PES, and one important question in this area concerns the sensitivity of the predicted results for the reaction dynamics on the choice of the functional form used for modeling the PES. Indeed, we present below a case in which different functional forms fit to the same set of data lead to profoundly different predictions for the rate constant. A related question involves the identification of the key properties of the PES that should be considered most carefully in the choice of a functional form. One such property is the barrier height, which strongly affects the reaction rate by determining the number of collisions that have sufficient energy to react. Another is the barrier shape or width, which controls the degree of reaction tunneling. As discussed below, the reaction path curvature arising from a particular functional form markedly
affects the results in some cases. Furthermore, we show below that requiring the PES have a particular bamer shape can strongly affect the degree of reaction path curvature obtained from a given functional form. A third important but more practical question deals with determining functional forms and fitting procedures that can provide efficient and accurate fits to a given set of property data. For example, we show below that the variation of surface parameters along the reaction path’-5 allows for an accurate and relatively systematic approach to the inversion of property data to yield a PES. In this paper we address the above questions within the context of the PES for the reaction OH + Hz HzO H (R1) +
+
which is important in combu~tion,~~’ and we present a new PES which better fits the available data than do previous (1) Blais, N. C.; Truhlar, D. G. J . Chem. Phys. 1974,61, 4186;1976,65, 3803E. (2)Truhlar, D.G.; Garrett, B. C.; Blais, N. C. J . Chem. Phys. q984,80, 232. (3)Steckler, R.;Truhlar, D. G.; Garrett, B. C. J . Chem. Phys. 1985.83, 2870. (4)Blais, N. C.;Truhlar, D. G. J. Chem. Phys. 1985,83, 5546. (5) Joseph, T.; Steckler, R.; Truhlar, D. G. J . Chem. Phys. 1987,87,7036. (6)Glassman, I. Combustion; Academic: New York, 1977. (7)Creighton, J. R.J . Phys. Chem. 1977,81, 2520.
0022-3654/92/2096-53 1%03.00/0 0 1992 American Chemical Society
Isaacson
532 The Journal of Physical Chemistry, Vol. 96, No. 2, 1992
TABLE I: Renetant, Product, pad % M e Point Roperties of OH + H2Potential Eaergy surface^"^^ ab initioC SEd no. 3c no. 4c
Figure 1. OH + H 2 saddle point geometry (drawn to scale) and definitions of the internal coordinates.
for this reaction. Based on the large-scale polarization configuration interaction (POL-CI) calculations of Walch and Dunning,’O Schatz and Elgersma8 obtained an analytic PES (herein called no. 1 or SE) for reaction R1 and performed quasiclassical trajectory calculations to determine product vibrational distributions. Isaacson and Truhlar” also used this surface to compute rate constants for reaction R1 and its isotopic analogue OH D2 HDO + D (R2)
+
2. Available Data This section discusses the theoretical and experimental data on reactions R1 and R2 to which the surfaces presented in this paper are fit and/or to which the predictions obtained from these surfaces are compared. The POL-CI calculations of Walch and Dunning’O lead to the values for the reactant, product, and saddle point properties given in the first column of Table I. (See Figure 1 for the definitions of the geometric quantities for this system.) Although more reliable values for the reactant and product properties are available,14J5for the sake of consistency the present fits are based on the Walch-Dunning (WD) values for these quantities. Owing to the limited nature of the a b initio data, the new surfaces include only harmonic terms in the bending internal coordinates, and the dynamical calculations discussed below treat the bound vibrational degrees of freedom only in the harmonic approximation. (8) Schatz, G. C.; Elgersma, H. Chem. Phys. Lett. 1980, 73, 21. (9) Rashed, 0.;Brown, N. J. J. Chem. Phys. 1985,82, 5506. (10) Walch, S.P.; Dunning, T. H., Jr. J. Chem. Phys. 1980, 72, 1303. (11) Isaacson, A. D.; Truhlar, D. G. J. Chem. Phys. 1982, 76, 1380. (12) Dunning, T. H., Jr.; Harding, L. B.; Kraka, E. In Supercomputer Algorithms for Reactivity, Dynamics, and Kinetics of Small Molecules; Lagan&, A., Ed.; Kluwer: Dordrecht, The Netherlands, 1989; pp 57-71. (13) Ravishankara, A. R.; Nicovich, J. M.;Thompson, R. L.; Tully, F. P. J . Phys. Chem. 1981,85, 2498. (14) Hcrzberg, G. Molecular Spectra and Molecular Structure. I . Spectra of Diatomic Molecules; Van Nostrand: Princeton, NJ, 1950. (15) Bartlett, R. J.; Shavitt, I.; Purvis, G. D. J. Chem. Phys. 1979,71,281.
1.429 426 1 95.1
1.428 4271 95.1
1.428 4271 95.1
1.864 3637 93.0
OH 1.863 3624 93.0
1.864 3637 93.0
1.864 3637 93.0
1.84268 104.168 40328 16128 402P -16.13
H20 1.8081’ 104.62‘ 3864’ 1687’ 3975’ -15.20
1A426 104.16 4032 1612 4020 -16.73
1.8426 104.16 4032 1612 4020 -16.73
1.857 2.523 1.631 91.8 14.9 0.0 3368 1945 1248 686 440 16553 6.19
1.852 2.521 1.631 97.9 15.0 0.0 3368 1945 1248 686 440 16551’ 6.19
Saddle Point
+
using variational transition state theory and adiabatic transmission coefficients. However, this surface does not reproduce accurately the a b initio data to which it was fit. In addition, as discussed below, harmonic canonical variational transition state theory rate constants including small-curvature semiclassical adiabatic ground-state transmission coefficients (CVT/SCSAG) computed from this surface are too high at low temperatures. Spurious wells in the asymptotic reactant region of this surface were removed by the modifications of Rashed and Brown? who used their modified surface (herein called no. 2) to carry out a quasiclassical trajectory study of the effects of the initial energy distribution on reactivity for reaction R1. We present below two new surfaces that accurately fit most of the Walch-Dunning ab initio data.I0 Surface 3 also reproduces the positive energy portion of the ab initio barrier shape determined by Dunning, Harding, and Kraka,12 while surface 4 yields CVT/SCSAG rate constants for reaction R1 that agree with the measured values of Ravishankara et al.I3 In section 2 we summarize the data to which the new surfaces have been fit, and in section 3 we present the form of the new surfaces and the fitting procedures. Section 4 describes the calculations of rate constants employing the new surfaces, and section 5 discusses the results. Section 6 presents a summary.
H2
1.428 4277 95.1
1.863 2.522 1.619 91.6 15.0 0.0 336gh 1945h 124gh 686h 440h 1655ih 6.19
1.856 2.308 1.640 116.7 -15.8 0.0 3545 1921 857 572 830 15261‘ 6.09
Frequencies in cm-I, distances in bohr, angles in degrees, and energies in kcal/mol relative to the reactants of reaction R1. bSee Figure 1 for the definitions of the geometrical quantities. Walch and Dunning, ref 10. dSchatz and Elgersma, ref 8. cPresentwork. /Energy of the products of reaction R1, i.e., the reaction endoergicity. 8Obtained by a least-squares fit of data in ref 10 to the potential function ‘ / ~ ~ R R-( %I2 R I + ‘ / 2 k ~ ~ ( -R Re)2 2 + l/2ksa(@+ ‘/~~R&RI,&)(e - 0,) + ‘/2k~#(R2 - R,)(B - 8,) followed by a harmonic analysis. *As reported in ref 32. ‘Fit to the ab initio data on H20in ref 15. j Out-of-plane bending angle. In addition to the data in ref 10, further a b initio studies12 indicate that the shape of the potential barrier (i.e., the variation in the potential energy along the reaction coordinate) closely resembles that of the SE surface for positive energies with respect to the reactants.16 Thus, surface 3 is simultaneously fit to the data of ref 10 and to the positive energy portion of the SE barrier shape. On the other hand, Ravishankara et al.13 published rate constants for reactions R1 and R 2 over the temperature range 250-1050 K. Surface 4 is thus also fit to the data of ref 10 but is additionally required to give CVT/SCSAG rate constants for reaction R1 that agree with the experimental values” within experimental error over the temperature range 298-1000 K. The measured kinetic isotope effect, Le., the ratio of the rate constant for reaction R1 to that for reaction R2, is then used as a measure of the reliability of the various surfaces discussed herein. 3. New Potential Energy Surfaces The form of the present surfaces is relatively simple; it consists of two valence bond three-body terms plus harmonic terms involving the bending and torsional angles:
where the coordinates are labeled as in Figure 1. The function V3(Rk,Rm,R,) in eq 1 is a three-body potential obtained from a London eq~ation:~’.‘~ V3(Rk,Rm,Rn) =
Qk
l[(Jk
(16) Dunning,
+ Qm + Qn - Jm12 + (Jm - Jn)’ + ( J k - Jn)21/2)”2 (2)
T. H., Jr. Private communication.
The Journal of Physical Chemistry, Vol. 96, No. 2, 1992 533
Global Potential Energy Surfaces
and
TABLE Ik Fit Constants for the New Surfaces@ no. 3 no. 4
F5=
*
D " = DI3) 01
bl c1 11
RP a2
b2 Cl
a3 b3 c3
a4 b4 c4
Rt OR
bR Ce
fe ge
a*
-f"-
Pa
R* W
Y
R' fR
PR
P
z
R*
1.428 1.072 80 0.151 548 0.102083 017 0.073 339 295 0.045 852978 1.265 37 0.057 93 14.498 7057 3.673 792 68 1.842 63 1.821 26 0.042 74 5.332 637 69 0.148201 0.015 003 4.505 0.058929417 0.053 625 769 3.450 2.22 2.064 9137 -0.494 1174 1.432 688 29 0.158 7 3.505 438 37 1.321 70694 0.018 OOO 61 1 0.20 4.107 26472 0.7709 40.0 3.50 0.0045 0.10 4.107265 1.oo 5.462 56 3.20
q = a2 + b2(tanh [c2(T - R,)] + l}/2
1.428 1.072 80 0.151 548 0.064 207 223 0.077891 618 0.050 221 465 1.265 37 0.057 93 15.358 8191 4.061 522 10 1.842 63 1.821 26 0.042 74 5.332637 69 0.148 201 0.015003 4.505 0.046 381 036 0.039423881 3.450 2.14 2.064 9137 -0.494 1174 1.432688 29 0.1587 3.586 778 18 0.838 779 182 0.028 826 226 0.20 4.107 264 72 0.7709 40.0 3.45 0.004 831 638 0.40 4.107265 1.oo 5.46246 3.20
where T = (R, + R2)/2 and R, is the WD value for the equilibrium 0-H separation in H20. The use of T as the independent variable in these equations is needed to ensure that the H 2 0 vibrations have the proper symmetries. Equation 7 allows for a smooth transition from the vibrational frequency in OH to those of H20,with somewhat lower frequency values in the saddle point region, consistent with the WD values in Table I. Thus, the value of a, was obtained from the WD value of we for OH,and the value of bl (together with the values of the c3andf, constants discussed below) was obtained from the WD frequencies for H 2 0 given in Table I; the determination of the constants cI and tl is discussed below. Equation 8 similarly allows for a smooth transition from the equilibrium 0-H distance in OH to that in H 2 0 the values of the constants a2,b2,and c2were obtained by requiring that eq 8 reproduce the WD values for R, in OH, in H20, and at the saddle point. For i = 1 and 2, we also take
D{') = a3
and the exoergicity of reaction R1 (Le., the energy of the products relative to that of the reactants), while c3was obtained (together with the constants 6 , andfe) from the WD frequencies for H 2 0 given in Table I. The advantage of such a form for D\') and 05') is that since the derivative of eq 9 with respect to either R I or R2 is zero for T = R,, the calculated equilibrium geometry of H 2 0 is not affected. Lastly, we take
Di3) = a4 + b4(exp {-[4(R2 - R,)/c4I4) - 1)
V, = D\')(R2 =
Qi(Ri) = ['EdRi) + 'Ei(Ri)l/2
(3)
- 3Ei(Ri)]/2
(4)
where 'Ei and 'Ei are singlet (attractive) and triplet (repulsive) interactions between the pair of atoms defining Ri. These interactions are in turn approximated by Morse and anti-Morse functions: 'E,(Ri) = D{')(exp[-2ai(Rj - R 3 ] - 2 exp[-ai(Ri - e ) ] ) ( 5 )
+
3Ei(Ri) = Di3){exp[-2ai(Ri- e ) ] 2 exp[-ai(Ri - Rl)]} (6) where Dj'),Dj3),ai, and are adjustable parameters. Table I1 lists the values of all constants in the PES. The values of the parameters q,ai, and Dj') for i = 5 and 6, which are taken to be constant, were obtained from the WD values of Re, w , and D for H . The values of the parameters Os'),D\3) = Di37, and DF) = DP) are also taken to be constant; their determination is discussed below. The remaining parameters in eqs 5 and 6 are assumed to vary with geometry, Le., along the reaction path. This strategy conveniently builds into the basic London equation the flexibility to make localized changes in a PES,as has already been shown in several studies.'-5 In particular, for i = 1,2, and 3 (i.e., the 0-H interactions), we take
(7)
a)
+ Db')
(11)
Finally, the term Vknd in eq 1 is given by
where Qi and Ji are the Coulomb and exchange integrals for the pair of atoms defining R,. These are approximated by
(17) London, F. Z . Electrochem. 1929, 35, 552. (18) Parr, C. A.; Truhlar, D.G . J . Phys. Chem. 1971, 75, 1844.
(10)
which allows the barrier along the reaction path to be made broader or narrower; the determination of the four constants a4, b4, c4, and R, is discussed below. The constant term (Vo)in eq 1 is chosen so that the zero of energy corresponds to infinitely separated reactants each at its own equilibrium geometry:
@
ai = al + [ b , - t,(T - R,)] exp[-cl(T - R,)2]
+ b3 exp[-c3(T - R,)2]
(9) where a3and b3were obtained from the WD values of De for OH
All values are in atomic units, except for angles, which are in radians.
Jj(Ri) = ['Ej(Rj)
(8)
Vbcnd
= '/2ke(B - 60)' + !@a(. - ao)' + '/k#
(12)
where 8 is the torsional angle and Bo, ao,ke, k,, and k, are adjustable parameters which also vary with geometry. For Bo, we obtain a value that smoothly varies from the reactants to the products by taking
Bo =
+ j/2bo{tanh[ce(T- R,)] + 1)
(13)
where a#,be, and Cg were determined from the WD values of B in H 2 0and at the saddle p i n t and by assuming that B = 90° for reactants.I0 Similarly, since the value of a. is expected to go to zero at reactants and products, we take a. = a* exp[-pa(R3
- R ')2]
(14) where R * is the WD value of R3 at the saddle point. In order to obtain a value for a at the saddle points of the fitted surfaces that agreed with the WD value, we found it necessary to treat both a* and pa as adjustable constants; their determination is discussed below. The parameter kois expected to be a maximum at products and go to zero at reactants; thus, we take ke = fs expt-go(T
- RJ'I
(1 5 )
wherefs was obtained (together with the constants b, and c3 discussed above) from the WD frequencies for H 2 0 given in Table I; the determination of the value of go is discussed below. On the other hand, since the k, and k, parameters are expected to go to zero at both reactants and products, we take
ka =fa exp[-pa(R3 - R *)']Sa and
(16)
534
The Journal of Physical Chemistry, Vol. 96, No. 2, 1992 k, = fs exp[-p,(R3 - R,)21s,
where Saand S, are switching functions given by 1, R22R‘ Sa = w exp[-y(R,- R’)’] + (1 - w ) , R, c R’
I
Isaacson
r
(17)
I
I
I
I
I
(19
(18)
and
- R*)2]+ (1 - h), Rz c R’
(19)
The purpose of these switching functions is to keep the behavior of the lowest frequency bending vibrations reasonable along the reaction path; in particular, they prevent the corresponding frequencies from becoming imaginary on the reactant side of the saddle point. The determination of the various constants in eqs 16-19 is discussed below. The strategy for determining the remaining constants in eqs 5-19 was as follows: (1) The values of b4, c4, R,, pa, w,y , and R’were adjusted as a group to force the lowest a’ frequency and the potential energy to behave as reasonably as possible along the reaction path and to satisfy an additional criterion. For surface 3, the additional criterion was that the classical energy along the reaction path fit the SE barrier shape as closely as possible at positive energies, where the SE and ab initio barriers are close.I6 For surface 4, the additional criterion was that CVT/SCSAG rate constants (discussed in the next section) for reaction R1 calculated from the PES agree with the measured values of Ravishankara et al.I3 (2) For each choice of the seven constants.in (l), the values of the five independent constants Oil),D{’) = Os’),D‘3’) = DL’), u4, and CY* were adjusted by a weighted steepest descent procedure to reproduce the WD saddle point geometry, barrier height, and imaginary frequency as closely as possible. (3) For each choice of the five constants in (2), the values of the five constants cl, tl, go,f,, and z were determined by a Newton search to reproduce the WD frequencies for the five bound vibrational modes at the saddle point. (4) As a final independent step, the values of the five constants&, pa, R,, h, and R* were chosen as a group to force the a” frequency to be as reasonable as possible along the reaction path. Table I compares the reactant, product, and saddle point properties obtained from the new surfaces 3 and 4 to those of the SE surface and to the WD a b initio data. Due to the manner in which the new surfaces were fit, as discussed above, the reactant and product properties obtained from both new surfaces match the WD data. (Since the H 2 0 region of the SE surface was fit to more reliable data,I5 we would not expect the H 2 0 properties obtained with the SE surface to agree with the WD data.) Both new surfaces also provide a much better fit to the WD data a t the saddle point than does the SE surface. In fact, the only significant differences between the WD saddle point geometry, frequencies, and energy and those obtained from surfaces 3 and 4 occur in the values of R I and R6 (Le., the OH1 and H a b distances), and even these differences are 60.01 bohr. This demonstrates that varying the surface parameters with geometry provides a relatively simple and efficient method for obtaining an accurate fit of a PES to a set of ab initio data. In addition, these results show that although a more complicated form for the PES than what is given in eq 1 [e.g., including a V3(R,,R2,R4) term] could have been employed in this study, such extra complexity is not needed in order to fit the PES to the available data. On the other hand, it is important to note that while the new surfaces discussed herein are defined globally (Le., for all nuclear geometries), their form would be unsuitable for quantal or classical dynamics studies, since the hydrogen atoms are not treated equivalently. Further comparisons between the two new surfaces and the SE surface can be seen in Figure 2, which shows the classical potential energy curves ( VMEp) and the ground-state vibrationally adiabatic along the corresponding reaction potential energy curves paths. The reaction coordinate s is defined as the distance along the minimum energy path (MEP) on each surface through mass-scaled coordinates from the saddle point (s = 0) to reactants
(e)
-3
-1.8
I
-1.5
I
-1.2
I
I
-0.9
-0.6 s(bohr)
1
I
-0.3
0.0
0.3
Figure 2. Classical potential energies, V, (left scale), and ground-state vibrationallyadiabatic potential energies, (right scale), as functions of the corresponding reaction coordinate s for the SE surface (solid curves), surface 3 (dashed curves), and surface 4 (dot-dashed curves). The arrows indicate which scale (left or right) applies to each set of curves.
(s = -a) and products (s = +m), where the coordinates are scaled to a reduced mass of mOHmH,/moH+H, = 1.802023 amu.I1 VMEp(s) is the Born-Oppenheimer potential energy along the reaction path, while is the sum of VMEp(s) and the total (harmonic) zero-point energy for the bound vibrational degrees of freedom orthogonal to the MEP at s. As discussed above, the constants for surface 3 were chosen so that the positive energy portion of VM& would be as close as possible to that of the SE surface, and Figure 2 shows that the fit is quite good. The well in VMEP(s)on the reactant side of the saddle point is even deeper for surface 3 (1.6 kcal/mol) than for the SE surface (1.3 kcal/ mol), in better agreement with the a b initio results.12 Despite the similar VMEp(s)curves for surface 3 and the SE surface, however, those surfaces exhibit markedly different dynamics for reaction R1. As discussed below, this is partly due to the fact that the variational transition state at 0 K [i.e., the highest maximum of Jf(s)]is shifted from the conventional transition state (s = 0.0) by 0.14 bohr for surface 3 but only by 0.10 bohr for the SE surface. In contrast, the VMEp(s) and curves for the surface 4, which was partly determined by requiring that calculated CVT/SCSAG rate constants for reaction R1 agree with experiment, are much broader on the reactant side of the saddle point than those for surface 3, and there is no reactant side minimum in VMEp(s). In addition, the variational transition state at 0 K for surface 4 is at s = -0.17 bohr, leading to an even larger variational effect on the rate constants discussed in the next section than for those of surface 3. For both surfaces 3 and 4, Figure 2 shows spurious local maxima and minima in which arise from artifacts in the variations of the bound vibrational frequencies along the reaction path. Since these anomalies are of limited extent and are quite removed from the global maximum in they are not expected to have a significant effect on the reaction rate constants discussed in the next section, especially at low temper atures.
e(s)
e(s)
e(s),
c(s),
4. Rate Constant Calculations
The calculations of the rate constants for reactions R1 and R2 from the potential energy surfaces discussed herein were carried out by canonical variational transition-state theory (CVT),I9-*’ (19) Truhlar, D. G.; Garrett, B. C. Acc. Chem. Res. 1980, 13, 440. (20) Garrett, B. C.; Truhlar, D. G. J. Chem. Phys. 1979,70, 1593; J . Phys. Chem. 1979,83, 1052, 1079; 1983,87, 4553E. (21) Truhlar, D. G.; Isaacson, A. D.; Garrett, B. C. The Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; Vol. 4, p 65.
The Journal of Physical Chemistry, Vol. 96, No. 2, I992 535
Global Potential Energy Surfaces
TABLE III: Calculated and Experimental Rate Constants for Reaction R1" T. K TST CVT CVTIMEPSAG 200 298 400 600 1000 1500 2400
4000
2.08(-17)' 2.80(-17)d 2.79(-17)' 1.75(-15) 2.44(-15) 2.44(-15) 1.76(-14) 2.52(-14) 2.5 1 (-1 4) 1.82(-13) 2.67(-13) 2.66(-13) 1.58(-1 2) 2.32(-12) 2.3 1 (-1 2) 6.39(-12) 9.36(-12) 9.32(-12) 2.69(-11) 3.92(-11) 3.91(-11) 1.09(-10) 1.59(-10) 1.59(-10)
9.97(-18) 3.80(-18) 2.69(-18) 1.11(-15) 6.60(-16) 5.16(-16) 1.30(-14) 9.86(-15) 8.10(-15) 1.57(-13) 1.52(-13) 1.30(-13) 1.50(-12) 1.75(-12) 1.56(-12) 6.26(-12) 7.81(-12) 7.1 1(-12) 2.66(-11) 3.42(-11) 3.14(-11) 1.08(-10) 1.40(-10) 1.28(-10)
4.91 (-16) 1.72(-16) 2.98(-17) 5.08(-15) 3.08(-15) 1.54(-15) 2.86(-14) 2.25(-14) 1.48(-14) 2.15(-13) 1.70(-13) 1.62(-12) 1.94(-12) 1.70(-12) 6.30(-12) 8.01 (-1 2) 7.3 3 (- 12) 3.40(-11) 3.16(-11) 1.06(-10) 1.39(-10) 1.28(-10)
CVTISCSAG 6.87(-15) 3.13(-14) 3.32(-16) 2.37(-14) 6.06(-14) 6.25(-15) 7.26(-14) 1.27(-13) 3.56(-14) 3.35(-13) 4.62(-13) 2.60(-13) 1.91(-12) 2.54(-12) 2.00(-12) 6.78(-12) 9.01(-12) 7.89(-12) 2.69(-11) 3.56(-11) 3.26(-11) 1.08(-10) 1.42(-10) 1.30(-10)
exDtb
6.08 f 0.37(-15)
3.7 f 1.1(-14) 2.9 f 0.9(-13) 2.4 f 0.7(-12)
"Units are cm3molecule-' s-I; power of 10 is given in parentheses. bFrom ref 13; for T = 298 K we quote the actual tabulated values, while for T > 298 K we employ those authors' three-parameter fit with 30% error bars. 'Upper entry in each case is for the SE surface. "Middle entry in each case is for surface 3. 'Lower entry in each case is for surface 4. rate constant at 298 K and less than a 5% difference in the exactly as described for a CVT study" of this system employing CVT/SCSAG rate constant at 200 K. Since there is no local the SE surface. In CVT, the generalized transition state is optimized for each temperature by finding the maximum of the minimum in VMEp(s)for surface 4, for convenience the MEP was generalized free energy of activation with respect to the reaction determined over the range -2.10 I s I 0.30 bohr, and the scoordinate. (In computing this generalized free energy of actidependent quantities were extrapolated exponentially to the vation for nonzero values of the reaction coordinate, the overall reactants with a range parameter of 0.528 bohr-I. In this case, translations and rotations as well as the unbound motion along however, the rate constants were converged to at least three the reaction coordinate are projected out, leaving 3 N - 7 gensignificant figures with respect to variations in all of the numerical parameters of the calculation. As noted in section 2, in keeping eralized vibrational normal modes orthogonal to the N-atom with the limited nature of the a b initio data to which the new nonlinear reaction path.) Quantum mechanical tunneling along the reaction coordinate was included by multiplying the CVT rate surfaces presented herein have been fit, the bound vibrational constant by one of two semiclassical adiabatic ground-state (SAG) degrees of freedom have been treated in the harmonic approxitransmission coefficient^.'^-^^*^^ The MEPSAG transmission mation throughout. All of the rate constant calculations reported c o e f f i ~ i e n t ~is' ~based ~ ~ . on ~ ~tunneling along the MEP, while the here as well as the characterizations of the potential energy small-curvature (SCSAG) transmission c o e f f i ~ i e n t ~ ~aP* ~ ~ - * ~surfaces discussed in the preceding section were carried out with a prerelease version of the POLYRATE program.z9 proximately accounts for curvature coupling of the reaction path to the generalized normal modes of the transition state. The latter 5. Results transmission coefficient thus allows for corner cutting during Harmonic rate constants for reaction R1 computed from the tunneling and is more accurate.26,28 For comparison, we also SE surface, surface 3, and surface 4 with various levels of appresent conventional transition-state theory (TST) rate constants, proximation are given in Table 111. We first note that since the which result from locating the transition state at the saddle point values of the saddle point and reactant properties of surfaces 3 (s = 0.0). and 4 are nearly the same as the a b initio values,1° TST rate In order to compute CVT rate constants that were converged constants computed from these surfaces are very close to those to at least three significant figures, the MEP was computed with using the a b initio data." As mentioned previously," the difti step size between gradient calculations of 1 X bohr, and ferences between these rate constants and those obtained from the distance between calculations of the Hessian was 0.01 bohr. the S E surface are mainly due to differences in the moments of (All calculations were carried out in coordinates scaled to a reinertia caused by the errors in the saddle point geometry of the duced mass of mOHmH,/mFH,.I1) Due to the local minimum in S E surface. We next note that optimizing the location of the VM&) on the reactant side of the saddle point, the MEP for transition state results in a substantial decrease in the computed surface 3 was determined over the range -1.54 I s I 0.30 bohr, rate constants (CVT vs TST), especially at lower temperatures. and the s-dependent quantities were extrapolated exponentially" In addition, this effect increases on going from the SE surface to the reactants with an exponential range parameter of 0.75 to surface 3 and from surface 3 to surface 4. Since a t low tembohr-I. Choosing other values of this range parameter from 0.50 peratures the CVT bottleneck occurs near the maximum in the to 1.OObohr-l yields less than a 3% difference in the CVT/SCSAG the greater shift of this maximum adiabatic potential, (22) Garrett, B. C.; Truhlar, D. G. J . Chem. Phys. 1980, 72, 3460. toward the reactants on going from the SE surface to surface 3 (23) Garrett, B. C.; Truhlar, D. G.; Grev, R. S.; Magnuson, A. W. J . Phys. to surface 4 (see section 3) leads primarily to a more pronounced Chem. 1980.84, 1730; 1983, 87, 4553E. tightening of the H2 bond (and, for surfaces 3 and 4, of the OH (24) Truhlar, D. G.; Kuppermann, A. J . Am. Chem. Soc. 1971,93, 1840. bond as well) between the saddle point and the CVT bottleneck (25) Skodje, R. T.; Truhlar, D. G.; Garrett, B. C. J . Phys. Chem. 1981, and hence to a larger variational effect. For example, at 298 K 85, 3019.
e@),
(26) Skodje, R. T.; Truhlar, D. G.; Garrett, B. C. J. Chem. Phys. 1982,
77. . , -5955. .- - .
(27) Truhlar, D. G.; Isaacson, A. D.; Skodje, R. T.; Garrett, B. C. J. Phys. Chem. 1982,86, 2252; I983,87,4554E. (28) Marcus, R. A. J. Chem. Phys. 1966, 45,4493; 1969, 49, 2617.
(29) Isaacson, A. D.; Truhlar, D. G.; Rai, S. N.; Steckler, R.; Hancock, G. C.; Garrett, B. C.; Redmon,M. J. POLYRATE. Compur. Phys. Commun. 1987, 47, 91.
536 The Journal of Physical Chemistry, Vol. 96, No. 2, 1992 TABLE IV: Calculated and Experimental Kinetic Isotope Effects" T. K TST CVT CVT/SCSAG exptb 200 298 400 600 1000 1500 2400
5.53c 6.52d 6.53e 4.09 4.40 4.40 3.34 3.47 3.47 2.58 2.62 2.62 1.95 1.97 1.97 1.66 1.68 1.68 1.49 1.50 1.50 1.41 1.42 1.42
3.42 2.61 2.19 2.98 2.40 2.13 2.69 2.24 2.04 2.31 2.02 1.90 1.88 1.78 1.70 1.64 1.61 1.57 1.48 1.49 1.47 1.40 1.43 1.42
24.0 44.2 6.84 8.92 17.8 4.83 5.16 7.93 3.58 3.12 3.64 2.54 2.06 2.17 1.90 1.68 1.73 1.64 1.48 1.51 1.49 1.39 1.43 1.43
3.3 f 0.9 2.5 f 0.7 1.7 f 0.5 -1.8
-1.5
1.6 f 0.5
-1.2
-0.9
-0.6
-0.3
0.0
0.3
0
Figure 3. Magnitudes of the reaction path curvature, K , in units of (mass-scaledbohr)-' as functions of the corresponding reaction coordinate s for the SE surface (solid curve), surface 3 (dashed curve), and surface 4 (dot-dashed curve).
periment at higher temperatures, so that they are in a better position to be improved by a tunneling correction. However, the SCSAG method predicts overly large tunneling effects on the 4000 kinetic isotope effect. Even though the CVT/MEPSAG rate constants for reaction R1 for surface 4 are too low (see Table 111), the CVT/MEPSAG values of the kinetic isotope effect for surface 'Ratio of the rate constant for reaction R1 to that for reaction R2. 4 decrease from 3.24 at 298 K to 1.77 at 1000 K, in good bFrom Table V of ref 13. N S e e Table 111. agreement with experiment. Thus, the effects of including reaction the CVT bottlenecks for the SE surface, surface 3, and surface path curvature by the CVT/SCSAG method appear to be 4 occur at -0.098, -0.139, and -0.161 bohr, respectively, leading somewhat too large for surface 4 and are even larger for surface to increases in the frequencies of the modes that primarily cor3. respond to the H2 (OH) stretches of 738 (-12), 1129 (224), and A more detailed comparison of the reaction path curvature for 1281 (255) cm-', respectively, compared to those of the saddle these surfaces can be seen in Figure 3, which shows the magnitude point. of the reaction path curvature K ( S ) 2 1 ' 2 5 along the corresponding Thirdly, we find that the inclusion of reaction coordinate reaction paths. Although the peak in the reaction path curvature tunneling along the MEP markedly increases the computed rate for surface 3 at s = -0.64 bohr is smaller than that for surface constants (CVT/MEPSAG vs CVT) for T 5 600 K for all three 4 at s = -0.76 bohr, the fact that it occurs where the adiabatic surfaces studied. Since the adiabatic potential curve provides the barrier is higher (cf. Figure 2) leads to a much more pronounced effective tunneling barrier in the MEPSAG a p p r o ~ i m a t i o n , 2 ~ - ~ ~ .effect ~ ~ on the computed transmission coefficient. Thus, in order the similarity of the curves for the SE surface and surface to understand the effects of reaction path curvature, it is important 3 (see Figure 2) lead to similar effects in this approximation, while to consider both its shape and magnitude along the reaction the much brceder curve for surface 4 yields less pronounced as well as its relationship to the adiabatic potential curve, increases. In addition, we note that incorporating the curvature 6. Summary of the reaction path in the tunneling calculation also greatly As discussed in section 3, the variation of the surface parameters increases the computed rate constants (CVT/SCSAG vs along the reaction path provides an efficient and accurate method CVT/MEPSAG) for T 6 lo00 K for all three surfaces and that for fitting a PES to a set of ab initio and/or experimental inthe effects are much larger for surface 3 than for the SE surface formation on reaction R1. But the results in Tables I11 and IV or surface 4. show that the dynamical results predicted from a given PES Finally, when we compare the most reliable treatment depend on both the functional form of the PES as well as the set (CVT/SCSAG) with the experimental measurements of Raof data to which it is fit. In particular, the results show that the vishankara and co-workersI3 listed in Table 111, we find that barrier shape and the reaction path curvature can strongly inalthough all three surfaces give rate constants that agree with fluence the rate constants predicted from a given PES. On the experiment within experimental error at 1000 K, both the SE basis of the CVT/SCSAG rate constant results for reaction R1 surface and surface 3 CVT/SCSAG rate constants are too high listed in Table I11 and the kinetic isotope effect values presented for T 5 400 K. (Surface 4 was determined such that the in Table IV, surface 4 appears to provide a better description of CVT/SCSAG rate constants obtained from it agree with exthe PES for reaction R1 than either the SE surface or surface periment within experimental error from 298 to lo00 K see section 3. In addition, V,,,(s) for surface 4 does not exhibit a reac3.) The fact that the CVT/MEPSAG rate constants are sometant-side well. Note, however, that only the reaction valley for what below the experimental ones for T 5 4 0 0 K suggests that reaction R1 (Le., the region of the PES along the MEP) has been the errors in the SE surface and surface 3 CVT/SCSAG rate considered here. constants at lower temperatures could be due either to an adiabatic The question raised above concerning the relative contributions barrier that is too narrow or to a contribution from reaction path of the narrowness of the adiabatic barrier and the degree of curvature that is too large. In fact, the effects of including reaction reaction path curvature to the prediction of rate constants for path curvature are quite similar for the SE surface and surface reaction R1 and kinetic isotope effects that are much larger than 4, but the much broader adiabatic barrier of surface 4 leads to the corresponding experimental values at lower temperatures (see CVT/SCSAG rate constants that agree with experiment. section 5 ) can be further investigated in two ways. First, we note Further indication that the effects of including reaction path that the OH + H2geometry for which VMEp(s) for surface 3 curvature are too large for these surfaces is presented in Table crosses zero on the reactant side of the saddle point (geometry IV, which compares calculated and experimental values of the A) is R , = 1.858 bohr, R2 = 3.452 bohr, R6 = 1.426 bohr, 9 = kinetic isotope effect, Le., the ratio of the rate constant for reaction 98.3O, and a = 67.3O, while the OH + H2 geometry along the R1 to that for reaction R2. Although the TST values are higher reaction path for surface 4 at approximately the same value of than experiment over the whole temperature range, the CVT R2 (geometry B) is R, = 1.858 bohr, R, = 3.449 bohr, & = 1.428 values are below experiment at low temperatures and above ex-
e(s) e(s)
8;
J . Phys. Chem. 1992, 96, 537-540 bohr, 8 = 9 5 . 6 O , and a = 50.1'. Thus, while the major difference between these two geometries is in the value of angle a,the value of V,,,(s) for surface 4 at geometry B is about 2 kcal/mol. To gain a better understanding of the correct shape of the potential energy along the reaction path, we therefore invite workers who perform very large scale ab initio calculations to carry out high-quality calculations at four geometries of the OH H2 system, namely, the saddle point geometries of surfaces 3 and 4 listed in Table I and geometries A and B, to determine whether the energies at the saddle point and at geometries A and B are roughly 6 and 0 kcal/mol, respectively, or 6 and 2 kcal/mol, respectively. Second, to study whether the SCSAG method exaggerates the effects of reaction path curvature in the results presented in section 5, calculations employing the surfaces discussed herein are currently in progress to compute transmission coefficients by the large-curvature ground-state approximation, version 3 (LCG3) m e t h ~ d , ~which ' , ~ ~ is better suited for systems with large reaction path curvature because it involves a wider
+
(30) Garrett, B. C.; Joseph, T.; Truong, T. N.;Truhlar, D. G. Chem. Phys. 1989, 136, 271.
537
region of the PES than can be represented in terms of coordinates based on the MEP. These calculations will thus provide further information on the sensitivity of the transmission coefficients to the form of the PES used. In addition, further refinements of the PES for reaction R1 are currently in progress to provide a surface that reproduces the experimental reactant and product properties as well as the more recent anharmonic ab initio information of Kraka and D ~ n n i n g . ~Such ' a surface can then be used to study the effects of including vibrational anharmonicity in the calculation of rate constants for reactions R1 and R2, which have already been shown" to be substantial for the SE surface.
Acknowledgment. The author is grateful to Donald Truhlar and Joseph Simpson for helpful suggestions and to Thom Dunning for unpublished results. All calculations reported here were camed out on the Miami University V A X computer, and the computer time is appreciated. Registry NO. 'OH, 3352-57-6; H2,1333-74-0. (31) Kraka, E.; Dunning, T. H., Jr. To be published. (32) Schatz, G. C.; Walch, S. P. J. Chem. Phys. 1980, 72, 776.
Absolute Free Energy of Solvation from Monte Carlo Simulations Using Combined Quantum and Molecular Mechanical Potentials Jiali Gao Department of Chemistry, State University of New York at Buffalo, Buffalo, New York 14214 (Received: September 3, 1991; In Final Form: October 22, 1991)
A new method for computing the absolute free energy of solvation using combined quantum mechanical (QM) and molecular mechanical (MM) potentials and the Monte Carlo method is presented. In this procedure, the Hamiltonian associated with the QM region is first annihilated to form a pure classical (empirical) potential. The free energy for the latter system can then be easily evaluated using traditional approaches. The method is illustrated by computing the solvation free energy of Cl- in water, treating the ion quantum mechanically during the free energy simulation. A good agreement was found between the calculated and experimental results.
1. htroduction The development of free energy simulation techniques has greatly enhanced our understanding of chemical processes in solution.' The progress is largely due to the availability of intermolecular potential functions that can properly describe the molecular interactions in solution.* These effective potential functions are typically pairwise additive and the parameters are determined by fitting the computed structural and energetic results against experimentally determined properties? However, empirical potentials are not suited for describing detailed polarization effects and chemical processes involving bond forming and breaking.3 To overcome these difficulties, several models have recently been proposed, in which quantum mechanical (QM) and molecular mechanical (MM) potentials are combined in molecular dynamics simulation and energy minimizati~n!,~ In these approaches, part
of the system, e.g., the solute molecule or the active site region of an enzyme, is treated quantum mechanically, while the environment or solvent molecules are approximated by molecular mechanical potentials. Among other alternatives, the FieldBash-Karplus approach is particularly suited for liquid simulations since semiempirical AM1 and MNDO methods6 were adopted to form the combined QM/MM force field! The major advantages of using semiempirical quantum mechanical theories include their highly computational efficiency (ca. lo00 times faster than ab initio calculations7at comparable levels) and excellent results for many systems that have been parametrized? This method, however, has not been implemented in Monte Carlo simulations. Applying the AM 1(MNDO)/CHARMM potential, the free energy of activation for an SN2reaction in water has been computed.8 In that study, a partitioning of the Hamiltonian, H(X)
(1) For recent reviews, see: (a) Kollman, P. A.; Merz, K. M. Acc. Chem. Res. 1990,23,246. (b) Jorgensen, W. L. Acc. Chem. Res. 1989,22, 184. (c) Beveridge, D. L.; DiCapua, F. M. Annu. Rev. Eiophys. Chem. 1989,18,431. (d) Straastsma, T. P.; McCammon, J. A. J . Chem. Phys. 1991, 95, 1175. (2) (a) Jorgensen, W. L.; Tirado-Rives, J. J . Am. Chem. SOC.1988,110, 1657. (b) Weiner, S.J.; Kollman, P. A.; Case., D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.Jr.; Weiner, P. J . Am. Chem. SOC.1984, 106, 765. (c) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.;Karplus, M. J . Comput. Chem. 1983, 4, 187. (3) (a) McCammon, J. A.; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University Press: Cambridge, MA, 1987. (b) Brooks, C. L., 111; Karplus, M.; Pettitt, B. M. Adv. Chem. Phys. 1988, 71. (4) Field, M. J.; Bash, P. A.; Karplus, M. J . Comput. Chem. 1990, 11,700.
( 5 ) (a) Singh, U. C.; Kollman, P. A. J . Comput. Chem. 1986,7,718. (b) Weiner, S.J.; Seibel, G. L.; Kollman, P. A. Proc. Natl. Acod. Sci. U.S.A. 1986,83,649. (c) Hwang, J.-K.; King, G.; Greighton, S.; Warshel, A. J . Am. Chem. SOC.1988, 110, 5297. (d) Rullmann, J. A. C.; Bellido, M. N.; van Duijnen, P. T. J . Mol. Eiol. 1989, 206, 101. (e) Bash, P. A.; Field, M. J.; Davenport, R. C.; Petsko, G. A.; Ringe, D.; Karplus, M. Biochemistry 1991, 30, 5826. (6) (a) Dewar, M. J. S.;Thiel, W. J. Am. Chem. SOC.1977,90,4899. (b) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J . Am. Chrm. SOC.1985, 107, 3902. (7) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (8) Bash, P. A.; Field, M. J.; Karplus, M. J . Am. Chem. Soc. 1987, 109, 8192.
0022-365419212096-537$03.00/0 0 1992 American Chemical Society
