Cite this Article
Peng Ya, Jiang Zhong-An, Chen Ju-Shi. Exploring the methane combustion reaction: A theoretical contribution. Chinese Physics B, 2018, 27(2): 023401  
Permissions
Exploring the methane combustion reaction: A theoretical contribution
Peng Ya †, Jiang Zhong-An ‡, Chen Ju-Shi
School of Civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083, China

 

† Corresponding author. E-mail: pengyaustb@sina.com jza1963@263.net

Abstract

This paper represents an attempt to extend the mechanisms of reactions and kinetics of a methane combustion reaction. Three saddle points (SPs) are identified in the reaction using the complete active space self-consistent field and the multireference configuration interaction methods with a proper active space. Our calculations give a fairly accurate description of the regions around the twin first-order SPs (3A′ and 3A″) along the direction of O(3P) attacking a near-collinear H–CH3. One second-order SP2nd (3E) between the above twin SPs is the result of the C3v symmetry Jahn–Teller coupling within the branching space. Further kinetic calculations are performed with the canonical unified statistical theory method with the temperature ranging from 298 K to 1000 K. The rate constants are also reported. The present work reveals the reaction mechanism of hydrogen-abstraction by the O(3P) from methane, and develops a better understanding for the role of SPs. In addition, a comparison of the reactions of O(3P) with methane through different channels allows a molecule-level discussion of the reactivity and mechanism of the title reaction.

PACS: 34.10.+x;31.15.xv;34.50.Lf
Keyword:combustion reaction;reaction mechanism;rate constant;saddle point
1. Introduction

The reaction of ground-state atomic oxygen with methane is an important initial step of oxidation in combustion processes.[1–5] Since methane is the simplest organic fuel and stored in not only various minerals but also atmosphere as the most abundant hydrocarbon, its combustion has attracted considerable interest, both from experimental and theoretical perspectives.[6–18] The reaction, which represents an important initial step in hydrocarbon combustion, provides an ideal platform for studying bimolecular reactions involving one transfer of a hydrogen atom.[16,17,19] The hydrogen-abstraction by O(3P) is commonly presumed to be the first step to generate a significant concentration of CH3 in flames, which could be oxidized to CO and CO2 in subsequent combustion steps. The exploration of reaction mechanisms and the next computation of related rate parameters have traditionally focused on searching for and studying properties of the saddle point on potential energy surface (PES) of the chemical species of interest.[19–21] A first-order saddle point (SP1st) is a critical point, and it represents a local minimum of PES from all directions/degrees of freedom except for the one that corresponds to the reaction coordinate, thus forming a valley on the PES with respect to which it is a maximum. In contrast to these numerous studies of SP1st, there have been relatively few studies of the second-order saddle point (SP2nd), in whose vicinity the PES was found to exhibit an unusual topology.[22,23] Does any reaction proceed through an SP2nd and give desired products? Or may a chemical transformation proceed through an SP2nd? The literature is rather contradictory concerning such cryptic SP2nd.[22,24,25]

The title reaction has been extensively investigated experimentally, and a considerable body of work was devoted to the experimental characterization of the OH molecule produced in related processes via molecular beams, like the reactions of O atom with saturated, unsaturated, cyclic and aromatic hydrocarbons, in which reactions a small amount of OH radical rotational excitation has been found.[26,27] It was interpreted as resulting from a direct reaction mode and a favored near-collinear O–H–C approach of the O atom to the C–H bond under attack. This point was also confirmed by Ausfelder et al.[11] They used isotope substitution to provide confirmation of the nuclear framework dynamics of this reaction and studied influences on the nascent product OH or OD radicals rotation with fine-structure state partitioning. The absolute rotational energy release is found to be essentially independent of the hydrogen isotope for all observed OH and OD product levels, which indicates a near-collinear abstraction mechanism for this hydrogen-abstraction with the product rotation determined primarily by the geometry at the point of repulsive energy release. As for the other product CH3 radical, the ν2 out-of-plane bending vibrational distribution was determined by Suzuki and Hirota[28] in an infrared diode laser absorption kinetic spectroscopy experiment. Their results suggested that the CH3 moiety was gradually relaxed to a planar structure during these reactions. Laser photolysis of NO2 has been combined with laser-induced fluorescence detection of the nascent OH product to investigate the dynamics of the reactions of O(3P) with methane.[10] The OH spin–orbit (SO) states were found to be nonstatistically populated, which was interpreted by a direct abstraction mechanism with a preferentially collinear C–H–O approach of the O(3P) atom to the C–H bond under attack.

