Prediction of Experimentally Unavailable Product Branching Ratios for

Article pubs.acs.org/JACS

Prediction of Experimentally Unavailable Product Branching Ratios for Biofuel Combustion: The Role of Anharmonicity in the Reaction of Isobutanol with OH Jingjing Zheng, Rubén Meana-Pañeda, and Donald G. Truhlar* Department of Chemistry, Chemical Theory Center, and Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455-0431, United States S Supporting Information *

ABSTRACT: Isobutanol is a prototype biofuel, and sorting out the mechanism of its combustion is an important objective where theoretical modeling can provide information that is unavailable and not easily obtained by experiment. In the present work the rate constants and branching ratios for the hydrogen abstraction reactions from isobutanol by hydroxyl radical have been calculated using multi-path variational transition-state theory with small-curvature tunneling. We use hybrid degeneracy-corrected vibrational perturbation theory to show that it is critical to consider the anharmonicity difference of high-frequency modes between reactants and transition states. To obtain accurate rate constants, we must apply different scaling factors to the calculated harmonic vibrational frequencies at the reactants and at the transition states. The factors determining the reaction rate constants have been analyzed in detail, including variational effects, tunneling contributions, the effect of multiple reaction paths on transmission coefficients, and anharmonicities of low- and high-frequency vibrational modes. The analysis quantifies the uncertainties in the rate calculations. A key result of the paper is a prediction for the site dependence of hydrogen abstraction from isobutanol by hydroxyl radical. This is very hard to measure experimentally, although it is critical for combustion mechanism modeling. The present prediction differs considerably from previous theoretical work.

1. INTRODUCTION Isobutanol is a promising biofuel or fuel component.1−4 The reaction of hydroxyl (OH) radical with the fuel components is an important elementary reaction in all combustion mechanisms. For isobutanol, this reaction has four possible products:

calculate the rate constants for these reactions; however, the calculated rate constants do not agree with experiments for the temperature dependence. In our previous work, we demonstrated the importance of including the torsional anharmonicity in reaction path calculations using reaction R1a as an example.12 In particular, we found that, without considering the variation of the torsional anharmonicity along the reaction path, one cannot obtain a reasonable free energy profile for this reaction, nor can one locate the variational transition state. In the present article, we will calculate the total rate constants of all four reactions and the branching ratios. The uncertainty of rate constant calculations is a contentious issue for the kinetics modeling community. One often sees very good agreement between calculated rate constants and experimental results in the literature. However, it would not be surprising to learn that many cases of such agreement are largely due to error cancellations. The rate constants calculated by variational transition-state theory (VTST)13−18 depend on many factors, such as barrier height, vibrational frequencies, anharmonicity, consideration of multiple structures, variational effects, and tunneling contributions. Sometimes improving one of these factors may even lead to calculated rate constants in worse agreement with experiments because some error

(CH3)2 CHCH 2OH + OH → (CH3)2 CH•CHOH + H 2O

(R1a)

→ (CH3)2 •CCH 2OH + H 2O

(R1b)

→ (CH3)CH(•CH 2)CH 2OH + H 2O

(R1c)

→ (CH3)2 CHCH 2O• + H 2O

(R1d)

The total rate constants of hydrogen abstraction by hydroxyl radical from isobutanol (i.e., the sums of the reaction rates for the above four processes) have been measured at both low and high temperatures.5−9 But, due to many secondary reactions and the difficulties of detecting product radicals, it is difficult to determine by experiment the branching ratios, i.e., the relative rates of the hydrogen abstractions from the four distinguishable sites of isobutanol. This set of reactions also provides great challenges for computation because many possible sources of error may prevent the calculated reaction rates from being accurate. Some theoretical studies10,11 have attempted to © 2014 American Chemical Society

Received: February 1, 2014 Published: March 20, 2014 5150

dx.doi.org/10.1021/ja5011288 | J. Am. Chem. Soc. 2014, 136, 5150−5160

Journal of the American Chemical Society

Article

reproduce accurate zero-point energies (ZPEs) of a set of small molecules. The parametrized scaling factors are expected to be transferrable for stable molecules, and they have been widely used. But when the ZPE of a molecule or a transition state is very large, a small percentage uncertainty in the scaling factor could lead to considerable errors in estimating ZPEs and in rate calculations, and the validity of the general scaling factor method has not ben demonstrated for transition states. Here we use hybrid28 degeneracy-corrected29 second-order30−32 vibrational perturbation theory (HDCVPT) to assess the anharmonicity at the transition state. The full details of HDCVPT for transition states are given in ref 28. A detailed analysis will be performed to discuss the effects of barrier heights, torsional anharmonicity, high-frequency-mode anharmonicity, variational effects, and tunneling contributions on accurate rate calculations for reactions R1a−R1d.

cancellation disappears. This phenomenon can be very discouraging for developing new methods and for investing effort to improve calculations. It is therefore critical to assess the sources of error to learn which ones limit the accuracy and by how much. This is especially important when one predicts detailed results that cannot easily be measured experimentally. The present study reports a prediction of such experimentally inaccessible and/or unavailable quantities, in particular, the sitespecific reaction rates of OH with various sites of isobutanol, and we supplement the calculations with a detailed discussion of possible sources of error in this kind of calculation. Among all the factors limiting the accuracy of rate calculations, the most attention has been given to classical barrier heights, which have indeed been a challenge to calculate and which have limited the attainable accuracy in rate constant predictions ever since transition-state theory was proposed in the 1930s. Nevertheless, since the classical barrier height V⧧ is a key factor determining chemical reaction rates, it is often used to determine the branching ratio (for a reaction that has more than one reactive channel) or reaction mechanism, especially for large, complex systems. When the difference ΔV⧧ of barrier heights is large, consideration of the Boltzmann factor, exp(−ΔV⧧/kBT), where kB is Boltzmann’s constant and T is temperature, may be sufficient to make a correct conclusion about which of two competitive reactions has a larger rate constant, but the uncertainty in the barrier height difference leads to uncertainty in the quantitative prediction of a branching ratio. The barrier height issue is widely appreciated, but the main focus of this article is on other sources of uncertainty, which are less widely studied. As barrier height calculations become more accurate, the uncertainties in the other factors become more crucial to estimate. According to transition-state theory, a chemical reaction rate is also proportional to the ratio of the transition state partition function to the partition function of the reactant or reactants. The single-structure, harmonic-oscillator (HO), rigid-rotor approximation is often used to calculate partition functions. Anharmonicities, multiple structures, and coupling of overall rotation to internal rotation all tend to become more important when the size of the molecule increases, and it is necessary to consider all three of these issues to make rate calculations quantitatively accurate, or sometimes even to make them qualitatively accurate. Furthermore, “large” molecules (e.g., molecules with alkyl chains having four or more carbons) often have more troublesome wide-amplitude motions, usually dominated by torsions. In a general sense, all four of these problematic considerations are manifestations of anharmonicity. Although anharmonicity becomes more important in complex systems, the HO approximation is nevertheless widely used because accurate treatment of anharmonicity is often computationally too expensive to be applied on an everyday basis to large molecules and transition states. In our recent studies19−24 we have shown the importance of including all the conformational structures and torsional potential anharmonicity in calculating partition functions and rate constants. Single-bond torsions are almost always lowfrequency modes, but we must also consider anharmonicity in high-frequency modes. Often the anharmonicity of highfrequency vibrational modes is taken into account by scaling harmonic vibrational frequencies.25−27 The scaling factor for a given model chemistry (a model chemistry is a combination of an electronic structure method or exchange-correlation functional and a basis set) can be determined by parametrization to

