Path-dependent variational effects and multidimensional tunneling in multi-path variational transition state theory: rate constants calculated for the reactions of HO2 with tert-butanol by including all 46 paths for abstraction at C and all six paths for abstraction at O

Multi-path variational transition state theory (MP-VTST) provides a conformationally complete framework for calculating gas-phase rate constants. For reactions in which the transition state has distinguishable torsional minima (which include most reactions), there are multiple possible reaction paths. In principle MP-VTST includes the contributions from all the reaction paths, and it should explicitly treat the variational and tunneling effects of each path, but in practice one may need to truncate the number of paths included in MP-VTST calculations in order to achieve a balance between computational cost and accuracy. In this work, we present calculations including all paths for two prototype combustion reactions, namely the two hydrogen abstraction reactions from tert-butanol by HO2 radical. For both reactions we included all the reaction paths. Since abstraction at C has 46 paths, it provided a good opportunity to carry out a case study in which we investigated the errors introduced by truncating the number of paths. For the reaction studied, we found that the variational and multidimensional tunneling transmission coefficients are very different for different reaction paths, which provides new evidence that MP-VTST is necessary for treating path-dependent variational effects and multidimensional tunneling. We found that tunneling transmission coefficients can be much larger for higher-energy paths than for lower-energy ones. Interestingly, the simple hypothesis that higher barriers are narrower does not explain this finding in the present case; we found instead that the effect is due to higher-energy barriers having the possibility of tunneling at energies farther below the barrier top. We also show that a previously applied criterion for judging convergence with respect to the number of paths may not be reliable at low temperature.


Introduction
The recently developed multi-path variational transition state theory 1,2 (MP-VTST) provides a framework for taking account of all reactant and transition state conformations and all conformationally different reaction paths for calculating thermal rate constants and predicting experimentally unavailable branching ratios for complex reactants. It should be especially useful in combustion chemistry, atmospheric chemistry, and all problems involving reactions of chain molecules. MP-VTST is an extension of multi-structural variational transition state theory 3 (MS-VTST), which has now been used to study the chemical kinetics of a variety of chemically important systems, including hydrogen abstraction reactions in combustion of biofuels [4][5][6][7][8] and silyl anion-silane polymerization reactions in nanodusty plasma chemistry. 9 The internal rotations (torsions) of a molecule can lead to multiple conformational structures, i.e., multiple local minima (for reactants and products) or multiple first-order saddle points (for transition structures) on the potential energy surface, and therefore multiple reaction paths. The contributions from all the reaction paths that correspond to the same chemical reaction are included in the MP-VTST calculations, and the variational effect and quantum tunneling contribution of each path is treated explicitly in MP-VTST. In principle, all the reaction paths should be included in the calculations, and one should perform variational transition state theory calculations with multi-dimensional tunneling calculations for all the paths. Such a calculation can be based on minimum-energy paths (MEPs) generated by direct dynamics. Practically, however, the computational cost for such calculations, which involve calculating the energies, gradients, and generalized normal mode analyses at large numbers of points along each reaction coordinate, can be demanding. Although density functional theory can greatly accelerate such calculations as compared to reliable wave-functionbased calculations, it is often unaffordable to accurately calculate all the reaction paths, which are generated from multiple conformers of the transition state, even for a medium-sized transition state (for instance, the transition state for hydrogen abstraction reaction from C-3 of 1-butanol by HO 2 radical consists of 262 distinguishable conformers 10 ).
In our previous work, 4,8 the number of paths included in MP-VTST was decided by ordering the paths in order of increasing barrier height and comparing calculations averaged over the first P paths and averaged over the first P À 1 paths (or P and P À 2 paths if the paths come in optically active pairs), and if the computed final rate constants have no significant differences (for instance, if the relative difference is smaller than 15%), then we stopped including more reaction paths. For a given reaction, it would be tempting to assume that various reaction paths are similar to each other except for the quantitative energetics, in which case this would be a reasonably reliable method. However, this is a dangerous assumption. In the current work, we will confirm that this kind of agreement does not necessarily mean that the MP-VTST rate constant has been converged using the first P paths because the transmission coefficients can be very different for different paths.
Hydrogen abstractions by small radicals (such as H, O, OH, HO 2 , and CH 3 ) are a class of important reactions in the combustion mechanisms of common fuels. 11,12 In Sarathy and co-workers' work, 11 rate constants for HO 2 + tert-butanol H-abstraction reaction are theoretically computed based on Evans-Polanyi approach, which not only ignores the temperature-dependence of the activation energy but also assumes a perfect correlation between bond energy and barrier height. Ab initio variational transition state calculations are needed in order to more accurately predict the thermal rate constants. In this work, we choose the tertbutanol + HO 2 reaction to study the effect of truncating the number of paths included in MP-VTST calculations, and we use the MP-VTST rate constants calculated based on all reaction paths as a reference. To be more specific, the reactions we consider here are: Reaction (R1) has 46 reaction paths (23 pairs of mirror images), and (R2) has 6 paths (3 pairs of mirror images). (R1) is the reaction of major interest in the current work since it provides a large number of reaction paths and also because it has a larger rate constant than reaction (R2).

