Critical tests of variational transition state theory and semiclassical

Mar 20, 1991 - (1) For an overview of the collinear tests see the following references: (a). Truhlar, D. G. ... 1982, 86, 2252. .... Kramers-Brillouin...
5 downloads 0 Views 821KB Size
Download PDF
J . Phys. Chem. 1991, 95, 10374-10379

10374

Critical Tests of Variational Transition State Theory and Semiclassical Tunneling Methods for Hydrogen and Deuterium Atom Transfer Reactions and Use of the Semiclassical Calculations To Interpret the Overbarrier and Tunneling Dynamics Bruce C. Garrett* Molecular Science Research Center, Pacific Northwest Laboratory, Richland, Washington 99352

and Donald G.Truhlar* Department of Chemistry and Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 55455-0431 (Received: March 20, 1991)

Variational transition state theory with semiclassical tunneling calculations is used to calculate rate constants which are compared to the best available quantum mechanical results for three cases involving H atom transfer and three cases involving D atom transfer. The mass combination is heavy-light-heavy, as appropriate for testing semiclassicalmethods for proton and deuteron (or hydride) transfer in organic and biochemical systems. In five of the six cases at 312.5 K and five (but not the same five) of the six cases at 423.2 K, semiclassical results based on an uncoupled treatment of the vibrational-rotational modes agree with the the best available quantum results within the reliability of the latter. In the other two cases the discrepancies are larger, a factor of 3.1 at 312.5 K and a factor of 2.3 at 423.2 K, and possible reasons for this are discussed. In particular, the effect of a coupled treatment of the bending degrees of freedom in these two cases is presented and compared to the uncoupled treatment. Coupling the bending modes reduces the maximum discrepancy at 312.5 K to a factor of 2.4 but increases the largest discrepancy at 423.2 K to a factor of 4.5. The approximate semiclassical calculations are used to interpret the dynamics, especially in cases where agreement with the fully quantal results is very good.

1. Introduction Variational transition state theory calculations with semiclassical tunneling coefficients (VTST/ST) with an uncoupled treatment of the vibrational-rotational modes has been widely tested against quantum mechanical results for many reactions, both against accurate quantum dynamics for collinear reactions and against best available quantum dynamics for three-dimensional reactions, by comparing VTST/ST and the accurate or best available quantum thermal rate constants calculated for the same potential energy When the semiclassical tunneling coefficient includes multidimensional effects by the least-action ground-state (LAG) method4-or by choosing between the small-curvature semiclassical adiabatic ground-state (SCSAGS)or large-curvature ground-state (e.g., LCG6 or LCG3’) methods-the agreement between the VTST/ST and the accurate or best available quantum rate constants has been uniformly excellent, typically within a factor of 1.5 or better at 300 K. Although we have made special (1) For an overview of the collinear tests see the following references: (a) Truhlar, D. G.; Isaacson, A. D.; Skodje, R. T.; Garrett, B. C. J . Phys. Chem. 1982,86,2252. (b) Truhlar, D. G.; Hase, W. L.; Hynes, J. T. J. Phys. Chem. 1983.87, 2664, 5523 (E). (c) Truhlar, D. G.; Garrett, B.C. Annu. Reu. Phys. Chem. 1984,35,159. (d) Garrett, B. C.; Truhlar, D. G. J . Chem. Phys. 1984, 81, 309. (e) Truhlar, D. G.; Garrett, B. C. J. Chim. Phys. 1987,84,365. (f) Lynch, G. C.; Truhlar, D. G.; Garrett, B. C. J . Chem. Phys. 1989,90, 3102. (2) For comparisons of VTST/ST and accurate quantum dynamics rate constants for three-dimensional calculations see refs. la and I f and: (a) Claw. D. C. J. Chem. Phys. 1985, 83, 1685. (b) Garrett, B. C.; Truhlar,‘ D. 6.1 Schatz. G. C. J. Am. Chem. Soc. 1986. 108. 2876. (3) Additional tests of semiclassical tunneling methods are based on comparing them to accurate quantum dynamics calculations for tunneling probabilities in the threshold region: (a) Haug, K.; Schwenke, D. W.; Truhlar, D. G.; Zhang, Y.; Zhang, J. 2.H.; Kouri, D. J. J . Chem. Phys. 1987,87, 1892. (b) Zhang, J. 2. H.; Zhang, Y.; Kouri, D. J.; Garrett, B. C.; Haug, K.; Schwenke, D. W.; Truhlar, D. G. Faraday Discuss. Chem. Soc. 1987,84,371. (c) Lynch, G. C.; Halvick, P.; Truhlar, D. G.; Garrett, B. C.; Schwenke, D. W.;and Kouri, D. J. Z . Naturforsch. 1989, 44a, 427. (4) Garrett, B. C.; Truhlar, D. G. J. Chem. Phys. 1983, 79, 4931. (5) (a) Skodje, R. T.; Truhlar, D. G.; Garrett, B. C. J . Phys. Chem. 1981, 85, 3019. (b) Skodje, R. T.; Truhlar, D. G.; Garrett, B. C. J. Chem. Phys. 1982, 77, 5955. (6) Garrett, B. C.; Truhlar, D. G.; Wagner, A. F.; Dunning, T. H. J . Chem. Phys. 1983, 78, 4400. (7) Truhlar, D. G.; Isaacson, A. D.; Garrett, B. C. In Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; Vol. 4, p 65.