The kinetics of the methane combustion reactions have been studied and the reaction rate coefficients were determined in the temperature range from 298 K to 2500 K. However, it is considerably less reliable due to the uncertainties in the reaction stoichiometry and tunneling effects.[29,30] The experimental data for the reaction was reexamined by Cohen[31] at temperatures below about 600 K. It was found that the values for rate coefficient were overestimated by factors of approximately 2 to 3 in the 250 K–400 K temperature range, with the error increasing as the temperature decreases, and a new expression of the experimental results was proposed as . This expression represents the curvature in Arrhenius plots of the experimental data but leading to higher values of the rate constant than the previous one at low-temperature. This reaction was studied in a crossed-beam experiment, in which the ground-state methyl products were probed using a time-sliced velocity imaging technique, and a direct rebound-type reaction pathway, for which the small impact parameter collisions dominate, was proposed.

Theoretical studies of this reaction have been carried out using various ab initio and dynamical methods;[1,32–34] however, all studies were based on the reaction channel with only one transition state (TS). The ground PES, constructed by Gonzaleza and coauthors[35] using the PMP4 calculations, was employed to obtain the barrier height of 13.3 kcal/mol for SP1st (3A″). The OHC bending angle is 180.0° and the imaginary frequency is 2259i cm−1. In 2012, high accuracy was achieved by using a high-level electron correlation method to obtain the barrier height of 14.08 kcal/mol and the OHC bending angle of 179.3° at the SP1st (3A″), which was taken as a benchmark value.[34] Later, a new full-dimensional analytical PES-2014[32] for O(3P) + CH4 was accurately constructed based on the CCSD(T) = FULL/aug-cc-pVQZ level, and the barrier height obtained was 14.0 kcal/mol with the corresponding OHC bending angle of 180.0°. The thermal rate coefficients of the reaction were also obtained with ring polymer molecular dynamics (RPMD)[15,36–38] and with variational transition state theory calculations on this new surface.[33] The TST (transition state theory) transmission coefficients were obtained by the least-action tunneling (LAT) approximation. In the classical high-temperature limit regime, the RPMD rate theory result agrees with the experiment perfectly, although it overestimates the experimental rate coefficients in the low-temperature regime;[15,36–38] however, there are still somewhat uncertain aspects in the previous theoretical studies. First, rough estimation of the barrier height extracted from experimental data by Arrhenius plots gives 8.7 kcal/mol to 11.7 kcal/mol,[2,33,39,40] which is much lower than previous theoretical reported values.[32,34,35] Second, since the approach of O(3P) along a C–H bond leads to Jahn–Teller conical intersections,[3,41–45] which are not isolated geometries but form a seam. They are widespread in polyatomic systems and make the SPs more and more elusive, representing a theoretical challenge.[16,17,34,46]

One of the reasons is the inadequacy of a single-configuration reference wave function for the post-Hartree–Fock calculations.[46] The geometries and roles of SPs in the reaction CH4 + O(3P) have not been clarified in previous studies. In order to better understand the reaction mechanism and clarify the role of the SPs, recently our group reported an SP1st (3A″) with the ab initio method and kinetic calculation to obtain a fairly accurate description of the regions around the 3A″ TS in the O(3P) attacking a near-collinear H–CH3 direction with a barrier height of 12.53 kcal/mol,[47] which is lower than those reported before.[32,34,35] Until now, there were few accurate ab initio works about the 3A′ state PES reported due to the existence of a conical intersection, which allows the system to relax to lower electronic states and poses a considerable theoretical challenge.[48–50] To accurately study the combustion reaction , we need to consider the first excited state 3A′. The purpose of this work is to investigate the mechanism and kinetics of the reaction CH4 + O(3P), including the SPs and the associated reaction on the excited state 3A′.

2. Methodology
2.1. Ab Initio electronic structure calculations

In this work, we present an extensive theoretical study on hydrogen-abstraction reaction of CH4 + O(3P) based upon the accurate ab initio electronic structure calculations, which are performed using the MOLPRO 2008.1 program package.[51] Since the Cs symmetry is the most relevant for all the stationary point geometries (reactants, products, complexes in the entrance and exit channels, and saddle points), for computational convenience, the reaction system is placed in the Cs symmetry and the internal coordinates are employed in most calculations. The coordinates are defined in Fig. S1 in Appendix A: Supplementary material, in which the bending angle θ is and is the dihedral angle.

Fig. S1. (color online) Coordinate definition geometry of CH4O.