2. COMPUTATIONAL METHODS 2.1. Electronic Structure Methods. The M08-HX33 exchange− correlation functional and the 6-311+G(2df,2p)34,35 basis set were used for conformational structure searching, geometry optimizations, and dynamics calculations. For H, C, and O atoms, the 6311+G(2df,2p) basis set is identical to the MG3S36 basis set; therefore, we use the MG3S notation through this article for its brevity. The CCSD(T)-F12a37−39 method with the jun-cc-pVTZ40 basis set was used to perform single-point calculations of barrier heights using the M08-HX/MG3S-optimized geometries. The CCSD(T)-F12a/jun-cc-pVTZ calculation can be viewed as an efficient way to yield (without extrapolation) results close to the complete basis set (CBS) limit of the standard CCSD(T) method, where CCSD(T) denotes coupled cluster theory with single and double excitations and a quasiperturbative treatment of connected triple excitations. Therefore, the CCSD(T)-F12a/jun-cc-pVTZ calculations will be abbreviated CCSD(T)/CBS, and they are taken as the benchmarks to validate the accuracy of the Kohn−Sham density functional calculations. Møller−Plesset second-order perturbation theory (MP2)41 and the MPW1K42 exchange−correlation functional with the 6-31+G(d,p) basis set are used for the calculation of cubic and quartic force constants needed for the anharmonicity calculations All electronic structure calculations were performed by using the Gaussian 0943 program, except for the CCSD(T) calculations, which were carried out by using the Molpro44 program. 2.2. Multi-structural Method with Torsional Anharmonicity. The multi-structural method with torsional anharmonicity (MST)19,20,45 is used for calculating conformational−rovibrational partition functions for reactants and transition structures. In particular, the MST method used here is based on a coupled torsional potential.20 A full explanation of the MS-T method based on a coupled torsional potential can be found in ref 20, but the method is briefly described here for the reader’s convenience and for discussion purposes. For a molecule or a transition state with J distinguishable conformational structures, the MS-T conformational−rotational−vibrational partition function is J MS‐T Q con = ‐rovib

∑ Q jT

(1)

j=1

Q jT ≡ Q rot, j exp(− βUj)Q jQHTj F

Tj = (2π ℏβ)t /2

∏m = 1 ωj , m F−t ∏m = 1 ωj , m

det Dj t ∏ τ = 1 Mj , τ

(2)

×

t

∏ exp(−βWj , η/2)I0(βWj , η/2) η=1

5151

(3) dx.doi.org/10.1021/ja5011288 | J. Am. Chem. Soc. 2014, 136, 5150−5160


Article

where QTj is the contribution of structure j to the total conformational−rotational−vibrational partition function with torsional anharmonicity, Qrot,j is the symmetry-corrected classical rigid-rotor rotational partition function for structure j (i.e., it is based on the classical rigid-rotor model but with the addition of quantum-mechanical symmetry numbers), Uj is the relative potential energy of structure j, is its quasiharmonic (QH) oscillator vibrational partition function QQH j partition function, and Tj is a factor to include coupled torsional potential anharmonicity for all torsional modes in the structure j. We use the notation quasiharmonic to refer to an approximation that uses the harmonic oscillator partition function formulas, but with effective frequencies that account for anharmonicity (and possibly also for systematic errors in the electronic model chemistry). In the factor Tj, t is the number of torsional modes, F is the number of bound vibrational modes (stretches, bends, and torsions), Dj is the moment of inertia matrix from Kilpatrick and Pitzer,46 Wj,η is the effective torsional barrier height based on a coupled torsional potential, and ωj,m and ωj , m are respectively the normal-mode quasiharmonic frequencies and the torsion-projected normal-mode frequencieswith both sets scaled. For reasons explained in section 2.3, the scaling factors we used are 0.973 for the reactants harmonic frequencies27 and are various reparametrized scaling factors (see Table 2 and section 3.2) for the transition-state harmonic frequencies. We define the factor

FXMS‐T =

where K is the number of distinguishable conformational structures of the transition state, Q⧧elec is the electronic partition function of the conventional transition state (which is denoted as ⧧), ΦMS‑T,R is the bimolecular reactants’ partition function per unit volume with torsional potential anharmonicity (with its zero of energy at the is the potential energy of the lowest-energy reactant structure), QT,⧧ k rovibrational partition function of saddle point k with torsional is the tunneling transmission potential anharmonicity (T), κSCT k is the CVT coefficient in the small curvature approximation, ΓCVT k recrossing transmission coefficient (see eq 8), and V⧧ is the classical barrier height (defined as the change in the potential energy from the zero of energy used for the reactants’ partition function to the zero of energy used for the transition-state partition function, where the zero of energy is defined as the lowest energy of all conformations of reactants and transition state, respectively). We define γ as a generalized transmission coefficient by P

γ=

(4)

to illustrate the effect of multi-structural torsional anharmonicity on a reactant (X = isobutanol) or a transition state (X = TS), and we define the factor MS‐T Frxn =

MS‐T FTS MS‐T Fisobutanol

(7a)

K

(5)

kSCT‐CVT/SCT = γ

