Biofuel Combustion. Energetics and Kinetics of Hydrogen Abstraction

Article pubs.acs.org/JPCA

Biofuel Combustion. Energetics and Kinetics of Hydrogen Abstraction from Carbon‑1 in n‑Butanol by the Hydroperoxyl Radical Calculated by Coupled Cluster and Density Functional Theories and Multistructural Variational Transition-State Theory with Multidimensional Tunneling I. M. Alecu,† Jingjing Zheng, Ewa Papajak, Tao Yu, and Donald G. Truhlar* Department of Chemistry and Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455-0431, United States S Supporting Information *

ABSTRACT: Multistructural canonical variational transition-state theory with small-curvature multidimensional tunneling (MS-CVT/SCT) is employed to calculate thermal rate constants for hydrogen-atom abstraction from carbon-1 of n-butanol by the hydroperoxyl radical over the temperature range 250−2000 K. The M08-SO hybrid meta-GGA density functional was validated against CCSD(T)F12a explicitly correlated wave function calculations with the jul-cc-pVTZ basis set. It was then used to compute the properties of all stationary points and the energies and Hessians of a few nonstationary points along the reaction path, which were then used to generate a potential energy surface by the multiconfiguration Shepard interpolation (MCSI) method. The internal rotations in the transition state for this reaction (like those in the reactant alcohol) are strongly coupled to each other and generate multiple stable conformations, which make important contributions to the partition functions. It is shown that neglecting to account for the multiple-structure effects and torsional potential anharmonicity effects that arise from the torsional modes would lead to order-of-magnitude errors in the calculated rate constants at temperatures of interest in combustion.

1. INTRODUCTION Over the past several years, n-butanol has received considerable attention as a potential alternative fuel source. There are reasons why n-butanol is a more optimal choice than the current most widely used liquid fuel blending component, ethanol; these include its higher energy density, decreased NOx emissions, lower vapor pressure, better blending properties, higher tolerance to water contamination, and nontoxicity to humans.1−3 Furthermore, it was also recently shown that, in a direct-injection spark-ignition engine, n-butanol can match the performance of ethanol while less of it is consumed.4 In light of such favorable properties, understanding the complex combustion chemistry of n-butanol under a wide range of conditions of practical interest has become a central challenge, and the development of an extensively validated kinetic model for nbutanol combustion is an active area of research.5−16 Most of the kinetic mechanisms that have been proposed for modeling n-butanol combustion can do so adequately only for high temperature conditions where the combustion behavior is governed by reactions involving the abstraction of H atoms from the fuel by small radicals such as H, OH, and CH3, the thermochemistry and kinetics of which are well characterized for many fuels. The chemistry of peroxyl radicals, which is considerably less well-established, is thought to drive the © 2012 American Chemical Society

combustion processes occurring in low- and intermediatetemperature alkane oxidation17 and is also expected to dominate under the same conditions in the combustion of oxygenated hydrocarbons such as n-butanol. Due to the complexity of peroxyl chemistry, to the best of our knowledge, only one attempt to model the combustion of n-butanol at the lower end of the practical temperature range has thus far been published.15 Sensitivity analyses performed during the available modeling studies5−16 have shown that the H-atom abstraction reactions from the various positions in n-butanol by the hydroperoxyl radical play a key role in the combustion of n-butanol at intermediate temperatures. Among these, the abstractions from the OH and carbon-2 sites of n-butanol are the slowest and contribute only marginally to the overall rate of reaction between n-butanol and HO2; abstractions from the carbon-3 and carbon-4 sites are more favorable and can each appreciably contribute to the overall rate of abstraction by HO2, but the process that is by far the most important at the temperatures of interest is the abstraction reaction from the carbon-1 position: Received: August 26, 2012 Revised: November 14, 2012 Published: November 14, 2012 12206

dx.doi.org/10.1021/jp308460y | J. Phys. Chem. A 2012, 116, 12206−12213

The Journal of Physical Chemistry A

Article

CH3(CH 2)3 OH + HOO → CH3(CH 2)2 CHOH + HOOH

(R1)

Thus, to unravel the rich combustion chemistry of n-butanol in this temperature regime, an accurate determination of the kinetics of this reaction system is critical. Unfortunately, such reactions are challenging to study experimentallyeven if the complicated issues of HO 2 generation and detection are overcome, these reactions are still too slow to measure much less isolate from secondary chemistry at low temperatures, whereas at the elevated temperatures where the rate constant would become appreciable, the HO2 radical and HOOH are no longer thermally stable. In light of these experimental limitations, dependable estimates for the kinetics for the H-atom abstractions by HO2 mustat least at presentcome from theory, preferably from reliable, well validated theory. Recently, Zhou et al.18 employed conventional transition-state theory with a one-dimensional asymmetric Eckart model for tunneling and a separable treatment of hindered rotations (1D-HR) to calculate rate constants for H-atom abstraction from all sites on n-butanol by HO2. The rate constants reported in the study of Zhou et al.18 were substantially different than the ones previously calculated by Vranckx et al.13 using transition-state theory with zero-curvature tunneling. Zhou et al. attributed the discrepancies between their calculations and those of Vranckx et al. to the different treatments of tunneling and, more importantly, to the different approximations used to treat the effects of hindered-rotationswith the latter study13 approximating the anharmonic effects of such vibrations by using a perturbative treatment to calculate effective anharmonic frequencies.19 Using the effective frequencies of an anharmonically perturbed oscillator to compute the contributions of torsional modes to the partition function can lead to significant errors because it does not account for the torsional potential anharmonicity and multiple structures generated by largeamplitude motion in torsions;20−23 thus, Zhou et al.18 suggested that a better treatment of torsional anharmonicity might be needed to reliably characterize these reaction systems. Seal et al.24 have used multistructural canonical variational transition-state theory with multidimensional small-curvature tunneling (MS-CVT/SCT)25 to investigate the H-atom abstraction by HO2 from the carbon-3 site on n-butanol. The MS-CVT/SCT formalism marks an improvement over the dynamical method used by Zhou et al.18 and Vranckx et al.13 because it includes a more physical multidimensional model for quantum mechanical tunneling and a more appropriate treatment of torsional anharmonicity, which captures the effects arising from the coupling of torsions to each other and to external rotation. The MS-CVT/SCT rate constants calculated by Seal et al. exhibit an appreciably different temperature dependence than those calculated in the previous two studies, resulting in substantial discrepancies at temperatures of interest: at 500 K, the rate constant from Zhou et al. is over a factor of 13 larger than that from Seal et al., while at 1400 K, the rate constant calculated by Vranckx et al. is about a factor of 4 smaller than the MS-CVT/SCT result (Figure 1). Such alarmingly large differences in the rate constants calculated for the carbon-3 abstraction raise concerns about whether or not the available literature estimates of the rate constant for the most important of these H-atom abstractions by HO2 from n-butanol, reaction R1, are reliable. In the present