These calculations are based upon the dual-level strategy.[21,52] The details of this computational scheme can be found elsewhere, thus only a brief outline will be given here. As the first step, extensive electronic structure determinations, including the geometry optimizations and vibrational frequencies are performed using the State-Averaged Complete Active Space SCF (SA-CASSCF),[53,54] in which the high level multireference SCF energies and wave functions for a series of ground and excited states are calculated by a state-averaged variational calculation. The selection of the active space is crucial in the CASSCF and icMRCI+Q calculations.[55,56] Here various active spaces in CASSCF calculations have been tested and finally the active space, which consists of twelve orbitals, referred to as (14e, 12o), is chosen for the present calculations. The active space (14e, 12o) is also used for the subsequent icMRCI calculations, and this choice is sufficient to describe the CH bond-breaking and OH bond-forming accurately. It should be noted that the 2s electrons of C and O are put into the closed shell, which means that both 2s orbitals are doubly occupied in all reference configuration state functions, but still are correlated in the icMRCI calculations to account for the core–valence correlations (called CV). The rest of the inner electrons are kept frozen and not correlated. A comparison with the frozen-core and all-electron results for the title reaction have been made, and the result of all-electron is better than that of the frozen-core.[32,34,47,57] The 1s orbitals of carbon and oxygen are not correlated but optimized at the SA-CASSCF level. A reasonable consideration of the basis set dependence is performed. To obtain the CV effects contribution to the energies, especially barriers,[32,58] the cc-pCVDZ basis sets are chosen but modified.[47,59–61]

Since the SA approach included each of the states in question with equal weight and the configuration space is equivalent for all states, SA-CASSCF calculation with 0.5, 0.5 weighting on one root in 3A′ and one in 3A″ symmetries, respectively, is performed using the orbitals obtained with the unrestricted Hartree–Fock method as the starting guess to obtain the orbitals in this work. At the same level, the corresponding frequency calculations of the optimized geometries are also done to prove the characters of minima without imaginary frequencies, the SP1st and SP2nd with only one and two imaginary frequencies, respectively.

The intrinsic reaction coordinate (IRC) calculations are performed by starting from the saddle point geometry at this level in order to identify a given transition state connecting the desired reactant and product.[62] The numerically gradient vector and hessian matrix including rotational partition functions were evaluated at each point along the IRC by harmonic vibrational frequencies calculations. As for the locating of SP2nd, which is also a minimum energy crossing point between 3A′ and 3A″ electronic states, three vectors were evaluated at the SA-CASSCF level: non-adiabatic derivative coupling, gradient of the lower state, and gradient of the upper state.[49,50,63]

In our earlier calculations for the title reaction with SP1st (3A″), only stationary points and each point in the IRC were calculated to evaluate relative energy by using the subsequent icMRCI method, which is for more sophisticated calculations of the single point energy. The calculated potential energy curve (PEC) of ab initio single point energies obtained from higher-level electronic structure calculations might not be as smooth as the counterpart from lower-level IRC calculations, and further rate calculations on the title reaction showed that this does not lead to significant variations of the maximum length of the reaction coordinate. We therefore decided to calculate a few symmetry unique geometries by the icMRCI method in the present work. These geometries were chosen carefully to represent the entrance and exit valley regions most accurately. The icMRCI is a direct MRCI procedure performed in an internally contracted configuration basis, and the MRCI wave functions include all single and double excitations relative to the CASSCF reference functions.[53,54] A large amount of correlation energy is described via single and double electron excitations, and further correlation energy due to higher excitations is approximated by the Davidson correction (+Q).[64] There are at least two reaction channels/transition states (3A′ and 3A″) involved in the title reaction.[25,46] Therefore, two reference states (3A′ and 3A″) are required and used for generating the internally contracted pairs in the icMRCI procedure.

2.2. Kinetic calculations

The canonical unified statistical (CUS) theory,[65–67] which considers a complex region between two bottlenecks and does not give a bound even in classical mechanics, and the canonical variational transition-state (CVT) theory is employed in our present work for the kinetic calculations using the general polyatomic rate constants code POLYRATE 2015 program.[68]

The general equation for rate constant calculations in CUS is: where is the rate constant of CVT, which locates the dividing surface between reactants and products at a point along the reaction path that minimizes the generalized transition-state theory rate coefficients for a given temperature T; is the Boltzmann constant, h is Planck’s constant, is the partition function per unit volume for reactants, and .

The CUS theory recrossing factor is where is calculated at the highest maximum of , at the second highest maximum, and at the minimum. Herein, SO coupling effects in the reactants are included by modifying the electronic partition function.[21,47]

The reaction path is calculated by the Euler steepest-descent method from −5.0 Bohr on the reactant side to +3.0 Bohr on the product side. The rotational partition functions of each point on the reaction path are obtained by ab initio frequency calculations classically, and the vibrational modes are treated as quantum mechanical separable harmonic oscillators, with the generalized normal-modes defined in redundant curvilinear coordinates. The geometries, gradient vector, and Hessian matrix of all points along IRC are obtained from CASSCF, while the energy is obtained from icMRCI +Q. The energy splitting between two electronic states for OH, between and are included in the CVT and CUS calculations. Since the title hydrogen-abstraction reaction shows that , the quantum-tunneling effect, which plays an important role at low temperature for a reaction of heavy–light–heavy mass combination,[69,70] is accounted for by the multidimensional small-curvature-tunneling (SCT) transmission coefficient . The rate constants reported are obtained from the CVT/SCT and CUS/SCT calculations.