efforts to include a wide variety of mass combinations and types of potential energy surfaces in these compari~ons,I-~ it is always possible that larger errors may be found in new cases. Thus further tests are of great interest, especially since the semiclassical methods are applicable to polyatomic reactionss for which accurate quantum dynamics calculations are not feasible. This paper is our contribution to the proceedings of the Conference on Photoinduced Proton Transfer Dynamics in Chemistry, Biology, and Physics in Honor of Michael Kasha. Testing semiclassical methods for the heavy-light-heavy mass combination is particularly relevant for proton-transfer processes in organic and biochemical systems that serve as the focus for the present symposium. With this background, the recent quantum studies of Schatz et al.9 of the three-dimensional C1 + HCl hydrogen-transfer reaction are of great interest and provide a stimulating challenge to semiclassical methods. Schatz et al. reported quantum mechanical rate constant calculations for both H and D transfer for three different potential energy surfaces, denoted BCMR,’O PK3,” (8) (a) Isaacson, A. D.; Truhlar, D. G. J . Chem. Phys. 1982.76, 1380. (b) Brown, F. B.; Tucker, S. C.; Truhlar, D. G. J . Chem. Phys. 1985,83,4451. (c) Hancock, G. C.; Rejto, P.; Steckler, R.; Brown, F. B.; Schwenke, D. W.; Truhlar, D. G. J. Chem. Phys. 1986,85,4997. (d) Isaacson, A. D.; Truhlar, D. G.; Rai, S. N.; Steckler, R.; Hancock, G. C.; Garrett, B. C.; Redmon, M. J. Computer Phys. Commun. 1987,47,91. (e) Steckler, R.; Dykema, K. J.; Brown, F. B.; Hancock, H. C.; Truhlar, D. G.; Valencich, T. J. Chem. Phys. 1987,87, 7024. (f) Joseph, T.; Steckler, R.; Truhlar, D. G. J . Chem. Phys. 1987,87, 7036. (g) Garrett, B. C.; Redmon, M. J.; Steckler, R.; Truhlar, D. G.; Baldridge, K. K.; Bartol, D.; Schmidt, M. W.; Gordon, M. S. J . Phys. Chem. 1988, 92, 1476. (h) Baldridge, K. K.; Gordon, M. S.; Steckler, R.; Truhlar, D. G. J . Phys. Chem. 1989.93.5107. (i) Garrett, B. C.; Joseph, T.; Truong, T. N.; Truhlar, D. G. Chem. Phys. 1989, 136, 271. 0’) Tucker, S. C.; Truhlar, D. G. J. Am. Chem. SOC.1990, 112, 3338, 3347. (k) Lu, D.-h.; Maurice, D.; Truhlar, D. G. J . Am. Chem. SOC.1990, 112, 6206. (1) Zhao, X. G.; Tucker, S. C.; Truhlar, D. G. J . Am. Chem. SOC.1991,113,826. See also (m) Kreevoy, M. M.; Ostovie, D.; Truhlar, D. G.; Garrett, B. C. J. Phys. Chem. 1986, 90, 3766. (n) Kim, Y.; Truhlar, D. G.; Kreevoy, M. M. J . Am. Chem. SOC.,in press. (0)Lu, D.-h.; Truong, T. N.; Melissas, V. S.; Lynch, G. C.; Liv, Y.-P.; Garrett, B. C.; Steckler, R.; Isaacson, A. D.; Rai, S.N.; Nancock, G. C.; Lauderdale, J. G.; Joseph, T.; Truhlar, D. G. Comput. Phys. Commun., in press. (9) Schatz, G. C.; Amaee, B.; Connor, J. N. L. J . Chem. Phys. 1990,93, 5544. (10) (a) Bondi, D. K.; Connor, J. N. L.; Manz, J.; Romelt, J. Mol. Phys. 1983, 50, 467; (b) Bondi, D. K.; Connor, J. N. L.; Garrett, B. C.; Truhlar, D. G. J . Chem. Phys. 1983, 78, 5981.

0022-3654191 12095-10374%02.50/0 , 0 199 1 American Chemical Society I

,

H and D Atom Transfer Reactions

The Journal of Physical Chemistry, Vol. 95, No. 25, 1991 10375