Theory
The thermal rate constant in multi-path variational transition state theory is given by 1,2 where k ConTST 1 is the rate constant computed by conventional transition state theory (ConTST) without tunneling using the lowest-energy transition structure, F MS-T act is the multistructural torsional anharmonicity factor of activation, which is defined as: where MS-R1, MS-R2, and MS-ConTS denote the multistructural reactant 1, the multi-structural reactant 2, and the multi-structural conventional transition state, and GM denotes ''global minimum''; Q MS-T is the multi-structure (MS) conformational-rovibrational partition function with torsional anharmonicity (T) computed based on a coupled torsional potential; 13 and Q SS-QH is the rovibrational partition function computed using the rigid-rotator-quasiharmonic-oscillator approximation based on a global minimum structure. When we say GM for a transition state we mean the lowest-energy first-order saddle point that connects reactants to products, i.e., the lowest-energy transition structure. Note also that a ''transition state'' is a dividing surface that passes through all the transition state structures, and the MS-ConTS is a dividing surface that passes through all the saddle points, which are also called transition structures. In a full calculation, F MS-T act includes the contributions from all the conformational structures, i.e., there is no truncation error introduced by limiting the sums of structures.
The variational effect and multidimensional tunneling (MT) contributions of each reaction path are explicitly treated using the path-averaged generalized transmission coefficient: where the subscript p denotes the index of a reaction path; P denotes the total number of reaction paths included in MP-VTST; Q ‡-SS-T p is the single-structural rovibrational partition function with torsional anharmonicity of the transition structure of path p; G CVT p is the variational recrossing transmission coefficient computed for path p; and k MT p is the multidimensional tunneling transmission coefficient of path p. The factor k MT p is equal to the ratio k CVT/MT p /k CVT p , where CVT/MT refers to the CVT rate constant including multidimensional tunneling. We should emphasize that the multidimensional tunneling transmission coefficient not only includes the effect of quantum mechanical tunneling (which increases the rate constant, especially at lower temperatures and for light particles), but also includes the contribution from nonclassical reflection (which decreases the rate constant) at energies above the barrier height. Thus a better technical name for k MT p would be the transmission coefficient that accounts for quantum mechanical effects on reaction coordinate motion, including its coupling to other degrees of freedom, but we call it tunneling transmission coefficient for brevity.
The multidimensional tunneling transmission coefficient is evaluated with the small-curvature tunneling approximation 14 (SCT) in the current work, but a few results are shown for the zero-curvature tunneling (ZCT) approximation [15][16][17] for comparison. In the ZCT approximation, the system tunnels through the vibrationally adiabatic ground-state potential curve defined for path p by where s is the signed distance (from the saddle point) along the MEP through mass-scaled coordinates with reduced mass m, V MEP, p (s) is the potential energy along MEP p, and e 0, p (s) is the local zero-point vibrational energy of modes transverse to this MEP. The vibrationally adiabatic ground-state potential curve is sometimes called the effective barrier for tunneling. In the SCT approximation, the tunneling is again governed by V G a, p (s), but ffiffiffi m p ds in the imaginary action integral (see next equation) is reduced to account for the shortening of the tunneling path by corner cutting in the small-curvature limit. 14,17 In practice, this reduction is accomplished by replacing m by m eff,p (s) o m, where m eff,p (s) depends on the reaction-path curvature at s. Then the SCT tunneling probability at energy E is where h is Planck's constant divided by 2p, s o is the turning point where V G a, p (s) = E on the left side of the maximum, and s o is the turning point where V G a, p (s) = E on the right side of the maximum.
In eqn (3), P is formally equal to the total number of distinguishable reaction paths (including paths that are mirror images), which is equal to the number N of distinguishable structures of the transition state (TS). Practically, due to the existence of a large number of reaction paths, we often have to limit the number of paths included in MP-VTST calculation by taking P o N. Notice that, if P = N, the denominator of eqn (3) is equal to Q MS-T TS ; and if P o N, and this introduces truncation errors. The possible sources of errors introduced in the final calculated rate constant are complex to analyze. In the present work, we focus on the errors introduced by truncating the number of reaction paths in MP-VTST calculations, regardless of other factors such as the accuracy of electronic structure calculations or the accuracy with which anharmonicity is treated.