to illustrate the effect of multi-structural torsional anharmonicity on the reaction rate. Note that FMS‑T OH for OH radical is equal to 1. 2.3. Hybrid Degeneracy-Corrected Second-Order Vibrational Perturbation Theory. The anharmonic ZPEs are calculated for the global minimum of isobutanol and for the lowest-energy saddle points of reactions R1a and R1b. These anharmonicity calculations are based on HDCVPT. These calculations require numerical third and fourth derivatives of the potential energy with respect to nuclear displacements, and these were calculated as first and second derivatives of analytic Hessians with step size 0.01 Å in MP2 calculations and step size 0.025 Å in MPW1K calculations. The step size used for the MPW1K calculations is larger due to the numerical integral noise in density functional calculations. Because high-order derivatives are very sensitive to the precision of analytic Hessians, a very fine grid was used in MPW1K calculations; this grid has 96 radial shells around each atom and a spherical product angular grid having 32 θ points and 64 φ points in each shell. The HDCVPT calculations were performed by Gaussian 0943 program (revision D.01). 2.4. Reaction Rate Calculations. Reaction rate constants are calculated under the assumption that the system is in the low-pressure plateau region, where the pressure is low enough for isolated binary collisions to be the reaction mechanism and there are no collisionstabilized pre-transition-state complexes (that could cause the rate constant to be pressure dependent) but high enough that reactant states are in local equilibrium with each other during the reaction. This plateau region is discussed elsewhere.23 We employ multi-path canonical variational transition-state theory (MP-CVT)22,23 with the small-curvature tunneling (SCT)47,48 approximation to calculate reaction rate constants. The MP-CVT/ SCT rate constant is given for a bimolecular reaction23 by

k MP‐CVT/SCT =

P

∑k = 1 Q kT,⧧

where P, which is less than or equal to K, is the number of distinguishable transition structures whose tunneling and CVT transmission coefficients are explicitly calculated. In the present work we will use P = 4, and we base the SCT tunneling and CVT transmission coefficients on the reaction paths through the four lowest-energy transition structures for each of the four reactions. Note that this requires calculating only two reaction paths for each reaction, since for each reaction the four lowest-energy reaction paths are two pairs of enantiomeric reaction paths. We will also compare the calculations with P = 4 to calculations with P = 2. Note that because the two lowest-energy paths are mirror images, one gets the same results with P = 2 as with P = 1. We then approximate eq 6 as

MS‐T Q con (X) ‐rovib

Q 1QH(X )Q rot,1(X)

T,⧧ ∑k = 1 κkSCT Γ CVT k Qk

⧧ T,⧧ kBT Q elec ∑k = 1 Q k exp(− βV ⧧) h Φ MS‐T,R

(7b)

Note that eq 7b is identical to eq 6 if P = K. The CVT recrossing transmission coefficient is calculated in the present work as12

ΓCVT = k

Q kT,CVT exp( − βV kCVT) Q kT,⧧ exp( − βV k⧧)

(8)

where V⧧k is the potential energy at saddle point k, VCVT is the potential k energy at the variational transition state along the reaction path is the rovibrational partition through saddle point k, and QT,CVT k function that includes torsional anharmonicity at the canonical is the variational transition state along the reaction path, and QT,⧧ k same quantity at the saddle point. Therefore, unlike the original version of MS-CVT/SCT21 and MP-CVT/SCT,22 torsional potential anharmonicity is included in the reaction path calculations, as described in a previous article.12 Note that potential energies denoted by V are measured with respect to a geometry-independent overall zero of energy (the equilibrium potential energy of reactants), whereas potential energies denoted by U are measured with respect to the lowest energy structure at the same value of the reaction path progress variable s, where s = 0 for the saddle points, s = −∞ for reactants, and s is finite for generalized transition states. To understand how the various factors affect the rate constant, we separate the multi-structural torsional anharmonicity in the rate calculations by rewriting eq 7b as ⧧ HO,⧧ ⧧ kBT Q elecQ 1 Q rot,1 h Φ1HO,R ⧧ exp(− βV )

MS‐T k MP‐CVT/SCT = γFrxn

⧧ K kBT Q elec T,⧧ exp(− βV ⧧) ∑ κkSCT Γ CVT k Qk MS‐T,R h Φ k=1

MS‐T TST ≈ γFrxn k1

(6) 5152

(9)



Article

where kTST is the rate constant of conventional transition state using 1 the quasiharmonic approximation and using only one lowest-energy structure of the reactants and one lowest-energy structure of the transition state. The reaction paths and reaction-path-dependent generalized normal-mode frequencies associated with the lowest-energy saddle points were calculated by direct dynamics. Generalized normal-mode frequencies along the reaction paths are calculated by using nonredundant internal coordinates49,50 with the reaction coordinate projected out, and these are used to evaluate the vibrationally adiabatic ground-state potential energy curves needed for the tunneling calculations. The conformational−rotational−vibrational partition functions are calculated by using the MSTor45,51,52 program. The POLYRATE53 and GAUSSRATE54 programs are used to calculate reaction paths and single-structure reaction rate constants; these two programs incorporate the MS-T torsional anharmonicity methods that are also available as a separate program.51

minimum structure due to a steric effect. Therefore, it is not reliable to assume that the conformation of the lowest-energy structure of the reactant remains unchanged in the lowestenergy conformation in the transition structure, as has sometimes been done.55 To find the lowest-energy transition structure in complex reactions and to calculate accurate reaction barrier heights, an exhaustive conformational search has to be performed. The energies of the transition structures relative to the lowest one for each reaction are shown in Figure 2. These energies are

3. RESULTS AND DISCUSSION 3.1. Conformational Structures: Geometries, Energies, and Partition Functions. An exhaustive conformational search was performed for isobutanol and for the transition structures of reactions R1a−R1d. Nine conformational structures are found for isobutanol. The numbers of conformational structures for the R1a−R1d transition states are 20, 18, 96, and 16, respectively. Reaction R1c has much a larger number of transition structures than the others because the symmetry of the methyl groups in isobutanol is broken when OH radical abstracts a hydrogen atom from the methyl group. The lowest-energy structures of isobutanol and each of the TSs are plotted in Figure 1, which shows values for the H−O−C−C and O−C−C−H dihedral angles.

Figure 2. Energies of the transition structures relative to the lowest one for each reaction, as calculated by the M08-HX/MG3S method.

calculated by the M08-HX/MG3S method. In this figure, each bar represents a pair of enantiomers that are mirror images of each other. The structures that form hydrogen bonds in R1a, R1b, and R1c have relative energies at least 0.7, 0.5, and 2 kcal/ mol, respectively, lower than those without a hydrogen bond. The classical barrier heights calculated by M08-HX/MG3S and CCSD(T)/CBS are listed in Table 1. Note that classical Table 1. Calculated Forward Barrier Heights (in kcal/mol)a method M08-HX/MG3S CCSD(T)/CBS a

Figure 1. Lowest-energy structures of isobutanol and R1a−R1d transition structures. The numbers are for two dihedral angles, H−O− C−C and O−C−C−H.

R1a

R1b

R1c

R1d

−0.44 −0.34

−1.19 −1.37

0.80 0.94

3.25 2.99

The spin−orbit energy −0.20 kcal/mol for OH radical is included.