and sf-POLCI.'2 The sf-POLCI results are a global fit to a previous semiglobal surface6 called "scaled" in ref 6; it results from semiempirical adjustment of a surface fit to polarization configuration interaction calculations augmented by dispersion terms. The quantum mechanical calculations were carried out by the coupled-states distorted wave (CSDW) method for all six cases ( H and D on three surfaces). For two of the cases, quantum mechanical calculations were presented earlier') using a coupledchannel-hyperspherical J-shifted (CCH/JS) algorithm, which may be more accurate. The CSDW rate coefficients are much larger than those predicted by trajectory calculation^,'^ in agreement with our earlier conclusions6 that tunneling dominates the thermal rate process for this reaction. Schatz et a1.9 compared VTST/LCG calculations6for the scaled surface to CSDW calculations for surface sf-POLCI and found discrepancies up to a factor of 2.5 at room temperature, much larger than we have found in previous tests.'-) There are several possible reasons why such a larger error might occur in this case: (i) As already mentioned, VTST/ST calculations with an uncoupled treatment of the vibrational-rotational modes may be less accurate for some new systems than for those previously studied. The scaled and sf-POLCI surfaces have very anharmonic bending potentials, and this may cause greater mode coupling and larger errors. (ii) We have improved our treatment of stretch anharmonicity since the calculations of ref 6 to which Schatz et al.9 compared. In particular, those calculations modeled the stretching vibrations as Morse oscillator^,'^^'^ whereas we have since then demonstrated that systematically more accurate results are obtained by calculating stretching energy levels by the WentzelKramers-Brillouin (WKB)IdJ8 method. (iii) The discrepancies in ref 9 may actually be an artifact of performing quantum dynamics calculations on a different surface than had been used for the semiclassical calculations. (iv) The discrepancies of ref 9 may be due more to errors in the CSDW calculations rather than to errors in the VTST/ST calculations. For example, for the two cases where both CSDW and CCH-JS results are available, the maximum factor between them is 2.0 and the average factor is 1.6. A recent comparison15 of CSDW tunneling probabilities for the 0 H D OH D and 0 H D OD + H reactions to converged quantum mechanical calculations for the same potential energy surface showed better agreement, but the maximum and average factors of deviation were found to be 1.4 and 1.2, respectively, still not insignificant. The purpose of the present paper is to present new VTST/ST calculations for all six cases for which Schatz et al. presented CSDW calculations, on precisely the same potential energy surfaces, and to use these to learn more about the accuracy of semiclassical tunneling methods for H and D atom transfers. Calculations are presented both for the uncoupled treatment of vibrational-rotational modes employed previously and also for a coupled treatment of the bending modes using the centrifugal oscillator19 formulation. The present tests are of special interest for those interested in organic and biochemical problems since proton and hydride transfers in those cases correspond to the heavy-light-heavy mass combination for which the present reaction serves as the most well studied prototype.

+

-

+

+

-

2. Calculations The calculations for uncoupled vibrational-rotational modes were performed using our standard methods presented in detail (11) Persky, A.; Kornweitz, H. J . Phys. Chem. 1987, 91, 5496. (12) Schatz, G. C.; Amaee, B.; Connor, J. N. L. J . Phys. Chem. 1988,92, 3190. ....

(13) Sun, Q.; Bowman, J. M.; Schatz, G. C.; Sharp, J. R.; Connor, J. N . L. J . Chem. Phys. 1990, 92, 1677. (14) Kornweitz, H.; Broida, M.; Persky, A. J . Phys. Chem. 1989, 93, 251. Schatz, G. C.; Amaee, B.; Connor, J. N. L. J . Chem. Phys. 1990,92,4893. (15) Halvick, P.; Zhao, M.; Truhlar, D. G.; Schwenke, D. W.: Kouri. D. J. J . Chem. SOC.,Faraday Trans. 1990,86, 1705. (16) Morse, P. M. Phys. Reo. 1929, 34, 57. (17) Garrett, B. C.;Truhlar, D. G. J . Phys. Chem. 1979, 83, 1052, 1079. (18) Laudau, L. D.; Lifshitz, E. M. Quontum Mechanics, 2nd ed.; Pergamon Press: Oxford, U.K., 1965; p 163. (19) Natanson, G. A. J . Chem. Phys. 1990, 93, 6589.

e l ~ e w h e r e . ' ~ * ~ *In ~ ,particular, ~ " - ~ ~ the VTST calculations were carried out by canonical variational theory (CVT),IS'*~'the ST calculations were carried out by both the least-action ground-state (LAG)4 and large-curvature ground-state version-3 (LCF3)7,8i,80,22 methods, the reference path for defining generalized transition states was the collinear minimum energy path through mass-scaled coordinate^,'^,^^ stretch anharmonicity was treated by the WKB method,Id and, in the first set of calculations, bending anharmonicity was included in an uncoupled-mode approximation in curvilinear coordinates.6,21 All calculations include the ratio of electronic partition functions a t the transition state and at reactants, where the former includes only the 22ground state and the latter includes the 2P3/2 and 2Pl/2states of C1. In the second set of calculations, the energy levels of the two bending degrees of freedom were calculated by a WKB approximation to the centrifugal oscillator form~lation,'~ and the partition functions for these modes were calculated by numerical summation of Boltzmann factors. The centrifugal oscillator includes the coupling between the two bending degrees of freedom. The WKB quantization equation isI9

where 0 = I@ - TI, O< and O> are classical turning points, @ is the Cl-H-Cl bond angle, eUfiis the quantized energy, u2 is the bending number, K is the vibrational angular momentum quantum number, AVb is the bending potential, Goo is the Wilson G matrix element defined and h is Planck's constant. Since 0 is intrinsically nonnegative, O< is 0 for K = 0 and is >O for K # 0. The partition function, with zero of energy at the potential for the linear geometry, is

where kB is Boltzmann's constant and T is temperature. This contrasts to

(3) where Q& is computed over the -O> to +e> interval in the way in the uncoupled treatment. In both sets of calculations the bending potential was fit to a quadratic-quartic potential c ~ r v e ~ ' , ~ ~

Avb(@,s)= 72Fds)(@- *)' + h A o ( $ ) ( @- T)4

(4)