Computational details
An exhaustive conformational search was carried out for reactant, product, and transition structures by using M08-HX 18 /MG3S 19 for the electronic structure and MSTor 20 for the vibrational anharmonicity and vibration-rotation coupling. In order to choose a reasonably accurate electronic structure method for carrying out VTST calculations, classical barrier heights for both forward and reverse reactions (V ‡ f and V ‡ r ), and energy of reaction (DE) were computed based on the lowest-energy structures. (Classical barrier heights and energies of reaction are difference in Born-Oppenheimer energies for stationary structures; they do not include vibrational zero-point energy or thermal energy.) Four hybrid density functionals, M05, 21 M06-2X, 22,23 M08-HX, 18 and M08-SO, 18 were used, combined with MG3S 19 and ma-TZVP 24 basis sets. The explicitly correlated coupled cluster theory CCSD(T)-F12a 25-27 with the jun-cc-pVTZ basis 28 was chosen as the benchmark. All the DFT calculations were performed with integration grids of 99 radial shells and 974 angular points per shell using a locally modified Gaussian 09 29,30 program; coupled cluster calculations were carried out with Molpro2010 31 software.
Based on the benchmark results discussed below, M08-HX/ MG3S was used for VTST calculations. Canonical variational transition state theory calculations were carried out in nonredundant internal coordinates 32,33 by using the Gaussrate 34 and Polyrate 35,36 programs. Multidimensional tunneling transmission coefficients were computed based on the smallcurvature tunneling approximation and reaction paths extending from s = À2.0 to +2.0 bohr, where s is the reaction coordinate taken as the distance from the saddle point along the isoinertial minimum-energy path scaled to a reduced mass of 1 amu. In calculating the rovibrational partition functions and performing generalized normal mode analysis, all the vibrational frequencies are scaled with a factor of 0.973. 37

Results and discussion
Benchmarking density functionals We will judge the accuracy of density functionals for the present reactions by estimating errors as the differences of their predictions of classical barrier heights and classical energies of reaction at stationary-point geometries optimized by M08-HX/MG3S to the results obtained by CCSD(T)-F12a/juncc-pVTZ with the same geometries. A classical barrier height is the difference in potential energy between the lowest-energy saddle point (of the transition state) and the lowest-energy structure of reactants (forward barrier height) or the products (reverse barrier height); it is the main determinant of the magnitude of the rate constant. A classical energy of reaction is the difference in potential energy between the lowest-energy structure of the products and the lowest-energy structure of the reactants; it is the main component of the free energy of reaction, which is related to the chemical equilibrium constant, that connects the forward and reverse rate constants according to the principle of detailed balance. Classical barrier heights and energies of reaction computed by various density functionals for reactions (R1) and (R2) are tabulated in Table S1 (tables labeled with ''S'' are in the ESI †). Forward barrier height V ‡ f , reverse barrier height V ‡ r , and energy of reaction DE are not independent because they are connected by the relation V ‡ f À V ‡ r = DE. One may argue that for some purposes it is reasonable to compute mean unsigned errors (MUEs) based only on forward barrier heights and energies of reaction, rather than based on all three of them, whereas for other purposes it is appropriate to consider all three. Therefore we report two kinds of MUEs in Fig. 1: MUE(6) includes both forward barrier heights, both reverse barrier heights, and both energies of reaction; MUE(4) only includes forward barriers and energies of reaction. MUE (6) and MUE(4) differs less than 0.4 kcal mol À1 in all the cases, and the use of either error criterion leads to similar conclusions about the relative accuracies of the method. The numerical results for the three most accurate of the eight methods tested are shown in Table 1. Among the model chemistries we tested, M08-HX/MG3S has the smallest MUEs (0.6 kcal mol À1 ) and therefore is chosen for rate constant calculations.