barrier height is the energy difference between the lowestenergy structure of reactants and the lowest-energy transition structure. The spin−orbit energy, −0.20 kcal/mol, for the OH radical is included in all barrier height calculations. (The spin− orbit energy is assumed negligible at the transition state because there is no low-lying excited electronic state there.) The barrier heights calculated by the M08-HX/MG3S method agree with those computed by the CCSD(T)/CBS method within 0.3 kcal/mol for H abstractions from all sites. For the O site abstraction, the CCSD(T)-F12a method may not be reliable because of the large multi-reference character of this transition state, but this channel should have the least contribution to the total rate, so we still use the coupled cluster method as a benchmark to get a rough estimate of the barrier height of this channel.

The conformation of the lowest-energy transition structure need not be the same as that of the reactant. For example, in the TSs for R1a and R1b, the C−O bond is rotated with respect to its position in the lowest-energy conformer of isobutanol to form a hydrogen bond with a hydroxyl radical. This hydrogen bond stabilizes the transition structures. If the isobutanol fragment in the TS for R1b keeps the same conformation as that in the reactant global minimum, the energy of the TS is 3.4 kcal/mol higher. In the lowest-energy transition structure for R1d, although there is no hydrogen bond, the H−O−C−C dihedral angle is different from that in the reactant global5153



Article

Figure 3. Percentage contribution of the lowest-energy and the highest-energy transition structures to total MS-T partition functions of transition states. The red lines are for the lowest-energy structures and blue lines for the highest-energy structures. The energy labels are relative to the zeropoint-exclusive lowest-energy transition structure in the corresponding reaction.

calculated harmonic frequencies by a scale factor parametrized either to reproduce the accurate ZPE or to reproduce the accurate harmonic frequency. These scaling factors are denoted respectively by λZPE and λH.27 One can interpret λH as correcting for the inexactness of the electronic structure calculation, and λZPE as including both a correction for the harmonic frequency that arises from the model chemistry and a correction for the anharmonicity of the ZPE. Therefore, the ratio λZPE/λH can be viewed as the correction purely due to anharmonicity in the zero-point level; we label this ratio as λAnh in this article. Scaling factors can be parametrized to yield correct results on average for a database of molecules;27 we call these general parameters if the database is diverse and specific range parameters if the database contains only a set of closely related molecules and is not intended to be representative of broader functionality. Alternatively, one can use scale factors optimized for a specific reaction or even for a specific reactant or specific transition state; either of these choices would be called specificreaction parameter (SRP) choices. The generally parametrized frequency scaling factor λZPE obtained from the F38/10 database for M08-HX/MG3S method is 0.973; this is the product of λAnh = 0.988 and λH = 0.984.27 The F38/10 database contains small, stable molecules, mainly with high-frequency and mid-frequency modes (stretches and bends); furthermore, the ZPE is sensitive mainly to high- and mid-frequency modes because they make the largest contributions. Thus, there is no reason to expect λZPE to be valid for low-frequency modes. Furthermore, one could question whether it applies to butanol, because butanol is larger than the molecules for which parametrization was done, or to the hydrogen abstraction transition states, because they too are large and they are also metastable. We will pursue this question in the following way. We will assume that the general

One general question is how important it is to include the high-energy conformational structures in the calculation of the partition functions. The simplest argument would say that lowenergy conformers are more important than the high-energy conformers because the Boltzmann factor decreases exponentially when the energy increases, and that factor is the dominant consideration at very low temperatures. However, high-energy conformers usually have more floppy modes that make larger contributions to partition functions than do the tight modes. The percentage contributions (to the total MS-T partition functions) of the lowest-energy structures are compared to those of the highest-energy structures in Figure 3. (Note that, for computing such percentages, the Boltzmann factors are included in the partition function for each structure.) Above 1200 K, the highest-energy structure contributes even more than the lowest-energy structure to the total partition f unction for R1a, R1b, and R1c, and the two structures have similar contributions for R1d. Even at lower temperatures, between 600 and 1200 K, the contributions of the highest-energy structures are not negligible. Figure 3 shows that, for medium and high temperatures, there is no direct correlation between a structure’s contributions to the total partition function and its potential energy. Therefore, it is clear that all conformational structures, regardless of their relative energies, have to be taken into account in order to obtain quantitatively accurate thermochemical and kinetic results at medium and high temperatures. 3.2. Specific-Reaction-Parametrized Scaling Factor. A widely used prescription, which we call the quasiharmonic approximation, is to use the harmonic oscillator formulas for ZPEs and partition functions, but with effective frequencies intended to give more accurate results than the accurate harmonic frequencies would yield. Two special cases are where the effective frequencies are obtained by multiplying a set of 5154



Article

(d,p) methods. We noticed that the λAnh value given by the MP2/6-31+G(d,p) method is considerably smaller than that given by the MPW1K/6-31+G(d,p) method for R1a and R1b TSs. Although the reaction-coordinate mode has no contribution to the ZPE of the transition structure, this difference between MP2 and MPW1K is due to anharmonicity arising from the coupling between reaction-coordinate modes and the bound modes,28 which is larger for MP2 than for MPW1K for these two transition states. We chose the SRP scaling factors obtained from the MPW1K/6-31+G(d,p) method in our kinetics calculations because MPW1K is more accurate than MP2 for barrier height calculations,56,57 and transition-state structures and vibrational frequencies are often correlated with barrier heights. 3.3. Multi-structural Torsional Anharmonicity. Figure 4 shows the multi-structural anharmonicity factors FMS‑T and X

parameter λH is indeed broadly applicable because it corrects the error that arises from the M08-HX/MG3S method itself and it is supposed to be independent of the system. That places the doubt on the generality of λAnh, and it is certainly possible that the anharmonicity factor has a structure dependence that makes it different for isobutanol and transition states than for the molecules in F38/10. To answer this question of the broadness of the applicability of λAnh, we performed anharmonic ZPE calculations using the HDCVPT method. Ideally, the anharmonic calculations should be performed using the M08-HX/MG3S potential energy surface, and then a direct comparison between harmonic ZPE and anharmonic ZPE can provide a specific-reaction parametrized scaling factor for the M08-HX/MG3S method to be used in the kinetics calculations. However, we found that the precision of the analytic Hessian calculated by meta-GGA exchange-correlation functionals such as M08-HX is very sensitive to the integral grid, and we were unable to converge the higher-order numerical derivatives even by using the largest affordable grid. Therefore, MP2 and the hybrid GGA functional MPW1K with 6-31+G(d,p) basis set are used for anharmonic ZPE calculations using the lowest-energy structures of isobutanol and the transition state of each reaction. Table 2 Table 2. Harmonic Zero-Point Energies and Anharmonic Zero-Point Energies Calculated by Using the HDCVPT Method (in kcal/mol) isobutanol ZPE(Harm) ZPE(Anh) λAnh λZPE ZPE(Harm) ZPE(Anh) λAnh λZPE

R1a TS

R1b TS