where s is distance along the minimum energy path, Fo is the harmonic force constant, and A# is the quartic force constant. In all cases the harmonic and quartic force constants were equated to a2V/a@2and a4V/aa4, respectively, where the derivatives are evaluated with H-Cl distances fmed at their values on the collinear minimum energy path. Two of the potential energy surfaces, the BCMR and PK3, have collinear minimum energy paths with positive bending force constants. In these cases the energy levels of the bend potential for the uncoupled treatment were determined by a perturbation-variation method that is explained in detail el~ewhere,~) and those for the coupled treatment were evaluated by Gauss-Legendre quadrature of the WKB quantization condition (1). For the sf-POLCI surface, we still base the calculations on a collinear reference path, as previously,6 even though F&) is negative. Since A&) is positive the bend potential for the uncoupled treatment has a double-minimum character. The (20) Truhlar, D. G.; Kuppermann A. J . Am. Chem. SOC.1970,93, 1840. (21) Garrett, B. C.; Truhlar, D. G.; Grev, R. S.; Magnuson, A. W. J . Phys. Chem. 1980, 84, 1730. (22) Garrett, B. C.; Abusalbi, N.; Kouri, D. J.; Truhlar, D. G. J . Chem. Phys. 1985,83, 2252. (23) (a) Truhlar, D. G.J . Mol. Spectrosc. 1971, 38, 415. (b) Garrett, B. C.; Truhlar, D. G.J . Phys. Chem. 1979,83, 1915. (c) Garrett, B. C.; Truhlar, D. G.J . Am. Chem. SOC.1979, 101, 4534.

Garrett and Truhlar

10376 The Journal of Physical Chemistry, Vol. 95, No. 25, I991

+

TABLE I: Ratios of Rate Constants for the CI HCI Reaction on the BCMR Surface X 312.5 K 368.2 K 423.2 K Reference = CSDW TST" 0.16 0.19 0.22 CVT" 0.16 0.19 0.22 CVT/LCG3" 0.60 0.53 0.52 CVT/LAG" 0.58 0.52 0.51 CCH-JS 0.77 0.62 0.62

Reference = CCI-I/JS TST" CVT" CVT/LCG3" CVT/LAG" CSDW

0.21 0.21 0.78 0.76 1.30

0.30 0.30 0.85 0.83 1.60

0.35 0.35 0.82 0.81 1.60

Uncoupled bending vibrations. TABLE 11: Ratios of Rate Constants for the CI + HCI Reaction on the PK3 Surface X 312.5 K 368.2 K 423.2 K

Reference = CSDW TST" CVT" CVT/LCG3" CVT/LAG" CCH-JS

0.08 0.08 0.50 0.49 0.73

0.10 0.10 0.44 0.43 0.62

0.12 0.12 0.40 0.39 0.50

Reference = CCH-JS TST" CVT" CVT/LCG3" CVT/LAG" CSDW

0.1 1 0.1 1 0.69 0.67 1.38

0.16 0.16 0.7 1 0.70 1.62

0.24 0.24 0.80 0.78 2.00