Multi-structural anharmonicity factors, rate constants, and branching fractions
The F MS-T act factors computed by M08-HX/MG3S are given for 12 temperatures in Table S2 (ESI †) and for six temperatures in Table 2. Computed MP-CVT/SCT rate constants for the forward and reverse reactions of (R1) and (R2) are given at 12 temperatures in Table S3 (ESI †) and for six temperatures in Table 3. The rate constants are fitted for 200-2400 K by the four-parameter functions proposed previously. 2 where R is the gas constant in kcal mol À1 , and T is temperature. The complete set of fitting parameters for both forward and reverse reactions are in Table S4 (ESI †). Branching fractions into the two possible products are shown in Fig. 2. We see that the yield of product R2 increases from a very small number at low temperature to about 20% at high temperature.
Activation energies and phenomenological Gibbs free energies of activation    The derivatives were computed analytically from the fits (see Table S4, ESI †). We see that, over the temperature range studied, the activation energy increases by B25 kcal mol À1 for reaction (R1) and by B15 kcal mol À1 for reaction (R2). This temperaturedependence of activation energy is significant and should not be ignored in theoretical combustion modeling. Since reaction (R1) dominates the overall reaction, the results for the overall reaction are very close to those for reaction (R1). Fig. S1 (ESI †) shows phenomenological Gibbs free energies of activation, which are given for bimolecular gas-phase reactions by 40,41 DG act;fwd ¼ RT ln where k B and h are the Boltzmann and Planck constants, k MP-CVT/SCT fwd is in the units of m 3 molecule À1 s À1 , and p 0 is the standard-state pressure (10 5 Pa). As for the activation energy, the result is a rapidly increasing function of temperature.

Path-dependent variational effects
The variational transition states are slightly away from the saddle points. For example, numbering the paths in order of increasing classical barrier height, we find that at 200 K, the canonical variational generalized transition state is 0.002, 0.012, 0.009 Å away from the saddle point for paths 1, 41, and 45, respectively; at 600 K, they are 0.003, 0.023, 0.013 Å away from the saddle points for these paths; and at 2400 K, they are 0.008, 0.081, and 0.039 Å away from the saddle points for these paths. Fig. 4 shows the computed CVT variational recrossing transmission coefficients of reaction (R1) for temperatures from 200 K to 2400 K for paths 1, 7, and 41 (the results for paths 2, 8, and 42 are the same as for 1, 7, and 41 because the paths come in enantiomeric pairs). Over the entire temperature range, the variational recrossing effect is negligible for path 1, for which the recrossing transmission coefficient varies from 0.998 at 200 K to 0.996 at 2400 K. For reaction path 7, which emanates from a transition structure 0.97 kcal mol À1 higher than that generating path 1, the behavior is significantly different from that for path 1. Variational recrossing transmission coefficients of path 7 increase from 0.913 at 200 K to a maximum 0.941 at 600 K and then gradually decrease to 0.930 at 1500 K, followed by a rapid decrease to 0.750 at 2400 K. For reaction path 41, there is an almost linear decrease of the variational recrossing transmission coefficient with temperature from 250 to 2400 K. At 2400 K, the variational transmission coefficient for the lowest-energy path is larger than the one for path 41 by a factor of 1.4.    tunneling transmission coefficients increase dramatically as the barrier heights increase; the tunneling transmission coefficients of the highest-energy paths are very much larger than the ones of the lowest-energy path, with the largest ratio being 17 and several other ratios being in the 9-12 range. As expected, the tunneling transmission coefficients are much closer to unity at 1000 K, because the Boltzmann factor for overbarrier reaction is much larger, but perhaps not so expected is the finding that the tunneling transmission coefficients cluster close to one another -with an average value of 1.35 and a standard deviation of only 0.04. We may now ask: why are the tunneling transmission coefficients larger for the paths with higher classical heights?