Figure 1. Comparison of literature computed rate constants for abstraction from carbon-3: Seal et al. (solid black line), Zhou et al. (dashed red line), and Vranckx et al. (dotted blue line).

article, we address this issue by computing accurate rate constants for (R1) over a wide range of temperatures using the MS-CVT/SCT formalism.

2. THEORY The thermal rate constants for (R1) are calculated using multistructural canonical variational transition-state theory with small-curvature multidimensional tunneling (MS-CVT/ SCT).25 This formulation of CVT is appropriate for the calculation of rate constants for complex reactions involving flexible species, i.e., reactants and/or transition states that exist in multiple conformations capable of interconverting among themselves by torsions. Standard transition-state treatments usually assume a model whereby the partition functions employed are calculated from the properties of just the lowest-energy conformation of the reactant or the transition state via the rigid-rotor harmonic (or anharmonic) oscillator model. Such treatments can generally be denoted as singlestructural (SS), and clearly, an SS-type of CVT (SS-CVT) is inadequate for quantitatively describing the complex kinetics of a reaction system composed of molecules with chain structures such as that under present investigation. To better capture the dynamical nature of flexible reaction systems, the SS partition functions entering CVT are substituted by more physically meaningful conformational−rotational−vibrational partition functions, denoted as MS-T26 partition functions, that account for torsional anharmonicity without assuming mode separability. The MS-CVT/SCT formalism for bimolecular reactions has been derived and discussed in detail in recent studies;24,27,28 therefore, we only briefly summarize the central concepts and equations here. It will be convenient to express and calculate the MS-CVT/SCT rate constant in terms of its more familiar single-structural analog kSS‑CVT/SCT: k MS ‐ CVT/SCT = F MS ‐ TkSS ‐ CVT/SCT

(1) MS‑T

In eq 1, the multistructural factor F adjusts the singlestructural CVT/SCT29−31 result to account for the anharmonic effects brought about by the existence of torsions in the transition state and/or reactants, the coupling between such modes, and the coupling of these modes to external rotation, by 12207



Article

substituting the corresponding MS-T26 partition functions (which capture these effects) for the single-structure rigid-rotor harmonic-oscillator (SS-RRHO) partition functions normally used in SS-CVT calculations: F

MS ‐ T

=

2

MS ‐ T FTS 2

(R1). These searches were conducted using the M08-SO34 hybrid meta-GGA density functional approximation and the MG3S35,36 basis set, and the same electronic model chemistry was then used to calculate the energies, gradients, and Hessians for all of the torsional conformers located for the TS. The optimized geometries and Hessians for all torsional conformers located for the TS in this study, and those previously located at the same level of theory for n-butanol,37 were then used to obtain the MS-T partition functions used in the MS-CVT/SCT rate constant calculation. All vibrational frequencies were scaled by 0.983, a previously determined38 empirical scale factor for vibrational frequencies obtained with M08-SO/MG3S. An integration grid consisting of 99 radial shells around each atom and 974 angular points in each shell was employed in all M08-SO/MG3S calculations to ensure that all calculated properties were converged. The MSTor program39,40 was used to calculate the MS-LQ and MS-T partition functions. The performance of the M08-SO/MG3S electronic model chemistry was validated against high-level single-point calculations employing the CCSD(T)-F12a41−43 explicitly correlated method with the jul-cc-pVTZ44 basis setcarried out on the M08-SO/MG3S geometries obtained for the lowest-energy conformers of all species along the reaction coordinate for (R1). The energetics obtained from this analysis are shown in Table 1, and they demonstrate that the M08-SO/MG3S

∏i = 1 FR,MSi ‐ T

=

MS ‐ T ∏i = 1 Q R,SSi ‐ RRHO Q con ‐ rovib,TS 2 SS ‐ RRHO MS ‐ T QTS ∏i = 1 Q con ‐ rovib,R, i

(2)

It is informative to express each anharmonicity correction factor, FMS‑T , as the product of a multistructural local-harmonic α component and a torsional potential component, F‑MS‑LH and α FTα , respectively;25,27,28 the former of these provides a measure of how much of the total anharmonicity correction can be attributed to the fact that, as the temperature is increased, the partition function receives increasingly more sizable contributions from more than just the one torsional conformer (i.e., structure) used in SS-RRHO treatments, whereas the latter factor captures the extent of the anharmonic correction stemming from the fact that the potentials for torsional motions are not quadratic. Then, with the subscript α denoting either the TS or one of the reactants, the above consideration yields FαMS ‐ T