Uncoupled bending vibrations. energy levels of the double-minimum potential were found by a uniform semiclassical m e t h ' ~ d . For ~ ~ the coupled treatment in this case the bend potential has a single minimum since 0 is intrinsically positive. Again we used Gauss-Legendre quadrature of the WKB quantization condition (1). We will report four rate constants: TST, conventional transition state theory without tunneling; CVT, canonical variational theory without tunneling; CVT/LCG3, canonical variational theory with large-curvature ground-state version-3 tunneling; CVT/LAG, canonical variational theory with least-action ground-state tunneling. One reason for presenting the CVT results is that they can be used to gauge the extent to which tunneling dominates the rate constant in each case. In particular, the CVT results are our best estimate (improved canonical variational transition state theory2' leads to identical results with at least three significant figures for all cases considered here) of what the rate constant would be if there were no tunneling. 3. Results and Discussion 3.1. Comparison of Uncoupled Semiclassical and Best Available Quantal Rate Constants. First we consider the H-transfer reaction on the BCMR and PK3 surfaces. We present these cases first because they are the two cases for which both CSDW and CCH-JS results are available. Tables I and I1 each present two sets of ratios, defined by RX(r ) = kX( 79 /Pf(7,) (5) where X may denote TST, CVT, CVT/LCG3, CVT/LAG, CSDW, or CCH-JS, and the reference result used in the denominator is taken in turn as CSDW and CCH-JS. If we take the deviation of the CSDW and CCH-JS results from each other as a measure of the uncertainty in the best available quantun dynamics values, which is quite reasonable, then we see that in both cases the CVT/LCG3 results agree with the CCH-JS results within their estimated reliability. This is very encouraging. (24) Connor, J . N. L. Chem. Phys. Lett. 1969, 4, 419.

TABLE 111: Ratios of Rate Constants for the CI the BCMR Surface (Reference = CSDW) X 312.5 K 368.2 K TST" 0.33 0.36 CVT" 0.33 0.36 CVT/LCG3" 1.04 0.88 CVT/LAG" 1.03 0.88

+ M31 Reaction on 423.2 K 0.39 0.39 0.81 0.81

Uncoupled bending vibrations. TABLE IV: Ratios of Rate Constants for the CI + DCI Reaction on the PK3 Surface (Reference = CSDW) X 312.5 K 368.2 K 423.2 K TST" 0.25 0.33 0.39 CVT" 0.25 0.33 0.39 CVT/LCG3" 1.07 1.03 0.99 CVT/LAG" 1.07 1.03 0.99 a

Uncoupled bending vibrations.

TABLE V Ratios of Rate Constants for the CI the sf-POLCI Surface (Reference = CDSW) X 312.5 K 368.2 K

+ HCI Reaction on 423.2 K

Uncoupled Bending Vibrations TST CVT CVT/LCG3 CVT/LAG

1.09 0.28 0.91 0.91

0.84 0.24 0.61 0.61

0.68 0.21 0.45 0.45

Centrifugal Oscillator TST CVT CVT/LCG3 CVT/LAG

0.47 0.15 0.44 0.42

0.39 0.14 0.31 0.3 1

TABLE M: Ratios of Rate Constants for the C1 the sf-POLCI Surface (Reference = CSDW)

X

312.5 K

0.33 0.13 0.24 0.24

+ DCI Reaction on

368.2 K

423.2 K

Uncoupled Bending Vibrations TST CVT CVT/LCG3 CVT/LAG

5.05 1.84 3.09 3.08

3.24 1.23 1.84 1.83

2.34 0.90 1.25 1.25

Centrifugal Oscillator TST CVT CVT/LCG3 CVT/LAG

1.45 0.86 1.35 1.32

1.08 0.61 0.87 0.86

0.87 0.47 0.63 0.62

Tables I11 and IV present similar results for D atom transfer on these two potential energy surfaces. These tables show extremely good agreement between the CVT/LCG3, CVT/LAG, and best available quantum dynamics calculations. The average ratio between the CVT/LAG and best available quantum dynamics results, always putting the smaller value in the denominator, is 1.08, which is excellent agreement for a comparison of semiclassical and quantal calculations in a case dominated by tunneling. Again we conclude that the semiclassical tunneling calculations agree with the quantum mechanical calculations within the reliability limits of the latter. Next we performed calculations with the sf-POLCI and scaled potential energy surfaces. The results agreed to about five significant figures. This eliminates possible problem (iii) of the list in the Introduction. Then we compared sf-POLCI results to the CSDW calculations which are the best available quantal results for this surface. The comparison is shown in the top half of Table V. First we see that the new results are about 35% larger than the previous ones6 due to the improved treatment of stretch anharmonicity (item ii in the list in the Introduction), and hence they agree better than the previously published ones6 with the CSDW calculations. In fact the agreement is to within an amazingly good 10% at 312.5 K, although there is still a factor of 2.2 discrepancy at 423.2 K. Since Tables I-IV show much better agreement for the first two surfaces considered, and the

The Journal of Physical Chemistry, Vol. 95, No. 25, 1991 10377

H and D Atom Transfer Reactions 10

I " " " : " " '

'

'

I I

8

" '

"

'

~

( b ) s = fo.24 a. .

i

I

e:",(

120 130 140 150 160 170 180 @ (degrees)

lo

-

*

8

,

3 :

: Ot -21 .

the potential energy surface over the whole well range. Since the uncoupled vibrational zero point energy is negative or only slightly positive at the critical geometries, and the ground state dominates 0 at the temperatures considered here, we conclude that the fit does not lead to any significant error. Next we examine another possible source of error in the treatment of the bend used in this section, namely the uncoupled approximation. 3.2. Semiclassical Calculations with the Centrifugal Oscillator Approximation. Calculations with the centrifugal oscillator formalism were carried out including 12 terms in eq 2, namely u2 = O , K = 0-5;v2 = 1, K = 0-3;andu2 = 2, K = 0,l. This was found sufficient to converge the rate constants to 4-5 significant figures. First we repeated the calculations for the BCMR and PK3 surfaces for which the bending force constant is positive. We did not expect a big change in these cases, and we did not find one. In all cases (both surfaces and both isotopes for the four methods and three temperatures in the tables), the calculated rate constants decreased but by only 1.3-2.5% for the BCMR surface and 0 5 1 . 4 % for the PK3 surface. The effect of coupling the bending vibrations was much larger for the sf-POLCI surface. Again all the calculated rate coefficients decreased, but now by the much more significant factors of 2.1-3.5 for conventional transition state theory, 1.6-2.1 for canonical variational theory, and 1.9-2.3 for calculations including tunneling by the LAG or LCG3 method. Since the effect of coupling the bending vibrations is sizeable, we recalculated the ratios to the best available quantum results and these are shown in the lower halves of Tables V and VI. Consider first Table V for the H transfer. The approximate theory now predicts too low a result by about a factor of 2 at the lowest temperature, which is consistent with our previous experience for difficult cases. The discrepancy from the best available quantum results is exacerbated though at 423.2 K. Next consider the results for D atom transfer in Table VI. The large discrepancy at the lowest temperature is ameliorated, and the CVT/LAG and CVT/LCG3 results now underestimate rather than overestimate the best available quantum ones at 423.2 K. The summary lesson of Tables V and VI, we believe, is that changes in the predictions due to coupling the bending vibrations are comparable in size to the discrepancies from the best available quantum results. This indicates that the difficulty of making accurate semiclassical predictions for this surface is associated more with the unusually large anharmonicity than with the nature of the tunneling process per se. 3.3. Interpretation of Semiclassical Results. One of the advantages of semiclassical calculations, such as the VTST/ST calculations presented here, is that they lead to a clear physical picture of the chemical reaction. In cases where the VTST/ST calculations are in excellent quantitative agreement with independent quantal calculations, as for the BCMR and PK3 potential energy functions considered in this paper, the semiclassical picture acquires added persuasiveness. First of all we define a "classical reaction percentage" as

-sf-POLCI A @ = 1Odeg

i, - -A@=20deg

1

I

I

\ t -\

'\\

' ' . ' . . . . ' , ,

,

, ' ,

,

C

'

,. ,

, ' .

, .

, I

120 130 140 150 160 170 180 Q, (degrees)

Figure 1. Bending potential as a function of CI-H-CI bond angle from original potential surface (solid curve) as compared to three harmonicquartic fit (dashed curves and diamonds): (a) at the saddle point; (b) at the variational transition state.

main difference between those cases and the present one is the bending anharmonicity, we think it is most likely that-if the fault in the comparison under discussion lies on the semiclassical side-the discrepancy is caused by the uncoupled treatment of the very anharmonic bend potential more than by the semiclassical tunneling method. Another possibility of course is that the discrepancy results mainly from errors in the quantum calculations. Finally, we consider the deuterium atom transfer on the sfPOLCI surface. Here the results shown in the top half of Table VI are again somewhat atypical of what we have observed in other cases, with the VTST calculations overestimating the best available quantum ones even at the lowest temperature when tunneling is not included. The error decreases with increasing temperature, which would be consistent with VTST overestimating the bending contribution to the zero point energy. In this case it appears that the discrepancy from the best available quantal results error at 312.5 K is too large to be attributed the approximations in those calculations. One possible source of error for the sf-POLCI surface is the quadraticquartic fit to the bend potential. As stated in section 2, the coefficients were obtained from the second and fourth derivatives at the collinear geometry, 9 = 180'. One way to evaluate these derivatives is, for each point along the collinear minimum energy path, to fit eq 1 to the bend potential at two angles, 180' - AQ and 180' - 2AQ, where AQ is small. We test the sensitivity to AQ for the C1 + HCl case. First we used AQ = 1' and compared the results to AQ = 2'; the rate coefficients for these two choices agree to better than three significant figures, confirming that these bend angles are indeed close enough to 180' to yield the true Taylor series. Figure 1 shows in fact that both AQ = 1' and AQ = 10' yield fits that agree quantitatively with

kCVT/CAG

kcla,,(%)

=

kCVr/LAC x

100%

(6)

The numerator in eq 2 is the rate constant we would calculate if we neglected quantum effects on the reaction coordinate (the CAG transmission coefficient is precisely defined by equations in previous work2'). Dividing this by our final semiclassical rate constant that includes these effects gives an upper bound to the fraction of the final rate calculation that comes from overbarrier processes. (It is a bound rather than a best estimateSbbecause, in the full semiclassical calculation, some of the tunneling is cancelled by nonclassical reflection.) The results are given in Table VII. Notice, for example that the percentage of reaction proceeding classically on the BCMR surface is only slightly higher for D than for H. This is a very clear indication of the inadequacy of expectations based on isotope-independent effective one-dimensional barriers and mass scaling, which would predict a bigger difference. We also see that in five of the six cases studied here

Garrett and Truhlar

10378 The Journal of Physical Chemistry, Vol. 95, No. 25, 1991 4

TABLE VII: Classical Reaction Percentages potential function BCMR

T, K 3 12.5 423.2 312.5 423.2 312.5 423.2

PK3 sf-POLCI

CI

kc,,,,, 72 CI

+ HCI

28 43 16 30 31 (37)" 46 (52)

h

i

+ DCI

8 5 3

32 49 24 40 60 (65) 72 (76)

a

F3:

a Where only one result is shown, the results for the uncoupled and coupled treatment of the bend agree to the number of figures shown; otherwise the result calculated by the uncoupled treatment is given first (followed by the result for the coupled treatment in parentheses).

less than half the reaction proceeds by the overbarrier mechanism, but the overbarrier contributions are not neligible either. Thus we examine the overbarrier dynamic bottlenecks in more detail, and various properties of these bottlenecks are shown in Table VIII. (All results in this table for the BCMR and PK3 surfaces are based on the uncoupled-mode treatment. Both sets of results are given for the sf-POLCI surface.) In canonical variational theory the dynamical bottleneck is located at the maximum of the free energy of activation profile, which is the generalized free energy of activation as a function of the reaction coordinate s.7,21-23c We have found previously23c that, even at finite temperatures, the factors controlling the location of the dynamical bottleneck may often be understood qualitatively in terms of its shape at absolute zero temperature, which is the same as the shape of the vibrationally adiabatic ground-state potential curve The latter is given for an atom-diatom reaction in the uncoupled-mode and coupled-mode approximations21.23c

e.

V?(s) = VMEP(s) + 4?r(s) + 2 t t e n d t ~ )

(7)

andI9 t t s ) = VMEPts)