3. Results and discussion
3.1. Stationary points

The ground state of the O atom (3P) has the electronic configuration . Ignoring for the moment the electrons of O atom, the remaining 2p4 electrons may occupy three orbitals 2px, 2py, and 2pz, which are represented in Fig. 1 to demonstrate the effective reaction directions and electronic states in the CH4 + O(3P) reaction system. For the Cs symmetry, this reaction correlates with different regions of the reaction: 3A′ + two 3A″ for reactant asymptote, for TS, for product asymptote. There are two surfaces that must be considered in the present work: the 3A′ state with the half-filled p orbital in the HCO plane and the 3A″ state with the half-filled p orbital perpendicular to the plane. The first triplet excited state (3A′) differs from the ground state (3A″) by the excitation of a single electron from the 2px orbital into the 2pz orbital, which weakens the influence of the 2pz orbital, causing the bond angle to slightly open as observed.

Fig. 1. (color online) Three orbitals , , and of O atom in CH4 + O(3P) reaction.

When the distance between C and O is smaller than 9.0 Bohr, the weakly-bound van der Waals (vdW) interaction of the ground-state O atom with methane along the H3C–H–O axis is helpful for orienting the system towards H3CH…O (WellR, and ) with a well depth of −0.46 kcal/mol relative to reactants. Czako et al. described this as a very shallow vdW minima with a depth of −0.18 kcal/mol, which could support few vibrational levels.[34] Gonzalez et al.[32,33] at the CCSD(T)=FULL/aug-cc-pVQZ//CCSD(T)=FC/cc-pVTZ level found the WellR depth of −0.5 kcal/mol at Bohr. All these calculations show that WellR does exist in the entrance channel and present a very weak long-range vdW interaction. The active orbitals in the reaction are the orbitals of the bond pair which transfers during the reaction (CH bond pair for reactants and OH bond pair for products, respectively) and the radical orbital (O atom for reactants and CH3 radical orbital for products, respectively).[25,32,33,47] The remaining CH bond pairs do not change qualitatively during the approach of highly symmetrical C3v geometry, which is an SP2nd ( and ) rather than an SP1st or minimum after WellR, except for a change in spatial orientation as the CH3 moiety goes from pyramidal for CH4 to planar for free CH3. Similarly, the rest orbitals of O atom remain localizing on the O during the reaction. At the C3v geometry, the molecule exists in a 3E electronic state and is therefore Jahn–Teller unstable. The O(2pπ) orbitals remain localizing on the O during the reaction. As a consequence of the electronic degeneracy, the 3E PES splits into two surfaces of symmetries 3A′ and 3A″, which are represented in Fig. 2. The PECs are very flat in the region of 175° to 185°, and the energy difference between the two electronic PESs is less than 0.01 kcal/mol at the CASSCF level. The singly occupied orbitals for the 3A′ state and 3A″ state change as the configuration of molecule deviates the linearity. Bending toward the lone pair will raise the energy more than bending toward the singly occupied orbital, due to the different electrostatic repulsions with the OH bond. Therefore, slight differences between singly occupied and lone pair orbitals make one state ( ) go higher, and another ( ) lower as the changes from 180°. One-dimensional stretching 3A′ and 3A″ PECs as a function of the OH and CH bonds, are represented in Fig. S2. Compared with bending PECs, stretching PECs change much faster but keep degeneracy occurring in orbitally degenerate.

Fig. 2. (color online) One-dimensional bending 3A′ and 3A″ PECs as a function of the O–H–C angle. For these calculations the other internal coordinates are fixed at the values of SP2nd (3E). The zero energy is taken at the Jahn–Teller conical intersection. The coordinates are defined in Fig. S1.
Fig. S2. (color online) One-dimensional stretching 3A′ and 3A″ potential energy curves as a function of the OH (a) and CH (b) bonds, respectively. For these calculations the other internal coordinates were fixed at SP1st (3A″) values. The zero energy is taken at the Jahn–Teller conical intersection.

Under Cs symmetry as the O(3P) atom approaches the methane along a C–H bond after WellR, two SP1st are found in the two respective surface of symmetries 3A″ and 3A′ with the length of breaking C–H bonds of 2.644 Bohr and 2.649 Bohr, respectively, and they are almost identical in energy but slightly different in the OHC bending angle and the tilt of the methyl group. The length of the breaking C–H bond in the two Cs symmetry SP1st is about 25.0% longer than the corresponding equilibrium distance in CH4 molecule, while the length of the forming O–H bond in SP1st geometry is about 20.0% longer than the corresponding equilibrium distance of OH. The geometry parameters and the relative energies of SP1st (3A′ and 3A″, respectively) and SP2nd (3E) are listed in the following Table 1.

