Kinetics of the Toluene Reaction with OH Radical

We calculated the kinetics of chemical activation reactions of toluene with hydroxyl radical in the temperature range from 213 K to 2500 K and the pressure range from 10 Torr to the high-pressure limit by using multistructural variational transition state theory with the small-curvature tunneling approximation (MS-CVT/SCT) and using the system-specific quantum Rice-Ramsperger-Kassel method. The reactions of OH with toluene are important elementary steps in both combustion and atmospheric chemistry, and thus it is valuable to understand the rate constants both in the high-pressure, high-temperature regime and in the low-pressure, low-temperature regime. Under the experimental pressure conditions, the theoretically calculated total reaction rate constants agree well with the limited experimental data, including the negative temperature dependence at low temperature. We find that the effect of multistructural anharmonicity on the partition functions usually increases with temperature, and it can change the calculated reaction rates by factors as small as 0.2 and as large as 4.2. We also find a large effect of anharmonicity on the zero-point energies of the transition states for the abstraction reactions. We report that abstraction of H from methyl should not be neglected in atmospheric chemistry, even though the low-temperature results are dominated by addition. We calculated the product distribution, which is usually not accessible to experiments, as a function of temperature and pressure.


Introduction
Aromatic hydrocarbons are important in both combustion chemistry and atmospheric chemistry; they are also major petrochemicals, and they are used in synthesis, as solvents, and as fuel additives. They are known to play a major role in emission and soot formation (particulate polycyclic aromatic hydrocarbons) during fossil fuel combustion and in the formation of photochemical smog in the urban air [1,2]. Therefore, the oxidation mechanisms and reaction kinetics of aromatic hydrocarbons are of great interest over wide pressure and temperature ranges. Toluene is one of the simplest aromatic hydrocarbons, and its reaction with the OH radical, a main oxidizing reactive species in both combustion and the atmosphere, has been extensively investigated both experimentally and theoretically [3][4][5][6][7][8][9][10][11][12][13]. The main reaction channels of toluene with OH radical are hydrogen abstraction reactions (R1-R4) and addition reactions (R5-R8), as shown in Figure 1. (Note that reactions R5-R8 in this figure involve the formation of bonded adducts, but the figure does not show the pre-reactive van der Waals complex of OH with toluene. We assume a fast pre-equilibrium between the reactants and this complex, and this has an effect on the tunneling calculations for the hydrogen abstraction reactions [11].) The first experimental study of the kinetics of the OH + toluene reaction was carried out by Davis et al. [3] at 298 K and 3-100 Torr total pressure; they found the rate constant to be pressure-dependent, which they interpreted as showing that a significant amount of reaction occurs by addition to the ring. Additional experiments were presented by Doyle et al. [4] (304 K, 1 atm), Hansen et al. [5] (298 K, 100-618 Torr), and Perry et al. [6] (298-473 K, 100-200 Torr), where the values in parentheses are ranges of temperature and total pressure. Neither Hansen nor Perry found significant pressure dependence. For temperatures between ∼325 K and 380 K, Perry et al. found nonexponential decay of the total rate constant, indicating reversible reaction; hence their analysis is not valid at those temperatures.
The most complete set of available experiments for OH + toluene kinetics are those of Tully et al. [ OH OH Figure 1: The reactions whose rate constants are calculated in this work.
assumption that the loss of OH followed pseudo-first-order kinetics without back reaction, but at temperatures in the interval 320-352 K, the loss of OH showed nonexponential decay, indicating that the assumption was not valid and that an accurate rate constant cannot be obtained from their experiments at those temperatures. Furthermore, at temperatures from 352 to 442 K, their extracted pseudo-firstorder rate constants depended linearly on the concentration of toluene only over a small concentration interval at low concentrations of toluene and the measured rate constants depended on the initial OH concentration, so the extracted rate constants have large uncertainties in this range. Tully et al. concluded that the observations of nonexponential decays and the nonlinear concentration plots "cannot satisfactorily be explained at this time." However, they did interpret some of their observations in terms of reversible reaction, i.e., decomposition of the newly formed toluene. . .OH adducts back to reactants. Based on these observations, only their experimental results in the temperature intervals not affected by these problems are meaningful for comparison with our theoretical estimates of the elementary-reaction rate constant for the forward reaction; these temperature intervals are 213-298 K and 504-1150 K. Seta et al. [8] measured the reaction rate at 919-1481 K (at 0.8-1.8 atm), where it is essentially all due to abstraction reactions. Vasudevan et al. [9] reported the reaction constant between 911 K and 1389 K at the pressure of ∼2.25 atm. Their rate constants range from 13% lower than those of Seta et al. at 919 K to 39% lower at 1389 K.
Uc et al. [10] used transition state theory (TST) with the Eckart approximation to zero-curvature tunneling to calculate the rate constants of the hydrogen abstraction reactions at 275-1000 K and the overall rate constants (taking into account all possible abstraction and addition reactions) at 200-600 K; they used a potential energy surface calculated by CBS/QB3//BHandHLYP/6-311++G(d,p). They included torsional anharmonicity by a one-dimensional approximation. Seta et al. [8] also estimated rate constants using TST with the Eckart approximation to zero-curvature tunneling, but in their case with B3LYP/6-31G(d) potential energy surfaces on which they adjusted the barrier heights to obtain agreement with experiment at high temperature, and they estimated the effects of anharmonicity originating from the -OH torsion and the C-H-O rocking. Li et al. [12] used TST with the Wigner tunneling approximation with G4//B3LYP/6-31G(2df,p) potential energy surfaces to estimate the rate constants for the abstraction paths; they do not mention including anharmonicity in the rate calculation. Recently, Pelucchi et al. [13] reported more reliable calculations of the rate constants for the abstraction reactions by employing canonical variational transition state theory (CVT) with the Eckart approximation to zero-curvature tunneling and with the potential energy surface calculated by M06-2X/6-311+G(d,p) and CCSD(T)/aug-cc-pVTZ with a correction for basis set effects. In their calculations, they considered the effect of torsional anharmonicity on the rate constants by using a two-dimensional hindered rotor method for hydrogen abstraction from the methyl group.
The above summary shows that the H-abstraction reactions have been studied much more than the pressuredependent addition reactions. The addition of OH to toluene is similar in many respects to the addition of OH to benzene, which was also studied in some of the papers cited above, as well as by several other workers. Lin et al. included the reversible reaction in an analysis of the OH + benzene problem and were able to obtain the equilibrium constant for the formation of the benzene...OH adduct as well as the forward and backward rate constants for its formation [14]. In later work, Uc et al. reanalyzed the experiments of Tully et al. on benzene and toluene with the help of theoretical equilibrium constants for adduct formation [11]. The calculations indicated that back reaction of the adduct becomes increasingly important at temperatures of 325 K and greater, and they were able to explain the signals observed by Tully et al. in terms of the contribution of the reverse reaction. This complication in the kinetics of the addition reaction does not affect the measured results at high temperature where the abstraction channel dominates the total reaction.
The work of Uc et al. stressed the role of the reverse dissociation reaction in leading to nonexponential decay of the OH concentration. This motivates explicit consideration of the pressure effect because the unimolecular dissociation reaction of the adduct is pressure-dependent. As far as we know, no theoretical study has been reported on the pressure effects on the toluene plus OH reaction. Experimentally, the pressure effect has mainly been studied at 298 K. As mentioned above, Davis et al. [3] found the rate constant to be pressure-dependent at 298 K and 3-100 Torr total pressure; but neither Hansen [5] nor Perry [6] found significant pressure dependence at 298 K and above 100 Torr. Tully et al. [7] reconfirmed their work by finding that the third-body dependence was only 9% at 298 K and 100 Torr in He, Ar, or SF 6 ; they concluded that the experiments at 298 K and over 100 Torr are close to the high-pressure limit. Lowering the pressure to 25 Torr in Ar lowered the rate by 10%, but the experimental uncertainty in the rate constants is 11% for that case. Lowering the pressure to 20 Torr in He lowered the rate by 21%, which is significant compared to the 10% experimental uncertainty in the rate constants for that case. With increasing temperature, the pressure fall-off of the addition reaction rate constant should be larger, although systematic study of the effect of pressure on rate constants for T > 298 K is missing. Therefore, it is interesting to explore the pressure dependence quantitatively.
In the present work, we carried out a theoretical study of kinetics of the toluene reaction with OH that includes the pressure effects.

Electronic Structure Calculations
An accurate potential energy surface is the foundation of reliable dynamics calculations. Although generally reliable wave function methods are available for small, weakly correlated systems, it is impractical to use them in direct dynamics calculations for medium-or large-size systems due to the high computational cost. A more practical scheme is to choose a density functional for a Kohn-Sham (KS) model chemistry for the specific system by benchmark calculations of the stationary points on the potential energy surface of that system.
We will use the following notation: Born-Oppenheimer potential energies are electronic energies (including nuclear repulsion) calculated for fixed nuclear coordinates. The classical energy of reaction ΔE is defined as the Born-Oppenheimer potential energy of the lowest-energy equilibrium structure of products minus the Born-Oppenheimer potential energy of the lowest-energy equilibrium structure of reactants; the classical barrier height V ‡ is defined as the Born-Oppenheimer potential energy of the lowest-energy saddle point minus the Born-Oppenheimer potential energy of the lowest-energy equilibrium structure of reactants. Adding the change in zero-point energy (products minus reactants) to the energy of reaction yields the enthalpy of reaction at 0 K. Adding to the classical barrier height the change in zero-point energy on proceeding to the transition state from the reactants yields the enthalpy of activation at 0 K. The temperature-dependent Arrhenius activation energy is obtained from the local slope of an Arrhenius plot of the logarithm of a rate constant versus the reciprocal of the temperature.
Gaussian 09 [15] and a locally modified version of Gaussian [16] were used for the density functional calculations.

Best Estimates of Reaction Energetics.
To find and validate a successful model chemistry for present system, we first calculated the geometrical structures of the stationary points (reactants, transition structures, and products) on the potential energy surface by using the M08-HX/MG3S method, which has been recommended for locating transition state structures in a previous benchmark study [17]. Subsequent frequency calculations using the same method confirmed that the structures we obtained are the desired minima and transition structures. Then we use these geometries to calculate the T 1 diagnostics by the CCSD [18] /jun-cc-pVTZ [19][20][21][22] method. The standard assumptions are that closedshell species with 1 > ∼0.02 and open-shell species with 1 > ∼0.045 are strongly correlated and need to be treated with a multireference method, whereas systems with lower T 1 diagnostics are weakly correlated and can be treated reliably by a single-reference method. As shown in Table 1, the obtained T 1 values are all smaller than these borders between small and moderate multireference character. This implies that we can use the single-reference CCSD(T)-F12a [23,24] /jun-cc-pVTZ method, a highly accurate method for weakly correlated systems, to obtain the best estimates of the classical energies of reaction and barrier heights and then use them to validate the KS model chemistry.
Molpro 2012.1 [25] was used for the coupled cluster calculations.

Selection of KS Model
Chemistry. In our experience, the M06-2X [26], M08-SO [27], and MN15 [28] density functionals combined with the MG3S [29] basis set perform well in calculations of potential energy surfaces for many reaction systems. Therefore, we tested the performance of these three KS model chemistries against the best estimates, and we selected the most accurate of them to be used for direct dynamics calculations for hydrogen abstraction reactions and addition reactions. Table 2 lists the calculated classical energies of reaction (Δ ) and the forward ( ‡ 푓 ) and reverse ( ‡ 푟 ) barrier heights using these KS model chemistries and compares them to the best estimates obtained by CCSD(T)-F12a/jun-cc-pVTZ//M08-HX/MG3S; the table also gives the mean unsigned derivations (MUDs) of these KS results from the best estimates.
As shown in Table 2, the M06-2X/MG3S method performs best for hydrogen abstraction reactions with a mean unsigned deviation (MUD) of 0.62 kcal/mol over all classical barrier heights and reaction energies of the abstraction reactions, and the M08-SO/MG3S model chemistry is the best choice for addition reactions with a MUD of 0.58 Channel

Pre-Complex -4.73
Pre-Complex -4.82 TS-add-m 0. 37 TS-add-p -0. 19 TS-add-ipso -0.87 TS-add-o -1.70 Figure 2: Enthalpy of activation profiles at 0 K for the toluene reaction with OH, as calculated by the M06-2X/MG3S method for hydrogen abstraction channels and by the M08-SO/MG3S method for OH-addition channels. Enthalpy of activation profiles at 0 K includes the zeropoint energy, which is calculated by the quasiharmonic approximation with frequencies scaled with the SRP scaling factors.
kcal/mol over all classical barrier heights and reaction energies of the addition reactions. Therefore, these two model chemistries were chosen for the direct dynamics calculations of the corresponding temperature-dependent reaction rate constants.
The enthalpy of activation profiles at 0 K are shown in Figure 2. [8,30] on radical addition to toluene has sometimes considered the reactive decomposition of adducts. In the present case we estimated the enthalpies of activation for decomposition of the adducts to make phenol + methyl radical are 4.7-8.5 kcal/mol (at 0 K) higher than the enthalpies of activation for nonreactive reversion to reactants. This agrees with the general trend found by Seta et al. [8]. Therefore these reactions are not major pathways, and we did not consider them further.

Dynamics Calculations
For dynamic calculations we mainly used the same methods as those used for toluene + H reactions by Bao, Zheng, and one of the current authors [30]. In particular, for the temperature-dependent high-pressure-limit (HPL) rate constants we used multistructural [31] canonical variational transition state theory (CVT) [32] with small-curvature tunneling [24,33]. As usual, we denote this method as MS-CVT/SCT. In CVT, the transition state is variationally located at the position of maximum free energy of activation along the minimum energy path [34]. In the later discussion we shall refer to the difference between CVT and conventional transition state theory (TST, in which the transition state passes through the saddle point) as a variationaltransition-state-location effect or sometimes simply as a variational effect.
In these calculations, we used the multistructural torsional anharmonicity method (MS-T) [35] to treat the multistructural anharmonic effects, by which we mean the torsional potential anharmonicity, the multiple conformers, and the torsion-rotation coupling. We used specific-reactionparameter (SRP) scaling factors [36,37] to scale all calculated vibrational frequencies of all species (reactants, pre-reactive van der Waals complex, transition structures, and products); these factors are explained in more detail below. When we use scaled frequencies with the harmonic oscillator formulas for partition coefficients, the result is called quasiharmonic; we also used scaled frequencies in the MS-T calculations. The ratio of the MS-T partition function to the quasiharmonic one for species is called MS-T 훼 .
The pressure-dependent rate constants of the addition reactions were calculated by the system-specific quantum Rice-Ramsperger-Kassel [30] (SS-QRRK) method.

MS-CVT/SCT Theory.
The temperature-dependent highpressure-limit rate constants calculated by the MS-CVT/SCT theory are expressed by where CVT is the CVT reaction rate constant obtained in the single-structure quasiharmonic oscillator approximation, Γ is the recrossing transmission coefficient (also called the variational-transition-state-location effect or the variational effect) that accounts for recrossing of the conventional transition state by calculating the flux through a dividing surface corresponding to the maximum free energy of activation, rather than a dividing surface corresponding to the conventional transition state, TST is the conventional transition state theory rate constant, SCT is the calculated tunneling transmission coefficient in the small-curvature tunneling approximation (although it is called the tunneling transmission coefficient, it also includes nonclassical reflection at energies over the barrier top), and MS-T act is defined [41] as the multistructural anharmonicity factor of the reaction. The MS-T act factor is calculated as the ratio of multistructural anharmonicity factor MS-T 훼 of the transition state ( = TS) to the multistructural anharmonicity factor of the reactant ( = R).
In the calculation of CVT , only the lowest-energy structures of the reactants and transition states are used, even when more than one conformer is located for the species. The higher-energy conformers are accounted for by MS-T act ( ). In the high-pressure limit the pre-reactive complex (an OH. . .toluene van der Waals complex) is fully equilibrated and thermalized, so we took the ground-state energy of this complex as the lowest-energy at which the tunneling occurs. This is called the pre-equilibrium model (PEM). This model was explained in detail in a recent study of the methanol reaction with OH radical [42].
The standard vibrational frequency scaling factor [36] ZPE is parametrized to reproduce the accurate zero-point energies (ZPEs) of a set of stable molecules, and this scaling factor is used for reactants, pre-reactive van der Waals complexes, and products, but previous work on transition states of some hydrogen abstraction reactions by OH radical [36,42] has indicated that this scaling factor is not suitable for such transition states. Therefore, we used a specificreaction-parameter scaling factor ZPE-SRP for each transition state to correct the anharmonicity and systematic errors in ZPE calculations for a given model chemistry; for a given model chemistry and a given species this can be factored as where Anh and H correct, respectively, the anharmonicity of the ZPE and the inaccuracy of the model chemistry. Note that H is characteristic of the model chemistry and Anh is characteristic of the species. For all species, we use the harmonic factor H that has been parametrized to obtain the accurate harmonic frequency in the F38/10 database [36]. We calculate Anh for a given species as the ratio of anharmonic ZPE estimated by hybrid [43], degeneracy-corrected [44], second-order vibrational perturbation theory (HDCVPT2) [45][46][47] to harmonic ZPE calculated for that species with the same model chemistry. In the calculations of Anh , the M06-2X/MG3S method was used for the hydrogen abstraction reactions, and the MPW1K/MG3S model chemistry was used for addition reactions. We set ZPE-SRP equal to the standard ZPE for reactants, pre-reactive van der Waals complexes, and products. Table 3 lists the calculated ZPE-SRP scaling factors of the transition states for all reactions in the present study, as well as the H factors and the standard ZPE factors. For hydrogen abstraction transition states calculated with M06-2X/MG3S, the table shows that the specific-reaction-parameter scaling factors ( ZPE-SRP = 0.960-0.966) are smaller than the standard scaling factor ( ZPE = 0.970). For the addition transition states the calculated ZPE-SRP scaling factors are the same as the standard factor ZPE (0.983). These results show that the hydrogen abstraction transition states are more anharmonic than the stable reactants or the transition states of the addition reactions.  Figure 1 are treated by the chemical activation mechanism:

The SS-QRRK Method. Addition reactions R5-R8 in
where T denotes the toluene molecule, T is temperature, and TOH * is the energized adduct. Applying the steadystate approximation to the energized adduct TOH * yields the phenomenological bimolecular stabilization rate constant, stab ( ) =

[T] [OH]
[TOH] where the square bracket denotes the concentration, M is the bath gas, and ( ) is the fraction of energized adduct TOH * with energy E. We use the ideal gas law, so [M] equals / , where p is the pressure and R is the gas constant.
The three elementary rate constants that need to be determined in (4) are the temperature-dependent highpressure-limit rate constant 1 ( ), the energy-dependent dissociation rate constant −1 ( ), and collisional deactivation rate constant 푐 . These are evaluated as follows: (a) 1 ( ) is evaluated by the MS-CVT/SCT method explained above.
(b) The high-pressure-limit dissociation rate constants (d) −1 ( ) is calculated by using the SS-QRRK method, which is like conventional QRRK theory [48,49] except that the frequency factor A is obtained from (6) and the threshold energy E 0 is obtained from (5). Because −1,∞ ( ) incorporates variational effects, multidimensional quantum mechanical tunneling, and multistructural anharmonicity, they are also inherited in −1 ( ).
(e) 푐 is estimated by collision theory as the product of the Lennard-Jones collision rate constant LJ and the collision efficiency 푐 . The collision efficiency 푐 is calculated by where 퐸 is the energy dependence factor of the density of states calculated by Troe's method [50] and the ⟨Δ ⟩ is the average vibrational energy transferred per collision during both energization and deenergization. In the present study-when not specified otherwise-we set ⟨Δ ⟩ equal to the experimental values [51]: 130 cm −1 for Ar and 75 cm −1 for He. In the calculations of LJ , we based the Lennard-Jones parameters on those used in reference [51] to derive the values of ⟨Δ ⟩; in particular, we used / B = 410 K and = 6.0Å for the adducts of OH and toluene [52], 120 K and 3.4Å for Ar [53], and 10 K and 2.55Å for He [51]. The Polyrate 2016 [54] and Gaussrate 17 [55] programs were used for the direct CVT/SCT calculations and SS-QRRK calculations, and the MSTor program [56] was used to calculate the conformational-rotational-vibrational partition functions.

Computational Details of Dynamics Calculations.
The calculations were converged with respect to the lengths of the calculated reaction paths; in isoinertial coordinates scaled to a mass of 1 amu, the ranges covered were −1.2 to 1.2Å for addition reactions, −1.5 to 1.5Å for hydrogen abstraction reactions, and −3.0 to 3.0Å for all dynamics calculations involving the pre-reactive van der Waals complex of OH with toluene. The step size was set to 0.005Å for dynamics calculations starting from isolated reactants of all channels except for H-abstraction from ortho-site (for which we used a step size of 0.01Å) and to 0.01Å for the dynamics calculations starting with the pre-reactive van der Waals complex. For addition reactions we set the ratio of the number of gradients to the number of Hessians to 9 or 10. For the hydrogen abstractions, we calculated more Hessians in order to obtain smooth generalized free energy of activation profiles; in particular we set the Polyrate variable INH equal to 1 or 2 for regions near the saddle point and to 5 for regions away from the saddle point.

Conformers.
The reactants, toluene and OH radical, and the products of hydrogen abstraction reactions (R1-R4) are simple molecules that have only one conformer. However, the other species are more complex; there are two torsions, the methyl torsion and the OH torsion, that can produce distinguishable conformers for the transition state of each reaction and the products of the addition reactions (R5-R8). The number of distinguishable conformers that we found for each of the complex species is specified in Table 4. As indicated in (1a), (1b), only the lowest-energy structure of each species is used in the direct dynamics calculations in MS-CVT/SCT calculations; the effect of the other structures on the rate constants is included by the multistructural anharmonicity factor of the reaction, MS-T act . We found only one conformer for the pre-reactive complex (see Figure 3(a).) Uc et al. [10] reported a s -symmetry structure for the pre-reactive complex by the BHandHLYP/6-311++G * * method. However, this s -symmetry structure was found to be a saddle point with a small imaginary frequency by the M06-2X/MG3S method. We also confirmed that the lowest-energy structures of transition states in all reaction channels share the same pre-reactive complex (Figure 3(b)) by the IRC calculations. Therefore, we used this complex structure in the tunneling calculations for the high-pressurelimit rate constants of the hydrogen abstraction reactions.    3 2 C 푎 The number in parentheses represents how many pairs of mirror structures are included in the conformers. b U denotes that the two defined torsions are considered as uncoupled, and C denotes that they are treated as coupled. Tables 5 and 6 give the calculated enthalpies of activation and enthalpies of reaction at 0 K as calculated by using the lowest-energy structures of all species. The tables show that the hydrogen abstraction reaction always has a larger enthalpy of activation than the addition reaction at the same site, and the ortho-addition reaction (R5) has the smallest enthalpy of activation. The hydrogen abstraction from methyl group of toluene has the lowest enthalpy of activation among the hydrogen abstraction reactions. Based on our best estimates, which were obtained by using the CCSD(T)-F12a/jun-cc-pVTZ method for electronic energy and by using the selected model chemistries and the SRP vibrational frequency scaling factors for zero point energy, the calculated enthalpies of activation are in turn -0.16, 3.19, 3.61 and 3.88 kcal/mol for hydrogen abstractions from the methyl group and ortho, meta, and para sites, and they are -0.94, -0.44, 0.26 and 0.56 kcal/mol for addition at ortho, ipso, para, and meta sites.

Energetics.
As shown in Table 5, the calculated enthalpy of activation is sensitive to the scaling factor used for vibrational   frequencies. Due to the greater anharmonicity of highfrequency modes of transition states of hydrogen abstraction reactions, the enthalpy of activation is overestimated with the standard scaling factor ZPE or harmonic scaling factor H as compared to that obtained with the SRP scaling factor ZPE-SRP . Several groups have investigated the hydrogen abstraction reactions of toluene by OH radical at various theoretical levels. As shown in Table 5, the CCSD(T)-CBS results of Pelucchi et al. [13] are in agreement with the current CCSD(T)-F12a results with standard frequency scaling factors within 0.2 kcal/mol; because they did not account for the greater anharmonicity at the transition states, their enthalpy of activation for abstraction at methyl is 0.8 kcal/mol higher than our best estimate. The G4 method used by Li et al. [12] and the G3 and CBS-QB3 methods used by Seta et al. [8] also overestimate the enthalpy of activation of the hydrogen abstraction reaction from methyl group as compared to our best estimates.
Wu et al. [57] used M06-2X, G3(MP2)-RAD//M06-2X and ROCBS-QB3//M06-2X calculations to calculate the 0 K enthalpies of addition reactions. As shown in Table 6, most of their results overestimate the enthalpy of activation as compared to our best estimates. Their M06-2X and G3(MP2)-RAD calculations predict the same order as we do for the enthalpies of activation for addition reactions, with the orthosite preferred. However, their ROCBS-QB3 results give a different order and predict the meta addition to have the lowest enthalpy of activation at 0 K.

Multistructural Anharmonicity.
We included the effect of multistructural anharmonicity due to molecular torsions by the multistructural torsion (MS-T) method, which corrects the rate constants obtained in the single-structure quasiharmonic approximation by multiplying by the multistructural anharmonic factor of the reaction. In the MS-T calculations, the coupled torsional potential is used for the transition states of ipso addition and ortho and methyl abstraction and for the products of ipso and ortho-addition because in these cases, the torsions that introduce the multistructural torsional anharmonicity are coupled with each other because of steric effects. For the other complex species, the uncoupled torsional potential is used because the torsions considered are only weakly coupled with each other. The calculated multistructural anharmonic factors MS-T act as functions of temperature are shown in Figure 4. The figure shows that the effect of multistructural anharmonicity usually increases with temperature, and it can change the calculated reaction rates by factors as small as 0.2 and as large as 4.2.

Transmission Coefficients.
The recrossing transmission coefficient Γ and the tunneling transmission coefficient 푆퐶푇 are important components of the dynamics calculations, and they are shown for each reaction (in the high-pressure limit) in Figures 5 and 6, respectively. Figure 5 shows that the most important recrossing coefficients are for H-abstraction from methyl and from the ortho site. These recrossing transmission coefficients can decrease the rate constants by more than a factor of two. Figure 6 shows that the tunneling transmission coefficient ranges from 0.35 to 1.9, where a value less than unity corresponds to nonclassical reflection dominating tunneling in the thermal average, and a value greater than unity corresponds to tunneling dominating nonclassical reflection in the thermal average. At 275 K, the tunneling transmission coefficient increases the ratio of abstraction from methyl to abstraction from the para position by a factor of 4.4, and the effect is even larger at lower temperatures.  . Figures 7(a) and 7(b) show the temperature-dependent HPL rate constants calculated by the MS-CVT/SCT method for the hydrogen abstraction reactions and addition reactions, respectively, and for comparison experimental data is also displayed in the figure. The hydrogen abstraction reactions dominate at high temperatures, whereas the addition reactions dominate at low temperatures. Therefore, we plot the experiment data at 213 K-298 K and 504 K-1481 K as rate constants of the addition reactions and the abstraction reactions, respectively. Figure 7(a) shows that our total abstraction rate constants agree well with the experimental results of Seta et al. [8] and Tully et al. [7] for the mid-temperature range (500 K < T < 1000 K) and are slightly higher (by a factor of 1.4) than the experimental data in the temperature interval 1000-1500 K. The calculations of Li et al. [12] and Uc et al. [10] are in good agreement with each other, but their results are lower than our results and those estimated by Pelucchi et al. [13] because both Li et al. and Uc et al. overestimated the barrier height for hydrogen abstraction from the methyl group as we saw in Table 5.

High-Pressure-Limit (HPL) Rate Constants: Abstraction.
The comparison of our abstraction rate constants with the work of Pelucchi et al. is particularly interesting because their work is in some respects similar to ours. They used conventional transition state theory with a correction for variationaltransition-state-location effects computed by M06-2X/6-311+G(d,p). They computed the classical barrier height by CCSD(T)/aug-cc-pVTZ//M06-2X/6-311+G(d,p) with a correction for basis set effects, and they calculated frequencies by M06-2X/6-311+G(d,p) with no mention of frequency scaling. (The 6-311+G(d,p) basis set they used is similar to the MG3S basis set used here for direct dynamics because for systems containing only C, H, and O, the MG3S basis set is equivalent to 6-311+G(2df,2p), which has more polarization functions than 6-311+G(d,p).) They included torsional anharmonicity by one-dimensional and two-dimensional models applied only at the transition state. They used an asymmetric Eckart model to evaluate the tunneling transmission coefficient. (The Eckart model is less reliable than the SCT method employed here because it does not include the true shape of the effective potential energy barrier, and it neglects cornercutting tunneling.) The net result of their calculation is in good agreement with our results (apparently due to some cancellation among the various factors that differ) at most temperatures, but there is a significant deviation for T < 450 K. Figure 8 shows the effects of frequency scaling factors. The figure clearly shows that using the standard scaling factors ZPE greatly underestimates the rate constants; for example, the HPL rate constant calculated with the standard ZPE factors is lower than that obtained with the SRP factor ZPE-SRP by a factor of 0.72 at 1600 K, by a factor of 0.40 at 500 K, and by a factor of 0.13 at 213 K. Perry et al. 1977 1000/T (＋ -1 ) Addition Figure 7: The reaction rate constants in the high-pressure limit: (a) sum of rate constants for abstraction reactions; (b) sum of rate constants for addition reactions. The SRP scaling factor is used for vibrational frequencies of transition states. For all reactions, tunneling is allowed to occur at all energies above the ZPE of the pre-reactive van der Waals complex. Some experimental data (symbols) and previous theoretical results (lines) are also given for comparison. Figure 8 also shows the variational-transition-statelocation effects. We find a variational-transition-statelocation effect increasing with temperature up to about 30% (i.e., decreasing the calculated rate constant by about 30%) at high T; this is in good agreement with the calculations of Pelucchi et al.
The most important conclusion from Figure 8 is that if we neglect variational-transition-state-location effects or do not account for the transition state anharmonicity being greater than predicted by the standard model, the agreement with experiment is noticeably degraded.

High-Pressure-Limit (HPL) Rate Constants: Addition.
The total HPL rate constants of addition reactions are plotted in Figure 7(b), and they agree reasonably with the experimental data of Tully et al. at T = 213-298 K. Figure 9 illustrates the calculated pressure dependence of the total rate constant for addition reaction with Ar as the bath gas. We see from (7) that the collision efficiency increases with increasing ⟨Δ ⟩, and consequently, the collision rate constant 푐 increases with increasing ⟨Δ ⟩. The addition rate constant therefore increases with increasing ⟨Δ ⟩ because the adducts are more readily stabilized. These expectations are consistent with the values in the figure, where the solid curves are for the experimental value of ⟨Δ ⟩, namely 130 cm −1 , and we also show what the results would be if ⟨Δ ⟩ were a factor-oftwo larger or a factor-of-two smaller. The effect of varying ⟨Δ ⟩ is relatively small; for example, at p = 100 Torr and T = 500 K, when ⟨Δ ⟩ is halved to 65 cm −1 , the rate constant is decreased by a factor of 0.90 as compared to that obtained by using ⟨Δ ⟩ = 130 cm −1 ; and using two times the ⟨Δ ⟩ value, that is 260 cm −1 , the rate constant is increased by a factor of 1.08.

Pressure-Dependent Rate Constants.
We see that the pressure dependence sets to an appreciable extent for T > ∼330 K, and the fall-off becomes steeper at temperatures higher than 500 K. The high-temperature  results in Figure 9 look nothing like the high-temperature results in Figure 7(b) because 1.2 atm (the highest pressure in Figure 9) is far from the high-pressure limit at high T. However, at low T, all the curves in Figure 9 are at the highpressure limit. Figure 10 illustrates the calculated pressure dependence of the total rate constant for addition reaction over a wider pressure range; this figure is for Ar as the bath gas, and it shows results only for the experimental value of ⟨Δ ⟩. Figure 10 shows that at T = 298.15 K, there is almost no pressure dependence for pressures larger than 100 Torr (log 10 p/bar = -0.9), and this result is in agreement with experimental observations [7]. The experiments of Tully in Figure 7(b) are for T ≤ 298 K and for log 10 p/bar in the range -1.6 to -0.6. Figure 10 shows that this is in the high-pressure limit, so the comparison of HPL rate constants to experiment in Figure 7(b) is warranted. Figure 11 shows the final total reaction rate constants calculated as the sum of the rate constants of R1-R8 reactions at several pressures, and it compares the calculations to the available reliable experimental data; Table 7 lists a few specific values.

Total Rate Constants.
In the high-temperature region (T > 1000 K), the calculated total reaction rate constant at the experimental pressures (20 Torr to ∼1 atm) is dominated by the hydrogen abstraction reaction (which is pressure-independent in our PEM model), and the current calculations agree well with the available experiments. As given in Table 7, we overestimate the rate constants at 500 K and 568 K by a factor of ∼1.5-1.7 as compared to those obtained by Tully in Ar, but only slightly underestimate the rate constants at 298.15 K. At 230-300 K, we reproduce the experimental data very well. At very low T, below 230 K, we overestimate the experimental results by as much as a factor of 2.5.
In the mid-temperature range (700-350 K), our calculations reveal the nonmonotonic temperature dependence and the strong pressure dependence of the total rate constants (at 500 K and 10 Torr, the rate constant is calculated to be 31% lower than the HPL). The nonmonotonic temperature dependence is a result of competition of the addition reactions and hydrogen abstraction reactions. The experimental results are very uncertain in this region.

Branching
Fractions. An advantage of theoretical studies is that they can yield rate constants for individual products,  whereas experimental studies are usually based on loss of reactant and so yield only total rate constants. Figure 12 shows the total addition fraction at various pressures in Ar, where this branching fraction is defined as the sum of the addition rate constants divided by the sum of all the rate constants (abstraction and addition). In the low-temperature region (213 K < T < 300 K) that is important for atmospheric chemistry, the reaction is mainly dominated by the addition reactions with a total branching fraction larger than ∼0.7. Figure 12 also shows that in the HPL the abstraction reaction dominates for T > 500 K. At lower pressure the switch in dominance occurs at lower temperatures, for example for 10 Torr it occurs at 400 K. Figure 13(a) shows HPL branching fractions for each of the eight reactions. It shows that the hydrogen abstraction reactions from the methyl group and meta site are the dominant reactions at high temperatures, and at low temperatures, the ortho-addition reaction has a significant dominance. We also notice the significant contribution of hydrogen abstraction from the methyl at low T. This is explained by the obtained negative (near-zero) enthalpy of activation at 0 K when the anharmonicity of the high-frequency modes of transition state is considered by using a SRP scaling factor of frequency, as shown in Table 5. At 298 K, we predict the contribution of methyl hydrogen abstraction to be ∼30%, which is larger than the branching ratio (7%) used in some atmospheric chemistry models [57], for example, MCM3.2 and SAPRC-11. The present results could be used to update those models. Figures 13(b)-13(d) show branching fractions at p = 1.2 atm, 100 Torr, and 10 Torr with Ar as bath gas. We observe significant pressure effects on the fractions of the two dominant reactions, addition at the ortho-site and abstraction from the methyl site. The main effect of decreasing the pressure from the HPL to 1.2 atm is that the addition reaction essentially turns off above 1000 K, and abstraction from the methyl increases significantly for temperatures above 600 K. The peak of the branching fraction curve for abstraction from the methyl site increases and moves to higher temperature; the maximum increases from 0.45 at ∼550 K (in HPL) to 0.78 at ∼650 K (at p = 10 Torr). With the decrease of pressure, the temperature where the dominant reaction shifts from orthoaddition and abstraction at methyl switches shifts from ∼400 K in HPL to 325 K at 10 Torr. At 10 Torr, the addition reactions have essentially turned off by 700 K.

Summary
In this work, we present a theoretical and computational study of the kinetics of toluene reaction with OH over a wide temperature range from 213 K to 2500 K. We focus on hydrogen abstraction and addition reaction channels (R1-R8), which dominate, respectively, at high and low temperatures. The high-pressure limit rate constants of all channels have been estimated with multistructural canonical variational transition state theory with the small-curvature tunneling approximation (MS-CVT/SCT) by using the multistructural torsional anharmonicity (MS-T) method to take account of anharmonicity resulting from torsional motions and multiple conformers. We also corrected the anharmonicity of highfrequency modes of the transition states by using the specificreaction-parameter (SRP) scaling factors for vibrational frequencies.
The pressure-dependent rate constants of addition reactions at several pressures have been estimated by the SS-QRRK method. The final total reaction rate constants calculated as a sum of the rate constants of reactions R1-R8 are compared to the available experimental data. We get good agreement with those experimental results in the temperature region in which the experiments appear to be most reliable. We predict the very weak pressure dependence of the rate constants for T < 300 K, but obtain strongly pressure-dependent reaction rates at higher temperatures. We also observe significant contributions of hydrogen abstraction from the methyl group even at low temperature (it accounts for ∼30% and ∼20% of the total rate at 298 K and 213 K, respectively), and that could be important for understanding the fate of alkyl aromatics in the atmosphere. energy 푎 ( ) (kcal/mol) of the high-pressure-limit rate constants of reverse reactions of R5-R8. These are obtained using the fitting parameters of k -1 . Table S4: energy dependence factor of the density of states ( 퐸 ) as a function of temperature as calculated by the Whitten-Rabinovitch method.