=

FαMS ‐ LQ FαT

⎞ ⎛ Q MS ‐ LQ ⎞⎛ Q MS ‐ T con ‐ rovib, α con ‐ rovib, α = ⎜⎜ SS ‐ RRHO ⎟⎟⎜⎜ MS ‐ LQ ⎟⎟ ⎠⎝ Q con ‐ rovib, α ⎠ ⎝ Qα

(3)

QMS‑LQ con‑rovib,α

Table 1. Forward and Reverse Barrier Heights and Reaction Energy (kcal/mol) for (R1)

where is the multistructural local-quasiharmonic (LQ) partition function obtained by keeping all the structures but treating each of them by local-harmonic formulas, as described in detail elsewhere.25−28 Note that LQ denotes using the harmonic oscillator formulas with frequencies based on a local normal-mode analysis, but because the frequencies are scaled29 to account for anharmonicity (as well as for systematic errors in the electronic structure model), the treatment is actually quasiharmonic. An important aspect of the MS-T formalism is that the result tends to a local harmonic approximation in the low-temperature limit and to the correct free-rotor limit in the high-temperature limit, with the transition through the hindered-rotor regime governed by effective local periodicities.26 The small-curvature tunneling30 (SCT) approximation is used to calculate multidimensional tunneling ground-state transmission coefficients for (R1) in this study. Following the ground-state approximation29,32 used in the original MS-CVT/ SCT paper,25 the reaction path used for the tunneling calculations is taken to be the minimum-energy path (MEP) through the saddle point with the lowest value of VGa (where VGa is the sum of the potential energy along the MEP and the local zero point vibrational energy of modes transverse to it); i.e., we employ the path passing through the transition structure with the lowest zero-point-inclusive energy. Ground-state transmission coefficients using the zerocurvature tunneling32,33 (ZCT) approximation were also computed to enable a more direct comparison with the calculations of Vranckx et al.13 and also to examine whether reaction-path curvature effects, which can be very large due to the phenomenon of corner cutting, are important in the case of (R1).

electronic model chemistry

V‡f

V‡r

ΔU

13.94

5.11

8.83

13.96

5.50

8.46

12.25

4.00

8.25

12.27

4.39

7.88

MUE

a

classical energetics: CCSD(T)-F12a/jul-cc-pVTZ// M08-SO/MG3S M08-SO/MG3S 0 K energetics:c CCSD(T)-F12a/jul-cc-pVTZ// M08-SO/MG3S M08-SO/MG3S a

0.26b

0.26b

b

ZPE exclusive. Mean unsigned errors in the three quantities with respect to CCSD(T)-F12a/jul-cc-pVTZ//M08-SO/MG3S values. c ZPE inclusive based on scaled M08-SO/MG3S frequencies.

electronic model chemistry is well-suited to characterizing the energetics of (R1), with a mean unsigned error with respect to that calculated by CCSD(T)-F12a/jul-cc-pVTZ of just 0.26 kcal/mol. We further note that the M08-SO/MG3S electronic model chemistry was recently used to explore the barrier heights and reaction energies for the analogous methanol + HO2 reaction system, where it was also found to give results in good agreement with experiment and high-level wave function theory (WFT) calculations.27 Two other studies45,46 also demonstrated that the M08-SO/MG3S electronic model chemistry can give accurate transition structures and transition-state energies. The minimum-energy path, the potential along this path (VMEP), and the vibrationally adiabatic ground-state potential (VGa ) used to obtain the SS-CVT/SCT rate constants for (R1) in this work were obtained using the multiconfiguration Shepard interpolation method47−50 (MCSI). The accuracy and efficiency that can be attained with MCSI make it a very appealing method for investigating the dynamics of large systems. This approach involves the construction of a semiglobal analytic multidimensional potential energy surface

3. CALCULATIONS The first step was carrying out structure searches along the degrees of freedom corresponding to torsions of the TS for 12208



Article

Table 2. Torsional Angles (deg), Local Periodicities (M), and Relative Zero-Point-Inclusive and Zero-Point-Exclusive Energies (kcal/mol) of 14 Transition Structures for (R1)a torsional anglesa

relative energies

structure (j)

τ1(8−7−6−5)

τ2(7−6−5−1)

τ3(6−5−1−2)

τ4(5−1−2−4)

τ5(2−1−9−12)

τ6(1−9−12−15)

Mjb

ZPE- exclusive

ZPE- inclusive

1 2 3 4 5 6 7 8 9 10 11 12 13 14

−95 −98 −97 98 −98 97 98 −98 −98 −101 98 98 97 −89

−13 −23 −19 22 −22 −1 1 −19 −18 −9 4 3 −1 27

36 54 48 −38 53 −17 −22 49 49 −41 −25 −23 −19 −13

−62 −64 −63 50 −65 53 53 −62 −64 50 54 54 55 −70

−54 −176 172 −44 −173 −58 178 65 58 −48 66 58 −174 −43

−178 180 60 −56 −63 180 59 −178 59 −60 −179 58 −179 −60

2.4912 2.4802 2.7027 2.5678 2.4368 2.5290 2.7816 2.3761 2.6508 3.0180 2.7711 2.6197 2.6618 2.3538