MP2/6-31+G(d,p) 92.17 92.56 89.36 88.67 0.969 0.958 0.956 0.945 MPW1K/6-31+G(d,p) 88.24 93.61 93.75 87.08 91.63 91.26 0.987 0.979 0.973 0.971 0.963 0.958 87.82 86.64 0.987 0.973

R1c TS

R1d TS

92.51 90.64 0.980 0.966

92.16 90.16 0.978 0.965

93.01 91.23 0.981 0.965

92.35 90.44 0.979 0.964

lists the calculated ZPEs. From these we calculated λAnh as the ratio of anharmonic ZPE to harmonic ZPE for the given model chemistry. Then we calculated the SRP scaling factor λZPE as the product of λAnh and the general parameter λH discussed in the previous paragraph. We also give λAnh and λZPE in Table 2. The values λAnh and λZPE for isobutanol as obtained by MP2 and MPW1K methods agree well with each other, and also they agree well with the general values for the M08-HX/MG3S method. We also verified that the scaling factor λZPE = 0.973 for OH radical reproduces the experimental ZPE with only a 0.003 kcal/mol discrepancy. Therefore, the scaling factor λZPE = 0.973 is assumed to be applicable to isobutanol and OH radical without specific-reaction modification. As found in a previous work,25 λAnh is almost independent of the model chemistry for the F38/10 database, and we have just seen that it is also almost independent of model chemistry for isobutanol. However, Table 2 shows that λAnh values of transition states are smaller than their counterparts for isobutanol; this result is found with both methods. This means that the studied transition states here need a larger correction for anharmonicity. The SRP scaling factors λZPE are listed in Table 2 for both methods; they are obtained by multiplying λAnh by 0.986 and 0.984 respectively for the MP2/6-31+G(d,p) and MPW1K/6-31+G-

Figure 4. (a) Multi-structural torsional anharmonicity factors FMS‑T for X isobutanol and all transition states. (b) Multi-structural torsional for the four reaction rates. The right anharmonicity factor FMS‑T rxn ordinate scale is for the R1c transition state and reaction rate, and the left ordinate scale is for the remaining channels, R1a, R1b, and R1d.

FMS‑T rxn . Because the transition state of R1c has many more conformational structures than the transion states for the other reactions, the factors for the R1c transition state and R1c reaction rate are much larger than for the remaining channels. The transition state of R1c has 96 conformational structures, and FMS‑T has a maximum value of 85 at 1400 K, which means X that each conformational structure contributes to the total partition function almost equally, on average, at this temperature, even though some structures have relatively high energies. The multi-structural anharmonicity factors of reaction, MS‑T FMS‑T for the rxn , become smaller than the corresponding FX transition state because of cancellation with FXMS‑T for 5155



Article

height are not populated, as discussed previously,23 and that nonclassical reflection from the barrier top is negligible. The tunneling transmission coefficient is greater than 2 for reactions R1c and R1d at temperatures lower than 300 K. The lowerenergy paths have slightly larger tunneling transmission coefficients, as shown in Figure 5. 3.5. Rate Constants. The calculated total MP-CVT/SCT rate constants are plotted in Figure 6, together with the

isobutanol. However, FMS‑T still has a considerable impact on rxn the reaction rate, especially for R1c, and it shows that one may make a large error in the calculations of the rate constants if multi-structural torsional anharmonicity is ignored. 3.4. Transmission Coefficients. Figure 5 shows the recrossing transmission coefficients ΓCVT and tunneling transmission coefficients κSCT for all calculated paths.

Figure 6. Calculated total MP-CVT/SCT rate constants (k = ka + kb + kc + kd) using various approaches: red dotted line is calculated using standard frequency scaling factor 0.973 for reactants and transition states (P = 4); blue dashed line is calculated with each transition state scaled by its own SRP scaling factor and reactants scaled by 0.973 (P = 4); black solid line is calculated by correcting the error of MS-T method for the blue dashed line (P = 4); black dashed line is the same as the black solid line except P = 2 is used for generalized transmission coefficient calculations; green dot-dash line is calculated by using Atkinson’s SAR method with updated substituent factors from refs 61 and 62.

Figure 5. Recrossing transmission coefficients and small-curvature tunneling transmission coefficients κSCT. The blue lines are R1a, red are R1b, green are R1c, and purple are R1d. Solid lines are for lowestenergy paths, and dashed lines are for the third lowest-energy paths. Note that the lowest-energy and the second lowest-energy paths are energetically degenerate, as are the third and fourth.

experimental data and rates calculated by the Atkinson’s structure−activity relationship (SAR) method59,60 with updated substituent factors.61,62 The total rate constants using the MPCVT/SCT formula, including torsional anharmonicity as a function of the reaction coordinate s and with all frequencies scaled by the standard λZPE (0.973) for the M08-HX/MG3S method, are denoted as kstd. They are much lower than experimental values and cannot correctly reproduce temperature dependence, although the M08-HX/MG3S barrier heights used in the rate calculations are in excellent agreement with the best estimated values. The total rate constants calculated with transition-state frequencies scaled by the SRP scaling factors λZPE obtained from the MPW1K/6-31+G(d,p) method (see Table 2) and with reactant frequencies scaled by the standard λZPE (0.973) are denoted as kSRP. Because smaller scaling factors are applied to transition states, the rate constants kSRP increase. Although the scaling factors are reduced by only 0.8−1.6%, the rate constants kSRP increase dramatically, e.g., by a factor of 22 at 200 K and by a factor of 2 at 2400 K, and they predict the correct temperature dependence. It is clear from the above that the negative temperature dependence at low temperature is caused mainly by the larger

First consider the recrossing transmission coefficients. Figure 5 shows that all of them deviate very significantly from unity at temperatures below 500 K, and the deviations for reactions R1a, R1b, and R1c remain large up to 2400 K. The deviation of the recrossing transmission coefficient from unity is called a variational effect. In previous work,12 it was demonstrated that the large variational effect of reaction R1a is caused by a significant change of one of the frequencies along the reaction path, which is a common occurrence.58 Here we see that the cause of the large variational effects for the isobutanol reactions other than R1a is the same as for R1a. For a given reaction, different paths have different variational effects, and this fact is accounted for in the MP-CVT/SCT formula by weighting the paths according to their partition functions and transmission coefficients. Since we assume that the system is in the low-pressure plateau region, the tunneling transmission coefficients for R1a and R1b are unity because their barrier heights are negative, and it is assumed that the energy levels lower than barrier 5156



Article

⎛ T + 31.91 ⎞2.424 ⎟ k final = 4.179 × 10−13⎜ × ⎝ 300 ⎠ ⎡ 1.658(T + 31.91) ⎤ exp⎢ ⎥ ⎣ R(T 2 + 1018) ⎦