Table 1.

Properties of saddle points (SPs) on the reaction path of in comparison with those from other ab initio calculations (bond length in unit Bohr, dihedral angle in unit degree, energy in units kcal/mol, and frequencies in unit cm−1).

.

The calculated barriers on both surfaces are around 14 kcal/mol, and the corresponding harmonic vibrational frequencies, performed at CASSCF levels, are 2284i and 2376i[47] in unit cm−1, respectively, for SP1st (3A′) and SP1st (3A″), while the single point calculations at the icMRCI+Q(14e,12o)/optCVDZ level give the barrier heights of 13.28 kcal/mol and 12.53 kcal/mol. The TST calculations in the next rate constant prediction lead to zero-point corrections of 1.46 kcal/mol for the 3A′ state and 1.56 kcal/mol for the 3A″ state. Therefore, including zero-point correction energy reduces the barrier to around 11.07 kcal/mol and 11.72 kcal/mol, respectively. Both in SP1st (3A′) and SP1st (3A″) are not zero, which indicates that it is a near-collinear H–CH3 direction as the O(3P) attacking methane with of 0.0° and 180.0°, for SP1st (3A′) and SP1st (3A″), respectively. This accounts for the steric effect in the reactivity since the barrier height of SP1st (3A′) is relatively higher and the corresponding transition state is relatively tighter than those of SP1st (3A″). Hence, we speculate that, everything else being equal, the shape of the excitation function in the post-threshold region in the crossed-beam experiment might be concave-up with higher collision energy. This conjecture awaits future experimental confirmation.

The second order saddle point SP2nd (3E) with linear configuration is found by the optimization of lower-energy points of degeneracy with lower symmetry. This structure presents two imaginary harmonic vibrational frequencies, which are presented in Fig. S3. Mode I: 2519i cm−1 corresponds to the transfer of a hydrogen atom from the minimum reactant side to that of the product side, and Mode II: 162i cm−1 to the OHC bending from one SP1st to another. The energy difference between SP1st (3A″) and SP2nd (3E) is 0.06 kcal/mol; however, it increases to 0.82 kcal/mol at single point calculations using icMRCI+Q(14e,12o)/optCVDZ. (It may be due to the fact that the total number of contracted configurations changes rapidly from linearity to non-linearity geometries in the icMRCI calculations.) An exhaustive geometry optimization for the SP2nd (3E) performed at the icMRCI level obtained the similar geometry with Bohr and Bohr, but this kind of calculation goes far beyond our currently reasonable computational resources, and neither frequency for SP2nd (3E) nor geometry optimizations for other stationary points could be achieved. The second order saddle point SP2nd (3E) is also confirmed by the topology of PESs in the vicinity of the conical intersection. The electronic degeneracy at SP2nd (3E) is lifted along two nuclear coordinates: the gradient difference vector (g) and the derivative coupling vector (h). These span a two-dimensional subspace of the nuclear coordinates called the branching space, the remaining dimensions forming the intersection space.[49,50,63] In the branching space, the adiabatic PESs form a double cone with a vertex at the origin, which are shown in Fig. 3 to illustrate the Jahn–Teller conical intersection. It may be important in chemical reaction since it serves as a funnel allowing radiationless electronic transitions from 3A′ to 3A″. The control of reactivity in the vicinity of conical intersections would allow one to control the possible products.

Fig. 3. (color online) Topography of two-dimensional PESs in the branching space of the Jahn–Teller conical intersection (3A′ and 3A″). The corresponding nuclear configurations are distributed in the branching space with derivative coupling :0.05:1.0 and gradient difference :0.05:1.0. The zero energy is taken at the Jahn–Teller conical intersection.
Fig. S3. (color online) The two imaginary frequency modes of the secondary order saddle point SP2nd (3E) between SP1st (3A′) and SP1st (3A″). Mode I: 2519 cm−1 for the transfer of a hydrogen atom; Mode II: 262 cm−1 for the OHC bending.

The minimum energy path for SP1st (3A′) is affirmed by an IRC calculation in Fig. S4, the single point energy of each point in the path is obtained at icMRCI+Q calculations. The pathways from the SP1st (3A″)[47] and SP1st (3A′) as the steepest descent paths WellR to WellP are defined and found in mass-weighted cartesian co-ordinates. The energy of SP2nd (3E) can easily be lower along the Mode II, leading to an SP1st or a minimum. It has a chemical significance, since there may be two minima/intermediates being connected dynamically by this SP2nd.

Fig. S4. (color online) O(3P)+CH4 ( ) reaction path through SP1st (3A′) obtained at the CASSCF(14e, 12o)/optCVDZ level.
3.2. Rate constant prediction