0.00 0.09 0.42 0.45 0.52 0.55 0.58 0.59 0.62 0.74 0.80 0.89 0.91 0.99

90.99 91.12 91.58 91.53 91.62 91.42 91.51 91.61 91.76 92.01 91.71 91.92 91.66 92.02

a

Each of these structures also has a nonsuperposable mirror image. bThese M values refer to torsions 1−6; the seventh torsion is the methyl torsion with Mj = 3 for all structures j.

for the reaction system of interest using a reactive force field trained against a small number of higher-level electronic structure calculations. A modified version of the Page−McIver algorithm51 in which the Hessians are reused at non-Hessian steps of the calculation (also known as the MPM algorithm52) was used to calculate the MEP. The Hessian was calculated at every tenth point, and the step size used between points was 0.001 Å in isoinertial coordinates scaled to a reduced mass of 1 amu. The MCSI method involves the interpolation of the offdiagonal elements of a diabatic potential matrix in which the diagonal elements are obtained from molecular mechanics. In the present case, we employ the MM3 method for the molecular mechanics portions of the calculations, and the M08SO/MG3S electronic model chemistry is used to calculate the energies, gradients, and Hessians of 6 nonstationary Shepard points along the MEP and three stationary points (reactant complex, product complex, and the saddle point). As mentioned above, a scale factor of 0.983 was used to scale the harmonic frequencies obtained with the M08-SO/MG3S electronic model chemistry.38 This factor accounts in part for systematic errors in the electronic structure model chemistry and also for the systematic difference between a harmonic ZPE and real anharmonic one. For the latter reason, even the local harmonic results include some anharmonicity, and all local harmonic results are actually quasiharmonic (they use the harmonic formulas but not the harmonic frequencies). The nonstationary Shepard points were placed at locations along the MEP that were strategically selected following a previously described procedure,48 at values along the reaction coordinate (s) of −0.09, +0.09, −0.21, +0.11, −0.62, and +0.15 Å. The MCSI calculations were carried out using the MCSI53 and MC-TINKERATE54 programs, and M08-SO/MG3S calculations were performed with a locally modified version55 of the Gaussian 09 program suite.56

image for a total of 188 transition structures. Six of the torsions of the transition structures are strongly coupled, and the strong coupling (SC) scheme previously described26 was used with the Voronoi algorithm40 to allocate torsional volumes for assigning an effective local periodicity to these modes in each structure for the calculation of MS-T partition functions. The methyl rotor in the transition state was assumed to be nearly separable (NS), so that the overall scheme used in this study corresponds, in the notation established previously,26 to NS:SC = 1:6. There are 14 pairs of transition structures within 1.0 kcal/ mol (ZPE exclusive) of the lowest-energy saddle point, and the torsional angles and local periodicities of the strongly coupled torsions, as well as the relative ZPE-inclusive and ZPE-exclusive energies of (one structure from each set of) these pairs ZPEinclusive and ZPE-exclusive energies of one each of these pairs of transition structures obtained using the M08-SO/MG3S electronic model chemistry are provided in Table 2. The atomnumbering scheme used for the torsions in Table 2 is illustrated in Figure 2. Information similar to that in Table 2 for the other 80 pairs of transition structures, along with the geometries in Cartesian coordinates of all structures, are given in the

4. RESULTS AND DISCUSSION The properties of the reactant n-butanol calculated using the M08-SO/MG3S electronic model chemistry were given in a previous publication.37 In the present study we located 94 conformers of the transition state, each possessing a nonsuperposable mirror

Figure 2. Lowest-energy transition structure for (R1). 12209



Article

so that the effect on the rate constant is actually an increase of a factor of 29; in particular, the multistructural factors for the reaction are shown in the last column, and they increase monotonically from 1.7 at 200 K to 29 at 2000 K. The VMEP and VGa curves obtained using the MCSI procedure are depicted in Figure 3. There is some noise that could

Supporting Information; all of the structures were used in the MS-T calculations. The classical barrier height in our work is defined as the energy difference between the lowest-energy structure of the saddle points and the lowest-energy structure of the reactant. In our lowest-energy saddle points, the conformations of the nbutanol portion are trans−gauche−gauche (T−G−G− and T+G+G+) although the lowest-energy structures of n-butanol are found to be trans−gauche−trans (T+G−T+ and T−G+T−). Note that T+ and T− stand for “trans” and correspond to (+150° to +180°) and (−150° to −180°), respectively; G+ and G− stand for “gauche” and correspond to (+30° to +75°) and (−30° to −75°), respectively. In the previous computational studies of reaction R1,13,18,57 comprehensive structure searches to characterize the numerous conformations of the TS were not conducted; instead, all of the kinetics calculations in the previous studies employed a transition structure connected to the global minimum structure of n-butanol (T+G−T+ or T−G+T−) on the MEP. The ZPE-inclusive barrier height calculated in the present work is about 1 kcal/mol lower than those in previous studies.13,18,57 This difference could be caused by several differences in the calculations, such as different electronic structure methods, different conformational searches, and different treatments of torsional anharmonicity. The 94 pairs of transition structures found in the present work span an energy range of ∼4 kcal/mol, so using an inappropriate (i.e., nonglobal-minimum) transition structure to calculate the barrier height could easily lead to overestimates of 1 kcal/ mol or even more. These findings highlight the danger of selecting the saddle point structure based on considerations of just the global minimum structure of the reactant, as this strategy can lead to the overestimation of barrier heights due to the selection of a TS structure that is in fact not the lowestenergy structure. Therefore, though attention must be paid to the accuracy of the electronic model chemistry used to calculate barrier heights, in the case of complex reactions such as (R1), it must also be ensured, through comprehensive conformational structure searches, that the barrier height is calculated from the properties of the appropriate reactant and transition structures. The species-specific F factors for all torsion-containing species whose partition functions are needed for determining the rate constants for (R1) are given in Table 3. The table shows that multistructural factors increase the transition-state partition function by up to a factor of 650. Part of this is canceled by the reactant F factor in the rate constant calculation