+ E$r(S) + t 2 o t ~ )

(8) respectively, where VMEpis the Bomappenheimer potential along the minimum energy path (MEP) or, for cases like the sf-POLCI potential function, where we use the collinear MEP as the reference function even though it is not the overall MEP, along the collinear reference path, is the ground-state (Le., zero point) vibrational energ of the bound stretching mode of the generalized transition state, tkndis the ground-state energy of one of the doubly degenerate bending modes relative to the value of the potential for the reference path, and tFo is the ground-state energy of the 2-D centrifugal oscillator, again relative to the values of the potential on the reference path. Columns 6, 8, and 9 of Table VI11 list two subsets of the terms in as well as itself for the maximum of the generalized free energy of activation curves at absolute zero (which is the same as the maximum of itself) and for the lowest and highest temperatures considered in this paper. Table VI11 also gives the two nearest-neighbor bond distances, RHa, and RHClb, and the harmonic stretching frequency

B

e

e

e

2

3

2

4

H-Cl distance (bohr)

Figure 2. Contour plots and interpretation of the dynamics for the C1 HCI and CI DCI reactions at 312.5 K. The plots correspond to collinear geometries. Axes are the H-CI bond distances, and the contour values correspond to 5, 10, 20, 40, and 80 kcal/mol. The straight lines at 45" angles to the axes are the representative tunneling paths for H and D. The point labeled % is the saddle point.

+

+