The TST, CVT, and CUS calculations were applied into this reaction with a single saddle point and two potential energy wells. Moreover, these calculations are based upon dual-level strategy extensive electronic structure determinations. The geometries, gradient vector, and Hessian matrix of all points along IRC are obtained from CASSCF, while the energy is obtained from icMRCI+Q//CASSCF(14e,12o)/optCVDZ. The SO coupling calculations for the O(3P) atom and CH4 molecular have been performed using SO pseudopotentials. The obtained energy splitting between two electronic states for OH is 118.2 cm−1, whereas 140.2 cm−1 and 259.9 cm−1 for and relative to the of O atom. These are included in the CVT and CUS calculations, as described before in Ref. [47]. SCT corrections were used to account for quantum-tunneling effects.

It should be noted that more symmetry unique geometries in the entrance and exit valley regions, calculated with the icMRCI method, were added in the present rate constant calculations of the 3A″ and 3A′ surfaces, respectively. The contribution to the rate from 3A″ and 3A′ surfaces have been calculated separately and summed to give the total rate. In Fig. 4, the thermal rate constants of the title reaction, obtained at different temperatures ranging from 298 K to 2500 K, are plotted as a typical Arrhenius behavior. The calculated reaction-rate constants (in units ) for the title reaction at 300, 600, 1000, and 2500 K, are listed in Table 2. At 300 K, the calculated rate constant CUS/SCT(3A′) with SO effects included is , which is lowered by 50% than without SO effects included. It shows that this kind of SO coupling effect in the entrance valley has an evident influence on the rate constant at low temperatures. However, with an increase of the temperature, this kind of influence becomes smaller. Moreover, all quantum mechanical effects become practically negligible at high temperature. The difference between CVT and conventional TST rate constants, the variational effect, is found at low temperature for this reaction, which is the expected behavior for a reaction with heavy–light–heavy mass combinations. Similarly, it becomes negligible at high temperature.

Fig. 4. (color online) Arrhenius plot for the reaction of . TST: presents the study using transition-state theory calculations; CVT/SCT: presents the study using canonical variational transition-state theory with small-curvature-tunneling (SCT) calculations; CUS/SCT: presents the study using canonical unified statistical theory with SCT calculations; RPMD(PES-2014): other study using ring polymer molecular dynamics approach calculations on PES-2014 from Ref. [33]; Exp.1: Experimental values from Ref. [31]; Exp.2: Experimental values from Ref. [30]. Arrhenius plot for the isotope effects on the rate constants. 3A′ presents the reaction path through SP1st (3A′) from [47], 3A″ presents the reaction path through SP1st (3A″).
Table 2.

Thermal rate coefficients (in unit ) for the O(3P) + CH4 reaction at 300, 600, 1000, and 2500 K using different dynamical methods.

.

The calculated rate constants of 600 K for the title reaction are and on the 3A″ and 3A′ PES, respectively, by using the CVT/SCT method. The total value of the rate constant is , which is in better agreement with the experimental result than previous theoretical values.[15] Although the values of the rate constant obtained from 3A″ PES at all temperatures are slightly larger than those from 3A′, the two surfaces make almost equal contributions to the reaction rate. The total rate constants of CVT and CUS are in good agreement with experiment over the temperature range from 298 K to 2500 K when the effect of tunneling is included using SCT transmission coefficient , while TST results are much higher and present noticeable difference in the Arrhenius plots, indicating the important role played by tunneling. Additionally, the values of the rate constant at 600 K and 1000 K accord better with the experimental results than those from Gonzalez’s calculation.[33] Therefore, the effect of SP1st (3A′) and SP2nd (3E) cannot be ignored.

The conventional TST in the present work locates the transition state dividing surface at the SP1st(3A′) and it seems to overestimate the rate constants, presumably due to the neglect of barrier recrossing. Since the harmonic approximation, which leads to an overestimation of the variational effects, is used in the present CUS/CVT calculations, their results could not perfectly match the experimental data at high temperature.[33,38] Our CUS/SCT results are larger at lower temperatures than theirs with CUS/LAT on PES-2014 and smaller at higher temperatures, but both results are in accordance with the available experimental data.

4. Conclusions