Figure 3. Calculated (a) VMEP and (b) vibrationally adiabatic groundstate potential, VGa , as a function of the reaction coordinate, s, scaled to a reduced mass of 1 amu for (R1) (abstraction at C-1). Points 2−7 are nonstationary Shepard points calculated using M08-SO/MG3S. Point 1 corresponds to the M08-SO/MG3S saddle point.

probably be eliminated by using straight direct dynamics or by using the MCSI method only to generate the MEP, followed by Hessian calculations on that MEP;58 however, the curves shown are smooth enough for variational TS calculations and SCT tunneling calculations, and so we did not need to recalculate them at a greater expense. Tunneling treatments going beyond the SCT approximation, such as large-curvature tunneling59−63 (LCT) and microcanonically optimized multidimensional tunneling 62,63 (μOMT), were found to be unnecessary when we explored the kinetics of the analogous H-atom abstraction reactions from methanol by HO2 at temperatures above 400 K,27 and so we limited our attention here to the less expensive SCT method for tunneling. The thermal rate constants we calculated for (R1) are summarized in Table 4. For use in practical applications, we fit the calculated multistructural rate constants for temperatures

Table 3. Species-Specific and Reaction Multistructural Factors for Reaction R1 T (K)

FMS‑LQ BuOH

FTBuOH

FMS‑T BuOH

FMS‑LQ TS

FTTS

FMS‑T TS

FMS‑T R1

250 298.15 300 400 500 600 800 1000 1300 1500 1800 2000

10.1 11.2 11.3 13.3 14.8 16.1 18.0 19.3 20.6 21.3 22.0 22.4

1.36 1.42 1.42 1.51 1.57 1.60 1.58 1.50 1.34 1.22 1.06 0.96

13.7 15.9 16.0 20.0 23.3 25.7 28.4 28.9 27.5 25.9 23.3 21.6

17.9 24.7 24.9 43.1 65.8 91.2 144.4 195.5 263.0 301.5 350.6 378.6

1.31 1.42 1.42 1.66 1.90 2.10 2.36 2.42 2.30 2.13 1.85 1.67

23.5 35.0 35.5 71.8 125.0 191.6 340.2 473.9 603.9 643.0 649.6 631.1

1.72 2.20 2.22 3.58 5.37 7.45 12.0 16.4 21.9 24.8 27.9 29.3 12210



Article

Table 4. Thermal Rate Constants (cm3 molecule−1 s−1) for Reaction R1 T (K) 250 298.15 300 400 500 600 800 1000 1300 1500 1800 2000

kSS‑TST 6.73 4.18 4.78 1.17 3.77 4.30 1.12 9.61 8.55 2.48 8.62 1.68

× × × × × × × × × × × ×

kSS‑CVT −26

10 10−24 10−24 10−21 10−20 10−19 10−17 10−17 10−16 10−15 10−15 10−14

6.17 3.86 4.42 1.09 3.52 4.01 1.05 8.91 7.58 2.10 6.59 1.23

× × × × × × × × × × × ×

kSS‑CVT/ZCT −26

10 10−24 10−24 10−21 10−20 10−19 10−17 10−17 10−16 10−15 10−15 10−14

7.11 2.13 2.38 2.82 6.46 6.11 1.32 1.03 6.76 1.87 5.45 8.76

between 250 and 2000 K to a four-parameter formula,64 and we obtained the expression k MS ‐ CVT/SCT = 3.097 × 10−15 ⎛ −8.606(T + 187.8) ⎞ ⎛ T ⎞3.643 ⎜ ⎟ exp⎜ ⎟ ⎝ 300 ⎠ ⎝ R(T 2 + 3.527 × 104) ⎠

× × × × × × × × × × × ×

−25

10 10−23 10−23 10−21 10−20 10−19 10−17 10−16 10−16 10−15 10−15 10−15

kSS‑CVT/SCT 3.25 6.81 7.53 5.61 1.01 8.39 1.58 1.15 7.24 1.97 5.65 9.02

× × × × × × × × × × × ×

−24

10 10−23 10−23 10−21 10−19 10−19 10−17 10−16 10−16 10−15 10−15 10−15

kMS‑CVT/SCT 5.59 1.50 1.67 2.01 5.43 6.25 1.89 1.88 1.59 4.88 1.57 2.64

× × × × × × × × × × × ×

10−24 10−22 10−22 10−20 10−19 10−18 10−16 10−15 10−14 10−14 10−13 10−13