anharmonicity in the transition state that lowers the zero-pointinclusive barrier heights. However, the rate constants kSRP are still lower than experimental values by a factor of 2.4 at 1000 K and by a factor of 1.5 at 250 K. Where do these errors originate? In our previous work,20 we found that the MS-T method, even with a coupled torsional potential, may have considerable uncertainty for a very strongly coupled system, e.g., the transition state of the hydrogen abstraction reaction of ethanol with hydroxyl radical. The MS-T method was evaluated by comparing the multi-structural reference classical (MS-RC) partition function to the partition function calculated by a fully coupled classical (FCC) phase space integral.20 Note that the MS-RC partition function is essentially the torsional part of the MS-T partition function (see ref 20 for details), and the error of the MS-RC partition function is considered to be an indication of the error of the torsional treatment in the MS-T method. We replot the ratio of the MS-RC partition function (QMS‑RC) to the FCC partition function (QFCC) in Figure 7 for the transition

(10)

where R is the gas constant in kcal/mol. Branching fractions of each hydrogen abstraction reaction from the MP-CVT/SCT calculations (the final ones with P = 4) and from McGillen et al.9 are plotted in Figure 8, together

Figure 8. Branching fraction of the hydrogen abstraction from all sites of isobutanol by hydroxyl radical. The blue lines are R1a, red are R1b, green are R1c, and purple are R1d. Solid lines are calculated in this work using the MP-CVT/SCT method (P = 4), and dashed lines are from McGillen et al.,9 obtained by global fitting with constraint by the acetone yield to a value closer to 40% obtained in their unpublished study, by the results from the shock tube study of Pang et al.8 over the temperature range 900−1150 K, and by SAR estimations above 1200 K. The branching fractions for R1b from Pang et al.8 (diamond symbol) and from Andersen et al.7 (square symbol) are also plotted for comparison. Figure 7. Plot of QMS‑RC/QFCC vs temperature for the transition state of hydrogen abstraction from the C-1 carbon of ethanol by hydroxyl radical. The data are taken from ref 20.

with the branching fractions of R1b reported by Pang et al.8 and by Andersen et al.7 It is difficult to determine the branching ratios experimentally because of the many side reactions and difficulties of detecting product radicals, and branching ratios, if reported at all, are usually inferred indirectly. Andersen et al. reported a 61 ± 4% branching fraction for R1b at 296 K, as inferred from assuming that acetone is produced only by abstraction from C-2 and that all abstraction from C-2 leads to acetone and from their own measurements of the end-product acetone yield at 296 K. McGillen et al.9 concluded that Andersen et al.’s branching ratio is too high due to a contribution of secondary acetone formation from isobutanal oxidation. McGillen et al. estimated the branching fractions by global fitting, constrained by assuming that an unpublished endproduct yield of 40% acetone means that 40% of the reaction occurs at C-2; they also used the results from the shock-tube study of Pang et al.8 over the temperature range 900−1150 K, and from SAR estimations by Bethel et al.62 above 1200 K, where the SAR estimation is used primarily to differentiate between abstraction at C-1 and abstraction at C-3. Our calculations give a 48% branching ratio for the R1b reaction, which is slightly closer to the result of McGillen et al. than to Andersen et al.’s result. Pang et al. combined their measured rate constants for nonR1b reactions with R1b rate constants used by Merchant and Green (personal communication between Pang et al. and Green) to get the R1b branching fractions around 1000 K, and these agree with our calculated values and with the SAR

state of the hydrogen abstraction reaction of ethanol with hydroxyl radical studied in ref 20; the ratio is denoted as g ≡ QMS‑RC/QFCC. For temperatures above 400 K, the MS-RC partition function is lower than the FCC partition function by a factor of 2 or more, and it is reasonable to assume that the same amount of error could be present in the MS-T partition functions for the transition states of the hydrogen abstraction reactions of isobutanol that are very similar to the ethanol reactions. Therefore, we divided the rate constants kSRP by the factor g in order to get the final rate constants kfinal. As shown in Figure 6, the final MP-CVT/SCT rate constants kfinal calculated using P = 2 and P = 4 are in excellent agreement with the available experimental measurements over the temperatures ranges of 250−370 and 900−1150 K. The differences between the P = 2 and P = 4 rates are about 10−20%. Our calculations also provide the value of the rate constants kfinal in the temperature region of 400−900 K, which is a very important range for engine ignition, but where there are no experimental data. We provide reliable values of the thermal rate constants from 200 to 2400 K. The final rate constant kfinal by the MPCVT/SCT method with P = 4 is fitted to our previously recommended four-parameter function,23,63 which gives 5157



Article

modes in R1b that significantly increases the relative R1b contributions. The site-specific reaction rate constants (the final ones calculated with P = 4) have been fit to the following forms:

method. We also plot our calculated rate constants and Pang et al.’s measured rate constants for non-R1b reactions in Figure 9 for direct comparison, and they are in good agreement.

Figure 9. Direct comparison with the experiment of Pang et al.

MP-CVT/SCT calculations and the global fits of McGillen et al. give very different branching fractions for the R1a and R1c reactions. In particular, the global fits, which are essentially SAR estimations at high temperature, predict that the dominant reaction channel is R1c above 1400 K, while MP-CVT/SCT predicts that R1a dominates for temperatures above 500 K. Note that the SAR method predicts the same rate for the R1c reaction if the OH group in isobutanol is replaced by a methyl group, i.e., if isobutanol is replaced by isopentane; however, there is a major difference between reaction with isobutanol and isopentane in that the hydrogen-bonding that is present for isobutanol is not present for isopentane. As shown in Figure 3, the high-energy transition structures without a hydrogen bond contribute more to the rate of reaction R1c than does the lowest-energy structure that has a hydrogen bond. The reason is that hydrogen-bonded structures tend to be more rigid and have less entropy. Because the classical barrier heights for transition structures with or without a hydrogen bond can be quite different, and because isopentane could have a quite different reaction rate for the R1c reaction than isobutanol, the SAR method based on the parametrization from alkane data might be inapplicable for estimating reaction rates for alcohols. The reaction rate difference between isobutanol and isopentane deserves further theoretical study for a better understanding of structure−reactivity relationships. Figure 8 shows that the hydrogen abstraction of R1b is not the dominant channel over the 200−2400 K range, although the classical barrier height of R1b is the lowest among the four reactions. Hydrogen abstraction from the OH site plays a minor role. One may expect that the R1c reaction becomes more important in the high-temperature region, and its branching fraction may cross over the branching fraction of R1b because its transition state has a much larger number of structures (this site has six hydrogen atoms to react). However, our calculated branching fraction for R1c is less than 14%. Although R1c has a much larger FMS‑T RXN factor than the other reactions, the effect of the large multi-structural torsional anharmonicity of R1c does not compete well with the anharmonicity of high-frequency