In the present work, we have presented two SPs1st and one SP2nd on an ab initio PES of reaction. These calculations are based upon dual-level strategy extensive electronic structure determinations at the CASSCF(14e,12o) level with a kind of optimized CVDZ basis set for the geometry optimizations, vibrational frequencies and reaction coordinate calculations, and at the icMRCI+Q(14e,12o) level for single point energies calculations. The core and core-valence correlation effects are included, which are necessary for an accurate quantum chemical description of the title reaction. Our calculations give a fairly accurate description of the regions around SP1st (3A′) and SP1st (3A″), which may play some role in the O(3P) attacking a collinear/near-collinear H–CH3 direction. For the first time, a second-order SP2nd (3E) between the above twin SPs with two imaginary vibrational frequencies at 162i cm−1 and 2519i cm−1 are found for this reaction. The corresponding 2519i cm−1 normal mode displays the beginnings of the transfer of a hydrogen atom from the minimum reactant side to that of the product side, and the other corresponding 162i cm−1 to the OHC bending from one SP1st to another. Our results indicate that the existence of SP2nd (3E) enriches the mechanism of reaction. Studying SP2nd (3E) can therefore help in understanding the nature of the transition states now recognised as being important for molecular reactions.

The thermal rate constants for this hydrogen-abstraction reaction are calculated with several dynamical methods including using the canonical unified statistical theory method over a large temperature range. These calculated rate constants of CVT and CUS are in agreement with experimental results; however, TST results are much higher and present noticeable difference in the Arrhenius plots, indicating the important role played by tunneling. The reasonably good agreement with available experimental results confirms the accuracy of our present work. The rate coefficient obtained from CUS/CVT is significantly lower than its previous RPMD result at low temperature below the crossover temperature, suggesting underestimation of the tunneling contribution. It may be attributed to the inaccuracies of CUS and CVT.

This theoretical study will play an important role in shaping our deeper understanding of SPs and bimolecular reactions with polyatomic molecules. Further work on the construction of the global PESs for reaction, suitable for the reactive scattering studies, is in progress, and various saddle points and PES intersections as revealed in our work, in particular conical intersections, will be included. A final judgment about the quality of our theoretical study on the mechanism and kinetics of the title reaction will only be possible once the further PES construction based on our present work and detailed dynamics calculations on the PES have been performed.