However, we must consider more than just the barrier height, and in particular it is interesting to compare the treatments of torsional anharmonicity, as these can heavily impact the calculated rate constants for complex reactions. In the work of Zhou et al., the torsional potential is obtained from scans around 6 internal rotations, which yielded 108 structures (3 × 3 × 3 × 2 × 2). These 108 structures were located by Zhou et al. by considering the abstraction of just one of the two available hydrogen atoms on carbon-1, and since each of these structures will possess a nonsuperposable mirror image structure that can only be generated from the abstraction of the other hydrogen atom, Zhou et al.’s calculations imply that there are actually 108 pairs of distinguishable transition structures. We note that this number exceeds the number of structures (94 pairs) obtained from the extensive conformational search in the present work, and we attribute this difference to Zhou et al.’s incorrect assumption that the torsions are separable. It is also interesting that the number of transition structures predicted by the individual 1-D scans for the C-3 reaction is the same (108) as that for the C-1 reaction in the work of Zhou et al.; however, the work by Seal et al.24 shows that the TS for the C-3 reaction has 131 pairs of structures, which is a further manifestation of the nonseparable nature of transition structure torsions in this kind of reaction. These discrepancies illustrate the need for a nonseparable method such as the MS-T method used in our work. Note that transition structures for both the C-1 reaction and the C-3 reactions have a chiral center, so that structures cannot be generated from their mirror images by internal rotations. We also note that in our work, the anharmonicities of seven torsions are treated by the MS-T method, but only six torsions are taken into account in the work of Zhou et al. (the torsion around the partial C−H bond was omitted); even though the methyl torsion does not generate new structures, it does contribute torsional potential anharmonicity. Another source of difference between our calculated rates and those obtained by Zhou et al. is the treatment of tunneling. The Eckart tunneling approximation was used in the work of Zhou et al., and the small curvature tunneling approximation with the reaction path calculated by the MCSI method was used in our calculations. It has been known for a long time that the Eckart tunneling approximation is oversimplified,33 and comparison of the ZCT and SCT columns of Table 3 confirms the importance of reaction-path curvature in the present reaction. At 500 K, Zhou et al.’s Eckart tunneling calculation increases the rate constant by a factor of 6.6,18 whereas the SCT transmission coefficient increases the rate constant by a factor of 2.9.

(4)

where R is the gas constant in kcal/mol, T is in K, and k is in cm3 molecule−1 s−1. The calculated multistructural rate constants are plotted in Arrhenius form along with those available from the literature in Figure 4. Figure 4 also shows the much lower rate constants

Figure 4. Comparison of computed rate constants for (R1) (abstraction at C-1): Zhou et al. (solid black line), Vranckx et al. (dashed black line), kMS‑CVT/SCT calculated in this work (blue triangles), kMS‑CVT/ZCT calculated in this work (red circles), and kSS‑CVT/SCT calculated in this work (green squares).

calculated with the single-structure local quasiharmonic method. The temperature range for which this reaction is most important in combustion is from about 700 to 1400 K (at higher temperatures, HO2 decomposes rapidly on the time scale of interest in combustion). At high temperature, we see better agreement of the present calculations with the previous work of Zhou et al.18 than with the work of Vranckx et al.13 At low and intermediate temperature our rate constants are lower than the calculations of Zhou et al.,18 which might at first be considered surprising because the barrier height used in our calculations is lower than that in the work of Zhou et al. 12211



Article

ing multidimensional tunneling. Comprehensive structure searches were performed to locate all of the torsional conformations of the transition state (in this work) and of nbutanol (in previous work),26 and the contributions of all of these structures were then explicitly accounted for in the final rate constant calculations in the form of conformational− rotational−vibrational partition functions that account for torsional anharmonicity. The implicit potential energy surface used to calculate the rate constants was validated by a high-level estimation of the barrier height, calculated as the energy difference between the lowest-energy conformations of the saddle point and the reactant; this validation employed the CCSD(T)-F12a/jul-cc-pVTZ//M08-SO/MG3S electronic model chemistry. Our calculated rate constants and the methods used to obtain them are compared with previous work. The final results of this paper are summarized in Table 1, Figure 4, and eq 4.

Another difference between the two sets of calculations is one that mainly affects the high-temperature results, namely that the reduced moment of inertia based on Pitzer’s formula is fixed at the most stable species for the torsional partition function calculations in the work of Zhou et al., but our MS-T calculations show that the product of Pitzer’s moments65 varies from 3.20 × 103 to 4.29 × 105 (amu Å2)7 and the determinant of the Kilpatrick and Pitzer moment-of-inertia matrix66 varies from 30.9 to 8.17 × 104 (amu Å2)7. Such a large variance of moments of inertia treated either by Pitzer’s independent torsion formula or by the Kilpatrick and Pitzer’s coupled torsion formula shows that it is insufficient to fix the moment of inertia as a constant in the torsional partition function calculations. Fixing the moments of inertia in all structures to the Pitzer moments of inertia of the global minimum changes the partition function of the transition state up to 10% at 2400 K using the current MS-T formula. Note that this fixed moment of inertia approximation is often used in solving 1-D Schrödinger equations for torsional potentials in literature, and not making this approximation is one of the improvements afforded by the MS-T treatment. The rate constants calculated by Vranckx et al. are considerably lower than ours and than those of Zhou et al. The potential energy surfaces used by Vranckx et al. and Zhou et al. are very similar. The much lower rates obtained by Vranckx et al. could be mainly due to their total neglect of multistructural effects (torsions are treated as small-amplitude vibrations with Hartree−Fock zero point energies and anharmonic fundamental frequencies). An underestimation of the rate when one uses a single-structural treatment is consistent with what one would expect from the results shown in Figure 4. The temperature-dependent Arrhenius activation energies computed from the four-parameter fit (eq 4) as

Ea = −R

■

ASSOCIATED CONTENT

S Supporting Information *

Graphical depiction of lowest-energy structure of the TS for (R1). Geometrical parameters and relative ZPE-inclusive and ZPE-exclusive energies of all transition structures for (R1) as optimized with M08-SO/MG3S. Complete ref 56. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address †

Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Ave., Cambridge, MA 02139. Notes

d ln k d(1/T )

The authors declare no competing financial interest.

