Transition state theory and Eckart's tunnelling factor: a good approximation for the calculation of bimolecular rate constants?

Chemical Physics ELSEVIER

Chemical Physics 183 (1994) 95-100

Transition state theory and Eckart’s tunneling factor: a good approximation for the calculation of bimolecular rate constants? J. Espinosa-Garcia

*, F.J. Olivares de1 Valle, J.C. Corchado

Departamento de Quimica Fisica, Universidad de fitremadura,

06071 Bahjoz,

Received 30 August 1993; in final form 14 February

Spain

1994

Abstract Bimolecular rate constants for polyatomic systems are commonly calculated by using non-variational transition state theory with Eckart’s tunneling correction. The reason is the high computational cost of a larger amount of information on the potential energy surface of polyatomic systems. The validity of such an approximation is here tested using the widely studied H + Hz, D + H2 and H + D2 reactions. With a large basis set, electron correlation, and spin projection for the free radicals to calculate the activation barrier, and the one-dimensional Eckart potential function fitted to the position and second derivative of the adiabatic curve at the saddle point to calculate Eckart’s tunneling factor, there was good agreement with the experimental and accurate quantum dynamics values for 7’> 300 K.

1. Introduction The complete construction of a potential energy surface (PES) is even today limited to a very few molecular systems, such as atom-atom or atom-diatom. Knowledge of the PES would allow the calculation of the rate constant, taking into account the effects of recrossing on the bottleneck to reaction [ I], the variational localization of the transition state [ 11, and the effects of the curvature of the reaction path [ 21 on the quantum mechanical tunnel effect. When one has to deal with larger polyatomic systems, however, the complete construction of the PES becomes prohibitive as the number of degrees of freedom rises. One has, therefore, to introduce approximations in order to calculate the rate constants. The conditions that will have to be satisfied by these approx-

* Corresponding

author.

0301-0104/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved SSDIO301-0104 (94)00069-M

imations are: ( 1) simplicity, so that only partial knowledge of the PES is necessary; (2) computational economy, so that it will be possible to use large basis sets; and (3) reproducibility of the results, to give an adequate approximation to the best theoretical and experimental values. The simplest level at which to study the rate constants in polyatomic systems is that of non-variational transition state theory (TST) [ 3-51, which only requires knowing the electronic energies and frequencies of the reactants, and the saddle point on the PES. The tunneling effect must be taken into account to correct the constant for quantum effects on the reaction coordinate. A straightforward way to include this effect is by numerical integration of the wave propagating in the one-dimensional Eckart potential [6] which has been fitted to the position and second derivative of the adiabatic curve at the saddle point. Therefore, the resulting tunneling factor depends not only on the

96

.I. E.rpinosa-Garth

et al. /Chemical

imaginary frequency that characterizes the saddle point, but also on the reaction energy and barrier height. This approach is not “exact”, but involves a onedimensional function widely used for hydrogen abstraction reactions [ 7-101 giving results that agree qualitatively with experiment [ I 11. These two approximations, TST and Eckart’s factor, satisfy conditions 1 and 2. One may consequently opt for the use of these two approximations to calculate rate constants in polyatomic systems. The calculations that will be presented here have the purpose of evaluating the lowest level of computation that must be included in the cited approximations to give a reasonable response to condition 3. As a test of our results, we use the most accurately known PES for any chemical reaction, H2 + H, and two isotopic variants, D+H2 and H +D,, for which there is a great deal of theoretical and experimental information available.

2. Method Molecular calculations were carried out using the GAUSSIAN 90 system of programs [ 12 1. Stationary point geometries were optimized to restricted (R) and unrestricted (U) second-order Moller-Plesset perturbation theory, MP2 [ 131, with full electron correlation using the 6-31G( p) basis set of Pople [ 141. We then used three levels of calculation. Level I. We make a single point calculation to fourthorder Moller-Plesset perturbation theory, MP4, with a frozen core (FC) approximation and single, double, triple and quadruple replacements with a somewhat larger basis set, 6-3 1 IG( p) [ 151. We may denote this level as: (R or U) MP4sdtq(FC)/6-31 lG(p)//(R or U) MP2(FULL)/6-31G(p). Level II. We use the same level, but with Dunning’s triple zeta basis set with polarization functions, p, (TZP) [ 161. We use the spin projection operator as implemented in GAUSSIAN 90 to annihilate, or, at least, to reduce, the spin contamination. The energy after spin decontamination will be PMP4, and we find expectation values of S* to be about 0.75 for the saddle point, indicating minor spin contamination. However, MP perturbation theory (especially the unrestricted version, UMP4) converges slowly, particularly in the case of free radicals [ 17,181. In the present work, we use

Physics 183 (1994) 95-100

two economical options that were put forward with the aim of taking the post-PMP4 energy into account: namely, the SAC4/Al empirical method [ 191, and the theoretical fourth-order invariant quantity of Feenberg [ 201 as an upper bound to the perturbation series for the correlation energy. Level III. Finally, for this level, we make a single point calculation to quadratic configuration interaction level, QCIsd( t) , with the basis set 6-3 11G ( p) We may denote this energy as: QCIsd(t)/6_311G(p)// MP2(FULL)/6-31G(p). The harmonic vibrational frequencies were calculated at the MP2/6-31G(p) level with a double objective: to confirm the nature of the saddle point and to provide the thermal corrections (TC) to the energy using the classical formalism of statistical thermodynamics. Such frequencies are known to be overestimated and so they were uniformly scaled by a factor of 0.95. For the reactants and products, we used experimental TC to avoid errors in the local shape. For the tunneling calculations, we use the standard symmetric one-dimensional Eckart potential fitted to the theoretically derived adiabatic barrier, using only the available information for reactants, products and saddle point energies (at the SAC4-A 1 level), and frequencies (at the MP2 level). This symmetric potential is given by [7] V,(s)

=4AE+Zl(

1 +Z)‘,

where AE * is the adiabatic barrier height and Z= exp( s/b), where s = - x,0, + x for reactants, saddle point and products, respectively, and b is the length scaling factor which is associated with the curvature of the potential at its maximum. There have certainly been lots of Eckart theory calculations done on the H + H2 reaction and its isotopic analogs [ 21-251. Thus, Pacey et al.‘s papers [ 21-231 involve use of the Eckart parameters as fitting functions to measured data; and Sato et al. [ 24,251 calculate the rate constants for the Hz + H reaction and its analogs at very low temperatures, using variational transitionstate theory including one-dimensional tunneling corrections.

3. Results and discussion Since we shall use the H, + H system exclusively as a test of our calculations, we shall restrict its study to

J. Espinosa-Garcia et al. /Chemical Physics 183 (1994) 95-100

the linear approximation. In this case, we found a length of 0.916 A at the MP2/6-31G(p) level, in excellent agreement with other values [ 26-371. The total energies with and without spin annihilation and with post-MP4 corrections are listed in Table 1 together with the energy changes relative to the reactants (zero level of energy). One of the best estimates of the barrier height is due to Liu and Siegbahn (LS) (9.8 kcal/mol) [30], and Liu [29] also gave upper and lower bounds for this barrier of 10.3 and 9.5 kcal/ mol, respectively. More recently, new ab initio calculations [ 31-331 have given a classical barrier of 9.65 kcal/mol, within Liu’s range. Levels II and III are outside this range, overestimating the activation barrier. At level I, the spin projection has a dramatic influence, lowering the barrier by 1.7 kcal/mol. The empirical post-PMP4 corrections have a comparatively minor influence, lowering the barrier by 0.5 kcal/mol. The differences between the two post-MP4 values, SAC4 and Feenberg, are negligible. We have investigated these effects on the barrier heights in the NH, + OH reaction [ 381, in the OH + H2 and HNO + H reactions [ 391 and in the NH3 + H reaction [ 403 with similar results. In the light of our experience with these reactions, we can conclude that, in order to obtain barrier heights that are comparable to experiment or to accurate theoretical results, one requires a large basis set, post-MP4 electron correlation, and spin projection in the case of free radicals. In the remainder of the work, therefore, we shall use the results from calculation level I. Table I Total energies ’ and-energy

UMP4

Hz

- 1.16769 - 1.16681

1

H,

AEb

Table 2 Unscaled harmonic vibrational frequencies (cm- ‘) and thermal corrections (LX) (300 K) (kcal/mol) for the saddle point This work d

0,

e

0,

r

w

TC

S’

WH’

PKb

H, DH, HD,

1037 993 793

951 911 727

947 907 724

879 841 671

H3 DHz HDz

2173 1867 1863

1945 1674 1673

2192 1881 1877

2130 1832 1829

H3 DH, HD,

2107i 2003i I.5871

1361i 129Oi 1021i

22953 21831 17291

19321 1836i 14561

HJ DH, HD,

8.20 7.38 6.89

“Ref. [58]. bRef. [26]. ‘Ref. [27]. d Bending vibration (doubly degenerate). eSymmetrical stretching vibration. f Imaginary frequency or reaction coordinate

vibration.

The harmonic vibrational frequencies and thermal corrections for the saddle point are listed in Table 2 together with other selected values from the bibliography. The linear saddle point was identified with one and only one negative eigenvalue of the Hessian matrix: that means one imaginary frequency. For the H3 system, the range of values of this imaginary frequency is from 1300 to 2300 cm - ’ (seeTable 2)) i.e., differences of around 1000 cm-‘, restricting the range somewhat to the more modern values, although a wider margin

changes ’ at various levels of theory

Level ’

II III

97

PMP4

SAC4

Feenberg

- 1.17390 - 1.17284

- 1.16826 - 1.16737

QCISD( T)

- 1.16828

I II III

- 1.64798 - 1.64581

I II 111

12.2 13.1

- 1.65074 - 164860

- 1.65777 - 1.65521

- 1.65231 - 1.65002 - 1.65030

’ In hartrees, H= -0.49981. bInkcal/mo1.0thervalues:9.8 [29],9.1 ’ I: MP4/6-311 +G(p)//MP2/6-3lG(p);

10.5 11.3

10.0 10.9

9.9 10.8 11.2

[26], 11 [28];9.5 [36]. II: MP4/TZP//MF2/6-31G(p);

III: QCISD(T)/6-311

+G(p)//MP2/6-31G(p)

J. Espinosa-Garcia ef al. /Chemical Physics 183 (1994) 95-100

98

'3.00 1

can be found in the literature [41-45]. Our value is in the upper zone of this range of values. In Figs. l-3, we show the function k(T) versus temperature over the range 200-l 000 K, for the Hz + H, D f H,, and H + Dz systems.

-

900

-

Y

3. I. H2 t H system

F -

The direct comparison of the theoretical and experimental rate constants must take into consideration the ortho-to-para and para-to-ortho conversion. In Fig. 1 by Kobs= #K, [ 42,461. we show Kobs approximated This approximation is correct for T> 300 K, while for T< 300 K, Schatz and Kuppermann’s accurate quantum calculations showed that Kobs is closer to K, + K_ , [47]. Firstly, if one compares the different theoretical values, one notes that the effects of the PES errors are appreciable. Thus, the rate constants calculated on Porter-Karplus surface no. 2 (PK2) [ 26) (lines 2) are greater than on the Liu-Siegbahnn surface [ 301 (with an accurate analytic representation by Truhlar and Horowitz [48]) (line 3) by a factor of between 5.3 (200 K) and 1.8 (600 K) [ 491. Secondly, this last rate constant is in excellent agreement with the experimental values which are available in the temperature range 1250

11 00

7

7.00 -

6 2 3

5.00 2;; 100

2.00

3.00

4.00

? 000/T(K)

Fig. 2. Arrhenius plot of log,& (cm’ mol-’ s-‘) against the reciprocal temperature (K) : The case of the H, + D reaction. Line 1: experimental value [27] ; line 2: experimental value [44]; line 3: theoretical value, CVT/MCPVAG, [ 521; line 4: theoretical value [53] ; line 5: our results at SAC4 level and Wigner tunneling; line 6: our results at SAC4 level and Eckari funneling.

300-444 K (line 1). Finally, our rate constant with Eckart’s tunneling factor (line 7) agrees with the experimental values and with the latest theoretical rate constant [ 50,5 11, with differences less than 30% with respect to the experimental values (range 300-400 K) and in the intermediate zone with respect to the theoretical values (range 300-I 000 K) . 3.2. D + H, system

10.50

Y _E?

8.50

2 7

6.50

I,, 1.00

I,, 2.00

,

,

/ 3.00

1 OOO/T(

, / 4.00

K)

Fig. 1. Arrhenius plot of log,& ( cm3 mol- ’ s- ‘) against the reciprocal temperature (K): The case of the H2 + H reaction. Line 1: experimental value [ 421; line 2: theoretical value on the PK2 surface [49 ] ; line 3: theoretical value on the LS surface [49] ; line 4: theoretical value on the LS surface [50]; line 5: theoretical value [ 5 I]; line 6: our results at SAC4 level and Wigner tunneling; line 7: our results at SAC4 level and Eckart tunneling.

One should note the excellent agreement in Fig. 2 between the experimental values (lines 1 and 2) and the theoretical values using the canonical variational theory with the MCPVAG transmission coefficient (CVT/MCPVAG) (line 3) [ 521 and using flux correlation function methods (line 4) [ 531. Bowman and Lee [ 541 using reduced dimensionality quantum methods, obtain a good agreement with the CVT/MCPVAG results [ 521. Our results, at the SAC4 level of calculation with the Eckart correction, present an excellent behaviour with respect to the forementioned values, with differences less than 30% in the range 300-1000 K. 3.3. H + D2 system The behaviour practically coincides with the other two systems studied (Fig. 3). Our calculation with

J. Espinosa-Garcia

et al. /Chemical

Physics 183 (I 994) 9.5-l 00

99

found that, as was to be expected, this factor underestimates the rate constants at low temperatures, but, for T> 300 K, gives a reasonable agreement with experimental values.

1 .oo

9.00

4. Conclusions 7.00

5.00

1.00

2.00

3.00

4.00

5.00

1000/T(K)

Fig. 3. Arrhenius plot of log,& ( cm3 mol- ’ s- ’ ) against the reciprocal temperature (K): The case of the H + D2 reaction. Line I : experimental value [27]; line 2: experimental value [45]; line 3: theoretical value, CVT/MCPVAG [ 521; line 4: our results at SAC4 level and Wigner tunneling; line 5: our results at SAC4 level and Eckart tunneling.

Eckart’s tunneling correction presents differences less than 20% in the range 300-1000 K. In the three reactions studied, the effects of the variational localization [ 551 and anharmonicities [ 491 are small, so that the most important factors will be the height of the activation barrier and the tunneling effect. For T> 300 K, the agreement between our theoretical values with the Eckart correction and the available experimental or theoretical values can be explained by the quality of the activation barrier obtained in the present work, within Liu’s proposed range of values, and because, for T> 300 K, the tunneling effect is similar whether obtained through the Eckart approximation or through other methods taking the PES curvature into account [ 491. For hydrogen abstraction reactions, it has been shown (see, for example, some recent papers [ 9,10,56 ] ) that the Eckart potential is narrower in the “tail” region than the true adiabatic potential, leading to transmission coefficients larger than those obtained with the true adiabatic potential at low temperatures. This behaviour would explain the poor agreement at low temperatures (T-C 300 K) in the three reactions studied here. Likewise, when we tested the behaviour of Wigner’s tunneling correction [57] in the three reactions, we

In the three reactions studied, we found good agreement with the experimental values and other more accurate theoretical values in the range 300-1000 K, with the agreement being poorer for T < 300 K. To attain these results with the conventional TST it is necessary to calculate the barrier height with spin projection, post-PMP4 energies and large basis set, and Eckart’s tunneling factor with a potential fitted to the available theoretical values. These results confirm the conclusions drawn in two previous studies with larger molecular systems [ 38-401. Hence, the use of this loworder approximation to obtain qualitative estimates of thermal rate constants may be justifiable in polyatomic molecular systems when a better knowledge of the PES and a correct description of the quantum-mechanical tunneling effect are not feasible.

Acknowledgement All calculations were carried out on the CONVEX 220 at CCUNEX (Badajoz, Spain). The authors would like to thank the DirecciQ General de Investigacidn Cientifica y Tecnica de1 Ministerio de Educaci6n Y Ciencia (Spain) for the partial support of this work (Project No. 930190).

References [I] B.C. Garret and D.G. Truhlar, J. Phys. Chem. 83 ( 1979) 1052; Accounts Chem. Res. 13 (1980) 440. [2] D.G. Truhlar, A.D. lsaacson and B.C. Garrett, in: Theory of chemical reaction dynamics, Vol. IV, ed. M. Baer (CRC Press, Boca Raton, FL, 1985) p. 65. [3] S. Glasstone, K.J. Laidler and H. Eyring, Theory of rate processes (McGraw-Hill, New York, 1941). [4] H.S. Johnston, Gas phase reaction rate theory (Ronald Press, New York, 1966). [5] K.J. Laidler, Theories of chemical kinetics (McGraw-Hill, New York, 1969). [6] C. Eckart, Phys. Rev. 35 (1930) 1303.

100

.I. Espinosa-Garcia et al. /Chemical Physics 183 (1994) 95-100

[7] R.P. Bell, The tunnel effect in chemistry (Chapman and Hall, London, 1980). [8] S. Scheiner and 2. Latajka, J. Phys. Chem. 91 (1987) 724. [9] A.K. Chandra, E.J. Padma Malar and D.S. Gupta, Intern. J. Quantum Chem. 41 (1992) 371. [ IO] A.K. Chandraand S. Rao, Intern. J. Quantum Chem. 47 ( 1993) 437. [ II] G. Leroy, M. Sana and A. Tinant, Can. J. Chem. 63 ( 1985) 1447. [ 121 M.J. Frisch, M. Head-Gordon, G.W. Trucks, J.B. Foresman, H.B. Schlegel, K. Raghavachari, M.A. Robb, J.S. Binkley, C. Gonzalez, D.J. DeFrees, D.J. Fox, R.A. Whiteside, R. Seeger, CF. Mehus, J. Baker, R.L. Martin, L.R. Kahn, J.J.P. Stewart, S. Topiol and J.A. Pople, GAUSSIAN 90 (Gaussian, Inc., Pittsburgh, 1990). [ 131 C. Moller and M.S. Plessett, Phys. Rev. 46 (1934) 618. [ 141 P.C. Hariharan and J.A. Pople, Theoret. Chim. Acta 28 (1973) 213. [15] R. Krishnan, J.S. Binkley and J.A. Pople, J. Chem. Phys. 72 ( 1980) 650. [ 161 T.H. Dunning Jr., J. Chem. Phys. 55 (1971) 716. [ 171 K.J. Raghavachari, Chem. Phys. 82 (1985) 4607. [ 181 R.H. Nobes, J.A. Pople, L. Radom, N.C. Handy and P.L. Knowles, Chem. Phys. Letters 138 (1987) 481. [ 191 M.S. Gordon and D.G. Truhlar. J. Am. Chem. Sot. 108 (1986) 5412. [20] E. Feenbcrg, Ann. Phys. (N.Y.) 3 (1958) 292; S. Wilson, Intern. J. Quantum Chem. 18 (1980) 905. [2l I H. Furue and P.D. Pacey. J. Phys. Chem. 91 (1987) 1584. [22] H. Fume and P.D. Pacey, J. Phys. Chem. 91 (1987) 4132. [23] IS. Jayaweeraand P.D. Pacey, J. Phys. Chem. 94 (1990) 3614. [24] T. Takayanagi, N. Masaki, K. Nakamura, M. Okamoto, S. Sato and G.C. Schatz, J. Chem. Phys. 86 (1987) 6133. [25] T. Takayanagi, K. Nakamura and S. Sato, J. Chem. Phys. 90 (1989) 1641. [26] R.N. Porter and M. Karplus, J. Chem. Phys. 40 (1964) 1105. [27] A.A. Westenberg and N. de Haas, J. Chem. Phys. 47 (1967) 1393. [28] 1. Shavitt, R.M. Stevens, F.L. Minn and M. Karplus, J. Chem. Phys. 48 ( 1968) 2700. [29] B. Liu, J. Chem. Phys. 58 (1973) 1925. [30] B. Liu and P. Siegbahn, J. Chem. Phys. 68 (1978) 2457.

[31] B. Liu, J. Chem. Phys. 80 (1984) 581. [32] D.M. Ceperley and B.J. Alder, J. Chem. Phys. 81 (1984) 5833. [33] R.N. Barnett, P.J. Reynolds and W.A. LesterJr., J. Chem. Phys. 82 (1985) 2700. [34] J.M. Bowman, Advan. Chem. Phys. 61 ( 1985) 115. 1351 A.J.C. Varandas, F.B. Brown, C.A. Mead, D.G. Truhlar and N.C. Blais, J. Chem. Phys. 86 (1987) 6258. [36] C.W. Bauschlicher Jr., S.R. Langhoff and H. Partridge, Chem. Phys. Letters 170 ( 1990) 345. [37] A.1. Boothroyd, W.J. Keogh, P.G. Martin and M.R. Peterson, J. Chem. Phys. 95 ( 1991) 4343. [38] J. Espinosa-Garcia, J.C. Corchado and M. Sana, J. Chim. Phys. Phys. Chim. Biol. 90 (1993) 1181. [39] J. Espinosa-Garcia, E.A. Ojalvo and J.C. Corchado, J. Mol. Struct. THEOCHEM 109 (1994) 131. [40] J. Espinosa-Garcia, S. Tolosa and J.C. Corchado, J. Phys. Chem. 98 (1994) 2337. [41] J.K. Cashionand D.R. Herschbach, J. Chem. Phys. 40 (1964) 2358. [42] W.R. Schultzand D.J. LeRoy, J. Chem. Phys. 42 ( 1965) 3869. [43] B.C. Garrett. D.G. Truhlar, A.J.C. Vamndas and NC. Btais, Intern. J. Chem. Kinet. 18 (1986) 1065. [44] J.V. Michael and J.R. Fisher, J. Phys. Chem. 94 (1990) 3318. 1451 J.V. Michael, J. Chem. Phys. 92 (1990) 3394. [46] D.G. Truhlar, J. Chem. Phys. 65 (1976) 1008. [47] G.C. Schatz and A. Kuppermann, J. Chem. Phys. 65 (1976) 4668. [48] D.G. Truhlar and C.J. Horowitz, J. Chem. Phys. 68 (1978) 2466; J. Chem. Phys. 71 (1979) (E) 1514. [49] B.C. Garrett and D.G. Truhlar. Proc. Natl. Acad. Sci. USA 76 (1979) 4755. [50] M.C. Colton and G.C. Schatz, Intern. J. Chem. Kinet. 18 (1986) 961. [51] T.J. Park and J.C. Light, J. Chem. Phys. 91 (1989) 974. [52] B.C. Garrett and D.G. Truhlar, J. Chem. Phys. 72 (1980) 3460. [53] T.J. Park and J.C. Light, J. Chem. Phys. 94 ( 1991) 2946. [54] J.M. Bowman and K.T. Lee, J. Chem. Phys. 79 (1983) 3742. [SS] B.C. Garrett,D.G. Truh1arandG.C. Schatz. J. Am. Chem. Sot. 108 (1986) 2876. [56] A.D. Isaacson, L. Wang and S. Scheiner, J. Phys. Chem. 97 (1993) 1765. [57] E.P. Wigner, Z. Physik. Chem. B 19 ( 1932) 203. [58] 1. Shavitt, J. Chem. Phys. 31 (1959) 1359.