⎛ T + 30.79 ⎞2.001 ⎟ ka = 6.587 × 10−13⎜ × ⎝ 300 ⎠ ⎡ 0.9924(T + 30.79) ⎤ exp⎢ ⎥ ⎣ R(T 2 + 948.1) ⎦

(11)

⎛ T + 36.32 ⎞2.548 ⎟ k b = 6.090 × 10−14⎜ × ⎝ 300 ⎠ ⎡ 2.213(T + 36.32) ⎤ exp⎢ ⎥ ⎣ R(T 2 + 1319) ⎦

(12)

⎛ T + 29.75 ⎞2.364 ⎟ kc = 7.055 × 10−14⎜ × ⎝ 300 ⎠ ⎡ 1.088(T + 29.75) ⎤ exp⎢ ⎥ ⎣ R(T 2 + 885.1) ⎦

(13)

⎛ T + 430.6 ⎞2.706 ⎟ kd = 1.719 × 10−14⎜ × ⎝ 300 ⎠ ⎡ 1.488(T + 430.6) ⎤ exp⎢ − ⎥ ⎣ R(T 2 + 1.854 × 105) ⎦

(14)

Note that the individual fitting parameters in eqs 11−14 should not be over-interpreted. Fitting parameters in such expressions do not have a precise physical interpretation. The sum of eqs 11−14 agrees with eq 10 within 7%; the small difference is due to small least-squares fitting errors in each expression.

4. CONCLUDING REMARKS We present reaction rate calculations for hydrogen abstraction reactions from isobutanol by hydroxyl radical using the multipath and multi-structural canonical variational transition-state theory with small-curvature tunneling. This work was only made possible by combining several recent developments in VTST and the methods available for treating vibrational anharmonicity, although we still had to make an approximation that needs follow-up in the future, namely that corrections to our treatment of torsional coupling effects in the reaction of OH with isobutanol are approximately the same as for the reaction of OH with ethanol. Despite our recent work in making the treatment of torsional anharmonicity more reliable, it might still be the largest uncertainty. The present calculations include the accurate classical barrier height, the variational effect, the tunneling contribution evaluated with multiple reaction paths and internal-coordinate-based vibrationally adiabatic potentials, and the anharmonicity of torsional vibration modes and high-frequency modes. Some effects decrease the value of the reaction rate, and some effects increase it. We found that multi-structural torsional anharmonicity is a determining factor for branching ratios at medium and high temperatures. We demonstrated that the anharmonicities of high-frequency modes are different in reactants and transition states, and we provided a practical way to take account this difference by using different frequency scaling factors for reactants and transition states. The specific-reaction parametrized scaling factors were obtained by calculating the 5158



Article

(11) Sarathy, S. M.; Vranckx, S.; Yasunaga, K.; Mehl, M.; Wald, P. O.; Westbook, C. K.; Pitz, W. J.; Kohse-Höinghaus, K.; Fernandes, R. X.; Curran, H. J. Combust. Flame 2012, 159, 2028−2055. (12) Zheng, J.; Truhlar, D. G. J. Chem. Theory Comput. 2013, 9, 2875−2881. (13) Wigner, E. Trans. Faraday Soc. 1938, 34, 29−41. (14) Eyring, H. Trans. Faraday Soc. 1938, 34, 41−48. (15) Garrett, B. C.; Truhlar, D. G. J. Chem. Phys. 1979, 70, 1593− 1598. (16) Truhlar, D. G.; Garrett, B. C. Acc. Chem. Res. 1980, 13, 440− 448. (17) Truhlar, D. G.; Isaacson, A. D.; Skodje, R. T.; Garrett, B. C. J. Phys. Chem. 1982, 86, 2252−2261; 1983, 87, 4554 (Erratum). (18) Fernandez-Ramos, A.; Ellingson, B. A.; Garrett, B. C.; Truhlar, D. G. In Reviews in Computational Chemistry; Cundari, T. R., Lipkowitz, K. B., Eds.; Wiley-VCH: Hoboken, NJ, 2007; Vol. 23, p 125. (19) Zheng, J.; Yu, T.; Papajak, E.; Alecu, I. M.; Mielke, S. L.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2011, 13, 10885−10907. (20) Zheng, J.; Truhlar, D. G. J. Chem. Theory Comput. 2013, 9, 1356−1367. (21) Yu, T.; Zheng, J.; Truhlar, D. G. Chem. Sci. 2011, 2, 2199−2213. (22) Yu, T.; Zheng, J.; Truhlar, D. G. J. Phys. Chem. A 2012, 116, 297−308. (23) Zheng, J.; Truhlar, D. G. Faraday Discuss. 2012, 157, 59−88. (24) Zheng, J.; Seal, P.; Truhlar, D. G. Chem. Sci. 2012, 4, 200−212. (25) Pople, J. A.; Scott, A. P.; Wong, M. W.; Radom, L. Isr. J. Chem. 1993, 33, 345−350. (26) Scott, A. P.; Radom, L. J. Phys. Chem. 1996, 100, 16502−16513. (27) Alecu, I. M.; Zheng, J.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2010, 6, 2872−2887. (28) Bloino, J.; Biczysko, M.; Barone, V. J. Chem. Theory Comput. 2012, 8, 1015−1036. (29) Kuhler, K. M.; Truhlar, D. G.; Isaacson, A. D. J. Chem. Phys. 1996, 104, 4664−4670. (30) Nielsen, H. H. Encycl. Phys. 1959, 37, 173−313. (31) Mills, I. M. Molecular Spectroscopy: Modern Research; Academic: New York, 1972. (32) Zhang, Q.; Day, P. N.; Truhlar, D. G. J. Chem. Phys. 1993, 98, 4948−4958. (33) Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2008, 4, 1849−1868. (34) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. v. R. J. Comput. Chem. 1983, 4, 294−301. (35) Lynch, B. J.; Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2003, 107, 1384−1388. (36) Lynch, B. J.; Truhlar, D. G. J. Phys. Chem. A 2003, 107, 3898− 3906. (37) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479. (38) Adler, T. B.; Knizia, G.; Werner, H.-J. J. Chem. Phys. 2007, 127, 221106/1−4. (39) Knizia, G.; Adler, T. B.; Werner, H.-J. J. Chem. Phys. 2009, 130, 054104/1−20. (40) Papajak, E.; Truhlar, D. G. J. Chem. Theory Comput. 2011, 7, 10−18. (41) Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618. (42) Lynch, B. J.; Fast, P. L.; Harris, M.; Truhlar, D. G. J. Phys. Chem. A 2000, 104, 4811−4815. (43) 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.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.;