■

(5)

ACKNOWLEDGMENTS We are grateful to Prof. William H. Green, Dr. Michael R. Harper, and Shamel S. Merchant for helpful discussions. This material is based on work supported as part of the Combustion Energy Frontier Research Center, funded by the Office of Basic Energy Sciences of the U.S. Department of Energy under Award No. DE-SC 0001198. Additional support was provided by the Office of Basic Energy Sciences under Grant No. DEFG02-86ER13579.

are listed in Table 5. As expected, the Arrhenius activation energy increases strongly as the temperature increases, and the Table 5. Activation Energy (kcal/mol) for Reaction R1 T (K)

activation energy (kcal/mol)

250 298.15 300 400 500 600 800 1000 1300 1500 1800 2000

10.6 12.1 12.2 14.1 15.2 15.95 17.15 18.2 20.1 21.3 23.2 24.5

■

REFERENCES

(1) Golombok, M.; Tierney, S. Ind. Eng. Chem. Res. 1997, 36, 5023− 5027. (2) Jacobson, M. Z. Environ. Sci. Technol. 2007, 41, 4150−4157. (3) Fargione, J.; Hill, J.; Tilman, D.; Polasky, S.; Hawthorne, P. Science 2008, 319, 1235−1238. (4) Wallner, T.; Miers, S. A.; McConnell, S. A Comparison of Ethanol and Butanol as Oxygenates Using Direct-Injection, SparkIgnition (DISI) Engine; ASME 2008 Internal Combustion Engine Division Spring Technical Conference, 2008, Chicago, Illinois. (5) Moss, J. T.; Berkowitz, A. M.; Oehlschlaeger, M. A.; Biet, J.; Warth, V.; Glaude, P.-A.; Battin-Leclerc, F. J. Phys. Chem. A 2008, 112, 10843−10855. (6) Sarathy, S. M.; Thomson, M. J.; Togbe, C.; Dagaut, P.; Halter, F.; Mounaim-Rousselle, C. Combust. Flame 2009, 156, 852−864. (7) Van Geem, K. M.; Pyl, S. P.; Marin, G. B.; Harper, M. R.; Green, W. H., Jr. Ind. Eng. Chem. Res. 2010, 49, 10399−10420.

old-fashioned way of characterizing a reaction by its activation energy as if it were a constant is shown to be very misleading.

5. SUMMARY In the present work, we calculated the rate constants of H abstraction from carbon-1 in n-butanol by the HO2 radical by using multistructural variational transition-state theory includ12212



Article

(8) Dagaut, P.; Sarathy, S. M.; Thomson, M. J. Proc. Combust. Inst. 2009, 32, 229−237. (9) Black, G.; Curran, H. J.; Pichon, S.; Simmie, J. M.; Zhukov, V. Combust. Flame 2010, 157, 363−373. (10) Grana, R.; Frassoldati, A.; Faravelli, T.; Nieman, U.; Ranzi, E.; Seiser, R.; Cattolica, R.; Seshadri, K. Combust. Flame 2010, 157, 2137− 2154. (11) Veloo, P. S.; Wang, Y. L.; Egolfopoulos, F. N.; Westbrook, C. K. Combust. Flame 2010, 157, 1989−2004. (12) Harper, M. R.; Van Geem, K. M.; Pyl, S. P.; Marin, G. B.; Green, W. H., Jr. Combust. Flame 2011, 158, 16−41. (13) Vranckx, S.; Heufer, K. A.; Lee, C.; Olivier, H.; Schill, L.; Kopp, W. A.; Leonhard, K.; Taatjes, C. A.; Fernandes, R. X. Combust. Flame 2011, 158, 1444−1455. (14) Hansen, N.; Harper, M. R.; Green, W. H. Phys. Chem. Chem. Phys. 2011, 13, 20262−20274. (15) Sarathy, S. M.; Vranckx, S.; Yasunaga, K.; Mehl, M.; Oβwald, P.; Metcalfe, W. K.; Westbrook, C. K.; Pitz, W. J.; Kohse-Hoinghaus, K.; Fernandes, R. X.; Curran, H. J. Combust. Flame 2012, 159, 2028−2055. (16) Yasunaga, K.; Mikajiri, T.; Sarathy, S. M.; Koike, T.; Gillespie, F.; Nagy, T.; Simmie, J. M.; Curran, H. J. Combust. Flame 2012, 159, 2009−2027. (17) Zador, J.; Taatjes, C. A.; Fernandes, R. X. Prog. Energy Combust. Sci. 2011, 37, 371−421. (18) Zhou, C.-W.; Simmie, J. M.; Curran, H. J. Int. J. Chem. Kinet. 2012, 44, 155−164. (19) Barone, V. J. Chem. Phys. 2004, 120, 3059−3065. (20) East, A. L.; Radom, L. J. Chem. Phys. 1997, 106, 6655−6674. (21) Pitzer, K. S. Quantum Chemistry; Prentice-Hall: Englewood Cliff, NJ, 1953. (22) Truhlar, D. G. J. Comput. Chem. 1991, 12, 266−270. (23) Sharma, S.; Raman, S.; Green, W. H. J. Phys. Chem. A 2010, 114, 5689−5701. (24) Seal, P.; Papajak, E.; Truhlar, D. G. J. Phys. Chem. Lett. 2012, 3, 264−271. (25) Yu, T.; Zheng, J.; Truhlar, D. G. Chem. Sci. 2011, 2, 2199−2213. (26) Zheng, J.; Yu, T.; Papajak, E.; Alecu, I. M.; Mielke, S. L.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2011, 13, 10885−10907. (27) Alecu, I. M.; Truhlar, D. G. J. Phys. Chem. A 2011, 115, 2811− 2829. (28) Zheng, J.; Truhlar, D. G. Faraday Discuss. 2012, 157, 59−88. (29) Truhlar, D. G.; Garrett, B. C. Acc. Chem. Res. 1980, 13, 440− 448. (30) Liu, Y.-P.; Lynch, G. C.; Truong, T. N.; Lu, D.-h.; Truhlar, D. G. J. Am. Chem. Soc. 1993, 115, 2408−2415. (31) Fernandez-Ramos, A.; Ellingson, B. A.; Garrett, B. C.; Truhlar, D. G. Rev. Comput. Chem. 2007, 23, 125−232. (32) Garrett, B. C.; Truhlar, D. G.; Grev, R. S.; Magnuson, A. W. J. Phys. Chem. 1980, 84, 1730−1748; Erratum. J. Phys. Chem. 1983, 87, 4554. (33) Truhlar, D. G.; Kuppermann, A. J. Am. Chem. Soc. 1971, 93, 1840−1851. (34) Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2008, 4, 1849−1868. (35) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. v. R. J. Comput. Chem. 1983, 4, 294−301. (36) Lynch, B. J.; Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2003, 107, 1384−1388. (37) Seal, P.; Papajak, E.; Truhlar, D. G. J. Chem. Phys. 2012, 136, 034306/1−10. (38) Alecu, I. M.; Zheng, J.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2010, 6, 2872−2887. (39) Zheng, J.; Mielke, S. L.; Clarkson, K. L.; Truhlar, D. G. MSTor, version 2011−3; University of Minnesota: Minneapolis, MN, 2012. (40) Zheng, J.; Mielke, S. L.; Clarkson, K. L.; Truhlar, D. G. Comput. Phys. Commun. 2012, 183, 1803−1812. (41) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479−483.

