The Journal of
Physical Chemistry
Registered in US.Patent Office 0 Copyright, 1981, by the American Chemical Society
VOLUME 85, NUMBER 21
OCTOBER 15, 1981
LETTERS A General Small-Curvature Approximation for Transition-State-Theory Transmission Coefficients Rex T. Skodje, Donald G. Truhlar," Chemical Physics Program and Department of Chemlstty, University of Minnesota, Minneapolis, Minnesota 55455
and Bruce C. Garrett Chemical Dynamics Corporation, Columbus, Ohio 43220 (Received: Ju& 13, 198 1)
We derive an expression for the effective mass along the reaction coordinate in a chemically reacting system. This expression is used in two convenient algorithms for including quantum mechanical effects on the reaction coordinate into transition-statetheory or generalized transition-statetheory. One algorithm is based on a uniform semiclassicalapproximation;the other is based on a parabolic approximation to the effectivepotential. Examples are presented to show that the theory accurately accounts for reaction-path-curvature effects on the tunneling contributions to the reactions H + H2,D + Dz,and D + H p Although the derivation, and hence expected reliability, of the new methods is most valid for systems with small reaction-path curvature, the methods are generally applicable and should provide useful accuracy even for medium-curvatureand some large-curvature systems. In the last few years it has become very clear that reaction-path curvature must be considered as an integral part of one-dimensional models for tunneling in chemical reactions; in particular one-dimensional transmission coefficients based on zero-curvature approximations for the reaction coordinate seriously underestimate the tunneling for many systems, and models that do include effects of reaction-path curvature are much more successful.l-' In the present paper we suggest and we test a new (1)D. G. Truhlar and A. Kuppermann, Chem. Phys. Lett., 9, 269 (1971);J. Chem. Phys., 56, 2232 (1972). (2)R.A. Marcus and M. E. Coltrin, J. Chem. Phys., 67,2609(1977). (3)B.C. Garrett and D. G. Truhlar, J.Phys. Chern., 83,200,3058(E), 1079 (1979);84,682(E)(1980). 0022-3654/81 /2085-3019$01.2510
semiclassical approximation for tunneling probabilities in systems with curved reaction paths. We also give a further approximation (for systems with a single well-defined maximum in their vibrationally adiabatic ground-state potential energy barrier) which reduces the tunneling probability to a simple but reasonably accurate analytic formula with the minimum number of parameters. (4)W.H.Miller, N. C. Handy, and J. E. Adams, J. Chem. Phys., 72, 99 (1980). (5) B. C. Garrett and D. G. Truhlar, J. Chem. Phys., 72,3460 (1980). (6)B. C. Garrett, D. G. Truhlar, R. S. Grev, and A. W. Magnuson, J. Phys. Chern., 84,1730(1980). (7)B. C. Garrett, D. G. Truhlar, R. S. Grev, and R. B. Walker, J. Chem. Phys., 73,235 (1980).
@ 1981 American Chemical Society
3020
The Journal of Physical Chemistry, Vol. 85, No. 2 1, 198 1
The most consistently successful calculations performed so far2v3t6’ for atom-diatom systems have been based on the tunneling path suggested for H + Hz by Marcus and Coltrin.2 In this method reaction-path curvature causes the system to tunnel on the concave side of the minimum-energy path (MEP). This is qualitatively similar to a negative centrifugal effect.“l0 Systems with positive kinetic energies show a positive internal centrifugal effect or “bobsled effect: i.e., the reaction-path curvature forces them to the convex side of the MEP. Tunneling systems, in contrast, show a negative internal centrifugal effect; they “cut the corner” to minimize the amount of exponential damping of the wave function in the passage through the classically forbidden region. This argument can be made more rigorous by analytic continuation of the least-action variational principle of classical mechanics to imaginary momenta.2 The new method presented here is a generalization and extension of this approach. We begin with an anharmonic generalization of the harmonic, zero-angular-momentum polyatomic reactionpath Hamiltonian of Miller et al.4 [their eq 4.21 and 4.23~1, which we divide into the “principal” terms Hoand the “interaction” terms H1and Hz % = Ho + HI+ Hz (1)
Letters
momenta to u, and s. For convenience we have combined the componpts u,, Pa,and B d s ) into (F - 1)-dimensional vectors 2, P, and Z(s) in eq 1-6. As discussed elsewhere: a consistent transmission coefficient K ~ ” ( T )for the adiabatic theory of reactions12 is J
m
~
~
exp(-E/kn ~ ( ~ ) d~
~
=
KVAG(T)
where PvAG(E) is the quantal transmission probability at total energy E for the one-dimensional vibrationally adiabatic ground-state (VAG) potential curve vaG(s)= vMEP(s) + fG(s)
(8) with maximum value VAG. The quantity &(s) is the vibrational zero-point energy of the fixed-s Hamiltonian, Le., + Hl- V m ( s ) ]evaluated at the lowest eigenvalue of [Ho s with p s = 0. Vibrational1 adiabatic ground-state transmission coefficients K*/V* (T) for conventional transition-state theory may than be computed bf
2
jVIaexp(-E/kT)
KvAG(T)
K*/VAG(T)
a
H1 =
cc
Ck,fJr(S)U,UflU,
~ > f l @ > -aa
A(s,Z) Hz = -{[P$ 2M
+ ...
(3)
- C C~J’$a&s)l~ - ~ , 2 1 (4) @#a a
A(s,Z) = [ l - ~‘-Z(S)]-~
(5)
K(s)Ia = -BaF(S)
(6)
Here we work in mass-scaled coordinates3J obtained by scaling each atomic Cartesian coordinate by (mass of atom/p)1/2. This makes the reduced mass be the same arbitrary quantity for every degree of freedom; for concreteness we usually choose p as the reduced mass of relative translational motion of reactants. Then s is defined as the reaction coordinate, and its magnitude measures arc length along the path of steepest descents3J1 though the mass-scaled coordinate system, The F - 1 coordinates u, are the quasinormal coordinates generated by the projected-mode analysis of Miller et al.;4 F is 3N - 5 or 3N - 6, N is the number of atoms, and all sums go up to F - 1. All u, are defined to be zero on the path of steepest descents and to be positive on the concave side of this path [this sign convention determines the sign in eq 51. The vibrational force constants are called k,,, k,,,, k,, etc. and the classical (i.e., Born-Oppenheimer) potential energy along the reaction coordinate is called Vmp(s) since the path of steepest descents is the minimum-energy path connecting reactants to produds. B,&) is the nonadiabatic coupling matrix: Bkds) is a component of the reaction-path curvature, and Paand p s are conjugate (8) R. A. Marcus, J . Chem. Phys., 45,4493 (1966).
(9)R.A. Marcus, J. Chem. Phys., 49,2617 (1969). (10)W.H.Miller, J. Chem. Phys., 61,1823(1974);62,1899(1975);63, 1166 (1975);Acc. Chern. Res., 9,306 (1976);Faraday DMCUSS. Chem. SOC., 62.40 (1977): S. ChaDman. B. C. Garrett. and W. H. Miller, J . Chern. P ~ Y S .63,27io , (wsj. (11)B. C.Garrett and D. G. W. . J. Phvs. . Chern., 83,1062,305NE) (i979j. ’
d~ (9)
When the barrier position of the VAG potential curve ~ (KT~) ~ ~are( the ! ! ‘ ) occurs at the saddle point,lZK * / ~ ~and same. When the vibrational coordinates are treated adiabatically the equation for the reaction-coordinate wave function becomes
I -
2~
d2 ds2
- + V,G(s) - E
where (A)(s) is determined by an appropriate average of A(s,ii) over the vibrational degrees of the freedom along the tunneling path. There are several ways to approximate (A)(s) and these lead to several versions of K ~ ~ ~ This (T). is the heart of the problem addressed in this paper. The simplest approximation is to assume zero curvature and that the system follows the minimum-energypath.l8 Then (A)(qAG=1. This is called the MEP approximation [ (T) and K ’ / ~ T)] ~ ;unfortunately ~ ~ ( it often leads to serious underestimates of the tunneling contribution^.^-' Attempts to evaluate ( A )(s) more rigorously suffer from the breakdown of the curvilinear (s,{u,}fi\) coordinate system for large positive u,14J5and from the breakdown of the vibrationally adiabatic separation in this coordinate system at low E.l0J6 Thus one must make physically motivated approximations and be careful to avoid singularities. The classical phase-average method of Miller et al.4 provides an improved estimate of (A)(s) by averaging over trajectories on both the concave and convex sides of the MEP, and the second-order method of Gray et al.lS provides improved estimates of the transmission probability without assuming the approximate s,u,separability that is implicit in eq 10 or in the use of (A)(s). In this (12)D.G. Truhlar, J. Chem. Phys., 63,2041(1970). (13)D.G. Truhlar and A. Kuppermann, J. Am. Chern. SOC.,93,1840 (1971). (14)J. M. Bowman, Ph.D. Thesis, California Institute of Technoloav. __ Pasadena, 1974. (15)R. T. Skodje, D. G. Truhlar, and B. C. Garrett, to be published. (16)S. K. Gray, W. H. Miller, Y. Yamamchi, and H. F. Schaefer. 111. J. Chem. Phys., 73,2733 (1980).
The Journal of Physical Chemistty, Vol, 85, NO. 21, 1981 5021
Letters
TABLE I: -
Transmission Coefficients for Collinear H t H, + H, t H
T.K
MEPSAG
MCPSAG
MCPPAG
SCSAG
SCPAG
accurate
200 300 400 600 1000 1500
2.08 1.25 1.12 1.05 1.02 1.01
25.5 2.81 1.67 1.23 1.08 1.03
31.1 2.81 1.50 1.14 1.04 1.02
22.5 2.67 1.63 1.22 1.07 1.03
28.3 2.65 1.43 1.12 1.03 1.01
29.8 3.32 2.01 1.44 1.20 1.08
article we attempt to estimate (A)(s) by considering a corner-cutting path similar to that suggested by Marcus and Coltrin2 for the collinear triatomic system H + H2. However, our derivation is general for any triatomic or polyatomic system and it uses a small-curvature limit and an exponential approximation to avoid the singularities associated with the breakdown of the coordinate system for high-curvature cases. I t is useful, for what follows, to reconsider the effect of curvature as a change in reduced mass and to approximate PVAG(E)semiclassically. The primitive semiclassical approximation to the tunneling probability at total energy E is’’ P(E)= exp[-28(E)] (11) where 0 represents the amount of exponential decay. For tunneling along the MEP this is given by
OmP(E) = h-11(2p[VaG(s)- E])’/2 ds
(12)
The integration limits for eq 12 are those s values for which the radicand is positive, i.e., the limits of the tunneling region. If s represents the distance along the MEP but the tunneling path is located at iit(s) on the concave side of the MEP, then it is easily shown by analytic geometry that an element of length along the tunneling path, d[, differs from ds by a factor of d[/ds = ( [ l - it(s).u’t(~)]~ + (dZt/ds)2}1/2
(13)
[Notice that it(s)-iit(s)> 0 for any point s,iit(s)on the concave side of the MEP.] For the tunneling path the action integral is
Combining (13) and (14) yields
@ ( E )= h-‘J(2/~eff(~)[Va~(~) - E])’/2 ds
(15)
where the effective mass associated with the length element ds is defined by pe&)/F [ l - it(s).iit(s)l2+ (diit/ds)2 (16) Comparison to eq 10 shows that the right-hand side of (16) may be considered as an approximation to l/(A)(s). Thus the metric factor for the tunneling path yields the effective mass for motion measured along the MEP. In the above exposition we have assumed a single tunneling path as motivated by classical S matrix theoryls or by the variational argument of Marcus and ColtrinS2To specify the optimum path, i.e., the path of least exponential decay, we consider the small-curvature limit. Marcus and Coltrin argued that within the classical turning points of the ground-state local vibrational motion the effective potential is indeed VaG(s),but outside this region the Born-Oppenheimer .potential V exceeds VaG(s) and one (17) L. D. Landau and E. M. Lifshiftz, “Quantum Mechanics”,2nd ed, Pergamon Press, New York, 1965, p 173. (18) T. F. George and W. H. Miller, J. Chem. Phys., 56,5722 (1972); 57, 2468 (1972).
should use V instead of VaG(s)in the integrand of (14). For this ansatz, in the limit of small curvature, and in the limit that the first square term in (16) dominates the second one, the path of classical turning points on the concave side of the MEP provides the least-action path (ie., least exponential damping). This generalizes the numerical result of Marcus and Coltrin who showed by actual calculation that this path has the lowest action from among a sequence of parametrized paths for collinear H + Hz.Notice that turning points are not strictly defined for coupled-mode systems so we neglect Hlto define the turning points [but not necessarily in the calculation of VaG(s)].Notice also that for the path of turning points all Pa are zero so the contribution of H2 vanishes. When, however, we substitute these turning points in (16) the right-hand side of (16) shows spurious behavior if it(s)dt(s) > 1. This spurious behavior is associated with the breakdown of the coordination system on the concave side of the MEP. (This is especially clear for the case F = 2, corresponding to collinear A BC. For this case K is the reciprocal of the ~ radius of curvature of the MEP, and ut > 1 / corresponds to a triple-valued region of natural collision coordinates!) For small enough curvature though, it(s)dt(s) < 1 and (16) is a valid result. Thus we treat (16) as a small-curvature expansion and we “un-expand” it to remove the singularity to infinity where it causes no trouble. This yields
+
peff(s)/p =
exp[-2a - u2 + (diit/dsI2]
(17)
where u = it(s).iit(s)
Notice that through quadratic terms in iit(s),eq 17 is identical with eq 16, but it avoids the spurious behavior associated with breakdown of (16) for large-curvature systems. In particular, unlike (16), eq 17 never yields pefp = 0 [which would cause a singularity in (10); see the comment following (16)], and, also unlike (16), it predicts that peff(s)always decreases monotonically if any component of the curvature is increased with the other components and iit(s)held constant. Equation 17 still has a serious defect. If shows a contribution even from modes for which the second square term of (16) dominates the first. Such would be the case, e.g., if the path of turning points shows a “bubble”, i.e., has its concave side facing the concave side of the MEP. In such a case a shorter, and hence better, tunneling path is achieved by cutting across parallel to the MEP rather than by following the path of turning points. An even more dramatic example is the case where K, = 0. In this case u, should be zero for the best tunneling path. These defects in (17) can be remedied by writing the exponential of a sum as a product of exponentials and retaining only the factors that lead to peff< p. This yields
(19) where f, = min
{
y p ( - 2a, - a,z
’+ W a Z )
3022
The Journal of Physical Chemistry, Vol. 85, No. 21, 1981
aa =
Letters
(21) w, = du,/ds (22) To finish the calculation we substitute (19) into (15) and use numerical quadrature. Then, instead of substituting P ={exp[(P- a)(AVAG- V2)l - 1) a IP (24) Ot(E)into the primitive semiclassical expression (ll),we P . use a parabolically uniformized algorithm including nonwhere classical reflection, as presented e1~ewhere.l~The final transmission coefficients obtained by this procedure are V2 = max(AEo,O) (25) called the small-curvature-approximationsemiclassical and adiabatic ground-state approximation, which, to be consistent with the systematic abbreviations we have introa = (27r/h)(p/f)'/2 (26) duced reviously? we abbreviate as SCSAG, i.e., K ~ ~ ~ ~ ( T ) and K * ~ ~ ~ ~ Notice ( T ' )that . although the SCSAG algor= ( 2 ~ h/ Iweffl) (27) ithm can be applied, with no spurious coordinate problems where weff is the effective parabolic frequency and = or singularities,to systems with any curvature of the MEP, 1/(kBT). it is expected to be most accurate for small-curvature and To include the effects of reaction-path curvature into medium-curvature systems for which the identification of the parabolic approximation we must use the new measure &(s) as the path of turning points is best justified physlength provided by ( A ) ( s )or peff(s)to compute a. In ically. For high-curvature systems the components of iit(s) particular, if the tunneling-path distance 4 and the reaction are expected to be greater than the distances to the turning coordinate s have the same origin, they are related by points, and thus ~,fi(s) provides a rough upper bound. This means we expect to calculate reasonably accurate tunneling = S , ' [ ~ e f f ( s ) / ~ I ' / ~ ds (28) probabilities for small-curvature systems and to underestimate the tunneling probabilities for high-curvature The effective force constant f is then found by a parabolic systems. To significantly improve the accuracy of the fit to the function [VAG(s),[(s)]. The effective frequency, tunneling probabilities for high-curvature systems one which is computed only for interpretative purposes, is still probably has to abandon the concept of peff(s)since the given by MEP does not provide a good reference path for the true tunneling path in such systems.22 IweffI = (f/~)'/' (29) For three-dimensional (or collinear) atom-diatom reResults calculated from eq 23-29 will be called SCPAG actions with collinear MEP's, only one component of s t ( ~ ) (small-curvature approximation parabolic adiabatic is nonzero. For such systems if the path of classical turning ground-state). To test the new algorithms we applied them points on the concave side of the MEP never strays into to collinear H + Hz, D + D2, and D + H2 using a realistic the multivalued region and if the first square term of (16) potential energy surface' for which accurate quantal results always dominates the second one, then results computed are available.'r21 The results are given in Tables 1-111. by eq 16 would be the same as the results computed by The tables also show the MCPSAG and MCPPAG (parthe Marcus-Coltrin-path semiclassical vibrationally adiaabolic approximation to MCPSAG) results, the zero-curbatic ground-state (MCPSAG) algorithm3 that we have vature semiclassical approximation (MEPSAG), and the used previously. The re-formulation as a small-curvature accurate (acc) transmission coefficient. The latter is deapproximation to the effective mass and the incorporation fined as the value that makes the adiabatic theory of reof the minimum criterion of eq 20 have allowed generalactions exact, i.e. ization to noncollinear MEP's and polyatomic systems, including those with rapidly changing zero point energies kaCC(T)= K a w k A ( T ) (30) [and hence large w(s)]and extension [eq 171 to remove the where kacc(T) is the accurate quantum mechanical rate singularities. The theory also allows for a simple parabolic coefficient, based on converged multidimensionalquantum approximation, which we consider next. K"( 5") is the accurate transmechanical For systems with he&) a constant, excellent results can mission coefficient, and k A ( T ) is the adiabatic-theory rate often be obtained with a parabolic approximation to V,G(s) coefficient for classical reaction-coordinatemotion. Notice if the parabolic force constant is chosen to give a good fit that kA(T) may also be computed by microcanonical varto the original VaG(s)over the significant tunneling region iational transition-state theory since, with classical reac(STR).20 This approximation is especially interesting tion-coordinate motion, that theory gives the same results because the thermally averaged tunneling probability is as the adiabatic theory of Notice, however, given with good accuracy by a simple analytic formula. that consistent transmission coefficients for conventional This formula requires only the barrier height AVAG(equal transition-state theory should contain another factor, as to VAGminus the reactants' zero point energy eRG), the explained above [this factor is called in ref 61. endoergicity AEo,and the force constant f for the parabolic A few additional details of the results in the tables: For approximation to VaG(s)as given by collinear atom-diatom collisions 2(s) and ii&> become VaG(s)= cRG s < -(2AVAG/f)'/' scalars. VaG(s), K ( s ) , and ut(s) were computed by the methods of ref 3,6, and 11 (Morse approximation I), and = VAG - 1/2fs2 -(2AVAG/f)'/2 5 s 5 f was computed by the method recommended in ref 20. [2(AVAG- L\E,)/f11/2 (23) Examination of Tables 1-111 shows that the SCSAG method agrees well with the MCPSAG method for the eRG + AB, [2(AVAG- h E ~ ) / f l '