anharmonic ZPE using hybrid degeneracy-corrected secondorder vibrational perturbation theory. We found that the negative temperature dependence of the total reaction rate at lower temperatures is controlled by anharmonicity of highfrequency modes in the transition-state structures. The calculated total rate constants are in excellent agreement with all available experimental values at low temperatures, 250− 370 K, and at high temperatures, 900−1150 K. Thermal rate constants and branching ratios for all hydrogen abstraction sites are given from 200 to 2400 K. Despite the excellent agreement with experiment, to the extent that experimental results are available, a detailed analysis of all the various factors affecting the calculated rate constants and branching ratios shows that large uncertainties are still present in reaction rate calculations, even when using state-of-art methods, and more work needs to be done to estimate these uncertainties accurately.

■

ASSOCIATED CONTENT

S Supporting Information *

Cartesian coordinate optimized geometries of the reactants and transition states and the Gaussian input files for HDCVPT calculations (.txt). This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

[email protected] Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors are grateful to Bill Green, Ron Hanson, and Genny Pang for helpful correspondence. The work of D.G.T. was supported in part by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences (BES), under grant no. DE-FG02-86ER13579, and the work of J.Z. was supported by the Combustion Energy Frontier Research Center, an Energy Frontier Research Center funded by DOE, Office of Science, BES, under award no. DE-SC0001198.

■

REFERENCES

(1) Connor, M. R.; Liao, J. C. Curr. Opin. Biotechnol. 2009, 20, 307− 315. (2) Atsumi, S.; Wu, T.-Y.; Eckl, E.-M.; Hawkins, S. D.; Buelter, T.; Liao, J. C. Appl. Microbiol. Biotechnol. 2010, 85, 651−657. (3) Peralta-Yahya, P. P.; Zhang, F.; Del Cardayre, S. B.; Keasling, J. D. Nature 2012, 488, 320−328. (4) Welz, O.; Savee, J. D.; Eskola, A. J.; Sheps, L.; Osborn, D. L.; Taatjes, C. A. Proc. Combust. Inst. 2013, 34, 493−500. (5) Wu, H.; Mu, Y.; Zhang, X.; Jiang, G. Int. J. Chem. Kinet. 2003, 25, 81−87. (6) Mellouki, A.; Oussar, F.; Lun, X.; Chakir, A. Phys. Chem. Chem. Phys. 2004, 6, 2951−2955. (7) Andersen, V. F.; Wallington, T. J.; Nielson, O. J. J. Phys. Chem. A 2010, 114, 12462−12469. (8) Pang, G. A.; Hanson, R. K.; Golden, D. M.; Bowman, G. T. J. Phys. Chem. A 2012, 116, 4720−4725. (9) McGillen, M. R.; Baasandorj, M.; Burkholder, J. B. J. Phys. Chem. A 2013, 117, 4636−4656. (10) Grana, R.; Frassoldati, A.; Faravelli, T.; Niemann, U.; Ranzi, E.; Seiser, R.; Cattolica, R.; Seshadri, K. Combust. Flame 2010, 157, 2137− 2154. 5159



Article

Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (44) Werner, H.-J.; Knowles, P. J.; Manby, F. R.; M. Schütz, Celani, P.; Knizia, G.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; C. Köppl, Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklaß, A.; Palmieri, P.; K. Pflüger, Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M.; Wolf, A. MOLPRO, 2010.1; University of Birmingham: Birmingham, UK, 2010. (45) Zheng, J.; Mielke, S. L.; Clarkson, K. L.; Truhlar, D. G. Comput. Phys. Commun. 2012, 183, 1803−1812. (46) Kilpatrick, J. E.; Pitzer, K. S. J. Chem. Phys. 1949, 17, 1064− 1075. (47) Lu, D.-h.; Truong, T. N.; Melissas, V. S.; Lynch, G. C.; Liu, Y.P.; Garrett, B. C.; Steckler, R.; Isaacson, A. D.; Rai, S. N.; Hancock, G. C.; Lauderdale, J. G.; Joseph, T.; Truhlar, D. G. Comput. Phys. Commum. 1992, 71, 235−262. (48) Liu, Y.-P.; Lynch, G. C.; Truong, T. N.; Lu, D.-h.; Truhlar, D. G.; Garrett, B. C. J. Am. Chem. Soc. 1993, 115, 2408−2415. (49) Jackels, C. F.; Gu, Z.; Truhlar, D. G. J. Chem. Phys. 1995, 201, 3188−3201. (50) Nguyen, K. A.; Jackels, C. F.; Truhlar, D. G. J. Chem. Phys. 1996, 104, 6491−6496. (51) Zheng, J.; Mielke, S. L.; Clarkson, K. L.; Meana-Pañeda, R.; Truhlar, D. G. MSTor, version 2013; University of Minnesota: Minneapolis, MN, 2013. (52) Zheng, J.; Meana-Pañeda, R.; Truhlar, D. G. Comput. Phys. Commun. 2013, 184, 2032−2033. (53) Zheng, J.; Zhang, S.; Lynch, B. J.; Corchado, J. C.; Chuang, Y.Y.; Fast, P. L.; Hu, W.-P.; Liu, Y.-P.; Lynch, G. C.; Nguyen, K. A.; Jackels, C. F.; Fernandez-Ramos, A.; Ellingson, B. A.; Melissas, V. S.; Villa, J.; Rossi, I.; Coitino, L.; Pu, J.; Albu, T. V.; Steckler, R.; Garrett, B. C.; Issacson, A. D.; Truhlar, D. G. POLYRATE; University of Minnesota: Minneapolis, 2013; locally modified version based on version 2010-A. (54) Zheng, J.; Zhang, S.; Corchado, J. C.; Chuang, Y.-Y.; Coitiño, E. L.; Ellingson, B. A.; Truhlar, D. G. GAUSSRATE, version 2009-A; University of Minnesota: Minneapolis, 2009. (55) Black, G.; Simmie, J. M. J. Comput. Chem. 2010, 31, 1236−1248. (56) Zheng, J.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2007, 3, 569−582. (57) Zheng, J.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2009, 5, 808−821. (58) Garrett, B. C.; Truhlar, D. G. J. Am. Chem. Soc. 1979, 101, 4534−4548. (59) Atkinson, R. Int. J. Chem. Kinet. 1986, 18, 555−568. (60) Atkinson, R. Int. J. Chem. Kinet. 1987, 19, 799−828. (61) Kwok, E. S. C.; Atkinson, R. Atmos. Environ. 1995, 29, 1685− 1695. (62) Bethel, H. L.; Atkinson, R.; Arey, J. Int. J. Chem. Kinet. 2001, 33, 310−316. (63) Zheng, J.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2010, 12, 7782−7793.

5160


Prediction of Experimentally Unavailable Product Branching Ratios for

Recommend Documents