(42) Adler, T. B.; Knizia, G.; Werner, H.-J. J. Chem. Phys. 2007, 127, 221106/1−4. (43) Knizia, G.; Adler, T. B.; Werner, H.-J. J. Chem. Phys. 2009, 130, 054104/1−20. (44) Papajak, E.; Truhlar, D. G. J. Chem. Theory Comput. 2010, 6, 597−601. (45) Zheng, J.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2009, 5, 808−821. (46) Xu, X.; Alecu, I. M.; Truhlar, D. G. J. Chem. Theory Comput. 2011, 7, 1667−1676. (47) Kim, Y.; Corchado, J. C.; Villa, J.; Xing, J.; Truhlar, D. G. J. Chem. Phys. 2000, 112, 2718−2735. (48) Albu, T.; Corchado, J. C.; Truhlar, D. G. J. Phys. Chem. A 2001, 105, 8465−8487. (49) Tishchenko, O.; Truhlar, D. G. J. Chem. Theory Comput. 2009, 5, 1454−1461. (50) Tishchenko, O.; Truhlar, D. G. J. Chem. Phys. 2010, 132, 084109/1−6. (51) Page, M.; McIver, J. W., Jr. J. Chem. Phys. 1988, 88, 922−935. (52) Melissas, V. S.; Truhlar, D. G.; Garrett, B. C. J. Chem. Phys. 1992, 96, 5758−5772. (53) Tischenko, O.; Higashi, M.; Albu, T. V.; Corchado, J. C.; Kim, Y.; Villà, J.; Xing, J.; Lin, H.; Truhlar, D. G. MCSI 2010; University of Minnesota, Minneapolis, 2010. (54) Albu, T. V.; Tischenko, O.; Corchado, J. C.; Kim, Y.; Villà, J.; Xing, J.; Lin, H.; Higashi, M.; Truhlar, D. G. MC-TINKERATE 2010; University of Minnesota, Minneapolis, 2010. (55) Zhao, Y.; Truhlar, D. G. MN-GFM: Minnesota Gaussian Functional Module, version 4.3; University of Minnesota: Minneapolis, MN, 2009. (56) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (57) Black, G.; Simmie, J. M. J. Comput. Chem. 2010, 31, 1236−1248. (58) Zheng, J. J.; Seal, P.; Truhlar, D. G. Chem. Sci. 2012, DOI: dx.doi.org/10.1039/C2SC21090H. (59) Garrett, B. C.; Truhlar, D. G.; Wagner, A. F.; Dunning, T. H., Jr. J. Chem. Phys. 1983, 78, 4400−4413. (60) Kreevoy, M. M.; Ostovi, D.; Truhlar, D. G.; Garrett, B. C. J. Phys. Chem. 1986, 90, 3766−3774. (61) Garrett, B. C.; Joseph, T.; Truong, T. N.; Truhlar, D. G. Chem. Phys. 1989, 136, 271−284. (62) Liu, Y.-P.; Lu, D.-h.; Gonzalez-Lafont, A.; Truhlar, D. G.; Garrett, B. C. J. Am. Chem. Soc. 1993, 115, 7806−7817. (63) Fernandez-Ramos, A.; Truhlar, D. G. J. Chem. Phys. 2001, 114, 1491−1496. (64) Zheng, J.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2010, 12, 7782−7793. (65) Pitzer, K. S. J. Chem. Phys. 1946, 14, 239−243. (66) Kilpatrick, J. E.; Pitzer, K. S. J. Chem. Phys. 1949, 17, 1064− 1075.

12213


Biofuel Combustion. Energetics and Kinetics of Hydrogen Abstraction

Recommend Documents