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.
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–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 .)
The first experimental study of the kinetics of the OH + toluene reaction was carried out by Davis et al.  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.  (304 K, 1 atm), Hansen et al.  (298 K, 100–618 Torr), and Perry et al.  (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.  (213–1150 K, 20–200 Torr). They analyzed their experiments under the 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-first-order 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.  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.  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.  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.  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.  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.  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 pressure-dependent 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 . 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 . 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.  found the rate constant to be pressure-dependent at 298 K and 3-100 Torr total pressure; but neither Hansen  nor Perry  found significant pressure dependence at 298 K and above 100 Torr. Tully et al.  reconfirmed their work by finding that the third-body dependence was only 9% at 298 K and 100 Torr in He, Ar, or SF6; 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.
2. 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.
2.1. 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 . 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  /jun-cc-pVTZ [19–22] method. The standard assumptions are that closed-shell species with > ~0.02 and open-shell species with > ~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  was used for the coupled cluster calculations.
2.2. Selection of KS Model Chemistry
In our experience, the M06-2X , M08-SO , and MN15  density functionals combined with the MG3S  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 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.
2.3. Decomposition of Adducts
Previous work [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. . Therefore these reactions are not major pathways, and we did not consider them further.
3. 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 . In particular, for the temperature-dependent high-pressure-limit (HPL) rate constants we used multistructural  canonical variational transition state theory (CVT)  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 . 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 variational-transition-state-location effect or sometimes simply as a variational effect.
In these calculations, we used the multistructural torsional anharmonicity method (MS-T)  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-reaction-parameter (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 .
The pressure-dependent rate constants of the addition reactions were calculated by the system-specific quantum Rice-Ramsperger-Kassel  (SS-QRRK) method.
3.1. MS-CVT/SCT Theory
The temperature-dependent high-pressure-limit rate constants calculated by the MS-CVT/SCT theory are expressed by where 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, is the conventional transition state theory rate constant, 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 is defined  as the multistructural anharmonicity factor of the reaction. The factor is calculated as the ratio of multistructural anharmonicity factor of the transition state (α = TS) to the multistructural anharmonicity factor of the reactant (α = R).
In the calculation of , 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 .
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 .
The standard vibrational frequency scaling factor  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 specific-reaction-parameter scaling factor 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 and correct, respectively, the anharmonicity of the ZPE and the inaccuracy of the model chemistry. Note that is characteristic of the model chemistry and is characteristic of the species. For all species, we use the harmonic factor that has been parametrized to obtain the accurate harmonic frequency in the F38/10 database . We calculate for a given species as the ratio of anharmonic ZPE estimated by hybrid , degeneracy-corrected , second-order vibrational perturbation theory (HDCVPT2) [45–47] to harmonic ZPE calculated for that species with the same model chemistry. In the calculations of , 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 equal to the standard for reactants, pre-reactive van der Waals complexes, and products.
Table 3 lists the calculated scaling factors of the transition states for all reactions in the present study, as well as the factors and the standard factors. For hydrogen abstraction transition states calculated with M06-2X/MG3S, the table shows that the specific-reaction-parameter scaling factors ( = 0.960–0.966) are smaller than the standard scaling factor ( = 0.970). For the addition transition states the calculated scaling factors are the same as the standard factor (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.
3.2. The SS-QRRK Method
Addition reactions R5-R8 in Figure 1 are treated by the chemical activation mechanism: where T denotes the toluene molecule, T is temperature, and TOH is the energized adduct. Applying the steady-state approximation to the energized adduct TOH yields the phenomenological bimolecular stabilization rate constant, 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 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 high-pressure-limit rate constant , the energy-dependent dissociation rate constant , and collisional deactivation rate constant . These are evaluated as follows:
(a) is evaluated by the MS-CVT/SCT method explained above.
(b) The high-pressure-limit dissociation rate constants is also evaluated by the MS-CVT/SCT method.
(c) We calculate a temperature-dependent Arrhenius preexponential factor and Arrhenius activation energy (T) by
(d) 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 incorporates variational effects, multidimensional quantum mechanical tunneling, and multistructural anharmonicity, they are also inherited in .
(e) is estimated by collision theory as the product of the Lennard-Jones collision rate constant 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  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 : 130 cm−1 for Ar and 75 cm−1 for He. In the calculations of , we based the Lennard-Jones parameters on those used in reference  to derive the values of ; in particular, we used ε/ = 410 K and σ = 6.0 Å for the adducts of OH and toluene , 120 K and 3.4 Å for Ar , and 10 K and 2.55 Å for He .
The Polyrate 2016  and Gaussrate 17  programs were used for the direct CVT/SCT calculations and SS-QRRK calculations, and the MSTor program  was used to calculate the conformational-rotational-vibrational partition functions.
3.3. 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.
4. Results and Discussion
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, .
We found only one conformer for the pre-reactive complex (see Figure 3(a).) Uc et al.  reported a -symmetry structure for the pre-reactive complex by the BHandHLYP/6-311++ method. However, this -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-pressure-limit rate constants of the hydrogen abstraction reactions.
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.
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 high-frequency modes of transition states of hydrogen abstraction reactions, the enthalpy of activation is overestimated with the standard scaling factor or harmonic scaling factor as compared to that obtained with the SRP scaling factor .
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.  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.  and the G3 and CBS-QB3 methods used by Seta et al.  also overestimate the enthalpy of activation of the hydrogen abstraction reaction from methyl group as compared to our best estimates.
Wu et al.  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 ortho-site 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.
4.3. 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 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.
4.4. 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.
4.5. High-Pressure-Limit (HPL) Rate Constants
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.
4.5.1. High-Pressure-Limit (HPL) Rate Constants: Abstraction
Figure 7(a) shows that our total abstraction rate constants agree well with the experimental results of Seta et al.  and Tully et al.  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.  and Uc et al.  are in good agreement with each other, but their results are lower than our results and those estimated by Pelucchi et al.  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.
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 variational-transition-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 corner-cutting 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 greatly underestimates the rate constants; for example, the HPL rate constant calculated with the standard factors is lower than that obtained with the SRP factor 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.
Figure 8 also shows the variational-transition-state-location effects. We find a variational-transition-state-location 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.
4.5.2. 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.
4.6. Pressure-Dependent Rate Constants
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-of-two 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.
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 high-pressure 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 (log10 p/bar = –0.9), and this result is in agreement with experimental observations . The experiments of Tully in Figure 7(b) are for T ≤ 298 K and for log10 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.
4.7. Total Rate Constants
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.
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.
4.8. 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 , 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 ortho-addition 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.
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 high-frequency modes of the transition states by using the specific-reaction-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.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.
This work was supported in part by the National Natural Science Foundation of China (91641127 and 91841301) and the U.S. Department of Energy (Grant DE-SC0015997).
Table S1: fitting parameters for the high-pressure-limit rate constants of the forward reactions of R1-R8. Table S2: fitting parameters for the high-pressure-limit rate constants (k–1) of the reverse reactions of R5-R8. Table S3: the Arrhenius pre-exponential factor (s−1) and Arrhenius activation 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. Table S5: collision efficiency and collisional deactivation rate constant ( ) as functions of temperature. (Supplementary Materials)
- B. T. Brem, L. Durdina, F. Siegerist et al., “Effects of fuel aromatic content on nonvolatile particulate emissions of an in-production aircraft gas turbine,” Environmental Science & Technology, vol. 49, no. 22, pp. 13149–13157, 2015.
- J. Noda, R. Volkamer, and M. J. Molina, “Dealkylation of alkylbenzenes: a significant pathway in the toluene, o-, m-, p-xylene + OH reaction,” The Journal of Physical Chemistry A, vol. 113, no. 35, pp. 9658–9666, 2009.
- D. D. Davis, W. Bollinger, and S. Fischer, “A kinetics study of the reaction of the OH free radical with aromatic compounds. I. Absolute rate constants for reaction with benzene and toluene at 300°K,” The Journal of Physical Chemistry C, vol. 79, no. 3, pp. 293-294, 1975.
- G. J. Doyle, A. C. Lloyd, K. R. Darnall, A. M. Winer, and J. N. Pitts, “Gas phase kinetic study of relative rates of reaction of selected aromatic compounds with hydroxyl radicals in an environmental chamber,” Environmental Science & Technology, vol. 9, no. 3, pp. 237–241, 1975.
- D. A. Hansen, R. Atkinson, and J. N. Pitts Jr., “Rate constants for the reaction of hydroxyl radicals with a series of aromatic hydrocarbons,” The Journal of Physical Chemistry, vol. 79, pp. 1763–1766, 1975.
- R. A. Perry, R. Atkinson, and J. N. Pitts Jr., “Kinetics and mechanism of the gas phase reaction of hydroxyl radicals with aromatic hydrocarbons over the temperature range 296-473 K,” The Journal of Physical Chemistry, vol. 81, pp. 296–304, 1977.
- F. P. Tully, A. R. Ravishanka, R. L. Thompson et al., “Kinetics of the reactions of hydroxyl radical with benzene and toluene,” The Journal of Physical Chemistry, vol. 85, pp. 2262–2269, 1981.
- T. Seta, M. Nakajima, and A. Miyoshi, “High-temperature reactions of OH radicals with benzene and toluene,” The Journal of Physical Chemistry A, vol. 110, no. 15, pp. 5081–5090, 2006.
- V. Vasudevan, D. F. Davidson, and K. Hanson, “High-temperature measurements of the reactions of OH with toluene and acetone,” The Journal of Physical Chemistry A, vol. 109, pp. 3352–3359, 2005.
- V. H. Uc, J. R. Alvarez-Idaboy, A. Galano, I. García-Cruz, and A. Vivier-Bunge, “Theoretical determination of the rate constant for OH hydrogen abstraction from toluene,” The Journal of Physical Chemistry A, vol. 110, no. 33, pp. 10155–10162, 2006.
- V. H. Uc, J. R. Alvarez-Idaboy, A. Galano, and A. Vivier-Bunge, “Theoretical explanation of nonexponential OH decay in reactions with benzene and toluene under pseudo-first-order conditions,” The Journal of Physical Chemistry A, vol. 112, no. 33, pp. 7608–7615, 2008.
- S.-H. Li, J.-J. Guo, R. Li, F. Wang, and X.-Y. Li, “Theoretical prediction of rate constants for hydrogen abstraction by OH, H, O, CH3, and HO2 radicals from toluene,” The Journal of Physical Chemistry A, vol. 120, no. 20, pp. 3424–3432, 2016.
- M. Pelucchi, C. Cavallotti, T. Faravelli, and S. J. Klippenstein, “H-Abstraction reactions by OH, HO2, O, O2 and benzyl radical addition to O2 and their implications for kinetic modelling of toluene oxidation,” Physical Chemistry Chemical Physics, vol. 20, no. 16, pp. 10607–10627, 2018.
- S.-C. Lin, T.-C. Kuo, and Y.-P. Lee, “Detailed rate coefficients and the enthalpy change of the equilibrium reaction OH + C6H6 HOC6H6 over the temperature range 345-385 K,” The Journal of Chemical Physics, vol. 101, pp. 2098–2105, 1994.
- M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., Gaussian 09, Revision C.01, Gaussian, Inc, 2010.
- Y. Zhao, R. Peverati, S. Luo et al., MNGFM6.7: Minnesota Gaussian Functional Module, Version 6.7, University of Minnesota, Minneapolis, Menn, USA.
- X. Xu, I. M. Alecu, and D. G. Truhlar, “How well can modern density functionals predict internuclear distances at transition states?” Journal of Chemical Theory and Computation, vol. 7, no. 6, pp. 1667–1676, 2011.
- P. J. Knowles, C. Hampel, and H. J. Werner, “Coupled cluster theory for high spin, open shell reference wave functions,” The Journal of Chemical Physics, vol. 99, pp. 5219–5227, 1993, Erratum: The Journal of Chemical Physics, vol. 112, pp. 3106-3107, 2000.
- T. H. Dunning Jr., “Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen,” The Journal of Chemical Physics, vol. 90, no. 2, pp. 1007–1023, 1989.
- E. Papajak and D. G. Truhlar, “Convergent partially augmented basis sets for post-Hartree-Fock calculations of molecular properties and reaction barrier heights,” Journal of Chemical Theory and Computation, vol. 7, no. 1, pp. 10–18, 2011.
- R. A. Kendall, T. H. Dunning Jr., and R. J. Harrison, “Electron affinities of the first-row atoms revisited: Systematic basis sets and wave functions,” The Journal of Chemical Physics, vol. 96, no. 9, pp. 6796–6806, 1992.
- E. Papajak, J. Zheng, X. Xu, H. R. Leverentz, and D. G. Truhlar, “Perspectives on basis sets beautiful: Seasonal plantings of diffuse basis functions,” Journal of Chemical Theory and Computation, vol. 7, no. 10, pp. 3027–3034, 2011.
- T. B. Adler, G. Knizia, and H. J Werner, “A simple and efficient CCSD(T)-F12 approximation,” The Journal of Chemical Physics, vol. 127, Article ID 221106, 2007.
- Y.-P. Liu, G. C. Lynch, T. N. Truong, D.-H. Lu, D. G. Truhlar, and B. C. Garrett, “Molecular modeling of the kinetic isotope effect for the [1,5] sigmatropic rearrangement of cis-1,3-pentadiene,” Journal of the American Chemical Society, vol. 115, no. 6, pp. 2408–2415, 1993.
- H.-J. Werner, P. J. Knowles, G. Knizia et al., Molpro, Version 2012.1, University of Birmingham, Birmingham, UK, 2012.
- Y. Zhao and D. G. Truhlar, “Density functionals with broad applicability in chemistry,” Accounts of Chemical Research, vol. 41, no. 2, pp. 157–167, 2008.
- Y. Zhao and D. G. Truhlar, “Exploring the limit of accuracy of the global hybrid meta density functional for main-group thermochemistry, kinetics, and noncovalent interactions,” Journal of Chemical Theory and Computation, vol. 4, no. 11, pp. 1849–1868, 2008.
- H. S. Yu, X. He, S. L. Li, and D. G. Truhlar, “MN15: A Kohn-Sham global-hybrid exchange-correlation density functional with broad accuracy for multi-reference and single-reference systems and noncovalent interactions,” Chemical Science, vol. 7, no. 8, pp. 5032–5051, 2016.
- B. J. Lynch, Y. Zhao, and D. G. Truhlar, “Effectiveness of diffuse basis functions for calculating relative energies by density functional theory,” The Journal of Physical Chemistry A, vol. 107, no. 9, pp. 1384–1388, 2003.
- J. L. Bao, J. Zheng, and D. G. Truhlar, “Kinetics of hydrogen radical reactions with toluene including chemical activation theory employing system-specific quantum RRK theory calibrated by variational transition state theory,” Journal of the American Chemical Society, vol. 138, no. 8, pp. 2690–2704, 2016.
- T. Yu, J. Zheng, and D. G. Truhlar, “Multi-structural variational transition state theory. Kinetics of the 1, 4-hydrogen shift the pentyl radical with torsional anharmonicity,” Chemical Science, vol. 2, pp. 2199–2213, 2011.
- D. G. Truhlar and B. C. Garrett, “Variational transition state theory,” Annual Review of Physical Chemistry, vol. 35, pp. 159–189, 1984.
- D.-H. Lu, T. N. Truong, V. S. Melissas et al., “POLYRATE 4: A new version of a computer program for the calculation of chemical reaction rates for polyatomics,” Computer Physics Communications, vol. 71, no. 3, pp. 235–262, 1992.
- B. C. Garrett and D. G. Truhlar, “Criterion of minimum state density in the transition state theory of bimolecular reactions,” The Journal of Chemical Physics, vol. 70, no. 4, pp. 1593–1598, 1979.
- J. Zheng and D. G. Truhlar, “Quantum thermochemistry: Multistructural method with torsional anharmonicity based on a coupled torsional potential,” Journal of Chemical Theory and Computation, vol. 9, no. 3, pp. 1356–1367, 2013.
- I. M. Alecu, J. Zheng, Y. Zhao, and D. G. Truhlar, “Computational thermochemistry: Scale factor databases and scale factors for vibrational frequencies obtained from electronic model chemistries,” Journal of Chemical Theory and Computation, vol. 6, no. 9, pp. 2872–2887, 2010.
- J. Zheng, R. Meana-Pañeda, and D. G. Truhlar, “Prediction of experimentally unavailable product branching ratios for biofuel combustion: The role of anharmonicity in the reaction of isobutanol with OH,” Journal of the American Chemical Society, vol. 136, no. 13, pp. 5150–5160, 2014.
- A. Fernandez-Ramos, B. A. Ellingson, B. C. Garrett, and D. G. Truhlar, “Variational transition state theory with multidimensional tunneling,” Reviews in Computational Chemistry, vol. 23, pp. 125–232, 2007.
- H. Zhang, X. Zhang, D. G. Truhlar, and X. Xu, “Nonmonotonic temperature dependence of the pressure-dependent reaction rate constant and kinetic isotope effect of hydrogen radical reaction with benzene calculated by variational transition-state theory,” The Journal of Physical Chemistry A, vol. 121, no. 47, pp. 9033–9044, 2017.
- J. L. Bao and D. G. Truhlar, “Variational transition state theory: Theoretical framework and recent developments,” Chemical Society Reviews, vol. 46, no. 24, pp. 7548–7596, 2017.
- J. L. Bao, R. Meana-Pañeda, and D. G. Truhlar, “Multi-path variational transition state theory for chiral molecules: The site-dependent kinetics for abstraction of hydrogen from 2-butanol by hydroperoxyl radical, analysis of hydrogen bonding in the transition state, and dramatic temperature dependence of the activation energy,” Chemical Science, vol. 6, no. 10, pp. 5866–5881, 2015.
- L. G. Gao, J. Zheng, A. Fernández-Ramos, D. G. Truhlar, and X. Xu, “Kinetics of the methanol reaction with OH at interstellar, atmospheric, and combustion temperatures,” Journal of the American Chemical Society, vol. 140, no. 8, pp. 2906–2918, 2018.
- J. Bloino, M. Biczysko, and V. Barone, “General perturbative approach for spectroscopy, thermodynamics, and kinetics: Methodological background and benchmark studies,” Journal of Chemical Theory and Computation, vol. 8, no. 3, pp. 1015–1036, 2012.
- K. M. Kuhler, D. G. Truhlar, and A. D. Isaacson, “General method for removing resonance singularities in quantum mechanical perturbation theory,” The Journal of Chemical Physics, vol. 104, no. 12, pp. 4664–4671, 1996.
- H. H. Nielsen, “The vibration-rotation energies of molecules and their spectra in the infra-red,” in Atoms III-Molecules I/Atome III-Moleküle I, S. Flügge, Ed., Encyclopedia of Physics, pp. 173–313, Springer, Berlin, Germany, 1959.
- I. M. Mills, “Vibration-rotation structure in asymmetric- and symmetric-top molecules,” in Molecular Spectroscopy: Modern Research, K. N. Rao and C. W. Mathews, Eds., pp. 115–140, Academic, New York, NY, USA, 1972.
- Q. Zhang, P. N. Day, and D. G. Truhlar, “The accuracy of second order perturbation theory for multiply excited vibrational energy levels and partition functions for a symmetric top molecular ion,” The Journal of Chemical Physics, vol. 98, no. 6, pp. 4948–4958, 1993.
- L. S. Kassel, “Studies in homogeneous gas reactions. I,” The Journal of Physical Chemistry C, vol. 32, no. 2, pp. 225–242, 1927.
- A. M. Dean, “Predictions of pressure and temperature effects upon radical addition and recombination reactions,” The Journal of Physical Chemistry C, vol. 89, no. 21, pp. 4600–4608, 1985.
- J. Troe, “Theory of thermal unimolecular reactions at low pressures. II. Strong collision rate constants. Applications,” The Journal of Chemical Physics, vol. 66, no. 11, pp. 4758–4775, 1977.
- H. Hippler, J. Troe, and H. J. Wendelken, “Collisional deactivation of vibrationally highly excited polyatomic molecules. II. Direct observations for excited toluene,” The Journal of Chemical Physics, vol. 78, no. 11, pp. 6709–6717, 1983.
- R. A. Eng, A. Gebert, E. Goos, H. Hippler, and C. Kachiani, “Incubation times, fall-off and branching ratios in the thermal decomposition of toluene: experiments and theory,” Physical Chemistry Chemical Physics, vol. 4, no. 16, pp. 3989–3996, 2002.
- A. Rahman, “Correlations in the motion of atoms in liquid argon,” Physical Review A: Atomic, Molecular and Optical Physics, vol. 136, no. 2A, pp. A405–A411, 1964.
- J. Zheng, J. L. Bao, R. Meana-Pañeda et al., Polyrate, Version 2016-2A, University of Minnesota, Minneapolis, Menn, USA, 2016.
- J. Zheng, J. L. Bao, S. Zhang et al., Gaussrate 17, University of Minnesota, Minneapolis, Menn, USA, 2017.
- J. Zheng, R. Meana-Pañeda, and D. G. Truhlar, “MSTor version 2013: a new version of the computer code for the multi-structural torsional anharmonicity, now with a coupled torsional potential,” Computer Physics Communications, vol. 184, no. 8, pp. 2032-2033, 2013.
- R. Wu, S. Pan, Y. Li, and L. Wang, “Atmospheric oxidation mechanism of toluene,” The Journal of Physical Chemistry A, vol. 118, no. 25, pp. 4533–4547, 2014.