Reference
[1] Espinosa-Garcia J Rangel C Garcia-Bernaldez J C 2015 Phys. Chem. Chem. Phys. 17 6009
[2] Espinosa-Garcia J Garcia-Bernaldez J C 2000 Phys. Chem. Chem. Phys. 2 2345
[3] Zhang J Liu K 2011 Chem. Asian J. 6 3132
[4] Yang J Y Shao K J Zhang D Shuai Q A Fu B N Zhang D H Yang X M 2014 J. Phys. Chem. Lett. 5 3106
[5] Dianat A Seriani N Ciacchi L C Bobeth M Cuniberti G 2014 Chem. Phys. 443 53
[6] Zhang L Luo W L Ruan W Jiang G Zhu Z H 2008 Chin. Phys. 17 2023
[7] Hou X W Wan M F Ma Z Q 2012 Chin. Phys. 21 103301
[8] Galashev A Y 2013 Chin. Phys. 22 123602
[9] Troya D Schatz G C Garton D J Brunsvold A L Minton T K 2004 J. Chem. Phys. 120 731
[10] Sweeney G M Watson A McKendrick K G 1997 J. Chem. Phys. 106 9172
[11] Ausfelder F Kelso H McKendrick K G 2002 Phys. Chem. Chem. Phys. 4 473
[12] Troya D Garcia-Molina E 2005 J. Phys. Chem. 109 3015
[13] Zhang J M Lahankar S A Garton D J Minton T K Zhang W Q Yang X M 2011 J. Phys. Chem. 115 10894
[14] Garton D J Minton T K Troya D Pascual R Schatz G C 2003 J. Phys. Chem. 107 4583
[15] Li Y L Suleimanov Y V Yang M H Green W H Guo H 2013 J. Phys. Chem. Lett. 4 48
[16] Liu R Yang M H Gabor C Bowman J M Li J Guo H 2012 J. Phys. Chem. Lett. 3 3776
[17] Jiang B Guo H 2014 J. Chin. Chem. Soc. 61 847
[18] Thanh L N Stanton J F 2017 J. Chem. Phys. 147 152704
[19] Jing F Q Cao J W Liu X J Hu Y F Ma H T Bian W S 2016 Chin. J. Chem. Phys. 29 430
[20] Orlando R N Francisco M B C Truhlar D G 1999 J. Chem. Phys. 111 10046
[21] Ma H T Liu X J Bian W S Meng L P Zheng S J 2016 Chem Phys Chem 7 1786
[22] Minyaev R M Getmanskii I V Quapp W 2004 Russ. J Phys. Chem. A 78 1494
[23] Yang C L Wang M S Liu W W Zhang Z H Ma X G 2013 Chin. Phys. 22 63102
[24] Xie T Wang D Y Bowman J M Manolopoulos D E 2002 J. Chem. Phys. 116 7461
[25] Walch S P Dunning T H 1980 J. Chem. Phys. 72 3221
[26] Andresen P Luntz A C 1980 J. Chem. Phys. 72 5842
[27] Kleinermanns K Luntz A C 1981 J. Phys. Chem. 85 1966
[28] Suzuki T Hirota E 1993 J. Chem. Phys. 98 2387
[29] Westenberg A A Haas N D 1967 J. Chem. Phys. 46 490
[30] Baulch D L Cobos C J Cox R A Esser C Frank P Just T Kerr J A Pilling M J Troe J Walker R W Warnatz J 1992 J. Phys. Chem. Ref. Data. 21 411
[31] Cohen N 1986 Int. J. Chem. Kinet. 18 59
[32] Gonzalez-Lavado E Corchado J C Espinosa-Garcia J 2014 J. Chem. Phys. 140 064310
[33] Gonzalez-Lavado E Corchado J C Suleimanov Y V Green W H Espinosa-Garcia J 2014 J. Phys. Chem. 118 3243
[34] Czako G Bowman J M 2012 Proc. Natl. Acad Sci. USA 109 7997
[35] Gonzaclez M Hernando J Millacn J Sayocs R 1999 J. Chem. Phys. 110 7326
[36] Suleimanov Y V Allen J W Green W H 2013 Comput. Phys. Commun. 184 833
[37] Suleimanov Y V Aoiz F J Guo H 2016 J. Phys. Chem. 120 8488
[38] Suleimanov Y V Espinosa-Garcia J 2016 J. Phys. Chem. 120 1418
[39] Corchado J C Espinosa-Garcia J Roberto-Neto O Chuang Y Y Truhlar D G 1998 J. Phys. Chem. 102 4899
[40] Zhao H L Wang W J Zhao Y 2016 J. Phys. Chem. 120 7589
[41] Chen F Xu T F Yan D D Li W Z 1995 Chin. Phys. B [Acta Phys. Sin. (Overseas Edition) 4 125
[42] Chen F Li W Z 1996 Chin. Phys. B [Acta Phys. Sin. (Overseas Edition) 5 321
[43] Cederbaum L S Domcke W Koppel H 1978 Chem. Phys. 33 319
[44] Opalka D Segado M Poluyanov L V Domcke W 2010 Phys. Rev. 81 042501
[45] Domcke W Mishra S Poluyanov L V 2006 Chem. Phys. 322 405
[46] Gonzalez C McDouall J J W Schlegel H B 1990 J. Phys. Chem. 94 7467
[47] Peng Y Jiang Z A Chen J S 2017 J. Phys. Chem. 121 2209
[48] Shu Y N Parker K A Truhlar D G 2017 J. Phys. Chem. Lett. 8 2107
[49] Wu Y N Zhang C F Ma H T 2017 RSC Adv. 7 12074
[50] Yarkony D R 1998 Acc. Chem. Res. 31 511
[51] Werner H J Knowles P J Knizia G et al. 2008 MOLPRO, version 2008.1, a package of ab initio programs, see http://www.molpro.net
[52] Liu J Y Li Z S Dai Z W Zhang G Sun C C 2004 Chem. Phys. 296 43
[53] Werner H J Knowles P J 1985 J. Chem. Phys. 82 5053
[54] Werner H J Knowles P J 1985 Chem. Phys. Lett. 115 259
[55] Werner H J Knowles P J 1988 J. Chem. Phys. 89 5803
[56] Knowles P J Werner H J 1988 Chem. Phys. Lett. 145 514
[57] Schirmer J Cederbaum L S Domcke W Niessen W V 1977 Chem. Phys. 26 149
[58] Zheng J J Zhao Y Truhlar D G 2009 J. Chem. Theor. Comput. 5 808
[59] Matteoli E Mansoori G A 1995 J. Chem. Phys. 103 4672
[60] Peterson K A Dunning T H 2002 J. Chem. Phys. 117 10548
[61] Partridge H Schwenke D W 1997 J. Chem. Phys. 106 4618
[62] Fukui K 1970 J. Phys. Chem. 74 4161
[63] Yarkony D R 1996 Rev. Mod. Phys. 68 985
[64] Langhoff S R Davidson E R 1974 Int. J. Quantum Chem. 8 61
[65] Hu W P Truhlar D G 1995 J. Am. Chem. Soc. 117 10726
[66] Hu W P Truhlar D G 1996 J. Am. Chem. Soc. 118 860
[67] Garrett B C Truhlar D G 1979 J. Phys. Chem. 83 1079
[68] Zheng J Zhang S Lynch B J et al. POLYRATE, version 2015, Computer Program for the Calculation of Chemical Reaction Rates for Polyatomics, see https://comp.chem.umn.edu/polyrate/
[69] Bruce C G Donaldf G T 1983 J. Chem. Phys. 79 4931
[70] Li Y L Suleimanov Y V Green W H Guo H 2014 J. Phys. Chem. 118 1989