- Email: [email protected]
Theoretical study on kinetics of the H2CO + O2 → HCO + HO2 reaction
Chemical Physics Letters 469 (2009) 81–84
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Theoretical study on kinetics of the H2CO + O2 ? HCO + HO2 reaction Lam K. Huynh a,1, Max Tirtowidjojo b, Thanh N. Truong a,* a b
Department of Chemistry, University of Utah, 315 S 1400 E, Rm. 2020, Salt Lake City, UT 84112, United States Dow Chemical Company, 2301 N Brazosport Blvd., B-1226, Freeport, TX 77451, United States
a r t i c l e
i n f o
Article history: Received 21 October 2008 In final form 9 December 2008 Available online 24 December 2008
a b s t r a c t A theoretical study on the kinetics of the H2CO + O2 ? HCO + HO2 reaction has been carried out from first principles using canonical variational transition state theory incorporating corrections from small curvature tunneling and hindered rotation treatments. Rate constants were calculated using the geometries, gradient and Hessian obtained at the BH&HLYP/cc-pVDZ level with the energy correction at the CCSD(T)/aug-cc-pVTZ//BH&HLYP/cc-pVDZ level. The rate constant expression is obtained as kðtÞ ¼ 2:95 1021 T 3:152 exp 18 T352 (cm3 molecule1 s1) for the temperature range 300–3000 K. The result shows an excellent agreement with available scattered experimental data over this temperature range and thus the expression can be used confidently for the whole range. Ó 2008 Published by Elsevier B.V.
1. Introduction Formaldehyde has recently generated serious environmental concerns because of its large accumulation in the atmosphere due to its thermodynamic stability [1,2]. It has also long been recognized as an important kinetic intermediate in the combustion of many hydrocarbon fuels, particularly of methane and methanol [3–6]. Recently, Michael et al. [7] has shown that an oxidation reaction, namely, H2CO + O2 ? HCO + HO2 is an essential component for the kinetic modeling of hydrocarbon combustion. As a result, it is included in various kinetic combustion models [8,9]. Unfortunately, there have been only seven recorded data of thermal rate constants for this reaction. Only one record from direct kinetic measurement data [10] is available on this elementary oxidation reaction H2CO + O2 ? HCO + HO2, which is in the small temperature range of 713–816 K and at a pressure of 0.07–0.08 bar and dates back from 1974; others are from theoretical estimates: three derived indirectly from complex mechanisms fit to experimental observations from the combustion of formaldehyde [11–13] and one from the CH3 + O2 reaction [7], two from the extensive literature review [14,15], one from a preliminary transition state theory (TST) calculation [13]. In the present study, we performed accurate theoretical study on rate constants for the H2CO + O2 ? HCO + HO2 reaction. Particularly, the rate constants are calculated using the CVT [16–19] with the inclusion of corrections of tunneling and hindered rotation treatments and the potential energy surface was corrected at the CCSD(T)/aug-cc-PVTZ level. The calculated rate constants are then
* Corresponding author. Fax: +1 801 581 4353. E-mail address: [email protected] (T.N. Truong). 1 Current address: Chemical Engineering Department, Colorado School of Mines, Golden CO 80401, Unites States. 0009-2614/$ - see front matter Ó 2008 Published by Elsevier B.V. doi:10.1016/j.cplett.2008.12.050
validated by comparison to those available in the literature so that rate constant expression can be used to obtain rate constants in the temperature range where there is no reliable rate constant available. 2. Electronic structure calculation All electronic structure calculations were carried out using the program package GAUSSIAN 03 [20]. It has been shown earlier that the combination of Becke’s Half-and-Half (BH&H) non-local exchange with Lee–Yang–Parr (LYP) correlation functionals can be used cost-effectively to predict the properties, e.g. geometries and frequencies of reactants, products, as well as transition states for hydrogen abstraction reactions by a radical [21–26]. The Dunning’s correlation-consistent polarized valence double-zeta basis set [3s2p1d/2s1p] denoted as cc-pVDZ [27] was employed with the DFT method to determine structural and vibrational frequencies for the stationary points. The quadratic configuration interaction method with inclusion of singles and doubles excitations (QCISD) using the same cc-pVDZ basis set is also used as a reference source for the DFT method in finding stationary structures as well as calculating vibrational frequencies. The minimum energy path (MEP) connecting the reactants and products is also obtained at the same BH&HLYP/cc-pVDZ level by tracing the steepest descent path from the transition state to the corresponding reactants and products using the Gonzalez–Schlegel method [28] with a step size of 0.1 (amu)1/2 bohr for a total of 200 steps. This step size is found to be sufficient for the convergence of the variational rate constant as mentioned in our previous study. Force constants at 22 selected points (12 points in the reactant channel and 10 points in the product channel) along the MEP were carried out to obtain the necessary potential energy surface information for CVT calculations. These points were found using our
82
L.KK. Huynh et al. / Chemical Physics Letters 469 (2009) 81–84
autofocusing technique [29], which uses a combination of the barrier shape, and the second derivatives of the Z-matrix geometrical parameters to predict the regions along the MEP that are sensitive in rate determination. In order to improve the energy information along the MEP, single point calculation for these points at a higher level of theory such as [CCSD(T)] is needed, in particular the CCSD(T) with aug-cc-pVTZ basis set was used. 3. Results and discussion Potential energy surface BH&HLYP optimized geometries of the reactants (H2CO and O2), products (HCO and HO2), the transition state as well as the van-der-Waal complex in the product channel are shown in Fig. 1. The experimental data available in the litera-
1.108 1.119 (1.116) 1.182 1.208 (1.208)
115.4° 115.2° (116.5°)
1.191 1.212 (1.208)
REACTANTS
1.306 1.342
105.6° 103.9°
124.6° 124.4°
1.126 1.137
1.166 1.188
0.967 0.977
PRODUCTS
90.7° 92.0°
144.5° 145.3°
2.064 2.172
0.979 0.984
COMPLEX
ture are also given in parentheses. Compared to the experimental data [30], the BH&HLYP/cc-pVDZ yields a shorter bond length for oxygen and formaldehyde molecules with the largest difference of 0.026 Å for O2 (0.008 Å and 0.017 Å for C–H and C@O bonds in H2CO, respectively). The transition state has the product-like structure with the 1C–4H bond distance of 1.440 Å, which is much larger than the C–H bond length of 1.108 Å in H2CO and the 4H–5O bond distance is very close to the O-H bond length of the OOH product (1.137 Å compared to 0.967 Å). The geometry of the complex is very close to that of the transition state except the 1C–4H bond length is longer and the 5O–4H is shorter. In addition, the angle (3H–1C–4H) decreases by around 10° when the transition state becomes the complex in the exit channel. BH&HLYP harmonic vibrational frequencies for the reactants, products, transitions state and complex are listed in Table 1. It can be seen that the BH&HLYP frequencies are very close to the available experimental data. For formaldehyde, the average difference is 136 cm1. The difference mainly comes from the symmetric and asymmetric H–O stretching. For oxygen molecule, the predicted stretching frequency is 238 cm1 larger than the experimental value from Ref. [30]. It is noted that the scaling factor of 0.97 was applied for these frequencies; this will reduce the gap between the calculated and experimental data. By comparison to the QCISD results, which are also given in Figs. 1 and 2, BH&HLYP predicts similar geometries of the stationary structures. In particular, for the reactants and products, the geometries at the two levels are almost the same with the maximum bond length difference of 0.035 Å for the O–O bond in OOH radical. The bond angles between the two levels are almost identical. For the transition state, the geometries at the two levels are still very close to each other except for the bond length between the two oxygen atoms for which the BH&HLYP bond length is 0.03 Å shorter than that using the QCISD level of theory. Such a difference was also found for the OOH radical. Since the O–O bond is only a spectator in the H2CO + O2 reaction, we can expect this difference to have only a minor effect on the calculated rate constants. Regarding frequencies, as shown in Table 1, BH&HLYP yields almost identical results for reactants, products and even for transition state compared to those from QCISD calculations. Specifically, BH&HLYP/cc-pVDZ tends to give slightly higher frequencies. The largest differences are found for the O–O stretches in the oxygen molecule and in the transition state, which are around 150 cm1. However, the contribution of the differences to the rate constants, especially the vibrational partition function, is small due to the large magnitude of these frequencies. Thus, it can be concluded that the BH&HLYP is a cost-effective method for structural and frequency information for this reaction. The schematic energy profile for the title reaction is shown in Fig. 2. The most accurate energetic results from this study are those
Table 1 Harmonic vibrational frequencies (in cm1) for the transition state (TS), reactants, products, and complex.
1.266 1.299
105.1° 103.9°
1.137 1.155
153.7° 152.1°
1.123 1.135
101.4° 100.2°
O2 H2CO
1.440 1.431
1.165 1.186
HCO OOH TS
TRANSITION STATE Complex Fig. 1. Optimized geometries (distances in Å and angles in degrees) of the reactants, products, complex and transition state at the BH&HLYP/cc-pVDZ and QCISD/ccpVDZ levels (italic numbers). The numbers in the parentheses are the experimental values from Ref. [30].
a
BHLYP/cc-pVDZ
QCISD/cc-pVDZ
Experimentala
1818 1247, 1302, 1570, 1923, 2999, 3068 1140, 2043, 2745 1250, 1502, 3781 2033i, 107, 124, 319, 429, 613, 1108, 1157, 1406, 1673, 2045, 2832 80, 109, 140, 216, 266, 574, 1095, 1272, 1560, 2057, 2807, 3509
1645 1193, 1275, 1544, 1816, 2944, 3004 1135, 1894, 2679 1086, 1435, 3675 2286i, 112, 149, 323, 419, 595, 1088, 1154, 1226, 1655, 1919, 2769 67, 101, 125, 184, 230, 492, 1095, 1102, 1486, 1922, 2747, 3531
1580 1167, 1249, 1500, 1746, 2783, 2843 1081, 1868, 2485 1098, 1392, 3436
From Ref. [30].
N/A
N/A
83
L.KK. Huynh et al. / Chemical Physics Letters 469 (2009) 81–84
TRANSITION STATE
5.60 4.75 4.68
-10 HCO+HO2
Baulch et al. (1992)
COMPLEX 40.47 BH&HLYP/cc-pVDZ 44.03 QCISD/cc-pVDZ 40.83 CCSD(T)/aug-cc-pVTZ
41.10 42.05 40.19
H2CO+O2 Fig. 2. Schematic energy profile of the reaction H2CO + O2, at the different levels of theory as indicated. Numbers are in kcal/mol.
Log{k(T)/cm3molecule-1s-1}
-15
Tsang et al. (1986) Baldwin et al. (1974) Michael et al. (2005)
-20
This work: CVT/SCT/HR
-25
-30
-35
-40
from single-point energy CCSD(T)/aug-cc-pVTZ calculations at the BH&HLYP/cc-pVDZ optimized geometries. The H2CO + O2 reaction is endothermic with the reaction energy of 40.19 kcal/mol. The calculated classical barrier for this reaction is 40.83 kcal/mol which is 0.64 kcal/mol above the exit channel. The weak binding complex in the exit channel is only 4.68 kcal/mol below the dissociative products. Since the reaction would occur at rather high temperatures due to its high barrier, one can expect that if the trajectory has sufficient energy to cross over the dynamical bottleneck, i.e. the maximum of the free energy curve along the reaction coordinate within the CVT formalism, it will not become trapped in this shallow well. Thus, such a complex would not play a significant role in the kinetics of this reaction. In this case, we can reasonably assume that the dynamical bottleneck is in the vicinity of the transition state, i.e. in the CVT calculations the reaction path information for formation of the complex is not included. It is interesting to note that the QCISD level of theory with the smaller cc-pVDZ basis set predicts very similar trend to that of the CCSD(T)/aug-cc-pVTZ level. BH&HLYP/ cc-pVDZ, on the other hand, yields a classical barrier height noticeably lower by 0.72 kcal/mol compared to the CCSD(T) value. 3.1. Rate constant The calculated rate constants using the TST and CVT methods in the temperature range of 300–3000 K are given in Table 2. The ratios between the rate constants kCVT/kTST illustrate the magnitude of the re-crossing effect. It can be seen that the ratio kCVT/kTST increases as the temperature increases; thus the re-crossing effect is less important at high temperature range. At room temperature, it accounts for 55% of the total rate. However at 1000 K, the recrossing effect only lowers the rate constants by 19%.
Table 2 Calculated TST and CVT rate constants (cm3 molecule1 s1). CVT=SCT=HR
T (K)
kTST
kCVT/kTST
HR factora
SCT
k
300 400 500 600 700 800 900 1000 1500 2000 2500 3000
5.01E40 5.64E33 1.12E28 9.20E26 1.21E23 4.96E22 9.38E21 1.02E19 1.77E16 9.54E15 1.21E13 7.19E13
0.45 0.56 0.64 0.68 0.72 0.74 0.76 0.77 0.80 0.81 0.81 0.81
1.07 1.10 1.13 1.15 1.16 1.17 1.18 1.18 1.15 1.10 1.05 1.01
2.25 1.59 1.33 1.19 1.12 1.07 1.04 1.02 0.98 0.97 0.97 0.97
5.45E40 5.54E33 1.06E28 8.63E26 1.12E23 4.62E22 8.72E21 9.48E20 1.60E16 8.26E15 9.99E14 5.68E13
a
kliterature
2.34E23b 7.76E22b
2.21E16c 1.95E14c
Correction factor to the rate constant due to the hindered rotation treatment of the HCO internal rotation motion. b Suggested values from Ref. [10]. c Suggested values from Ref. [7].
0.0
1.0
2.0 1000/T (1/K)
3.0
Fig. 3. Arrhenius plot of the CVT/SCT/HR forward rate constants and available data from the literature.
Due to the relatively flat shape of the VMEP curve in the transition state region and the barrier being quite low (4.68 kcal/mol in the reverse direction), we can expect the tunneling contribution to be small despite there is a large motion of hydrogen atom along the reaction coordinate. This is confirmed by the small curvature tunneling (SCT) calculation. The effect of tunneling is the most at room temperature where it contributes a factor of only 2.25 to the total rate. At higher temperature, e.g. larger than 1000 K, the SCT tunneling factor is negligible. As mentioned earlier internal rotational motion of the HCO group is treated as a hindered rotation (HR). It can be seen that when the temperature increases, the hindered rotation factor, i.e. the correction factor to the rate constants, increases until it reaches its maximum value at 950 K and then decreases. However, its contribution to the rate constants is small, within 20% for the whole temperature range. The Arrhenius expression for the rate constants calculated at the CVT method, including the hinder rotation treatment and the quantum effects, in the temperature range of 300–3000 K given in as:
kðtÞ ¼ 2:95 1021 T 3:152 exp
18 352 1 ; cm3 molecule s1 T ð1Þ
Compared to the latest values of the rate constants available in the literature (also given in Table 2), our results give lower values. The CVT/SCT/HR rate constants and available rate constants in the literature are plotted in Fig. 3. At the low temperature range, this study gives a lower values of rate constants when compared to the available data whereas in the high temperature range (>1500 K) it gives values between the lower values reported by Tsang and Hampson [15] and higher values introduced by Micheal et al. [7]. It can be seen that the CVT/SCT/HR rate constants are within a factor of two when compared to the available data over the wide range of temperature 300–3000 K. Particularly, the deviation of the calculated rate constants are less than 40% when compared to the Baulch et al. data [14] from 700–1000 K, less than 55% compared to the direct measured kinetic data by Baldwin et al. [10] from 713– 813 K; and less than 85% of the Michael et al. data [7] from 1600 to 2100 K. This is noted that the data suggested by Baulch et al. are the only reported data having the uncertainty factor of 3.16. 4. Conclusion In this Letter, the kinetics of the hydrogen abstraction reaction H2CO + O2 ? HCO + HO2 is studied based on the variational
84
L.KK. Huynh et al. / Chemical Physics Letters 469 (2009) 81–84
canonical transition state theory (CVT). The potential energy information was calculated from a sufficiently accurate level of theory. In particular, structural and frequency information along the reaction coordinate were calculated at the BH&HLYP/cc-pVDZ level, which was shown to have a similar level of accuracy as the QCISD level for these properties. Potential energy along the reaction coordinate was calculated at the CCSD(T)/aug-cc-pVDZ//BH&HLYP/ccpVDZ level. We found that the variational/re-crossing effect is important to this reaction, especially at low temperature range where the rate constants are lowered by a factor of up to 2.2; while the correction factor for treating the internal rotation of the HCO group as a hindered rotor is largest at 950 K but less than 20% for the whole temperature range. The calculated rate constants with the corrections agree well with those scattered available data in the literature. Therefore, the rate constant expression can be used confidently for the whole temperature range 300–3000 K. References [1] J.R. Barker, A.A. Herstrom, D.T. Tingey, Water Air Soil Pollut. 86 (1996) 71. [2] O. Hernandez et al., J. Hazard. Mater. 39 (1994) 161. [3] J. Warnatz, Rate Coefficients in the C/H/O System, Springer-Verlag, New York, 1984, p. 197. [4] C.K. Westbrook, F.L. Dryer, Prog. Energy Combust. Sci. 10 (1984) 1. [5] T.S. Norton, F.L. Dryer, Combust. Sci. Technol. 63 (1989) 107. [6] J. Peeters, G. Mahnen, Proc. Combust. Inst. 14 (1973) 133.
[7] N.K. Srinivasan, M.-C. Su, J.W. Sutherland, J.V. Michael, J. Phys. Chem. A 109 (2005) 7902. [8] G.P. Smith et al.,
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Theoretical study on kinetics of the H2CO + O2 ? HCO + HO2 reaction Lam K. Huynh a,1, Max Tirtowidjojo b, Thanh N. Truong a,* a b
Department of Chemistry, University of Utah, 315 S 1400 E, Rm. 2020, Salt Lake City, UT 84112, United States Dow Chemical Company, 2301 N Brazosport Blvd., B-1226, Freeport, TX 77451, United States
a r t i c l e
i n f o
Article history: Received 21 October 2008 In final form 9 December 2008 Available online 24 December 2008
a b s t r a c t A theoretical study on the kinetics of the H2CO + O2 ? HCO + HO2 reaction has been carried out from first principles using canonical variational transition state theory incorporating corrections from small curvature tunneling and hindered rotation treatments. Rate constants were calculated using the geometries, gradient and Hessian obtained at the BH&HLYP/cc-pVDZ level with the energy correction at the CCSD(T)/aug-cc-pVTZ//BH&HLYP/cc-pVDZ level. The rate constant expression is obtained as kðtÞ ¼ 2:95 1021 T 3:152 exp 18 T352 (cm3 molecule1 s1) for the temperature range 300–3000 K. The result shows an excellent agreement with available scattered experimental data over this temperature range and thus the expression can be used confidently for the whole range. Ó 2008 Published by Elsevier B.V.
1. Introduction Formaldehyde has recently generated serious environmental concerns because of its large accumulation in the atmosphere due to its thermodynamic stability [1,2]. It has also long been recognized as an important kinetic intermediate in the combustion of many hydrocarbon fuels, particularly of methane and methanol [3–6]. Recently, Michael et al. [7] has shown that an oxidation reaction, namely, H2CO + O2 ? HCO + HO2 is an essential component for the kinetic modeling of hydrocarbon combustion. As a result, it is included in various kinetic combustion models [8,9]. Unfortunately, there have been only seven recorded data of thermal rate constants for this reaction. Only one record from direct kinetic measurement data [10] is available on this elementary oxidation reaction H2CO + O2 ? HCO + HO2, which is in the small temperature range of 713–816 K and at a pressure of 0.07–0.08 bar and dates back from 1974; others are from theoretical estimates: three derived indirectly from complex mechanisms fit to experimental observations from the combustion of formaldehyde [11–13] and one from the CH3 + O2 reaction [7], two from the extensive literature review [14,15], one from a preliminary transition state theory (TST) calculation [13]. In the present study, we performed accurate theoretical study on rate constants for the H2CO + O2 ? HCO + HO2 reaction. Particularly, the rate constants are calculated using the CVT [16–19] with the inclusion of corrections of tunneling and hindered rotation treatments and the potential energy surface was corrected at the CCSD(T)/aug-cc-PVTZ level. The calculated rate constants are then
* Corresponding author. Fax: +1 801 581 4353. E-mail address: [email protected] (T.N. Truong). 1 Current address: Chemical Engineering Department, Colorado School of Mines, Golden CO 80401, Unites States. 0009-2614/$ - see front matter Ó 2008 Published by Elsevier B.V. doi:10.1016/j.cplett.2008.12.050
validated by comparison to those available in the literature so that rate constant expression can be used to obtain rate constants in the temperature range where there is no reliable rate constant available. 2. Electronic structure calculation All electronic structure calculations were carried out using the program package GAUSSIAN 03 [20]. It has been shown earlier that the combination of Becke’s Half-and-Half (BH&H) non-local exchange with Lee–Yang–Parr (LYP) correlation functionals can be used cost-effectively to predict the properties, e.g. geometries and frequencies of reactants, products, as well as transition states for hydrogen abstraction reactions by a radical [21–26]. The Dunning’s correlation-consistent polarized valence double-zeta basis set [3s2p1d/2s1p] denoted as cc-pVDZ [27] was employed with the DFT method to determine structural and vibrational frequencies for the stationary points. The quadratic configuration interaction method with inclusion of singles and doubles excitations (QCISD) using the same cc-pVDZ basis set is also used as a reference source for the DFT method in finding stationary structures as well as calculating vibrational frequencies. The minimum energy path (MEP) connecting the reactants and products is also obtained at the same BH&HLYP/cc-pVDZ level by tracing the steepest descent path from the transition state to the corresponding reactants and products using the Gonzalez–Schlegel method [28] with a step size of 0.1 (amu)1/2 bohr for a total of 200 steps. This step size is found to be sufficient for the convergence of the variational rate constant as mentioned in our previous study. Force constants at 22 selected points (12 points in the reactant channel and 10 points in the product channel) along the MEP were carried out to obtain the necessary potential energy surface information for CVT calculations. These points were found using our
82
L.KK. Huynh et al. / Chemical Physics Letters 469 (2009) 81–84
autofocusing technique [29], which uses a combination of the barrier shape, and the second derivatives of the Z-matrix geometrical parameters to predict the regions along the MEP that are sensitive in rate determination. In order to improve the energy information along the MEP, single point calculation for these points at a higher level of theory such as [CCSD(T)] is needed, in particular the CCSD(T) with aug-cc-pVTZ basis set was used. 3. Results and discussion Potential energy surface BH&HLYP optimized geometries of the reactants (H2CO and O2), products (HCO and HO2), the transition state as well as the van-der-Waal complex in the product channel are shown in Fig. 1. The experimental data available in the litera-
1.108 1.119 (1.116) 1.182 1.208 (1.208)
115.4° 115.2° (116.5°)
1.191 1.212 (1.208)
REACTANTS
1.306 1.342
105.6° 103.9°
124.6° 124.4°
1.126 1.137
1.166 1.188
0.967 0.977
PRODUCTS
90.7° 92.0°
144.5° 145.3°
2.064 2.172
0.979 0.984
COMPLEX
ture are also given in parentheses. Compared to the experimental data [30], the BH&HLYP/cc-pVDZ yields a shorter bond length for oxygen and formaldehyde molecules with the largest difference of 0.026 Å for O2 (0.008 Å and 0.017 Å for C–H and C@O bonds in H2CO, respectively). The transition state has the product-like structure with the 1C–4H bond distance of 1.440 Å, which is much larger than the C–H bond length of 1.108 Å in H2CO and the 4H–5O bond distance is very close to the O-H bond length of the OOH product (1.137 Å compared to 0.967 Å). The geometry of the complex is very close to that of the transition state except the 1C–4H bond length is longer and the 5O–4H is shorter. In addition, the angle (3H–1C–4H) decreases by around 10° when the transition state becomes the complex in the exit channel. BH&HLYP harmonic vibrational frequencies for the reactants, products, transitions state and complex are listed in Table 1. It can be seen that the BH&HLYP frequencies are very close to the available experimental data. For formaldehyde, the average difference is 136 cm1. The difference mainly comes from the symmetric and asymmetric H–O stretching. For oxygen molecule, the predicted stretching frequency is 238 cm1 larger than the experimental value from Ref. [30]. It is noted that the scaling factor of 0.97 was applied for these frequencies; this will reduce the gap between the calculated and experimental data. By comparison to the QCISD results, which are also given in Figs. 1 and 2, BH&HLYP predicts similar geometries of the stationary structures. In particular, for the reactants and products, the geometries at the two levels are almost the same with the maximum bond length difference of 0.035 Å for the O–O bond in OOH radical. The bond angles between the two levels are almost identical. For the transition state, the geometries at the two levels are still very close to each other except for the bond length between the two oxygen atoms for which the BH&HLYP bond length is 0.03 Å shorter than that using the QCISD level of theory. Such a difference was also found for the OOH radical. Since the O–O bond is only a spectator in the H2CO + O2 reaction, we can expect this difference to have only a minor effect on the calculated rate constants. Regarding frequencies, as shown in Table 1, BH&HLYP yields almost identical results for reactants, products and even for transition state compared to those from QCISD calculations. Specifically, BH&HLYP/cc-pVDZ tends to give slightly higher frequencies. The largest differences are found for the O–O stretches in the oxygen molecule and in the transition state, which are around 150 cm1. However, the contribution of the differences to the rate constants, especially the vibrational partition function, is small due to the large magnitude of these frequencies. Thus, it can be concluded that the BH&HLYP is a cost-effective method for structural and frequency information for this reaction. The schematic energy profile for the title reaction is shown in Fig. 2. The most accurate energetic results from this study are those
Table 1 Harmonic vibrational frequencies (in cm1) for the transition state (TS), reactants, products, and complex.
1.266 1.299
105.1° 103.9°
1.137 1.155
153.7° 152.1°
1.123 1.135
101.4° 100.2°
O2 H2CO
1.440 1.431
1.165 1.186
HCO OOH TS
TRANSITION STATE Complex Fig. 1. Optimized geometries (distances in Å and angles in degrees) of the reactants, products, complex and transition state at the BH&HLYP/cc-pVDZ and QCISD/ccpVDZ levels (italic numbers). The numbers in the parentheses are the experimental values from Ref. [30].
a
BHLYP/cc-pVDZ
QCISD/cc-pVDZ
Experimentala
1818 1247, 1302, 1570, 1923, 2999, 3068 1140, 2043, 2745 1250, 1502, 3781 2033i, 107, 124, 319, 429, 613, 1108, 1157, 1406, 1673, 2045, 2832 80, 109, 140, 216, 266, 574, 1095, 1272, 1560, 2057, 2807, 3509
1645 1193, 1275, 1544, 1816, 2944, 3004 1135, 1894, 2679 1086, 1435, 3675 2286i, 112, 149, 323, 419, 595, 1088, 1154, 1226, 1655, 1919, 2769 67, 101, 125, 184, 230, 492, 1095, 1102, 1486, 1922, 2747, 3531
1580 1167, 1249, 1500, 1746, 2783, 2843 1081, 1868, 2485 1098, 1392, 3436
From Ref. [30].
N/A
N/A
83
L.KK. Huynh et al. / Chemical Physics Letters 469 (2009) 81–84
TRANSITION STATE
5.60 4.75 4.68
-10 HCO+HO2
Baulch et al. (1992)
COMPLEX 40.47 BH&HLYP/cc-pVDZ 44.03 QCISD/cc-pVDZ 40.83 CCSD(T)/aug-cc-pVTZ
41.10 42.05 40.19
H2CO+O2 Fig. 2. Schematic energy profile of the reaction H2CO + O2, at the different levels of theory as indicated. Numbers are in kcal/mol.
Log{k(T)/cm3molecule-1s-1}
-15
Tsang et al. (1986) Baldwin et al. (1974) Michael et al. (2005)
-20
This work: CVT/SCT/HR
-25
-30
-35
-40
from single-point energy CCSD(T)/aug-cc-pVTZ calculations at the BH&HLYP/cc-pVDZ optimized geometries. The H2CO + O2 reaction is endothermic with the reaction energy of 40.19 kcal/mol. The calculated classical barrier for this reaction is 40.83 kcal/mol which is 0.64 kcal/mol above the exit channel. The weak binding complex in the exit channel is only 4.68 kcal/mol below the dissociative products. Since the reaction would occur at rather high temperatures due to its high barrier, one can expect that if the trajectory has sufficient energy to cross over the dynamical bottleneck, i.e. the maximum of the free energy curve along the reaction coordinate within the CVT formalism, it will not become trapped in this shallow well. Thus, such a complex would not play a significant role in the kinetics of this reaction. In this case, we can reasonably assume that the dynamical bottleneck is in the vicinity of the transition state, i.e. in the CVT calculations the reaction path information for formation of the complex is not included. It is interesting to note that the QCISD level of theory with the smaller cc-pVDZ basis set predicts very similar trend to that of the CCSD(T)/aug-cc-pVTZ level. BH&HLYP/ cc-pVDZ, on the other hand, yields a classical barrier height noticeably lower by 0.72 kcal/mol compared to the CCSD(T) value. 3.1. Rate constant The calculated rate constants using the TST and CVT methods in the temperature range of 300–3000 K are given in Table 2. The ratios between the rate constants kCVT/kTST illustrate the magnitude of the re-crossing effect. It can be seen that the ratio kCVT/kTST increases as the temperature increases; thus the re-crossing effect is less important at high temperature range. At room temperature, it accounts for 55% of the total rate. However at 1000 K, the recrossing effect only lowers the rate constants by 19%.
Table 2 Calculated TST and CVT rate constants (cm3 molecule1 s1). CVT=SCT=HR
T (K)
kTST
kCVT/kTST
HR factora
SCT
k
300 400 500 600 700 800 900 1000 1500 2000 2500 3000
5.01E40 5.64E33 1.12E28 9.20E26 1.21E23 4.96E22 9.38E21 1.02E19 1.77E16 9.54E15 1.21E13 7.19E13
0.45 0.56 0.64 0.68 0.72 0.74 0.76 0.77 0.80 0.81 0.81 0.81
1.07 1.10 1.13 1.15 1.16 1.17 1.18 1.18 1.15 1.10 1.05 1.01
2.25 1.59 1.33 1.19 1.12 1.07 1.04 1.02 0.98 0.97 0.97 0.97
5.45E40 5.54E33 1.06E28 8.63E26 1.12E23 4.62E22 8.72E21 9.48E20 1.60E16 8.26E15 9.99E14 5.68E13
a
kliterature
2.34E23b 7.76E22b
2.21E16c 1.95E14c
Correction factor to the rate constant due to the hindered rotation treatment of the HCO internal rotation motion. b Suggested values from Ref. [10]. c Suggested values from Ref. [7].
0.0
1.0
2.0 1000/T (1/K)
3.0
Fig. 3. Arrhenius plot of the CVT/SCT/HR forward rate constants and available data from the literature.
Due to the relatively flat shape of the VMEP curve in the transition state region and the barrier being quite low (4.68 kcal/mol in the reverse direction), we can expect the tunneling contribution to be small despite there is a large motion of hydrogen atom along the reaction coordinate. This is confirmed by the small curvature tunneling (SCT) calculation. The effect of tunneling is the most at room temperature where it contributes a factor of only 2.25 to the total rate. At higher temperature, e.g. larger than 1000 K, the SCT tunneling factor is negligible. As mentioned earlier internal rotational motion of the HCO group is treated as a hindered rotation (HR). It can be seen that when the temperature increases, the hindered rotation factor, i.e. the correction factor to the rate constants, increases until it reaches its maximum value at 950 K and then decreases. However, its contribution to the rate constants is small, within 20% for the whole temperature range. The Arrhenius expression for the rate constants calculated at the CVT method, including the hinder rotation treatment and the quantum effects, in the temperature range of 300–3000 K given in as:
kðtÞ ¼ 2:95 1021 T 3:152 exp
18 352 1 ; cm3 molecule s1 T ð1Þ
Compared to the latest values of the rate constants available in the literature (also given in Table 2), our results give lower values. The CVT/SCT/HR rate constants and available rate constants in the literature are plotted in Fig. 3. At the low temperature range, this study gives a lower values of rate constants when compared to the available data whereas in the high temperature range (>1500 K) it gives values between the lower values reported by Tsang and Hampson [15] and higher values introduced by Micheal et al. [7]. It can be seen that the CVT/SCT/HR rate constants are within a factor of two when compared to the available data over the wide range of temperature 300–3000 K. Particularly, the deviation of the calculated rate constants are less than 40% when compared to the Baulch et al. data [14] from 700–1000 K, less than 55% compared to the direct measured kinetic data by Baldwin et al. [10] from 713– 813 K; and less than 85% of the Michael et al. data [7] from 1600 to 2100 K. This is noted that the data suggested by Baulch et al. are the only reported data having the uncertainty factor of 3.16. 4. Conclusion In this Letter, the kinetics of the hydrogen abstraction reaction H2CO + O2 ? HCO + HO2 is studied based on the variational
84
L.KK. Huynh et al. / Chemical Physics Letters 469 (2009) 81–84
canonical transition state theory (CVT). The potential energy information was calculated from a sufficiently accurate level of theory. In particular, structural and frequency information along the reaction coordinate were calculated at the BH&HLYP/cc-pVDZ level, which was shown to have a similar level of accuracy as the QCISD level for these properties. Potential energy along the reaction coordinate was calculated at the CCSD(T)/aug-cc-pVDZ//BH&HLYP/ccpVDZ level. We found that the variational/re-crossing effect is important to this reaction, especially at low temperature range where the rate constants are lowered by a factor of up to 2.2; while the correction factor for treating the internal rotation of the HCO group as a hindered rotor is largest at 950 K but less than 20% for the whole temperature range. The calculated rate constants with the corrections agree well with those scattered available data in the literature. Therefore, the rate constant expression can be used confidently for the whole temperature range 300–3000 K. References [1] J.R. Barker, A.A. Herstrom, D.T. Tingey, Water Air Soil Pollut. 86 (1996) 71. [2] O. Hernandez et al., J. Hazard. Mater. 39 (1994) 161. [3] J. Warnatz, Rate Coefficients in the C/H/O System, Springer-Verlag, New York, 1984, p. 197. [4] C.K. Westbrook, F.L. Dryer, Prog. Energy Combust. Sci. 10 (1984) 1. [5] T.S. Norton, F.L. Dryer, Combust. Sci. Technol. 63 (1989) 107. [6] J. Peeters, G. Mahnen, Proc. Combust. Inst. 14 (1973) 133.
[7] N.K. Srinivasan, M.-C. Su, J.W. Sutherland, J.V. Michael, J. Phys. Chem. A 109 (2005) 7902. [8] G.P. Smith et al.,