w , , ~for the reference collinear geometry at the maximum of the free energy of activation profile. First consider the results for the BCMR and PK3 potentials. In both these cases the variational transition states for both isotopic transfers are located at the saddle point, and the difference between the bending zero point energies for the Cl-H-C1 and Cl-D-C1 transition states is about 0.6 kcal/mol. For the sf-POLCI surface, the saddle point has equal H-Cl distances of 2.78 ao, a classical potential (VMEp)of 9.01 kcal/mol, and a harmonic stretching frequency of 336 cm-I (all these quantities are the same for H and D transfer at the saddle point). Table VI11 shows that the dynamical bottlenecks occur quite a bit earlier for both H and D where the classical potential is much lower, but the stretching frequency is much higher. For the D-transfer case, when the vibrations are uncoupled the last term of eq 7 is negative at the dynamical bottlenecks, and is lower than VMEp+ .e,: The resulting vibrationally adiabatic barrier height, relative to reactants at classical equilibrium, is 0.9 kcal/mol lower for Cl-D-Cl than for Cl-H-Cl. On the basis of only the uncoupled results one might have concluded that the discrepancy from the best available quantum results in the top half of Table VI is due to our uncoupled mode treatment underestimating the bending zero point requirements for Cl-D-Cl. However, the coupled-mode treatment leads to the conclusion that is 1.8 kcal/mol lower for D transfer than for H transfer, which casts doubt on this interpretation. Further work on bending anharmonicity, and perhaps on bend-stretch and vibration-rotation coupling, will be required before these points are completely understood. The other point of interpretation that we wish to emphasize is the nature of the tunneling process. The LAG approximation optimizes the shape of the tunneling path for each pair of tunneling path termini (one in the reactant valley and one in the product valley) that is considered, whereas the LCG3 approximation treats it as a straight line in isoinertial coordinatesz5for each such pair.

e

e

TABLE VIII: Properties of Variational Transition States

BCMR PK3 sf-POLCI

+

CI HCI CI + DC1 C1 + HCI CI DC1 CI HCI

+ +

CI

+ DCI

0.0-423.2" 0.0-423.2" 0.0-423.2" 0.0-423.2" O.Ob

0.0' 423.2b 423.2c

O.Ob 0.0' 423.2b 423.2'

2.77 2.77 2.8 1 2.81 3.26 3.26 3.26 3.26 3.17 3.17 3.16 3.16

2.77 2.77 2.81 2.81 2.49 2.49 2.49 2.49 2.52 2.52 2.52 2.52

8.55 8.55 8.54 8.54 7.26 7.26 7.26 7.26 7.64 7.64 7.65 7.65

344 344 347 347 2090 2090 2089 2089 1422 1422 1421 1421

9.04 9.04 9.03 9.03 9.89 9.89 9.89 9.89 9.43 9.43 9.43 9.43

10.57 10.12 11.09 10.49 10.10 10.37 10.10 10.37 9.21 9.53 9.21 9.53

"In these cases the variational transition states are located at the saddle point over this whole temperature range. buncoupled treatment of bending vibrations. Coupled bending vibrations.

H and D Atom Transfer Reactions

The Journal of Physical Chemistry, Vol. 95, No. 25, I991

+

TABLE IX: Ratios of Rate Constants for the Collinear C1 HCI and CI + DCI Reactions to Accurate Quantal Ones on the BCMR Surface 300 K 400 K 600 K CI HCI CVT/LCG3 0.80 0.97 1.19 CVT/LAG 0.76 0.93 1.15

+

CI CVT/LCG3 CVT/LAG

0.87 0.86

+ DCI 0.81 0.80

1.01

1.oo

10379

whether the ‘good agreement in Table I [is] fortuitous due to [possible] cancellation of errors in the VTST treatment of the collinear and bending dynamics”. Such checks are presented in Table IX where we give ratios of CVT/LCG3 and CVT/LAG rate constants to the accurate quantum collinear rate constants reported in ref lob. The VTST calculations with large-curvature-tunneling transmission coefficients are accurate within 25% for all cases in the table. We conclude that good results in three dimensions do not result from cancellation of errors in treating the bend with errors in the dynamics.

