Artificial Intelligence and Machine Learning Approaches for Predicting Hydrogen Abstraction Reactivity of Toluene by the OH/NH2/CH3 Radicals

Hydrogen abstraction reactions constitute a fundamental class of elementary processes in combustion chemistry, atmospheric oxidation, and pyrolysis of hydrocarbon fuels. Aromatic compounds such as toluene are widely used as model systems due to their structural simplicity and relevance to real fuels. In particular, hydrogen abstraction from different sites of the toluene molecule, including benzylic and ring positions, plays a critical role in determining the formation of key intermediates such as benzyl radicals and subsequent oxidation products. Accurate prediction of reaction rate constants for such processes is therefore essential for constructing reliable kinetic models that can describe complex reacting systems.^1,2 Traditionally, the investigation of gas-phase reaction mechanisms relies on high-level quantum chemical calculations combined with statistical rate theories. Methods such as coupled-cluster theory, variational transition state theory (VTST),³ and Rice–Ramsperger–Kassel–Marcus (RRKM) theory⁴ have been successfully employed to determine potential energy surfaces (PESs) and temperature-dependent rate constants. However, these approaches are computationally expensive, especially for systems involving multiple reaction pathways and large molecular structures. Moreover, the construction of detailed kinetic models often requires the estimation of thousands of rate coefficients, making the application of purely ab initio approaches impractical for large-scale reaction networks.

In recent years, artificial intelligence (AI) and machine learning (ML) have emerged as transformative tools in chemical kinetics and reaction mechanism studies.^5,6 These approaches enable the rapid prediction of key chemical properties, including activation energies, rate constants, and even entire potential energy surfaces, by learning from existing datasets. Unlike traditional methods that rely solely on solving the Schrödinger equation, ML models can identify complex, nonlinear relationships between molecular descriptors and reactivity, thereby significantly reducing computational cost. It has been demonstrated that ML-based methods can accurately estimate reaction rate coefficients across diverse reaction families and can be scaled efficiently to large datasets.^7,8

The application of AI/ML in gas-phase reaction kinetics is particularly promising. Recent studies have shown that machine learning models can be trained to predict activation barriers and reaction rates with near ab initio accuracy, enabling rapid screening of reaction pathways and exploration of chemical space.^9,10 In addition, ML techniques have been successfully integrated with kinetic modeling frameworks to predict reactivity trends in networks of coupled gas-phase reactions, providing insights into the propagation of uncertainty and the influence of individual rate constants on overall system behavior. Furthermore, the use of ML-based interatomic potentials and active learning strategies has enabled the automated construction of high-quality PESs, allowing efficient simulation of reaction dynamics without the need for extensive quantum chemical calculations.

Beyond kinetics prediction, AI has also demonstrated strong capabilities in uncovering reaction mechanisms. By analyzing large datasets of chemical reactions, ML models can identify dominant reaction pathways, classify reaction types, and even suggest new mechanistic hypotheses. In complex systems such as combustion or atmospheric chemistry, where numerous competing pathways exist, these capabilities are particularly valuable. AI-driven approaches can thus complement traditional theoretical methods by providing rapid, data-driven insights into reaction mechanisms while maintaining reasonable accuracy.¹¹

Despite these advances, the application of AI/ML to hydrogen abstraction reactions involving aromatic systems remains relatively limited. In particular, comparative studies of different radical species, such as amino and methyl radicals, using a unified framework that combines high-level quantum chemistry and machine learning are still scarce. Understanding the differences in reactivity between these radicals is essential for improving kinetic models relevant to combustion and atmospheric processes.

In this work, we present an integrated approach combining quantum chemical calculations and machine learning techniques to investigate hydrogen abstraction reactions of toluene with NH₂ and CH₃ radicals in the gas phase. By constructing a comprehensive dataset of reaction descriptors and employing advanced ML models, we aim to develop accurate and computationally efficient predictive models for reaction rate constants. This study not only provides detailed mechanistic and kinetic insights into the reactivity of toluene but also demonstrates the potential of AI-driven approaches as a powerful tool for accelerating research in gas-phase chemical kinetics.

COMPUTATIONAL METHODS

Electronic structure calculations

