875
J. Phys. Chem. 1994,98, 875-886
Interpolated Variational Transition-State Theory and Semiclassical Tunneling Calculations of the Rate Constant of the Reaction OH + at 200-3000 K Vasilios S. Melissas and Donald C. Trublar' Department of Chemistry and Supercomputer Institute. University of Minnesota, Minneapolis; Minnesota, 55455-0431 Received: September 7, 1993; In Final Form: October 29, 1993"
-
+
+
A b initio molecular orbital calculations were carried out for the reaction O H C2H6 H2O C2H5 using second-order Maller-Plesset perturbation theory and employing a very large basis set. Correlation energy is found to play an important role in determining the barrier height of the saddle point and also the geometry and the vibrational frequencies of both the transition state and the equilibrium points. The final calculated values for the forward and reverse classical barrier heights are 4.0 and 21.6 kcal/mol, respectively. These have been used to treat the kinetics of the reaction in the temperature range from 200 to 3000 K by using interpolated canonical variational transition-state theory and the centrifugal-dominant, small-curvature tunneling approximation, including information at the reactants, products, transition state, and one extra point along the minimum energy path. The calculated rate constants agree well with experiment over a wide temperature range.
1. Introduction Of the numerous radicals of relevance to both atmospheric and combustion chemistry, the hydroxyl radical plays a particularly important role. For example, reaction with the hydroxyl radical is the primary loss channel for alkanes in the stratosphere, a loss mechanism which receives some competition from photolytic decomposition of alkanes by solar UV light.' The hydroxyl radical is also a principalchain carrier in hydrocarbon combustion. Thus, the reaction of OH with ethane, the simplest hydrocarbon after methane, is of particular importance in both combustion and atmospheric applications, and its theoretical investigation is very interesting. Many experimental studies of the kinetics of this reaction have been reported, using both absolute and, especially at elevated temperatures, relative rate The results have been evaluated and combined into four recommended expres~ i o n s * ~ J ~ J ~ , first ~ ~ -three t h e being three-parameter fits and the last a two-parameter one. The most recent recommended expression2' covers the temperature range from 210 to 380 K, and the most reliable three-parameter oneIgcovers the temperature range from 226 to 800 K. For the overlapping temperature range of the two fits, the low-temperature fit2' is considered more reliable. An Arrhenius plot of the three-parameter recommended expressionl9 shows significant curvature at temperatures below about 700 K, and this may be an indication that tunneling plays an important role in this temperature range; however, no curvature is included in the Arrhenius fit of the low-temperature data.21 The lack of an analytical potential energy function or of any ab initio calculationshas impeded theoretical progress in dynamical studies of reaction (Rl). In the present study, we report accurate ab initio calculations for the structures, energetics, and vibrational frequencies of the reactants, products, and saddle point of reaction ( R l ) and we also report results at the same level for an extra point along the minimum energy path (MEP23-Z5).The extra point is displaced from the saddle point a small distance toward the products. Then, we use these calculationswith a recently proposed approximation schemG7forquantizedvariationaltransition-state theory (VTST) and semiclassical tunneling calculations to predict reaction rate constants. This scheme is called interpolated VTST or IVTST, *Abstract published in Advance ACS Abstracts, December IS, 1993.
0022-3654/94/2098-0875$04.50/0
-
-
and it was recently applied to the reactionsof OH + CHI H2O CHsz8and OH CD4 HDO CD329with ab initio surface information. The resulting rate constants and non-Arrhenius behavior for those reactions are in good agreement with experiments. In the present work, we use the first-order global version of interpolated canonicalvariational transition-state theory with small-curvaturetunneling (SCT) corrections to predict reaction rate constants and activation energies. The SCT calculations are based on the relatively new centrifugal-dominant, smallcurvature semiclassicaladiabatic g r o u n d - ~ t a t e ~(CD-SCSAG, ~,~ or, for short, small-curvature tunneling, SCT) method, and we also compare them to calculations with a zero-curvature tunneling24,26,31*32 (ZCT) calculation.
+
+
+
'
2. Computational Methodology
For the electronic structure calculations, we employ restricted Hartree-Fock (RHF) theory33for closed-shell systems (H20, C2H6) and unrestricted Hartree-Fock (UHF) theory33for openshell systems (OH, C2H5, the saddle point, and the extra point along the MEP). Second-order Mraller-Plesset perturbation theory33(MP2), unrestricted MP233(UMP2), and UMPZ with spin p r ~ j e c t i o n(PMP2IIUMP2) ~~-~~ are employed, based on these RHF and UHF reference wave functions. The motivation of the spin-projection method arises from previous studies where the transition structures calculated at the UHF level are affected by a considerable amount of spin contaminati~n,~~ and as a result, the reaction barriers can be overestimated by up to 10 kcal/mol when correlationcorrections are calculated by unrestricted Mraller-Plesset perturbationtheory. In this rest of this paper, we omit the U when referring to UHF and UMPZ calculations. Thus, H F denotes H F or UHF depending on whether the number of electrons is even or odd, respectively, and likewise for MP2. All correlation energy in MP2 calculationsis ~ c a l e d (scaling 3~ all correlation energy is abbreviated SAC) to yield an improved calculation of the barrier height at the MP2 optimized structure. This is called the MP-SACZ//MPZ level. The total electronic energy in the MP-SACZ//MPZ calculations is written as36
where 3 2 is a scaling factor which is assumed to be independent of geometry for a given system and HF denotes RHF or UHF, 0 1994 American Chemical Society
876
The Journal of Physical Chemistry, Vol. 98, No. 3, 1994
depending on the system. The motivation of the SAC method for large basis sets, such as used for the calculations here, is to make the calculated exoergicity of the reaction agree perfectly with experiment so that the treatment of reactants and products is balanced and a set of forward and reverse rate constants consistent with the equilibrium constant can be obtained. The 3 2 factor can be computed from the accurate, MP2, and H F binding energies for each bond type. Further details of the SAC method may be found in the original work36and our previous application.28 A correlation-consistentpolarized-valence triple-{(cc-pVTZ) basis set of Dunning3' is the starting point for our large-scale calculations. For H, this is a (5s2pld) primitive Gaussian basis set contracted to [3s2pld] (18 primitive functions/l4 contracted functions), whereas for C and 0,it is a (10sSp2dlf) primitive set contracted to [4s3p2dlfl (60 primitive functions/30 contracted functions). Pure d and f functions are used in the d and f shells, i.e., the d functions have five components and the f functions seven. We used this basis unchanged for some calculations, although for others, we used an adjusted version of this basis set where the same 3 2 scaling factor applies to the computed dissociation energies of both the C-H and 0-H bonds. This adjusted basis, called second-adjustedcc-pVTZ (adj2-cc-pVTZ), is explained in detail in section 4. The basis set was adjusted in order to treat electron correlation in the forming and breaking bonds of reaction (Rl) on an equal footing as required to obtain an accurate barrier height for both directions of reaction. The geometries of all stationary points, Le., minimum energy structures and the saddle point, were fully optimized at the H F and MP2 levels of theory with the adj2-cc-pVTZ basis set. In addition, the geometries of all equilibrium points, i.e., reactants and products,were fully optimized at the HF/cc-pVTZ and MP2/ cc-pVTZ levels of theory. We used C1 point group symmetry for each species as the initial geometry of the optimization procedure, and we fully optimized3sall the Z - m a t r i ~ parameters. ~~q~ Internal coordinates4Iwere selected for the geometry optimization, since the total energy is invariant to overall translation and rotation. The resulting point group symmetries at all available levels are C,, for OH, CZ,for HzO, C, for CZHS, and D3d for CZH6, whereas the symmetry is C, for the CzHs-H-OH saddle point at the HF/ adj2-cc-pVTZ level and C1 at the MPZ/adjZcc-pVTZ level. The symmetry is CIfor the extra point along the reaction path on the product side of the saddle point (no optimization is required for this point). Vibrational frequencies were calculated at four levels for all equilibrium points, in particular at the H F and MP2 levels for both the cc-pVTZ and adj2-cc-pVTZ basis sets. However, only HF/adjZ-cc-pVTZ and MPZ/adjZ-cc-pVTZ frequency calculations were carried out for the saddle point, and only MP2/ adj2-cc-pVTZfrequencieswere calculated at the extra point along the MEP. Since unrestricted H F wave functions are not spin eigenfunctions, we have monitored the expectation value of Sz in all unrestricted H F and MP2 wave functions (where S denotes electron spin angular momentum in units of f i ) . For doublets, S2wasalwaysin therange0.756-0.796, where 0.75Oisthecorrect value. (For completeness, we note that the largest deviation of S2,0.046, from the correct value corresponds to the saddle point calculation at the HF/adjZ-cc-pVTZ level.) Although this amount of spin contamination is not considered severe, we mention that the forward barrier height of reaction (Rl) at the MP2/ adj2-cc-pVTZlevel-for which spin projection is not included-is larger by 1.80 kcal/mol than its value at the PMP2//MP2/ adj2-cc-pVTZ level. Various amounts of memory and disk space were required for the optimization and frequenciescalculationsof thevarious species under study. Most of the optimization and frequencies calculations required no more than 12 megawords of memory and 2
Melissas and Truhlar TABLE 1: Optimized Minimum Energy Structures (bond lengths in 4 bond angles in degrees) for OH and HzO level/basis set HF/a-pVTZ HF/adjZ-cc-pVTZ MPZ/cc-pVTZ MPZ/adjZ-cc-pVTZ experimental
OH (Cw) 0-H 0.9512 0.9512 0.9659 0.9667 0.9706'
H20 (Cd 0-H 0.9406 0.9407 0.9575 0.9586 0.9576b
HOH 106.00 105.96 103.66 103.57 1O4.4gb
Reference 48. Reference 49.
TABLE 2 Optimized Minimum Energy Structures (bond lengths in A, bond and dihedral angles in degrees) for c2H5 (c8)
level/ basis set RC(IMC(II .. .. RC(IFH(I)
Rc(ltH(3 RC(2kH(5)
a
B Y
61' 626
HF/cc-pVTZ and HF/adjZ-cc-pVTZ . ~
MPZ/cc-pVTZ and Mp2/adiZcc-~VTZ , .
1.4948 1.0893 1.0838 1.0734 111.57 111.34 120.47 119.66 82.57
1.4809 1.0918 1.0851 1.0728 111.60 111.55 120.67 119.42 84.54
61 is the dihedral angle between the H ( 2 ) 4 ( 1 ) 4 ( 2 ) and C(1)C(Z)-H(l) planes. 62 is the dihedral angle between the H(4)4(2)C( 1) and C(Z)-C(l)-H( 1) planes.
gigabytes of disk space; however, the MP2 frequencycalculations in the saddle point region required 18 megawords of memory and 8 gigabytes of disk space. D i r e ~ tSCF ~ ~ methods , ~ ~ were used in all our calculations, recomputing the two-electron integrals as needed. Moreover, analytic Hessian calculationswere performed at both the HF44and MP245levels of theory. A mass-scaled coordinate s y ~ t e mwas ~ ~defined , ~ ~ as
where mi is the mass of atom i, 1.1 is the reduced mass of the relative translational motion of reactants, which in this case equals 10.858 amu, and Ri, is a Cartesian coordinate with y = x, y , or z. The geometry of a nonstationary point along the MEP on the product side of the saddle point was determined by taking a step of length 0.005 a0 in the direction of the mass-scaled eigenvector corresponding to the saddle point imaginary frequency. All ab initio calculations were carried out using the GAUSSIAN9247program on the Cray X-MP/EA and Cray Y-MP8/ 864 computers at the Minnesota and San Diego Supercomputer Centers, respectively. 3. Structures and Vibrational Frequencies
The optimized geometries for the minimum energy structures are given in Tables 1-3. We found that the H F bond lengths of both closed- and open-shell species are generally shorter than the experimental value^^^-^^ (which are available only for OH, HzO, and C Z H ~ )but , by at most 2%. However, when electron correlation is included at the MP2 level, bond lengths agree with the experimental data within 0.5%. Vibrational analysis was performed at all points, and the vibrational frequencies for the equilibrium species calculated at
VTST and Semiclassical Tunneling Calculations
The Journal of Physical Chemistry, Vo1. 98. No. 3, 1994 877
TABLE 3: Optimized Minimum Energy Structures (bond length in & bond and dihedral angles in degrees) for CZH6 ( a d )
TABLE 5
Optimized Staggered (C,)Transition-State
Structures at the HF/adjZ-cc-pVTZ and MPZ/adjZ-cc-pVTZ Levels, respectively (bond lengths in A, bond and dihedral angles in degrees) O(1) RH(4).0(l) /
level/ basis set
HF/cc-pVTZ and HF/adj2-cc-pVTZ
Rc-c Rc-H
MPZ/cc-pVTZ and MPZ/adjZ-cc-pVTZ
1S245 1.0841
&Ib
experimental‘
1.5177 1.0850
111.21 120.00
Q
111.17 120.00
1S243 1.0895
111.95 120.00
a Reference 50. &I is the dihedral angle between the H(2)-C(1)C(2) and C(l)-C(2)-H(l) planes.
TABLE 4
Vibrational Frequencies (cm-’)of the Equilibrium
structures HF/ MP2/ HF/ sym- cc- adj2-cc- ccmode metry pVTZ pVTZ pVTZ 1 1 2 3
ug
al al b2
4039 4127 1753 4227
4127 1755 4226
H20 3872 3847 1650 1646 3993 3963 (13.52)” 3239 3127 3049 1520 1505 1419 1098 998 488 3330 3169 1516 1218 825 135(38.08)“
n/aa n/a nja n/a n/a n/a nja n/a n/a n/a nja n/a n/a n/a n/a
3033e 292Oe 2842“ 1440” n/a 1366e 1138e n/a 54Oe 3112c 2987e 144W 117Y nja n/a
3101 1441 1038 343 3099 1418 3158 1535 1239 3180 1536 834(47.75)O
3043f 1449’ 1016‘ 303’ 30611 1438f 3175’ 1559 1246‘ 314Cf 1526‘ 829
2954 13888 9958 2898 289@ 13798 296gh 14688 11908 29858 14698 8228
a‘ a’ a’ a’ a‘ a’ a’ a’ a’ a” a” a” a” a” a”
3264 3180 3113 1603 1582 1523 1100 1068 465 3362 3214 1600 1294 861 162
3264 3180 3113 1603 1582 1523 1100 1068 465 3362 3214 1600 1294 861 162
1 2 3 4
al, al, al, al, a2u azu e, e, e, e,, e,, e,,
3158 1548 1049 328 3152 1522 3194 1614 1323 3222 1619 881
3158 1548 1049 328 3152 1522 3194 1614 1323 3222 1619 881
C2H6 3101 1441 1038 343 3099 1418 3158 1535 1239 3180 1536 i834
H(6)
level/basis set, transition-stateorder RC(I)-W Rc( I )-H( I ) Rc(I)-H(I) RC(ZFH(4) Rc(2FH(9 RH(4)-0(l)
Ro(IFH(~) a
B
4040
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
5
har- fundamonic mental (exp) (exp)
OH 3837 3814 (5.45)’
C2H5 3239 3127 3049 1520 1505 1419 1098 998 488 3330 3169 1516 1218 825 135
6 7 8 9 10 11 12
MP2/ adj2-ccpVTZ
~l)-H(7l
y \ W 7 )
373Sb 3570b
Y
6 c
r
f &d
Values in parentheses are total harmoniczero-point energies in kcal/ mol. Reference 48. Reference 51. n/a: not available. Reference 52. /Reference 53. g Reference 54. Reference 55. the HF and M P 2 levels for both the cc-pVTZ and adj2-cc-pVTZ basis sets are listed in Table 4. Since the adjustment to the ccpVTZ basis affected only the oxygen basis functions, calculated geometries and vibrational frequencies for C2H5 and C & species with the cc-pVTZ basis set are the same as the ones calculated with the adj2-cc-pVTZ basis set. Moreover, Table 4 lists some e~perimenta1~8-5l-55 values. It is worthwhile to note that the calculated frequencies are harmonic values, and the differences
;wei;
HF/adjZ-cc-pVTZ, first‘ 1.5110 1.0862 1.0825 1.2785 1.0801 1.2349 0.9472 110.51 111.11 106.25 114.59 178.91 99.03 119.77 180.00 114.65 180.00 0.00
MPZ/adjZ-cc-pVTZ, secondb 1.508 1.0866 1.0841 1.1670 1.0834 1.3601 0.9675 110.65 111.01 107.01 113.45 175.42 95.90 119.87 180.00 116.51 180.00 0.00
‘Saddle point. Hill top. T is the dihedral angle between the H(2)C(l)-C(2) andC(l)-C(2)-H(1) planes. &isthedihedralanglebetwecn the H(4)-C(2)-C( 1) and C(2)-C( 1)-H(l) planes. e x is the dihedral angle between the H(S)-C(2)-C(l) and C(2)-C(l)-H(4) planes. f # is the dihedral angle between the 0(1)-H(4)-C(2) and H ( 4 ) 4 ( 2 ) - C ( l ) planes. w is the dihedral angle between the H(7)-0( 1)-H(4) and O(1)H(4)-C(2) planes. between experimental harmonic and fundamental frequencies may be as large as the differences between calculated and fundamental frequencies in Table 4. Harmonic frequencieshave been estimated from experiment only for 16 out of 31 cases; the available harmonic frequencies are given in Table 4 as well, and the calculated frequencies are in much better agreement with these harmonicvalues than with the fundamentals. Discrepancies between calculated and experimental harmonic frequencies are largest for the bond-stretching modes. Notice that, in all stable species for which harmonic frequencies are available in Table 4, using the adj2-cc-pVTZ basis a t the MP2 level usually improves the frequencies as compared to using the cc-pVTZ basis, with the largest error in the former case relative to the available experimental harmonic values equal to 79 cm-l, whereas in the latter case, the largest error is 102 cm-*. In general, the good agreement of the MPZ/adjZ-cc-pVTZ frequencies in Table 4 with the available experimental harmonic values provides confidence that this is a reasonable level for studying the transitionstate region. The optimized transition-state structures calculated for the OH C2H6 reaction a t both the HF and M P 2 levels using the adj2-cc-pVTZ basis set are given in Tables 5 and 6, respectively. As mentioned above, the transition state was located under the assumption of Cl symmetry, resulting in a C, structure a t the HF level and in a C l one a t the M P 2 level. The calculated Hessian matrix for these structures was found to possess only one imaginary eigenvalue, confirming that the transition-state geometry is
+
Melissas and Truhlar
878 The Journal of Physical Chemistry, Vol. 98, NO. 3, 1994
TABLE 6
Optimized Eclipsed (C,) Transition-State Structure and Geometry of the Nonstatio~ryPoint along the OH Cfi Reaction Path on the Product Side of the Saddle Point at the MPZ/adjZ-cc-pVTZ Level (bond lengths in & bond and dihedral angles in degrees)
+
TABLE 7: Transition-State Harmonic Vibrational Freqwnch (cm-l) and Harmonic Zero-Point Energies (kcal/mol) HF/adja’ 4089,3212,3203,3150, 1611,1586, 1539, 1524, CC-PVTZ 1321,1117, 1014,786,523,149,3058i a”
3277,3229,1604, 1349, 1272,884,466,185,32
ZPE 53.07 MP2/adj-
CC-PVTZ
S ( d
saddle point 0.000
nonstationary point 0.005
~
RC(I)-C(2) Rc(I)-H(I) RC(I)-H(Z)
Rc(I)-H(J)
Rc(z)-H(~) Rc(~)-H(~) Rc(z)-H(~) R~94)4(1) RO(ltH(7)
a
B Y
6 c
r
7
9 TIa
726
4= XId
9‘
og
1.5075 1.0866 1.0841 1.0841 1.1678 1.0835 1.0836 1.3595 0.9673 110.71 110.99 111.01 107.65 113.34 1 1 3.79 171.65 96.10 119.94 119.90 179.32 117.72 115.46 98.76 32.21
1.5075 1.087 1 1.0839 1.0843 1.1681 1.0843 1.0823 1.3581 0.9683 110.72 1 1 1.00 111.00 107.60 113.30 113.82 172.33 95.87 119.96 119.89 179.69 117.32 115.91 100.14 30.64
T I is the dihedral angle between the H(2)-C(l)-C(2) and C(1)C(2)-H(1) planes. 7 2 is the dihedral angle between the H(l)-C(2)C(1) and C(2)-C(l)-H(3) planes. 4 is the dihedral angle between the H(4)-C(2)-C( 1) and C(2)-C( 1)-H( 1) planes. XI is the dihedral angle between the H(5)-C(2)-C(1) and C(2)-C(l)-H(4) planes. exz is the dihedral anglebetweentheH(4)-C( 1)-C(2) andC(1)-C(2)-H(6) planes. f + is the dihedral angle between the 0(1)-H(4)-C(2) and H(4)-C(2)C(1) planes. * o is the dihedral angle between the H(7)-0( 1)-H(4) and 0(1)-H(4)