Multidimensional tunneling
The answer to this question will illuminate an interesting feature of the way that tunneling effects enter a multi-path calculation. In order to give this answer we need to consider some aspects of the tunneling in more detail, and we now turn to this. Fig. 6 shows both the ZCT and SCT tunneling transmission probabilities of various paths of reaction (R1) as functions of tunneling energy; the figure includes the lowest-and highestenergy paths and some intermediate-energy paths. SCT tunneling transmission probabilities are greater than ZCT ones because reaction-path curvature increases the tunneling probability by allowing corner cutting. The zero of energy for the tunneling calculations is the potential energy of the equilibrium structure of reactants. For each path the semiclassical tunneling probability rises approximately exponentially with energy until it reaches 0.5 at the energy of the maximum of V G a, p (s) for that path; this maximum is called V AG a, p , and it is usually higher for paths with a higher ordinal number, i.e., for paths with a higher classical barrier height.
The V G a, p (s) curves for reaction for path 1 (which, along with its mirror image, path 2, is the lowest-energy path), for path 41, and for path 45 (which, along with its mirror image, path 46, is the highest-energy path) are shown in Fig. 7, and the potential energies along the minimum energy paths (MEPs through the isoinertial coordinate system, which are sometimes called IRCs (intrinsic reaction coordinates)) for these three reaction paths are shown in Fig. 8. Fig. 9 shows the forming O-H bond distance as a function of the breaking C-H bond distance along the MEPs for paths 1, 41, and 45. The MEPs are very similar for these two critical bond distances, but the shapes of the vibrationally adiabatic potential energy curves are different because of the variations of the zero-point vibrational energy along the reaction paths. We see that V G a, p (s) peaks close to, but not precisely at s = 0, which, by convention, is placed at the saddle point, which is the maximum of V MEP (s).  The semiclassical tunneling transmission coefficient for path p is given by where E 0 is the maximum of the ground-state vibrational energy of reactants and the ground-state vibrational energy of products (this is the lowest energy at which reaction can occur), Prob p (E) is the transmission probability of eqn (5), and Prob C p (E) is the classical transmission probability along path p at energy E. Since the classical transmission probability is a unit step function at E = V AG a , we can carry out the integral in the denominator, which yields the combination of eqn (5) and (12) shows that the different k p for various paths may result from different effective reduced masses or different V G a, p (s). However, since the higher-energy paths also have larger transmission coefficients in the ZCT approximation, and the reduced masses are all the same in the SCT approximation, we look to the shape of the V G a, p (s) for an explanation.
The energy at which the integrand of eqn (12) peaks is called the representative tunneling energy, E rep; p T ð Þ; some examples (in the SCT approximation) are shown in Table 4. Table 4 shows that representative tunneling energies are 2.6-5.2 kcal mol À1 below the effective barrier maxima. A simple hypothesis might be that the higher-energy paths show larger tunneling transmission coefficients because the effective barriers are thinner near their tops. However, examination of the V G a, p (s) curves shows that, in this energy range, these curves have about the same widths (see Table S7, ESI †). Consistent with this, Table 4 shows that the tunneling probabilities are about the same for all paths at energies 2 to 3 kcal mol À1 below the maximum of the effective potential energy curve for tunneling. So this hypothesis Fig. 7 Vibrationally adiabatic potential energy (in kcal mol À1 ) along the minimum energy paths (MEP) for paths 1, 41, 45 of reaction (R1), where s is the signed distance from the saddle point (i.e., transition structure) along each curved MEP in isoinertial coordinates scaled to a reduced mass of 1 amu. Note that each conformer of the transition state generates its own reaction path, and so each curve in this figure corresponds to a different reaction path, where the paths are shown in Fig. 9. The transition structures are shown next to the vibrationally adiabatic ground-state potential energy curve for the path they generate. , where s is the signed distance from the saddle point (i.e., transition structure) along the curved MEP in isoinertial coordinates scaled to a reduced mass of 1 amu. Note that each conformer of the transition state generates its own reaction path, and so each curve in this figure corresponds to a different reaction path, where the paths are shown in Fig. 9. The transition structures are shown next to the vibrationally adiabatic ground-state potential energy curve for the path they generate. does not explain the observation of larger tunneling transmission coefficients for the higher-energy paths.
We plot the integrand of eqn (12) for the SCT tunneling approximation at 200 K in Fig. 10 for various paths of reaction (R1). The curves in this figure start at an energy of 107.7 kcal mol À1 because that is the energy of the lowest-energy state of the products (i.e. CH 2 (CH 3 ) 2 COH radical and H 2 O 2 ) of this endothermic reaction (the zero of energy is the potential energy of separated reactants at their equilibrium geometries; on this scale the potential energy of products is 17.1 kcal mol À1 , and their zero-point energy is 90.6 kcal mol À1 ). The argument of the exponential suggests that it would be instructive to replot the results in Fig. 10 as functions of EÀV AG a,p , and this is done in Fig. 11. Fig. 11 shows that the curves are very similar as functions of this variable, but curves extend lower for the higher-energy curves, giving more area under the curve. The reason for this is that the curves stop at 107.1 kcal mol À1 , as explained in the previous paragraph. Since eqn (12) shows that the area under the curve is proportional to the tunneling transmission coefficient, we see that the reason that higherenergy paths have larger tunneling transmission coefficients is that they can tunnel at energies farther below the maximum of the effective barrier for tunneling. Thus, for an endothermic reaction, the higher-energy paths have more to gain by tunneling because they can tunnel at energies farther below the barrier top.

