Article pubs.acs.org/JPCA
Quantum Chemistry, Reaction Kinetics, and Tunneling Effects in the Reaction of Methoxy Radicals with O2 Hongyi Hu and Theodore S. Dibble* Chemistry Department, State University of New York-Environmental Science and Forestry, Syracuse, New York 13210, United States S Supporting Information *
ABSTRACT: The reaction of the methoxy radical with O2 is the prototype for the reaction of a range of larger alkoxy radicals with O2 in the lower atmosphere. This reaction presents major challenges to quantum chemistry, with CCSD(T) overpredicting the barrier height by about 7 kcal/mol in the complete basis set limit. CCSD(T) calculations also indicate that the CH3OOO• analog of the HOOO• radical is energetically unstable with respect to CH3O• + O2, a finding that seems unlikely. The previous successful prediction of the barrier height using CCSD(T)/cc-pVTZ energies at CASSCF/6-311G(d,p) geometries is shown to rely on the use of a metastable Hartree− Fock reference wave function. The performance of several density functionals is explored and B3LYP is selected to examine the role of tunneling, including the competition between small curvature tunneling (SCT) and large curvature tunneling (LCT). SCT is found to be sufficient to describe tunneling, in contrast to the typical findings for bimolecular hydrogen-abstraction reactions. The previously proposed mechanism of a cyclic transition state yields rate constants for CH3O• + O2 that faithfully reproduces the experimentally derived Arrhenius pre-exponential term. Predictions of the branching ratios for the competing reactions CH2DO• + O2 → CHDO + HO2 (1a) and CH2DO• + O2→ CH2O + DO2 (1b) are also in good agreement with experiment. reaction.9 The rate constant for the alkoxy + O2 reaction can then be used to calculate a rate constant of the dominant unimolecular reaction channel. Unfortunately, rate constants are not known for alkoxy radicals derived from nonalkane hydrocarbons, let alone, e.g., oxygenated volatile organic compounds (VOCs).8 The fate of alkoxy radicals strongly influences the production of ozone and secondary organic aerosol in polluted air,10−13 and atmospheric chemists concerned with the mechanism of VOC degradation would appreciate a reliable theoretical approach for determining these rate constants.14−18 This study was initiated with the ultimate goal of validating theoretical methods that could be used to predict rate constants for a wide range of alkoxy radicals reacting with O2. The methoxy radical itself is an important intermediate of atmospheric oxidation of methane, a process that eventually results in formation of molecular hydrogen. Measured HD/H2 ratios are used to quantify the sources and sinks of molecular hydrogen on a global basis.19 The change in deuterium content during the conversion of CH3D to molecular hydrogen has been found to be greatly affected by the branching ratio, k1a/k1b, of singly deuterated methoxy radical (CH2DO•) reacting with O2:20
I. INTRODUCTION The rate constant for methoxy radical reacting with O2 has been measured multiple times1−7 because this reaction is the prototype for a wide range of alkoxy radical + O2 reactions of atmospheric interest.8,9 For larger alkoxy radicals, such as the 2pentoxy radical shown in Scheme 1, the reaction with O2 in the atmosphere competes with multiple unimolecular reactions.8,9 In many kinetics experiments on alkoxy radicals, the only quantity that can be determined is a ratio of rate constants, kO2/ kuni, between the O2 reaction and the dominant unimolecular Scheme 1. Reaction Pathways of 2-Pentoxy Radicals, Including Reaction with O2, β C−C Scission (Decomposition), and Isomerization (1,n H-Shift)
Received: September 11, 2013 Revised: November 23, 2013
© XXXX American Chemical Society
A
dx.doi.org/10.1021/jp409105q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
CH 2DO• + O2 → CHDO + HO2
(R1a)
CH 2DO• + O2 → CH 2O + DO2
(R1b)
Scheme 2. Structure of Saddle Points (Left, TS1; Right, TS1′) of Methoxy Radical + O2 Found by Bofill et al.27 a
Therefore, there is a strong need to explicitly determine the branching ratio of CH2DO• + O2 reaction (k1a/k1b) to provide data for atmospheric modelers. This branching ratio has been measured at room temperature by Nilsson et al.21 and in the temperature range 250−333 K by a collaboration that included the present authors.22 However, for application to the upper troposphere and lower stratosphere the branching ratio needs to be determined for temperatures down to about 200 K. Therefore, a secondary goal of this research was to use theory to reliably extend the branching ratio, k1a/k1b, to temperatures below the range in which measurements have been made. The rate constant for CH3O• reacting with O2 has been determined experimentally at room temperature and above,1−7 leading to a recommended rate constant of (1.6−1.9) × 10−15 cm3 molecule−1 s−1 at 298 K.8,23,24 Three research groups1−3 determined the absolute rate constant at temperatures spanning the range 298−973 K. Over the range 290 K ≤ T ≤ 610 K, these results are well fit by the expression k(CH3O+O2) = −14 7.8+4.7 exp[−(1150 ± 190)/T] cm3 molecule−1 s−1, −2.9 × 10 which corresponds to 1.6 × 10−15 cm3 molecule−1 s−1 at 298 K.8 The relative rate constant technique has also been applied to study the rate constants of this reaction, giving similar results at 298 K.4−7 The pre-exponential factor in the Arrhenius equation is rather low compared with typical hydrogenabstraction reactions, which led Jungkamp and Seinfeld to suggest that this reaction may not proceed via direct H-atom abstraction.25 Another oddity is that the rate constants obtained at T > 610 K greatly exceed that obtained by an extrapolation of an Arrhenius plot of the data at lower temperatures.3 Quantum chemical calculations have been carried out to explore the reaction mechanism as well as the importance of tunneling. Jungkamp and Seinfeld25 suggested that the reaction occurs through O2 addition to form methyl trioxy radical, CH3OOO•,26 followed by HO2 elimination via a five-member ring saddle point: CH3O• + O2 → CH3OOO•
(R2)
CH3OOO• → CH 2O + HO2•
(R3)
a The dashed line represents the breaking C−H bond and forming H− O bond, and the dotted line indicates the noncovalent O−O interaction in TS1.
For computing rate constants, Bofill et al. took tunneling into account using the asymmetric Eckart potential, resulting in a tunneling coefficient, κ(T), of ∼9 at 298 K. Setokuchi and Sato30 adopted Bofill et al.’s conclusion that hydrogen abstraction is the main channel of the alkoxy radical with O2 reaction and calculated the potential energy surface at the G2M(RCC1)31 level of theory. Setokuchi and Sato considered tunneling with a multidimensional approach: small-curvature tunneling (SCT). Their calculations suggested that tunneling is more modest at room temperature (κ ≈ 2) than found by Bofill et al. They did obtain a much larger tunneling coefficient at lower temperatures of atmospheric interest (∼8 at 200 K). Both papers obtained rate constants in reasonable agreement with experiment; however, their results did not explain the significant non-Arrhenius behavior observed by Wantuck et al., at T > 600 K.3 The disagreement between the two theoretical studies regarding the importance of tunneling at 298 K casts doubt on our ability to reliably predict the branching ratio of the two product channels of the CH2DO• + O2 reaction at T < 250 K. Both the existing data on this branching ratio (for 250 K < T < 333 K) and ongoing experiments on the rate constant for CH3O•/CD3O• + O2 reaction will provide stringent tests of the results of theoretical calculations of barrier heights and tunneling. Bofill et al. suggested, on the basis of large spin contamination in the UHF wave function, that the wave function of TS1 was strongly influenced by nondynamical correlation. However, the fact that both their UCCSD(T) single point calculations and the G2M(RCC1) calculations of Setokuchi and Sato obtained reasonable barrier heights might suggest that the multireference character is not so bad as to preclude the (cautious) use of single-reference methods. In this paper, we have employed single-reference quantum chemical methods to study tunneling and variational effects in the methoxy + O2 reaction on the doublet spin surface. Tunneling is studied using the multidimensional (semiclassical) tunneling treatments. Because of the need to use efficient methods of electronic structure calculation to compute multidimensional tunneling, we carefully evaluate the performance of various single-reference methods against each other and the experimental activation energy before proceeding with tunneling calculations.
Bofill et al.27 identified an error that undermined this conclusion. Moreover, they recalculated the potential energy surface at UCCSD(T)/cc-pVTZ//CASSCF(9,7)/6-311G(d,p) and found the cis-methyltrioxy radical to be 4.8 kcal/mol above reactants. Because they calculated a 50 kcal/mol barrier on the path leading from the trioxy radical to product CH2O and HO2, they concluded that the trioxy radical was not relevant to this reaction. Note that the simplest trioxy radical, HOOO•, is very challenging for quantum chemistry28,29 and that CH3OOO• has never been characterized by any experiment. Bofill et al. found a saddle point (TS1) for direct abstraction of a hydrogen atom from CH3O by O2. An oddity of the structure of the saddle point is a noncovalent interaction between the radical center of the methoxy radical and the oxygen atom of O2 that is not abstracting the hydrogen atom (Scheme 2) and a weakly bound (0.3 kcal/mol) prereactive complex. Bofill et al. found a second (acyclic) saddle point (TS1′) that was computed to be ∼8 kcal/mol higher in energy than TS1.
II. METHODS Electronic structures of all species were based on spinunrestricted wave functions, except formaldehyde and where noted. The geometries of relevant stationary points were optimized with several density functional theory (DFT) methods, quadratic configuration interaction (QCI),32 and B
dx.doi.org/10.1021/jp409105q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
coupled cluster (CC) methods.33 QCI and CC calculations were carried out with the frozen core approximation. The DFT methods include generalized gradient approximation (GGA) and meta-GGA methods (PBE1W,34 MPWLYP1W,34 M06-L,35 BLYP, 36,37 OLYP, 38 and B97D 39,40 ), hybrid methods (B3LYP, 36,37 X3LYP, 41 O3LYP, 38,42 M06-2X, 43 CAMB3LYP44), double hybrid DFT method (B2PLYP45), and long-range corrected46 variants of some of the above methods (LC-BLYP, LC-OLYP, LC-M06L). Quadratic configuration interaction with single and double excitations (QCISD) and with a perturbative estimate of the triple excitation (QCISD(T)) were used; the analogous CCSD and CCSD(T) methods were also used. Stability checks were performed to make sure the UHF and UDFT wave functions were stable (that is, the lowest energy roots of the UHF or UDFT equations). In several cases, the quadratically convergent Hartree−Fock method47 was applied instead of Pulay’s direct inversion in the iterative subspace (DIIS) extrapolation48 to ensure the UHF calculations converged on the stable wave function. Frequency analyses were carried out for all the stationary points to ensure that radicals and complexes had all real vibrational frequencies and that saddle points had a single imaginary frequency. Intrinsic reaction coordinate (IRC) calculations49 were done for the saddle points for DFT methods to confirm the identities of reactants and products. Most geometries are computed using Pople’s 6-311G(d,p) basis set, as was used in previous theoretical studies. Tests of basis set effects on geometry were carried out at B3LYP using a range of basis sets up to aug-cc-pVTZ in size. Tests of basis set effects on the activation barrier were carried out at UCCSD(T) using single point calculations with Dunning’s correlation consistent basis sets cc-pVxZ (x = D, T, Q). The UHF energy was extrapolated to the complete basis set (CBS) limit using a three-parameter exponential extrapolation: EXHF = ECBSHF + A exp(−αX), and the CCSD(T) correlation energy was extrapolated using a two-parameter extrapolation, EXCorr = ECBSCorr + BX−3.50 The sum of the two terms determines the CBS-extrapolated UCCSD(T) energy. As a variation of couple cluster methods, single point energies were also computed using Brueckner Doubles calculation with a triples contribution (denoted as BD(T)).51 All the above calculations were performed using Gaussian 09 program package.52 To test methods beyond CCSD(T), single point energy calculations were carried out using variant D of the completely renormalized single-reference coupled-cluster method with singles, doubles, and a noniterative treatment of triples, denoted as CR-CC(2,3).53−55 These calculations, based on restricted open shell HF wave functions, were carried out with the GAMESS program package.56,57 Furthermore, the single point energies were calculated using the coupled-cluster method with singles, doubles, and triples (denoted as CCSDT)58 using the CFOUR program.59 Rate constants were computed using the GAUSSRATE60 and POLYRATE61 codes. Conventional transition-state theory (TST), which was applied by Bofill et al.,27 assumes the saddle point is the transition state and neglects the effect of recrossings on the reaction rate; therefore, TST overestimates the rate constants. Variational transition-state theory (VTST) locates the transition state dividing surface that minimizes the rate constant. Canonical VTST (CVT) locates the variational transition state in a canonical ensemble. Improved canonical variational theory (ICVT) treats the recrossing effect more accurately than CVT while saving computational cost relative to
a fully microcanonical approach; the ICVT approach was used in this work. Variational effects can significantly change rate constants; for example, Setokuchi and Sato computed a 20% reduction in the rate constant for i-C3H7O + O2 due to the variational effect.30 The tunneling coefficient κ(T) is defined via the relationship k(T )QM = k(T )classical × κ(T )
where k(T)QM is the quantum mechanical rate constant (including tunneling) and k(T)classical is the TST or, here, the ICVT rate constant. Two semiclassical approaches are typically used to investigate multidimensional tunneling: the small curvature (SCT) approximation assumes the curvature of the reaction path is small to medium and the large curvature (LCT) approximation assumes the curvature of the reaction path is medium to large.62 The optimized microcanonical multidimensional tunneling (μOMT) computes κ(T) as a convolution of a Boltzmann distribution with the larger of κ(E)SCT and κ(E)LCT.63 Zero curvature tunneling (ZCT), which ignores the curvature of the reaction path and therefore tends to underestimate the extent of tunneling, is also calculated as a comparison.62 Because one of the goals of this work is to determine a method that can be easily applied to a series of alkoxy + O2 reactions, tunneling effects were also explored with the asymmetric Eckart potential.64
III. RESULT AND DISCUSSION III.A. Saddle Points and Reaction Mechanisms. A cyclic saddle point (TS1) with Cs symmetry and A″ electronic state was found for all ab initio methods and DFT methods. In addition to TS1, a second saddle point of C1 symmetry (TS1′) was found with all ab initio methods and most of the DFT methods. Selected geometric parameters of TS1 and TS1′ are shown in Figure 1, and Table 1 lists relative energies, imaginary
Figure 1. Geometry of TS1 and TS1′ at UCCSD(T)/6-311G(d,p), B3LYP/6-311g(d,p), and CASSCF(9,7)/6-311g(d,p) from Bofill et al.,27 and ROCCSD(T) for TS1 only. Selected internuclear distances are specified in Ångstroms and the dihedral angles (D) are in degrees.
frequencies, and the length of the breaking C−H and forming H−O bonds. Where both TS1 and TS1′ were found, TS1 is computed to be more stable by 2−9 kcal/mol than TS1′. In contrast to the result of Bofill et al. at CASSCF(9,7),27 all our results suggest that both saddle points are reactant-like, having a shorter C−H distance than O−H distance, as Table 1 shows. The energy barriers of TS1 computed using unrestricted quadratic configuration interaction (QCI) methods and coupled cluster (CC) methods are 13−20 kcal/mol relative C
dx.doi.org/10.1021/jp409105q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
Table 1. Activation barriers (E, kcal/mol, Including the Zero-Point Energy (ZPE) Correction), Imaginary Frequencies (ν*, cm−1), and Their Selected Geometric Parameters (RC−H and RH−O, Ångstroms) for TS1 and TS1′ Using a Range of Quantum Chemical Methods and the 6-311G(d,p) Basis Seta TS1 method
E
CASSCF(9,7)b UQCISD UQCISD(T) UCCSD UCCSD(T) ROCCSD(T) PBE1W MPWLYP1W BLYP OLYP M06L B97D B3LYP X3LYP O3LYP M062X B2PLYP LC-BLYP LC-OLYP LC-M06L CAM-B3LYP
39.6 19.1 c 19.6 13.8 8.7 −9.5 −8.5 −7.1 −0.4 −0.7 −2.7 6.2 6.4 6.2 11.3 12.3 11.5 12.1 6.1 9.6
TS1′
ν*
RC−H
RH−O
E
ν*
RC−H
RH−O
2301 1653 2327 1666 1309 507 621 611 473 1011 535 1043 1176 620 1860 1848 2000 2003 2290 1638
1.349 1.286 1.253 1.288 1.255 1.278 1.265 1.272 1.271 1.240 1.252 1.252 1.243 1.247 1.228 1.276 1.264 1.277 1.278 1.283 1.263
1.229 1.346 1.403 1.349 1.403 1.364 1.399 1.401 1.404 1.434 1.384 1.445 1.447 1.438 1.463 1.399 1.409 1.362 1.350 1.335 1.397
43.3 22.1 c 22.6 18.6
3221 2970 2774 3065 2769
1.321 1.292 1.276 1.295 1.275
1.267 1.338 1.372 1.337 1.372
−0.3 −0.2 0.7 3.6 4.9 1.7 9.1 9.3 8.0
378 484 510 655 1383 792 1769 1830 1395
1.190 1.199 1.204 1.210 1.228 1.220 1.263 1.265 1.245
1.624 1.597 1.588 1.535 1.433 1.527 1.400 1.392 1.424
N/A N/A N/A N/A N/A N/A
a
For several functionals where TS1′ was not found to be a saddle point, and this is indicated by N/A in the table. bResults of Bofill et al.27 cThis result is unavailable due to the inability to obtain CH3O that is a potential energy minimum at UQCISD(T).
to the reactants. These barrier heights are much larger than the experimental derived activation energy of 2.3−2.5 kcal/mol (at 298−610 K), and this difference is consistent with TS1 being a multireference problem. The fact that CASSCF(9,7) produces a barrier height of nearly 40 kcal/mol (see Supporting Information for ref 27) implies that, in addition to nondynamical electron correlation, dynamic electron correlation is critical to a proper treatment of this system. This should be taken into account when one considers whether the saddle point is as product-like as suggested by CASSCF calculations or more reactant-like, as indicated by all the single reference methods used here. The UQCISD and UCCSD methods agree with each other as to the barrier heights and imaginary frequencies. UQCISD(T) and UCCSD(T) agree with each other, also; however, the use of the perturbative estimate of the triples contribution decreases both the barrier height (by 5−6 kcal/mol) and the imaginary frequency (by ∼650 cm−1). This seems to imply the importance of dynamic electron correlation due to the connected triples; alternatively (or additionally), it may correspond to the connected triples calculation incorporating a significant amount of the nondynamical correlation of TS1. The barrier height calculated at ROCCSD(T) is in much better agreement with experiment than the UCCSD(T) results. The structure of TS1 calculated using ROCCSD(T) is slightly more product-like than the UCCSD(T) structure, but much more reactant-like than the CASSCF results. Due to the computational expense, geometry optimizations were not carried out at UQCISD(T) and UCCSD(T) with basis sets larger than 6-311G(d,p). However, basis set effects on UCCSD(T) activation barriers were explored by calculating the UCCSD(T) single point energies and extrapolating to the complete basis set (CBS) limit. In this exercise, geometries of the reactants and the saddle points were obtained at the three
levels of theory: UCCSD(T)/6-311G(d,p), UQCISD(T)/6311G(d,p), and B3LYP/6-311G(d,p). UCCSD(T)/CBS barrier heights vary by less than 0.3 kcal/mol (see Supporting Information) with variations in molecular structure. In contrast, using a larger basis set lowers the UCCSD(T) barrier height from 14.0 kcal/mol at 6-311G(d,p) to 10.0 kcal/mol at the CBS limit (without ZPE corrections). The ZPE correction is very small: −0.2 and +0.4 kcal/mol at UCCSD(T)/6311G(d,p) and UB3LYP/6-311G(d,p), respectively. This small correction still leaves the UCCSD(T)/CBS barrier height far higher than the activation energy of about 2.5 kcal/mol. As will be argued later in this paper, tunneling does not appear to enormously increase the rate constant for this reaction, so that this activation barrier should be close to the barrier height. As a result, it appears that UCCSD(T)/CBS badly overestimates the activation barrier for the methoxy + O2 reaction: by about 7 kcal/mol or a factor of ∼4. In an effort to improve the accuracy of the computed barrier height over the QCISD(T) and CCSD(T) results, the single point energies were computed using CR-CC(2,3), BD(T), and UCCSDT with geometries of reactants and TS1 optimized at both UCCSD(T)/6-311G(d,p) and UB3LYP/6-311G(d,p). CR-CC(2,3) uses a restricted open-shell Hartree−Fock wave function and is known to perform well for bond breaking processes and biradicals.54 BD(T) is based on Brueckner orbitals and can compensate to some degree for a poor Hartree−Fock wave function.65 CCSDT is the most complete of the standard coupled cluster methods we could afford. Results for all these methods, shown in Table 2, have barrier heights lower than UCCSD(T) values and similar to the ROCCSD(T) results, but still significantly higher than the experimentally derived activation energy of ∼2.5 kcal/mol. D
dx.doi.org/10.1021/jp409105q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
therefore, the wave function with lower CISD energy should be the stable one. The result, listed in Table 3S9 in the Supporting Information, confirms that the stable wave function is correct. As a further check, we computed CISD energies of TS1 using both UHF wave functions (stable and metastable) at the CCSD(T) geometries optimized using both wave functions. At both geometries, the UCISD (and UHF) energies were lower when the stable wave function was used. The UCCSD(T) correlation energy for the metastable wave function is larger than that for the stable wave function; this, rather than differences in geometries, accounts for most of the difference between our UCCSD(T) energies and those of Bofill et al. These results tend to confirm that the UCCSD(T) barrier height obtained by Bofill et al. is an artifact. This is despite the fact that the metastable wave function helps UCCSD(T) obtain good agreement with the experimental activation energy. Recall that the efficiency of DFT is needed for multidimensional tunneling calculations, which motivates our consideration of DFT methods despite that the failure of CCSD(T) implies DFT results will not generally be reliable. As can be seen from Table 1, results from the GGA and meta-GGA functionals (PBE1W, MPWLYP1W, M06L, BLYP, OLYP, and B97D) indicate that the energy of the cyclic saddle point (TS1) is below the reactant. Intrinsic reaction coordinate (IRC) calculations indicate that, in these cases, the minimum energy path (MEP) connects TS1 to the cis-methyltrioxy radical rather than isolated reactants or a weakly bound complex. The methyltrioxy radical is bound with respect to separated reactants at these levels of theory. In these cases, TS1 corresponds to a unimolecular elimination of HOO radical rather than a bimolecular hydrogen-abstraction reaction. In contrast to the GGA and meta-GGA functionals, Table 1 shows the energy of TS1 computed with hybrid and long-range corrected hybrid functionals is 6−12 kcal/mol higher in energy than reactants. The imaginary frequency of TS1 computed with these hybrid functionals varies enormously with the method (from 620 to 2290 cm−1). IRC analysis of TS1 using these functionals indicates the reactant is either separated CH3O + O2 or the prereactive complex noted by Bofill et al. IRC analyses were not carried out with QCI and CC methods. The effect of basis set size was explored above for CCSD(T) calculations, and here we briefly discuss basis set effects on geometries, relative energies, and the imaginary vibrational frequency. The results of an extensive set of calculations at B3LYP, shown in Table 4, suggest that the geometric parameters of TS1 are not sensitive to the basis set: the changes of the C−H and H−O bond span only 0.005 and 0.012 Å, respectively. As to the computed barrier height, adding a diffuse function caused a significant increase, whereas adding higher angular momentum functions changes the barrier height only modestly. Basis set effects were also tested with the O3LYP, LC-M06L, and M06-2X functionals using the 6311G(d,p) and 6-311G(3df,3pd) basis sets. The results, also shown in Table 4, are consistent with those obtained using B3LYP. The T1 diagnostic in QCI and CC methods is one measure of the importance of nondynamical electron correlation effects.68 Both saddle points have high values of T1 diagnostic: 0.050−0.053 for TS1 and 0.043−0.047 for TS1′. Lee and Scuseria suggested that a calculation of T1 diagnostic of 0.02 or higher should be viewed with caution.68 The value obtained here for the T1 diagnostic of TS1 is much higher than this, and comparable to that of the CN radical, which is notorious for
Table 2. Barrier Heights to TS1 (kcal/mol, without ZPE Corrections) of the CH3O + O2 Reaction Calculated at UCCSDT, BD(T), and CR-CC(2,3)a method for optimizing geometries barrier heights
UCCSD(T)
UB3LYP
CR-CC(2,3)/cc-pVDZ CR-CC(2,3)/cc-pVTZ BD(T)/cc-pVDZ BD(T)/cc-pVTZ CCSDT/cc-pVDZ
9.1 8.8 9.4 8.4 7.3
8.7 8.6 9.3 8.6 6.9
a
The geometries of reactants and saddle point were optimized by UCCSD(T)/6-311G(d,p) and UB3LYP/6-311G(d,p).
Clearly, this reaction is a challenging problem for modern electronic structure theory. Bofill et al. found the barrier height of the cyclic saddle point (TS1) to be 3.7 kcal/mol, using UCCSD(T)/cc-pVTZ single point energies with the CASSCF(9,7)/6-311G(d,p) optimized geometry.27 It seemed unusual that the UCCSD(T) single point energy of the CASSCF geometry results in a reasonable barrier height, whereas our UCCSD(T) activation barrier is far too high at UCCSD(T) and UB3LYP geometries. Therefore, we tried to reproduce the UCCSD(T) single point energy of Bofill’s CASSCF optimized geometry. Our results, listed in Table 3, show that the UCCSD(T)/cc-pVTZ single point Table 3. Relative Energies of Saddle Points (kcal/mol, without ZPE) Using Coupled Cluster Methods at CASSCF(9,7)/6-311G(d,p) Optimized Geometries from Bofill et al. method
TS1
TS1′
Bofill et al. UCCSD(T)/cc-pVTZ//CASSCF UCCSD(T)/cc-pVTZ(stable w.f.)//CASSCF UCCSD(T)/cc-pVTZ(metastable w.f.)//CASSCF
3.7 7.0 4.3
11.6 15.3 11.7
energy of TS1 at the CASSCF-optimized geometry is 7.0 kcal/ mol: 3.3 kcal/mol higher than Bofill et al.’s value. Because we had already verified the stability of the UHF wave function used in our CCSD(T) calculations, we suspected that the CCSD(T) results of Bofill et al. were not based on the lowest energy solution of the UHF equations, in other words, were the result of a metastable wave function. In fact, we very nearly reproduce their UCCSD(T)/cc-pVTZ relative energies with a Hartree− Fock wave function that is higher in energy than the stable wave function. Table 3 provides a comparison of barrier heights with both the stable and metastable UHF wave function for both TS1 and TS1′. The results above suggest that the lower barrier height obtained by Bofill et al. is based on this metastable UHF wave function. The expectation value of the total spin, ⟨S(S + 1)⟩, of TS1 is 0.87 for the metastable Hartree−Fock wave function whereas it is around 1.6 for the stable wave function. This spin contamination in a nominally doublet state appears to be typical of the presence of a nearly degenerate quartet state.66,67 Bofill et al. reported similar values of ⟨S(S + 1)⟩ at QCISD but did not specifically report the value of ⟨S(S + 1)⟩ of the UHF wave function obtained from UCCSD(T) calculations at their CASSCF geometry. To confirm the validity of the stable UHF wave function, we computed the energy of the CASSCF structure of TS1 at CISD/6-311G(d,p) using both wave functions. Unlike coupled cluster methods, CISD is variational; E
dx.doi.org/10.1021/jp409105q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
Table 4. Basis Set Effects on the Barrier Height (E, kcal/mol, Including the ZPE Correction), Imaginary Frequencies (ν*, cm−1), and Selected Geometric Parameters (RC−H and RH−O, Ångstrom) of the Cyclic Saddle Point (TS1) at Several Density Functionals B3LYP/6-31+G(d,p) B3LYP/6-311G(d,p) B3LYP/6-311+G(d,p) B3LYP/6-311G(3df,2p) B3LYP/6-311+G(3df,2p) B3LYP/cc-pVTZ B3LYP/aug-cc-pVTZ B3LYP/6-311G(3df,3pd) O3LYP/6-311G(d,p) O3LYP/6-311G(3df,3pd) LC-M06L/6-311G(d,p) LC-M06L/6-311G(3df,3pd) M06-2X/6-311G(d,p) M06-2X/6-311G(3df,3pd)
E
ν*
RC−H
RH−O
7.7 6.2 7.8 5.9 7.3 6.6 7.6 5.5 6.2 5.5 6.1 4.7 11.3 10.4
1060 1042 1108 993 1038 1001 1015 968 620 540 2290 2229 1860 1846
1.241 1.243 1.242 1.237 1.236 1.237 1.237 1.233 1.228 1.218 1.283 1.278 1.276 1.272
1.446 1.447 1.455 1.452 1.458 1.452 1.455 1.457 1.463 1.473 1.335 1.340 1.399 1.403
Scheme 3. Examples of Alkoxy Radicals Derived from Isoprene (Top) and Esters or Ketones (Bottom) Reacting with O2 and Forming a Conjugated π System
in Scheme 3. This is likely to have a significant effect on the electronic structure of the saddle point, and it is not clear that the success of G2M(RCC1) will extend to such systems. This point requires further investigation. III.B. Methyl Trioxy Radical and the Prereactive Complex. Cis and trans conformers of methyl trioxy radicals, both of which are Cs symmetry, were found as local minima for all methods and are shown in Figure 2. Some methods,
multireference effects.69 This confirms the point, suggested by Bofill et al., that a single reference wave function is a poor approximation of the true electronic wave function, and a multireference approach may be necessary to ensure accurate relative energies, frequencies, and geometries. This is certainly consistent with the failure of CCSD(T) to yield a reasonable activation barrier. However, as noted above, CASSCF grossly overestimates the activation barrier, which indicates that both dynamical and nondynamical correlations are critical to reliably characterizing TS1. In light of the multireference nature of TS1, it is worth examining how Setokuchi and Sato obtained rate constants close to the experimental values using the composite G2M(RCC1) approach. As noted previously, they found that tunneling increases the rate constant by only a factor of 2 at 298 K and the variational effect is very modest. If these conclusions are correct, they clearly computed a very accurate activation barrier. We repeated their calculations to better understand how the composite method succeeded where our calculations failed. Using the ROCCSD(T) approach 70 coded in Gaussian09 rather than the one71 used by Setokuchi and Sato, we reproduce the relative energies to within 0.1 kcal/mol. The complete set of G2M(RCC1) energies of reactants and TS1 are provided in the Supporting Information, and only the main result is discussed here. The ROCCSD(T)/6-311G(d,p) calculation obtains a barrier height of 8.3 kcal/mol (neglecting the small effect of ZPE), far lower than the UCCSD(T) value obtained using the same basis set. Note that Bofill et al., obtained ROCCSD(T) results in good agreement with their UCCSD(T) results at the same basis set. The success of the G2M(RCC1) approach for the CH3O + O2 reaction is clearly due, in large part, to the success of ROCCSD(T) in handling the nondynamical correlation in the wave function of TS1. It should be noted that this approach obtains reasonable rate constants for ethoxy and 2-propoxy radical reacting with O2, so this success might extend to studies of other alkane-derived alkoxy radicals, or for radicals with substituents that do not significantly affect the electronic structure of TS1. Unfortunately, many of the most relevant examples of alkoxy radicals (from isoprene,14 ketones, and esters) reacting with O2 lead to conjugated π bonds, as shown
Figure 2. Geometries of cis-methyl trioxy radicals, trans-methyl trioxy radicals, gauche-methyl trioxy radicals, and prereactant complex for various levels of theory using the 6-311G(d,p) basis set. Selected internuclear distances are listed in Ångstroms. Note the two different structures of the cis-methyl trioxy radical.
including UQCISD and UCCSD but not UQCISD(T) or UCCSD(T), also obtained a gauche conformer with C1 geometry as a local minimum (Table 5). Both cis and trans conformers have A″ electronic ground states. In contrast to the HOOO• radical, for which the trans conformer is favored over the cis,72 the cis-methyl trioxy radical is consistently computed to be more stable than the trans-methyl trioxy radical. As listed in Table 5, the cis-methyl trioxy radical has a longer O2−O3 bond than the trans-methyl trioxy radical, except at M06L. F
dx.doi.org/10.1021/jp409105q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
Table 5. Relative Energies (E, kcal/mol, Including the Zero-Point Energy (ZPE) Correction), O2−O3 Distances (R, Ångstroms) and Dihedral Angle C4HpO1O2 (ϕ, deg) of the Pre-reactive Complex and Methyl Trioxy Radicals Using Various Quantum Chemical Methodsa cis-trioxy
complex
trans-trioxy
gauche-trioxy
methods
E
R
E
R
ϕ
E
R
E
R
CASSCFb UQCISD UQCISD(T) UCCSD UCCSD(T) PBE1W MPWLYP1W BLYP OLYP M06L B97D B3LYP X3LYP O3LYP M062X B2PLYP LC-BLYP LC-OLYP LC-M06L CAM-B3LYP
−0.7 −0.2
3.579 3.204
8.9 15.6
1.515 1.556 1.564 1.549 1.564 1.797 1.803 1.804 1.798 1.872 1.813 1.614 1.600 1.662 1.481 1.541 1.453 1.442 1.448 1.523
180 180 180 180 180 0 0 0 0 0 0 180 180 180 180 180 180 180 180 180
10.0 16.3 c 17.1 11.7 −5.8 −5.6 −4.2 1.8 0.1 0.4 7.5 7.7 7.8 10.5 12.1 7.2 8.5 3.9 9.3
1.499 1.489 1.537 1.489 1.536 1.703 1.717 1.719 1.714 2.267 1.723 1.588 1.574 1.639 1.470 1.530 1.445 1.435 1.436 1.507
9.8 16.3
1.477 1.483
c
N/A −0.5
3.177 N/A N/A N/A N/A N/A N/A N/A
−0.9 −0.5 −0.4 −1.0 −0.5 −1.3 −1.3 −8.7 −0.7
2.736 2.669 3.366 2.587 3.149 2.574 2.598 2.347 2.699
17.0 10.9 −9.9 −9.4 −7.9 −0.8 −1.7 −3.0 6.4 6.6 6.9 9.9 10.9 7.1 8.4 1.8 8.7
N/A 17.1
1.483 N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A
7.2 8.5 2.3
1.436 1.426 1.430 N/A
a
The 6-311G(d,p) basis set is used throughout. N/A refers to cases where geometries could not be optimized at the particular level of theory. b CASSCF result of Bofill et al. cThis result is not available due to the unexpected failure of CH3O frequency calculation at UQCISD(T).
Trioxy radicals have energies ranging from −10 to +17 kcal/ mol with respect to CH3O + O2, depending on the choice of methods. We verified that the wave functions were stable in all cases. Furthermore, the orientations of the methyl group in the cis conformer differ among the methods: For methods that find trioxy radical are stable with respect to reactants (e.g., PBE1W, etc.), the hydrogen atom that is in the plane of symmetry (identified as Hp below and in Figure 2) is aligned to the terminal O (O1), making an HpC4O2O1dihedral angle of 0°. In contrast, methods indicating that the methyl trioxy radical is energetically unstable with respect to reactants (including QCI and CC methods and hybrid functionals) predict the hydrogen atom to be in the plane of symmetry orientated away from the terminal oxygen atom (dihedral angle of HpC4O2O1 = 180°) and the terminal oxygen atom to be equidistant from the two nearest hydrogen atoms. The latter result agrees with that of Bofill et al. at CASSCF(9,7)/6-311G(d,p) level. For some of the methods that predict an unstable methyl trioxy radical (hybrid DFTs, QCI, and CC), a prereactive complex was also located. Such a complex was also found by Bofill et al. at their CASSCF(9,7) calculations. We distinguish the complex from the cis-methyl trioxy radical primarily by the difference in the O2−O3 bond length: structures with shorter O2−O3 bond (2.6 Å) are referred to as complexes. The complex is only bound by between −0.2 and −0.9 kcal/mol relative to the isolated reactants, except at LC-M06L (−8.7 kcal/mol). The strong bonding computed with this method for long-range interactions (∼2.3 Å) appear unphysical, and the binding energy for the complex is almost certainly an artifact. Hybrid DFT calculations show the prereactive complex has Cs symmetry and an A″ electronic state. In contrast, UQCISD and UCCSD calculation predict that the complex has C1 symmetry.
Selected geometric parameters of the complex are shown in Figure 2 and Table 5, alongside those of methyl trioxy radicals. Note that Bofill et al. obtained a high barrier (50 kcal/mol) for HOO elimination from the methyl trioxy radical, which argues strongly against the trioxy radical being an intermediate in the CH3O + O2 reaction.27 It is interesting to compare the degree to which HOOO• and CH3OOO• challenge single reference methods, and if anything can be learned from the comparison. Figure 3 depicts the RO− O2 bond energy (R = H, CH3) versus the RO−O2 bond length using selected methods. CH3OOO• is computed to be less strongly bound than HOOO• at each level of theory. Given that HOOO• is only bound by 2.9 ± 0.1 kcal/mol (D0),73 it is reasonable to wonder if CH3O−OO• is bound by any chemically significant amount.
Figure 3. RO−O2 bond energy (kcal/mol, including ZPE) vs RO−O2 bond length (Ångstrom) for trans-HOOO• and cis-CH3OOO• using selected methods (labeled in italics for cis-CH3OOO•). G
dx.doi.org/10.1021/jp409105q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
Table 6. Variational Effect on the Location of Saddle Point TS1 (sICVT(T), Ångstroms) and the Rate Constants (Expressed as kICVT/kTST) of Methoxy Radicals Reacting with O2 at B3LYP/6-311G(d,p)a temperature (K) CH3O CD3O
a
CH2DO
H abs
CH2DO
D abs
sICVT(T) kICVT/kTST sICVT(T) kICVT/kTST sICVT(T) kICVT/kTST sICVT(T) kICVT/kTST
150
200
250
298
350
−0.0635 0.794 −0.0140 0.987 −0.0781 0.751 −0.0123 0.988
−0.0496 0.890 −0.0065 0.998 −0.0630 0.864 −0.0056 0.998
−0.0358 0.946 0.0006 1.000 −0.0455 0.932 0.0011 1.000
−0.0244 0.975 0.0069 0.998 −0.0309 0.968 0.0070 0.998
−0.0142 0.992 0.0128 0.994 −0.0186 0.988 0.0129 0.993
“H abs” and “D abs” refer to hydrogen abstraction and deuterium abstraction, respectively in the reaction of CH2DO + O2.
III.C. Variational and Tunneling Effects on Kinetics. Despite the wide range of computed barrier heights and the significant variation in the computed structures of TS1, we proceed to compute the effects of a variational transition state and multidimensional tunneling on rate constants. A DFT method is commonly used to investigate multidimensional tunneling effects, because higher-level approaches such as QCI and CC methods are too time-consuming for the repeated computation of Hessians that is needed in kinetic studies using POLYRATE. We follow this approach here. Given the failure of all the DFT methods used here to obtain a reasonable barrier height, we do not focus on absolute values of the rate constant. Instead, we will use the experimentally determined Arrhenius pre-exponential factor and deuterium kinetic isotope effects to evaluate whether calculations accurately treat tunneling for the methoxy + O2 reaction. The higher the barrier, the higher the imaginary frequency will be, so methods with a large barrier would seem to have artificially high values of the imaginary frequency. Compared with the experimental activation energy of ∼2.3 kcal/mol, B3LYP, X3LYP, and LC-M06L seem to be the most reasonable functionals (Table 1). X3LYP and B3LYP produce very similar barrier heights and imaginary frequencies; therefore, it is expected that these two functionals would yield very similar rate constants. Of these two functionals, we chose B3LYP to study the kinetics of methoxy + O2 reactions. Imaginary frequencies computed at LC-M06L agree well with CCSD and QCISD results but are much higher than CCSD(T) and QCISD(T) results. This is due to the unphysically strong binding energy of the complex (8.7 kcal/mol), which leads to an effective barrier height of ∼15 kcal/mol. Therefore, we did not investigate kinetics using the LC-M06L potential energy surface. As shown in Table 6, for the CH3O + O2 reaction, the B3LYP/6-311G(d,p) approach suggests that the variational effect is small at 298 K (
