Junwei Lucas
Bao
,
Rubén
Meana-Pañeda
and
Donald G.
Truhlar
*
Department of Chemistry, Chemical Theory Center and Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455-043, USA. E-mail: truhlar@umn.edu
First published on 16th June 2015
The goal of the present work is modeling the kinetics of a key reaction involved in the combustion of the biofuel 2-butanol. To accomplish this we extended multi-path variational transition state theory (MP-VTST) with the small curvature tunneling (SCT) approximation to include multistructural anharmonicity factors for molecules with chiral carbons. We use the resulting theory to predict the site-dependent rate constants of the hydrogen abstraction from 2-butanol by hydroperoxyl radical. The generalized transmission coefficients were averaged over the four lowest-energy reaction paths. The computed forward reaction rate constants indicate that hydrogen abstraction from the C-2 site has the largest contribution to the overall reaction from 200 K to 2400 K, with a contribution ranging from 99.9988% at 200 K to 88.9% at 800 K to 21.2% at 3000 K, while hydrogen abstraction from the oxygen site makes the lowest contribution at all temperatures, ranging from 2.5 × 10^{−9}% at 200 K to 0.65% at 800 K to 18% at 3000 K. This work highlights the importance of including the multiple-structure and torsional potential anharmonicity in the computation of the thermal rate constants. We also analyzed the role played by the hydrogen bond at the transition state, and we illustrated the risks of (a) considering only the lowest-energy conformations in the calculations of the rate constants or (b) ignoring the nonlinear temperature dependence of the activation energies. A hydrogen bond at the transition state can lower the enthalpy of activation, but raise the free energy of activation. We find an energy of activation that increases from 11 kcal mol^{−1} at 200 K to more than 36 kcal mol^{−1} at high temperature for this radical reaction with a biofuel molecule.
Hydroperoxyl radical (HO_{2}˙) plays a particularly important role in the ignition process at intermediate temperature (∼800–1200 K), partly because another important radical, HO˙, can be formed by the decomposition of H_{2}O_{2}, which is the product of the hydrogen abstraction reaction by hydroperoxyl radical. In the present study, we calculated the rate constants of hydrogen abstraction from different sites of 2-butanol by HO_{2}˙ radical using the recently developed multi-path variational transition state theory (MP-VTST)^{16,17} with the small-curvature tunneling (SCT) approximation,^{18} with the MP-VTST treatment being extended here to the treatment of a chiral reactant. Multi-structural variational transition state theory (MS-VTST) calculations, which have now been widely applied in previous work, have also been carried out for comparison. MP-VTST can be viewed as an extension of MS-VTST; first it accounts for the overbarrier contribution of all the reaction paths to the total one-way flux in a more accurate way, and secondly it treats the differences in the multi-dimensional tunneling paths and their contribution to the thermal rate constants more thoroughly.
The five possible hydrogen abstraction reactions that we considered in this work are the abstraction of a hydrogen atom from each of the four carbon sites of 2-butanol or from the oxygen atom. Carbon-2 of the 2-butanol molecule is a chiral center; therefore 2-butanol has two enantiomers: rectus (R) and sinister (S), i.e. (R)-2-butanol and (S)-2-butanol. When 2-butanol reacts with a nonchiral species, which in the present study is the HO_{2} radical, the thermal rate constants of the R enantiomer and the S enantiomer are the same. Therefore the rate constants of the studied reactions (R1)–(R5) can be computed based on a single enantiomer. In the present work, we use (S)-2-butanol to compute the rate constants. Therefore the reactions we consider are
S-CH_{3}CH(OH)CH_{2}CH_{3} + HO_{2}˙ → S-CH_{2}·CH(OH)CH_{2}CH_{3} + H_{2}O_{2} | (R1) |
S-CH_{3}CH(OH)CH_{2}CH_{3} + HO_{2}˙ → CH_{3}C·(OH)CH_{2}CH_{3} + H_{2}O_{2} | (R2) |
S-CH_{3}CH(OH)CH_{2}CH_{3} + HO_{2}˙ → S-CH_{3}CH(OH)C(H˙)CH_{3} + H_{2}O_{2} | (R3) |
S-CH_{3}CH(OH)CH_{2}CH_{3} + HO_{2}˙ → S-CH_{3}CH(OH)CH_{2}CH_{2}˙ + H_{2}O_{2} | (R4) |
S-CH_{3}CH(OH)CH_{2}CH_{3} + HO_{2}˙ → S-CH_{3}CH(O˙)CH_{2}CH_{3} + H_{2}O_{2} | (R5) |
The present calculations not only yield the total rate constant, which is the sum of the five elementary rate constants for reactions (R1)–(R5), but they also yield the five individual elementary rate constants; these site-specific quantities are seldom available experimentally, but they are very important for combustion mechanisms because each of the five products has a different further reactivity pattern after being produced.
Section 2 presents the Theory, and Section 3 has Computational details. Readers interested only in the results may skip to Section 4, which presents Results and discussion. To facilitate this, Section 4 reminds the reader of the definitions of all acronyms and mathematical symbols used in that section. Section 5 summarizes the main conclusions.
Each species can have several structures. For reactants and products, the structures are local minima on the potential energy surface and may be called conformers. For conventional transition states, the structures are saddle points with one imaginary frequency and they may be called transition structures. Only saddle points that are maxima on a minimum energy path are considered; that means that the fluxes through transition structures of a given transition state must be added like currents in parallel, not in series.^{17,19} (One could consider transition states in series by the canonical unified statistical model,^{20} but that is not the subject of the present study.) For generalized transition states, the structures are the points where the reaction paths pass through the dividing surface, and they may be called generalized transition structures.
For bimolecular reactions, we consider that “the partition function of the reactants” is the product of the partition functions for the two reactants, and “the energy of the reactants” is the sum of their energies. The electronic energy of a structure including nuclear repulsion is the potential energy for nuclear motion. The barrier height V^{‡}_{f} is the increase in potential energy from the lowest-energy equilibrium structure of the reactants to the lowest-energy conventional transition state. Adding zero point energy at these structures give the conventional-transition-state-theory enthalpy of activation at 0 K, which is denoted ΔH^{‡}_{0}. The energy of reaction ΔE is the potential energy of the lowest-energy equilibrium structure of the products minus that of the lowest-energy equilibrium structure of the reactants. Adding zero point energy at these structures gives the enthalpy of reaction at 0 K, which is denoted ΔH_{0}.
We are interested in species that have (in the general case) more than one structure, and these structures may be called conformers. We must replace the rovibrational partition function that is appropriate for a species with a single structure by a conformational–rovibrational partition function for species that have multiple structures. In many cases the conformers are connected by torsions; the local maxima along the torsional coordinates are not barriers along reaction paths connecting reactants to products and so they are not treated as transition structures or generalized transition structures in considering the kinetics; however, the heights of these barriers do have an effect on the conformational–rovibrational partition function when the barrier regions of configuration space are thermally accessible at the temperature of interest. The “classical barrier height” of the chemical reaction is the increase in potential energy from the lowest-energy equilibrium structure of the reactants to the lowest-energy conventional transition state. Adding zero point energy at these structures give the conventional-transition-state-theory enthalpy of activation at 0 K, which is denoted ΔH^{‡}_{0}.
With these concepts and terminologies as background, we can proceed to discuss the theory.
(1) |
(2) |
A key feature of the MS-T method is that it does not require the calculation of the barrier heights along the reaction paths for the interconversion of conformers; rather these barriers are implicitly approximated in terms of the local periodicities by using a Voronoi tessellation in the available torsional phase space,^{21–23} which are much easier to compute. This is a generalization of the relation valid for symmetric one-dimensional torsions by which the barrier between equivalent minima can be predicted approximately from the local frequency and the periodicity. Due to the coupling between the torsions, the number of conformational structures that one finds after the optimizations is less than the number of structures generated by ideal torsions, i.e., less than the ideal number of 3^{t} for t threefold torsions. Therefore, the periodicity is not well defined, but we replace it by a local periodicity as explained previously^{22,23} and extended to chiral reagents in Section 2.5.
The effect of multistructural torsional anharmonicity on the partition function of a species X can be characterized by the following factor
(3) |
To demonstrate the multistructural torsional anharmonic effect on the reaction rate, we define the following multistructural torsional anharmonicity factor
(4) |
k^{MS-CVT/SCT} = F^{MS-T}_{act}k^{SS-CVT/SCT}_{1} | (5) |
The variational transition state location needed for eqn (5) is obtained by maximizing the quasiclassical standard-state generalized free energy of activation ΔG^{GT,o}_{act} as a function of the reaction coordinate s, which is defined as the distance from the saddle point along the minimum energy path (MEP). For a bimolecular reaction without treating torsional anharmonicity as a function of reaction coordinate, we have
(6) |
k^{MP-CVT/SCT} = F^{MS-T}_{act}〈γ〉k^{TST}_{1} | (7) |
(8) |
In the MP-VTST method, the recrossing transmission coefficient of the reaction path p without treating torsional anharmonicity as a function of reaction coordinate s can be written as follows:
(9) |
In contrast, in MP(full)-VTST, we include the s dependence of the torsional anharmonicity, and so we obtain:^{24}
(10) |
Note that the location of the generalized transition state may differ (and usually will differ) from the one determined without treating torsional anharmonicity as a function of reaction coordinate, due to the change of the location of the maximum of the free energy of activation, which is now determined now by finding the maximum of
(11) |
Fig. 1 Schematic representation of the reactions and transition structures for the hydrogen abstraction reaction from 2-butanol by HO_{2} radical. |
The conformational–rotational–vibrational partition function of the transition state of (R3) includes both diastereomers, (2S,3R) and (2S,3S), (where the zero of energy is set at the lowest-energy structure, see eqn (1)), and their corresponding enantiomers, (2R,3S) and (2R,3R) are not included in the calculation because they are formed from (R)-2-butanol (whereas we computed the reaction rates only for (S)-2-butanol). For the other forward reactions (i.e., (R1), (R2), (R4) and (R5)), only the transition structures that can be generated from the S enantiomer of 2-butanol, should be included in the calculation of the partition function. Next we consider the computation of the local periodicities for all these cases.
The local periodicity M_{j,η} needs to be determined in order to compute the potential anharmonicity factor _{j,η} in eqn (1). When the torsions are strongly coupled, the Voronoi tessellation method^{23} is employed to assign the local periodicities. First we define the torsional space as a t-dimensional space described by the torsional coordinates; where the torsional coordinates are taken to be a set of dihedral angles φ_{1}, φ_{2},…, φ_{t}. The Voronoi method divides the torsional space into cells around a discrete set of points that correspond to structures. All structures used in Voronoi tessellation should be interconvertable by internal rotation, i.e. the entire torsional space should be accessible by internal rotation. Each geometry belongs to a specific cell, and each cell contains all the geometries with values of torsional angles that are closer to the primary structure, located at the center of the cell, than to any other structure. A primary structure is a conformer of a reactant, product, first-order saddle point (conventional transition state), or a point on the reaction path (generalized transition state).
As explained in the original MS-T paper,^{22} the Voronoi calculation must include both distinguishable and indistinguishable primary structures that can be generated from a single structure by internal rotation of one or more torsional coordinates (that is, dihedral angles). If the interconversion between distinguishable structures is impossible by internal rotation, these structures are not in the same torsional space and they are considered separately via independent Voronoi calculations. It is important to analyze the torsional space available for each structure before applying the Voronoi tessellation method to obtain the local periodicities.
We illustrate this by an example in Fig. 2. This figure shows the hydrogen abstraction reactions of butane, which is a nonchiral reactant, by atomic oxygen. A and A′ are the two enantiomeric transition structures formed during the hydrogen abstraction from C-2; B and B′ are obtained by rotating the torsional bonds of A and A′ respectively. Note that B and B′ cannot be superimposed; therefore, the effective periodicities of B and B′ should be computed separately, i.e. B′ should not be included in the Voronoi calculation for the effective periodicities of B. Generally speaking, the enantiomers of a molecule with a chiral carbon cannot be interconverted by torsions and when one computes the effective periodicities for conformers that belong to one specific configuration (R or S), the corresponding enantiomers should not be included in the Voronoi calculation. For example, C and C′ are the two enantiomeric transition structures formed during the hydrogen abstraction from C-1. These two structures should be included in the same Voronoi calculation because they can be interconverted by internal rotations (D and D′, which are obtained by rotating the torsional bonds of C and C′ respectively, are identical).
Fig. 2 A schematic illustration of the interconversion between enantiomers as illustrated by the simple case of butane reacting with oxygen atom. |
For a chiral species (chiral reactant, chiral transition state in the case where there are chiral reactants, or chiral transition state in the case where the reactants are not chiral), two kinds of the structures are used in the Voronoi calculations of the effective periodicities, namely indistinguishable primary structures and distinguishable primary structures. However, this does not mean we should include all the distinguishable primary structures in the Voronoi calculation. If the interconversion between some distinguishable structures and the structures we are considering is impossible by internal rotation, these structures are not in the torsional space we are considering and hence should not be included in the Voronoi calculation.
In the present study, since we only consider the transition states and products that can be formed from (S)-2-butanol in the Voronoi calculations of all the species except for P2 (the product radical of reaction (R2), which does not have a chiral carbon); the corresponding enantiomers are not included. For the diastereomeric transition structures of C-3 in the present case, i.e. (2S,3R) and (2S,3S), we compute the effective periodicities of these two classes of transition structures separately via the Voronoi tessellation method. For all species, the enantiomers did not need to be included in the conformational rovibrational partition function either because they are different molecules or because we here consider only the structures that can be formed from (S)-2-butanol. A schematic summary of the relations between various species for reaction (R3), along with an illustration, is provided in Fig. S1 in the ESI.†
In the present study, all the torsions except the methyl groups are assumed to be strongly coupled, and their effective periodicity values are assigned by the Voronoi method. Methyl groups are considered as nearly separable and the local periodicity value for internal rotation of each methyl group is 3.
Benchmark calculations were carried out in order to test and validate various density functionals for dynamic calculations. Various combinations of M08-HX, M08-SO,^{27} M06-2X^{35} density functionals with MG3S, jun-cc-pVTZ,^{36} jul-cc-pVTZ,^{37} maug-cc-pVTZ^{29} and aug-cc-pVTZ^{38} basis sets were used to compute the single-point energies based on M08-HX/MG3S geometries. Single-point energies calculated by the explicitly correlated electronic structure method CCSD(T)-F12a^{39,40} were chosen as the benchmark values. Three different basis sets, namely the orbital (AO) basis, the density fitting (DF) basis^{41} and the resolution of the identity (RI) basis,^{42} are required in F12 calculations. The basis set jun-cc-pVTZ was used as the orbital basis for constructing the molecular orbitals; the F12 integrals were evaluated using density fitting approximations with cc-pVTZ/MP2FIT basis;^{43} and the cc-pVTZ/JKFIT basis^{42} was chosen as the RI basis for approximating the many-electron integrals. All CCSD(T)-F12a calculations were performed by using the Molpro^{44} program.
The reaction paths were computed via direct dynamics and all the dynamics calculations were performed with the POLYRATE^{45} and GAUSSRATE^{46} programs. The electronic structure method, M08-HX/aug-cc-pVTZ, with the smallest mean unsigned error (MUE) in the validations discussed in section 4.3 was employed in all the dynamic calculations. For the dynamics calculations, all stationary-point geometries were re-optimized by M08-HX/aug-cc-pVTZ, and with this basis set a factor of 0.975 was used to scale the vibrational frequencies. The generalized normal mode analyses were performed using non-redundant internal coordinates. The tunneling transmission probabilities are computed using the small-curvature tunneling approximation. In the direct dynamics calculations, the minimum energy paths (MEPs) in isoinertial coordinates were computed by using the Euler steepest-descents (ESD) method with a step size of 0.005 bohr and the Hessian was updated at every 9th point along the path.
The lowest-energy structures for (S)-2-butanol, for the organic products, and for conventional transition states are shown in Fig. 3. Nine distinguishable structures were found for (S)-2-butanol. For products P1–P5, the numbers of distinguishable structures are 9, 12 (6 pairs of mirror images), 7, 9 and 3 respectively. Transition states TS1, TS2, TS3 (2S,3R), TS3 (2S,3S), TS4 and TS5 have 117, 35, 42, 39, 109 and 44 structures respectively. The numbers of distinguishable structures in each relative potential energy range for transition states TS1–TS5 are presented in Fig. 4. For each transition state, the zero of relative energy for the transition state structures is chosen as the energy of the global-minimum structure.
Fig. 3 Lowest-energy structures of (S)-2-butanol, products P1–P5 and transition structures TS1–TS5. Carbon, hydrogen and oxygen atoms are represented by grey, white, and red balls respectively. |
Fig. 4 Number of distinguishable structures in each relative potential energy range for transition states TS1–TS5. |
Fig. 5 Multistructural torsional anharmonicity factors F^{MS-T}_{X} for: (top) (S)-2-butanol, products P1–P5; (bottom) transition states TS1–TS5. |
Fig. 6 Multistructural torsional anharmonicity factors F^{MS-T}_{act} for: (top) Forward reactions (R1)–(R5); (bottom) Reverse reactions (R1)–(R5). |
The F^{MS-T}_{X} factor is not directly proportional to the number of distinguishable conformational structures, and it is interesting to notice that a large number of conformers for a species does not necessarily lead to a large value of F^{MS-T}_{X} factor. For instance, TS1 has 117 distinguishable conformers while TS2 only has 35 conformers, but their F^{MS-T}_{X} factors are very close to each other. From Table S1,† we see that TS1 has larger Q^{MS-T} than TS2 but smaller values of the rovibrational partition functions of the global minimum, which is caused by the higher values of the high-frequency vibrational modes.
Electronic model chemistry//M08-HX/MG3S | Forward barrier height V^{‡}_{f} | Reverse barrier height V^{‡}_{r} | Energy of reaction ΔE (kcal mol^{−1}) | MUE | MUE (w/o (R5)) | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
(R1) | (R2) | (R3) | (R4) | (R5) | (R1) | (R2) | (R3) | (R4) | (R5) | (R1) | (R2) | (R3) | (R4) | (R5) | |||
a For consistent comparisons, the same set of geometries (those obtained by M08-HX/MG3S) is used throughout this table. The lowest-energy structures for reactants, conventional transition states, and products are employed. Mean unsigned errors (MUEs) are calculated with respect to the CCSD(T)-F12a/jun-cc-pVTZ method. Barrier heights and energies of reaction are classical, i.e., zero-point-energy exclusive. | |||||||||||||||||
M08-HX/aug-cc-pVTZ | 20.89 | 12.09 | 17.66 | 17.61 | 22.07 | 4.46 | 4.64 | 4.37 | 2.83 | 2.23 | 16.44 | 7.45 | 13.28 | 14.78 | 19.84 | 0.22 | 0.21 |
M08-SO/aug-cc-pVTZ | 21.92 | 13.10 | 18.56 | 18.39 | 21.94 | 5.06 | 5.19 | 5.35 | 3.00 | 1.41 | 16.86 | 7.90 | 13.20 | 15.39 | 20.52 | 0.66 | 0.70 |
M06-2X/aug-cc-pVTZ | 19.02 | 10.61 | 15.87 | 15.71 | 19.83 | 3.11 | 3.36 | 3.20 | 1.40 | 0.43 | 15.90 | 7.25 | 12.67 | 14.31 | 19.40 | 1.14 | 1.11 |
M08-HX/maug-cc-pVTZ | 21.06 | 12.34 | 17.84 | 17.94 | 22.29 | 4.61 | 4.81 | 4.50 | 3.09 | 2.38 | 16.45 | 7.52 | 13.34 | 14.85 | 19.90 | 0.29 | 0.26 |
M08-HX/MG3S | 20.63 | 11.61 | 17.18 | 17.31 | 21.68 | 3.89 | 3.94 | 3.55 | 2.08 | 1.64 | 16.74 | 7.67 | 13.62 | 15.23 | 20.03 | 0.50 | 0.58 |
CCSD(T)-F12a/jun-cc-pVTZ | 20.32 | 12.38 | 17.57 | 17.47 | 21.72 | 4.52 | 4.76 | 4.43 | 2.91 | 1.90 | 15.80 | 7.62 | 13.14 | 14.56 | 19.82 | 0.00 | 0.00 |
The transition state of the oxygen-site abstraction reaction (R5) is a multi-reference system, and all the calculations carried out here are based on single-reference methods. The T1 diagnostic^{47} value of the lowest-energy structure of TS5 is 0.0928, which is considerably larger than the value of 0.045 that has been suggested as a criterion for multi-reference character in an open-shell molecule.^{48–52} However, since this reaction has the largest forward barrier height, it has the least contribution to the total reaction rate constant. Therefore, in this study, we are satisfied to use single-reference methods to estimate the barrier heights and energy of reaction for reaction (R5).
Mean unsigned errors (MUEs) are calculated with respect to the results obtained by the F12a explicitly correlated CCSD(T) electronic structure theory and are shown in Table 1 both including and excluding the least important reaction (R5). Based on these comparisons, we chose the M08-HX exchange–correlation functional with the aug-cc-pVTZ basis set to carry out the dynamics calculations.
T/K | (R1) | (R2) | (R3) | (R4) | (R5) | ||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
〈γ〉_{1} | 〈γ〉_{3} | 〈γ〉_{4} | URE^{b} (%) | 〈γ〉_{1} | 〈γ〉_{3} | 〈γ〉_{4} | URE (%) | 〈γ〉_{1} | 〈γ〉_{3} | 〈γ〉_{4} | URE (%) | 〈γ〉_{1} | 〈γ〉_{3} | 〈γ〉_{4} | URE (%) | 〈γ〉_{1} | 〈γ〉_{3} | 〈γ〉_{4} | URE (%) | ||
a Generalized transmission coefficients averaged over the first p lowest-energy path(s) are denoted as 〈γ〉_{p}. b The unsigned relative error (URE) is given by |〈γ〉_{4} − 〈γ〉_{3}|/〈γ〉_{4}. | |||||||||||||||||||||
200 | 196.5 | 231.8 | 243.5 | 4.8 | 20.84 | 26.94 | 26.88 | 0.2 | 207.0 | 214.7 | 218.6 | 1.8 | 45.43 | 59.42 | 67.68 | 12.2 | 1.42 | 1.78 | 1.95 | 8.4 | |
300 | 17.04 | 18.37 | 18.95 | 3.0 | 3.74 | 4.12 | 4.08 | 1.0 | 17.09 | 17.20 | 17.32 | 0.7 | 7.99 | 9.56 | 10.55 | 9.3 | 1.18 | 1.35 | 1.44 | 5.8 | |
400 | 5.79 | 6.00 | 6.13 | 2.2 | 2.02 | 2.12 | 2.11 | 0.8 | 5.62 | 5.63 | 5.66 | 0.4 | 3.62 | 4.08 | 4.37 | 6.6 | 1.10 | 1.20 | 1.25 | 4.2 | |
600 | 2.33 | 2.35 | 2.38 | 1.5 | 1.27 | 1.32 | 1.31 | 0.9 | 2.19 | 2.21 | 2.23 | 0.6 | 1.80 | 1.91 | 1.98 | 3.6 | 1.04 | 1.08 | 1.10 | 2.4 | |
1000 | 1.36 | 1.35 | 1.37 | 1.2 | 0.73 | 0.90 | 0.87 | 3.3 | 1.23 | 1.27 | 1.29 | 1.1 | 1.16 | 1.17 | 1.18 | 1.5 | 1.01 | 1.01 | 1.02 | 1.5 | |
1500 | 1.13 | 1.11 | 1.13 | 1.4 | 0.54 | 0.69 | 0.66 | 3.6 | 0.98 | 1.04 | 1.06 | 1.8 | 0.96 | 0.95 | 0.95 | 0.9 | 1.00 | 0.98 | 0.99 | 1.3 | |
2400 | 1.03 | 0.98 | 1.01 | 2.4 | 0.42 | 0.54 | 0.53 | 3.3 | 0.79 | 0.89 | 0.92 | 3.3 | 0.80 | 0.79 | 0.78 | 0.6 | 1.00 | 0.96 | 0.97 | 1.4 |
The reaction we considered for this test is reaction (R2), which has the largest contribution to the total reaction rate constant from 200 K to 2400 K; and the multistructural CVT/SCT rate constants are computed based on reaction path 1 (the lowest-energy path). As shown in Fig. 8, the differences between the MS-CVT/SCT and MS(full)-CVT/SCT rate constants are negligible at all temperatures (from 200 K to 3000 K). This result can be understood as follows. For path 1 of reaction (R2), the locations of the variational transition states determined by eqn (6) are −0.021 Å at 200 K, −0.022 Å at 298 K, and −0.21 Å at 1000 K respectively; and the locations determined by eqn (11) at these temperatures are −0.023 Å, −0.026 Å, and −0.22 Å respectively. We can see that the locations are all close to the saddle point (which is defined as s = 0), and the key point is that the multistructural anharmonicity does not differ much at the variational transition state and at the conventional TS. For other reaction paths, we also found that the locations of the variational transition state are all close to the saddle points (see Table S9 in ESI†). Therefore, including the torsional anharmonicity only at the saddle point is a good approximation.
Fig. 8 MS-CVT/SCT (denoted as “w/o Torsion” in the figure) and MS(full)–CVT/SCT (denoted as “with Torsion” in the figure) rate constant of reaction (R2), path1 computed with/without treating torsional anharmonicity as a function of reaction coordinate s at various temperatures. |
Thus, for the rate constants reported in the rest of the paper, the torsional anharmonicity is only included at the stationary points.
T/K | (R1) | (R2) | (R3) | (R4) | (R5) | Overall |
---|---|---|---|---|---|---|
200 | 2.00 × 10^{−33} | 1.76 × 10^{−25} | 4.09 × 10^{−30} | 8.42 × 10^{−32} | 4.42 × 10^{−36} | 1.76 × 10^{−25} |
300 | 1.66 × 10^{−27} | 2.78 × 10^{−22} | 2.54 × 10^{−25} | 9.68 × 10^{−27} | 5.93 × 10^{−29} | 2.78 × 10^{−22} |
400 | 2.48 × 10^{−24} | 2.13 × 10^{−20} | 1.09 × 10^{−22} | 6.72 × 10^{−24} | 3.28 × 10^{−25} | 2.14 × 10^{−20} |
600 | 8.21 × 10^{−21} | 3.44 × 10^{−18} | 1.06 × 10^{−19} | 1.47 × 10^{−20} | 3.25 × 10^{−21} | 3.57 × 10^{−18} |
1000 | 1.72 × 10^{−17} | 4.68 × 10^{−16} | 7.99 × 10^{−17} | 2.36 × 10^{−17} | 1.19 × 10^{−17} | 6.01 × 10^{−16} |
1500 | 1.62 × 10^{−15} | 9.25 × 10^{−15} | 4.34 × 10^{−15} | 1.79 × 10^{−15} | 1.32 × 10^{−15} | 1.83 × 10^{−14} |
2400 | 9.07 × 10^{−14} | 1.45 × 10^{−13} | 1.49 × 10^{−13} | 7.57 × 10^{−14} | 7.87 × 10^{−14} | 5.40 × 10^{−13} |
The product of the reverse reaction of (R2), P2, has no chiral carbon. However, the H_{2}O_{2} molecule can attack P2 from two different sides (top and bottom) and this leads to two possible products, (R)-2-butanol and (S)-2-butanol. In order to be consistent with the reactions we considered in this work as listed in the introduction, the reverse rate constants reported in Table S3† are based solely on (S)-2-butanol. If one wants to compute the reverse rate constants of (R2), which leads to both (R)- and (S)-2-butanol, one should multiply the rate constant we reported here by a factor of 2.
The MP-VTST rate constants for the forward reactions are fitted by the following equation^{54}
(12) |
(13) |
Forward reactions^{b} | Reverse reactions^{c} | |||||||
---|---|---|---|---|---|---|---|---|
lnA | n | T _{0} | E | lnA | n | T _{0} | E | |
a Rate constants are in the unit of cm^{3} molecule^{−1} s^{−1}, the parameters T_{0} and E are in the unit of K, kcal mol^{−1} respectively. b Forward reaction rate constants are fitted using eqn (12). c Reverse reaction rate constants are fitted using eqn (13). | ||||||||
(R1) | −39.1738 | 5.5754 | 104.563 | 11.230 | −38.0230 | 3.8710 | 527.398 | 3.305 |
(R2) | −36.4351 | 4.0045 | 139.974 | 6.581 | −39.8651 | 3.9670 | 303.680 | 2.459 |
(R3) | −38.0966 | 5.0931 | 117.587 | 9.240 | −37.9467 | 3.7753 | 490.295 | 3.544 |
(R4) | −38.3831 | 5.0739 | 138.408 | 10.766 | −31.7528 | 1.4836 | 576.178 | 7.653 |
(R5) | −37.3145 | 4.8844 | 91.719 | 13.875 | −34.0876 | 3.3045 | 483.334 | 2.014 |
Overall | −38.8559 | 5.6479 | 120.623 | 5.351 | −34.1251 | 3.3873 | 493.788 | 2.101 |
Fig. 9 Relative potential energy (in kcal mol^{−1}) of all the optimized structures for (S)-2-butanol and of transition structures of the hydrogen abstraction reactions (R1)–(R5). Green triangles correspond to structures that do not have hydrogen bond; blue squares are structures with a strongly bent hydrogen bond, and red circles are structures with a normal hydrogen bond, as defined in Section 4.8. |
The thermal rate constant is related to the quasiclassical standard-state generalized free energy of activation at the variational transition state (ΔG^{CVT,o}_{act}): the higher ΔG^{CVT,o}_{act} is, the lower the rate constant. Because we've shown that the positions of GT are close to the saddle points in the cases studied here, the quasiclassical standard-state generalized free energy of activation is approximately the quasiclassical free energy difference between the saddle point and the reactants. Fig. 10 shows the computed quasiclassical standard-state free energies, enthalpies and entropies at 800 K and 2000 K for all the transition structures of (R2).
Fig. 10 Relative energy (in kcal mol^{−1}), standard enthalpy (in kcal mol^{−1}), entropy (in cal K^{−1} mol^{−1}) and quasiclassical free energy (in kcal mol^{−1}) for all the transition structures of reaction (R2) at 800 and 2000 K. Red circles correspond to hydrogen–bonded structures. |
Fig. 10 shows that an H-bonded transition structure with lower energy or enthalpy (for instance, structure 4 in Fig. 10) is not necessarily lower in free energy than non-H-bonded structures; this is because the presence of the hydrogen bond can also reduce the entropy. The H-bonded structure 34 is energetically lower than all of the non-H-bonded structures, however its free energy is even higher than most of the non-H-bonded conformers, and therefore it contributes significantly less than other structures to the thermal rate constants. Therefore, we conclude here that a hydrogen-bonded transition structure does not necessarily contribute significantly to the thermal rate constants.
Many authors (we give only a couple of examples^{57,58}) have emphasized how a hydrogen bond can stabilize transition structures, and sometimes it is assumed that the stronger the interaction is, the lower the barrier and the faster the reaction. However, the results presented here show that basing such conclusions solely based on the energetic effects of the hydrogen bond not the entropic effects is not reliable. We have shown a hydrogen bond can reduce the entropy and thereby increase the free energy of the transition state. Hence, a strong hydrogen bond interaction may lead to a slower reaction rate; reliable conclusions must be based on free energies of activation, not barrier heights or enthalpies of activation.
(14) |
Fig. 11 Activation energies for the forward reactions (R1)–(R5) and the overall reaction. |
The temperature dependence of the activation energy is sometimes ignored in the empirical estimations of the rate constants,^{61} which is a convenient and reasonable choice for modeling a small range of temperatures, but this may cause significant errors over a wide range of temperatures, especially if one adopts the activation energy measured at a lower temperature to estimate the rate constant at very high temperatures in combustion chemistry or if one uses higher-temperature measurements of rate constants (where they are large enough to be measured) to estimate rate constants at 200–300 K for atmospheric applications.
The formula usually used for fitting rate constants when the temperature dependence of the activation energy is taken into account is:
(15) |
E_{a} = E + nRT | (16) |
However, we have found that this linear dependence is not observed in all cases, even at high temperature. Interestingly though, Fig. 10 shows a nearly linear behavior for T > 600 K, but eqn (16) is inapplicable at lower temperatures.
Activation energies of reactions are sometimes estimated using Evans–Polanyi-type correlations based on bond dissociation energies.^{15,62} This approach not only ignores the temperature dependence of the activation energies but also assumes a perfect correlation between the barriers and bond dissociation energies. The computational accuracy is sacrificed in such an approach in order to achieve efficiency for estimating rate constants for a large number of elementary reactions proposed^{63,64} in combustion modeling.
Fig. 12 Branching ratios of forward reactions (R1)–(R5) at various temperatures. |
Below 2400 K, the order of the contributions from reactions (R1)–(R5) to the total reaction rate constant is (R2) > (R3) > (R1) > (R5) > (R4); above 2400 K, the order is (R3) > (R2) > (R1) > (R5) > (R4). The computed forward reaction rate constants indicate that hydrogen abstraction from the C-2 site has the largest contribution to the overall reaction from 200 K to 2400 K, with a contribution ranging from 99.9988% at 200 K to 88.9% at 800 K to 21.2% at 3000 K, while hydrogen abstraction from the oxygen-site makes the lowest contribution at all temperatures, ranging from 2.5 × 10^{−9}% at 200 K to 0.65% at 800 K to 18% at 3000 K.
We also computed the branching ratios using multistructural CVT/SCT, and the relative errors with respect to the multi-path CVT/SCT branching ratios have been listed in Table 5. The relative error can be up to ∼50% at 3000 K for (R1) and for (R3)–(R5). We conclude that if one evaluates the rate constants without taking the multi-path effects into consideration, one cannot predict the branching ratios with high accuracy.
The convergence criterion we set for the multi-path variational transition state theory (MP-VTST) calculations in this work is that we add reaction paths, in order of increasing saddle point energy, until the difference in averaged generalized transmission coefficients upon adding one more reaction path is less than 15%. With this criterion, the calculations on all the five reactions are converged with the four lowest-energy paths for each case. We also found that for the present reactions it is unnecessary to include the torsional anharmonicity as a function of reaction coordinate s.
We found that the decreasing order of the forward barrier heights is (R5) > (R1) > (R3) > (R4) > (R2). The calculations yield the site dependence of the reaction rates, which is difficult to obtain experimentally. The predicted forward reaction rate constants and the branching ratios illustrate that the overall reaction rate constants are mainly contributed by the hydrogen abstraction from the C-2 site below 2400 K and by the hydrogen abstraction from the C-3 site above 2400 K, while the hydrogen abstraction from the oxygen site contributes the least.
Hydrogen-bonded structures are very common in chemical reactions, and one might think that hydrogen-bonded transition structures, because of their generally lower energy, would accelerate the reaction rate. In the current work, we investigated the role played by hydrogen-bonded transition structures and concluded that a hydrogen-bonded transition structure does not necessarily contribute significantly to the thermal rate constants due to enthalpy–entropy compensation. This interesting result corrects a false impression that – if used in estimating rate constants for mechanisms – can lead to errors in rate constants; its consequences are very broad, extending beyond combustion chemistry to fundamental kinetics modeling in general.
The temperature dependence of the activation energies is very dramatic, for example, the activation energy of reaction (R4) increases from 11.4 kcal mol^{−1} at 200 K to 19.7 kcal mol^{−1} at 600 K to 36.1 kcal mol^{−1} at 2400 K. We also found that this dependence is not linear, and this reinforces recent suggestions that eqn (15) is an inadequate representation of the temperature dependence of rate constants.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5sc01848j |
This journal is © The Royal Society of Chemistry 2015 |