Truncation of number of paths included in MP-VTST
The number of paths included in MP-VTST affects the values of the generalized transmission coefficients hgi P . Unsigned percentage errors as compared to using all N paths and unsigned percentage deviations between two successive approximations of the generalized transmission coefficients for reaction R1 are shown in Fig. 12. The unsigned percentage errors (UPEs) are computed as: where UPE P is the unsigned error of generalized transmission coefficient that includes P distinguishable paths (including paths that are mirror images), and hgi N is the exact generalized transmission coefficient that includes all 46 distinguishable paths. Unsigned percentage deviations (UPDs) are computed as: where UPD P is the unsigned deviation between the generalized transmission coefficient based on an average over P paths and the one that is based on PÀ2 paths. (We use PÀ2 rather than P À 1 because the paths occur in enantiomeric pairs.) Since the generalized transmission coefficient is the only P-dependent quantity in eqn (1), UPD P defined in eqn (14) is equivalent to the unsigned percentage deviation of the MP-VTST rate constants computed including P and P À 2 paths. a Evaluated at an energy 2 kcal mol À1 below V G a, p (s). b Evaluated at an energy 3 kcal mol À1 below V G a, p (s).

Fig. 10
Values of the integrand in eqn (12) at 200 K for various paths of reaction (R1) as functions of total energy E, with the zero of energy at the overall zero of energy, which is taken as the energy of the equilibrium structure of reactants. For small values of P, UPD P is fluctuating rather than converging as P increases; thus, even if UPD P is very small (o5%), the truncation error could still be quite large. For instance, at 200 K, UPD 8 is only 2.5%, however UPE 8 is 30%. Nevertheless, at high temperatures (for instance, 1500 K), including a few paths is accurate enough for the final generalized transmission coefficients.

Summary
In the current work we discussed the errors introduced in truncating the number of paths in MP-VTST calculations. For hydrogen abstraction from tert-butanol by HO 2 radical, we found that the variational effects and the contributions of multidimensional tunneling are very different among the various reaction paths. As a consequence, including only one or a few paths does not yield accurate rate constants at low temperature. Over the entire temperature range, the variational effect is negligible for the lowest-energy path; but for some of the higher-energy paths, the behavior is significantly different from that for the lowest-energy path; we gave an example where the variational recrossing transmission coefficients for different paths vary by up to a factor of 1.4. The differences are even greater -much greater -for the tunneling transmission coefficients. At 200 K, the highest-energy path has a smallcurvature tunneling transmission coefficient 17 times larger than the one for lowest-energy path. The higher-energy paths can have much larger tunneling transmission coefficients because the effective barrier for tunneling extends higher above the reactants and products, giving a greater energy range in which tunneling can operate. We conclude that the multi-path methodology is necessary for better understanding of the detailed kinetics and predicting the thermal rate constants. Fig. 12 Unsigned percent deviations of the generalized transmission coefficients for reaction (R1). The number of paths included (P) is the abscissa. The paths are indexed 1-46, with path numbering in order of increasing classical barrier height. The top figure is the deviation between including P paths and including all N paths. The bottom figure is the deviation between including P paths and including P À 2 paths. See eqn (13) and (14).