4. Summary The fact that the LAG and LCG3 cases agree so well in the We have tested semiclassical calculations of rate constants for present case is a direct consequence of the fact that the optimized hydrogen atom and deuterium atom transfer between two heavy tunneling paths and the straight-line ones are either identical or atoms against the best available quantum mechanical calculations almost identical at all important tunneling energies for all six cases. for six cases of three-dimensional reactions and also for two cases This is what we expect for this mass combination. Since the of collinear reactions, and in most cases we find excellent present reaction is symmetric, Le., it is a “degenerate” reaction agreement. The semiclassical calculations are used to obtain a in which reactants and products are the same, the straight-line physical picture of the dynamical bottlenecks for the overbarrier tunneling paths have a very simple form when transformed to processes as well as for the representative tunneling paths for the bond-length coordinates; in particular, they correspond to the H most dominant underbarrier processes. moving between fixed chlorines. Tables I-IV indicate that this Of course one of the purposes of tests like those reported here is an excellent approximation. is to anticipate the reliability of approximate methods for cases Figure 2 shows examples of the LCG3 tunneling paths for H where they cannot be tested against accurate quantum dynamics and D transfer on the PK3 surface. The axes of the plot are the because accurate calculations are too difficult. In such cases one RHCI,and RHClb nearest-neighbor bond distances. For each isotope may sometimes check against experimental results, but such checks we define a “representative tunneling energy” for a temperature are complicated by the possibility of errors in the potential energy T = 312.5 K as the energy at which the integrand F ( E ) exp(-@E) function, which may either add to or cancel errors in the dynamics. peaks in the e x p r e s ~ i o n ~for ~ *the ~ thermally averaged LCG3 Nevertheless there are a number of successful applications of transmission coefficient. Here is the ground-state tunneling variational transition state theory and semiclassical tunneling probability, E is total energy, and 0 is E / R T , where R is the gas calculations to reactions with four or more atoms,8J6and the reader constant. This yields E = 9.88 kcal/mol for H transfer and 9.54 kcal/mol for D transfer. In the LCG3 a p p r o x i m a t i ~ n , at ~ ~ ~ ” ~ may ~ ~ wish to consult these to learn more about the validity of the methods for more complicated cases. a given value of E, tunneling can initiate and terminate at any value of the reaction coordinate for which c ( s ) C E . We define Acknowledgment. We are grateful to Gillian Lynch for helpful the “representative tunneling path” as the one with largest tunassistance. This work was supported in part by the US. Deneling probability at the representative tunneling energy for the partment of Energy, Office of Basic Energy Sciences, through temperature considered, and this is the path we plot for each grant no. DE-FG02-86ER13579, Pacific Northwest Laboratory isotope in Figure 2. The figure shows that these tunneling paths is operated for the U S . Department of Energy by Battelle correspond to much larger Cl-Cl distances than the saddle point. Memorial Institute under contract DE-AC06-76RLO 1830. In particular, the saddle point Cl-Cl distance on the PK3 surface Registry No. D2, 7782-39-0. is 5.62 ao, whereas the Cl-Cl distances for the representative tunneling paths in Figure 1 are 5.90 a. for H transfer and 5.74 a. for D transfer. The highest energy points on the representative (26) See also, e.g.: (a) Truhlar, D. G.; Isaacson, A. D. J . Chem. Phys. paths for H and D are 10.22 and 8.88 kcal/mol, considerably 1982, 77, 3516. (b) Bicerano, J.; Schaefer, H. F. 111; Miller, W. H. J . Am. higher than the saddle point value of 8.54 kcal/mol. For the Chem. SOC.1983, 105, 2550. (c) Carrington, T. Jr.; Hubbard, L. M.; BCMR surface the corresponding values for the Cl-C1 distances Schaefer, H. F. 111; Miller, W. H. J. Chem. Phys. 1984,80, 4347. (d) Houk, K. N.; Rondan, N. G. J . Am. Chem. SOC.1984,106,4293. (e) Houk, K. N.; are 5.54 a. for the saddle point, 5.85 a. for the representative Rondan, N. G.;Mareda, J. Tefrahedron 1985,41, 1555. (0 Isaacson, A. D.; tunneling path for H at 312.5 K, and 5.69 a. for the representative Sund, M. T.; Rai, S.N.; Truhlar, D. G. J . Chem. Phys. 1985,82, 1338. (8) tunneling path for D at 3 12.5 K, and the corresponding values Viswanathan, R.; Raff, L. M. Thompson, D. L. J. Chem. Phys. 1985, 82, of the highest energies are 10.51 and 9.05 for H and D vs 8.55 3083. (h) Hase, W. L.; Duchovic, R. J. J . Chem. Phys. 1985,83, 3448. (i) Doubleday, C.; McIver, J.; Page, M.; Zielinski, T. J . Am. Chem. SOC.1985, for the saddle point. We expect similar deviations of the optimum 107, 5800. (j) Wardlaw, D. M.; Marcus, R. A. J . Chem. Phys. 1985, 83, tunneling path from the minimum energy path to occur for 3462. (k) Wardlaw, D. M.; Marcus, R. A. J . Phys. Chem. 1986, 90, 5383. proton-transfer reactions as well. (1) Doubleday, C.; McIver, J. W. Jr.; Page, M. J Phys. Chem. 1988,92,4367. 3.4. Collinear Collisions. One of the referees requested that (m) Klippenstein, S . J.; Khundkar, L. R.; Zewail, A. H.; Marcus, R. A. J. Chem. Phys. 1988,89,4761. (n) Page, M.; Lin, M. C.; He, Y.; Choudhury, we carry out further checks for collinear collisions to determine T. K. J. Phys. Chem. 1989, 93, 4404. (0)Blake, J. F.; Wierschke, S. G.; (25) As discussed in the second article of ref 17, there are an infinite number of coordinate systems in which the kinetic energy has no cross terms

and the reduced mass is the same for every square term, but they are all related to one another by axis rotations. Thus path shapes are the same in all such coordinate systems, and we refer to them collectively as isoinertial coordinate systems.

Jorgensen, W. L. J . Am. Chem. SOC.1989, 111, 1919. (p) Garrett, B. C.; Koszykowski, M. L.; Melius, C. F.; Page, M. J . Phys. Chem. 1990,9#, 7096. (9) Truong, T. N.; Truhlar, D. G. J . Chem. Phys. 1990,93, 1761. (r) Kaseki, S.; Gordon, M. S.J . Phys. Chem. 1989, 93, 118. (s) Truhlar, D. G.; Gordon, M. S. Science 1990, 249, 491. (t) Gonzalez-Lafont, A,; Truong, T. N.; Truhlar, D. G.J . Phys. Chem. 1991,95,4618. (u) Truong, T. N.; McCammon, J. A. J. Am. Chem. 1991, 113, 7504.