The geometries of all stationary points involved in the hydrogen abstraction reactions of toluene with hydroxyl (OH), amino (NH₂) and methyl (CH₃) radicals were determined using a hybrid computational protocol combining density functional theory (DFT) calculations and machine learning (ML) assisted optimization strategies. Initial structural guesses and conformational searches were accelerated using AI/ML-based algorithms to efficiently explore the potential energy landscape of the reacting system. The resulting structures were subsequently refined using quantum chemical calculations. All geometry optimizations and harmonic vibrational frequency calculations were performed at the M06-2X/aug-cc-pVTZ level of theory. The M06-2X functional is a hybrid meta-GGA functional specifically designed to provide reliable predictions for main-group thermochemistry, reaction kinetics, and noncovalent interactions, and it has been shown to perform well for reaction barrier heights and thermochemical properties.¹²

Frequency calculations at the same level were carried out to verify the nature of each stationary point on the potential energy surface (PES). Local minima were confirmed by the absence of imaginary frequencies, whereas transition states exhibited a single imaginary frequency corresponding to the reaction coordinate associated with hydrogen abstraction. Intrinsic reaction coordinate (IRC) calculations¹³ were further performed to ensure that each transition state properly connects the corresponding reactants and products. To obtain more accurate energetics, single-point energy calculations were performed at the CCSD(T)/aug-cc-pVTZ level using the optimized geometries obtained at the M06-2X level. The harmonic frequencies were scaled by a factor of 0.985 to account for anharmonic effects and systematic overestimation of vibrational frequencies.

The reliability of the single-reference coupled-cluster treatment was assessed using the T₁ diagnostic computed at the CCSD(T)/aug-cc-pVTZ level. The obtained values indicate that all closed-shell species exhibit T₁ values smaller than the commonly accepted threshold of 0.02, while open-shell structures remain within the acceptable limit of 0.045,¹⁴ confirming that the single-reference wavefunction description is adequate for the present system. In addition, spin contamination was examined for all open-shell species. The calculated ⟨S²⟩ values are close to the theoretical value of 0.75 for doublet states, indicating negligible spin contamination in the wavefunctions.

Kinetic modeling: Rate Constant Calculations

The rate constants for the hydrogen abstraction reactions were calculated using transition state theory (TST). The TST calculations were performed with the ChemRate program package,¹⁵ which allows the evaluation of temperature-dependent rate coefficients based on statistical rate theory and molecular properties obtained from electronic structure calculations.

All molecular partition functions were computed using the rigid-rotor harmonic-oscillator approximation. Vibrational frequencies obtained at the M06-2X level were used to evaluate vibrational partition functions and thermodynamic corrections. For the hydrogen abstraction transition states, the imaginary frequency associated with the reaction coordinate was excluded from the partition function calculation.

Quantum mechanical tunneling effects were accounted for using the Eckart tunneling correction, which provides a realistic description of barrier penetration for hydrogen transfer reactions. This correction is particularly important for reactions involving light atoms such as hydrogen, where tunneling can significantly enhance the reaction rate, especially at low temperatures. The rate constants were calculated over a temperature range relevant to combustion and atmospheric chemistry conditions. For the hydrogen abstraction reactions considered in this work, the rate coefficients primarily depend on the barrier heights and vibrational properties of the transition states, since these reactions proceed through a single transition state without the formation of intermediate complexes.

RESULTS AND DISCUSSION

The potential energy surface (PES) for the C₆H₅CH₃ + OH/NH₂/CH₃ systems has been characterized at the CCSD(T)//M06-2X/aug-cc-pVTZ level of theory to account for the possible H-abstraction pathways. The energy profile (in kcal/mol, relative to the reactants as zero) is shown in Figure 1.

The hydrogen abstraction reactions between toluene and the OH, NH₂, and CH₃ radicals may proceed through two distinct pathways: abstraction of a hydrogen atom from the methyl substituent (benzylic H-abstraction) or abstraction of a hydrogen atom directly from the aromatic ring (ring H-abstraction) at the ortho, meta, or para positions. These two pathways exhibit significantly different energetic characteristics on the potential energy surface (PES) and therefore lead to different kinetic and thermodynamic behaviors. In general, hydrogen abstraction from the benzylic position is expected to be more favorable because the resulting benzyl radical is strongly stabilized by resonance delocalization over the aromatic π system, which substantially lowers the reaction barrier.

Figure 1. Mechanisms of the C₆H₅CH₃ reactions with the OH/NH₂/CH₃ radicals. Energies (kcal/mol) were calculated at the CCSD(T)//M06-2X/aug-cc-pVTZ level of theory.