Mechanism and Kinetics of Methane Combustion, Part I: Thermal Rate

Subscriber access provided by Fudan University

Article

Mechanism and Kinetics of Methane Combustion, Part I: Thermal Rate Constants for the Hydrogen-abstraction Reaction of CH+ O(P) 4

3

Ya Peng, Zhong-an Jiang, and Jushi Chen J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b12125 • Publication Date (Web): 01 Mar 2017 Downloaded from http://pubs.acs.org on March 5, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Mechanism and Kinetics of Methane Combustion, Part I: Thermal Rate Constants for Hydrogen-abstraction Reaction of CH4 + O(3P) Ya Peng,∗ Zhong’an Jiang,∗ and Jushi Chen School of Civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083,China. E-mail: [email protected]; [email protected]

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract The mechanism and kinetics of gas-phase hydrogen-abstraction by the O(3 P) from methane are investigated using ab initio calculations and dynamical methods. Not only the electronic structure properties including the optimized geometries, relative energies and vibrational frequencies of all the stationary points are obtained from state-averaged complete active space self-consistent field calculations, but also the single-point energies for all points on the intrinsic reaction coordinate are evaluated using the internally contracted multireference configuration interaction approach with modified optimized cc-pCVDZ basis sets. Our calculations give a fairly accurate description of the regions around the 3A00 transition state in the O(3 P) attacking a near-collinear H–CH3 direction with a barrier height of 12.53 kcal/mol, which is lower than those reported before. Subsequently, thermal rate constants for this hydrogen-abstraction are calculated using the canonical unified statistical theory method with the temperature ranging from 298K to 1000K. These calculated rate constants are in agreement with experiments. The present work reveals the reaction mechanism of hydrogen-abstraction by the O(3 P) from methane, and it is helpful for the understanding of methane combustion.


Page 2 of 37

Page 3 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

Introduction

Methane is the main component of natural gas, and stored in not only various minerals 1–4 but also atmosphere as the most abundant hydrocarbon. 5,6 Consequently, it is one of the most important atmospheric species. The reactions of ground-state oxygen with hydrocarbons are key initial steps and of significant importance in the combustion processes. 7–10 However, even the reactions of O atom with the simplest hydrocarbon molecule are not well understood. It is commonly presumed that the first step is the hydrogen-abstraction by O(3 P) to generate a significant concentration of CH3 in flames, which could be oxidized to CO and CO2 in subsequent combustion steps. The slightly endothermic reaction O(3 P) + CH4 (X1 A1 ) → OH(X2 Π) + CH3 (X2 A002 ) is a prototype one for detailed investigations of the hydrocarbon combustion. 10–12 Thus, the topic of this hydrogen-abstraction by O(3 P) has attracted considerable interest, both from experimental and theoretical perspectives. 13–18 In the early 1970s, the reaction rate coefficients have been determined in the temperature range from 300K to 2500K. However, it is considerably less reliable due to the uncertainties in the reaction stoichiometry and tunneling effects, especially at low-temperature. 19,20 The experimental data for the reaction was reexamined by Cohen 21 at temperatures below about 600K. It was found that the values for rate coefficient were overestimated by factors of approximately 2 to 3 in the 250-400 K temperature range, with the error increasing as the temperature decreasing, and a new expression of the experimental results was proposed as k(T ) = 2.69 × 10−18 T 2.3 exp(−3570/T )cm3 molecule−1 s−1 . This expression represents the curvature in Arrhenius plots of the experimental data but leads to higher values of the rate constant than the previous one at low-temperature, in better agreement with the experimental values. Later, a wide variety of methods were employed to the experimental study of this hydrogen-abstraction reaction. The dynamics of the reactions O(3 P) + CH4 /CD4 were investigated experimentally, and the characterisation of the internal energy distributions of



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

products OH and CH3 were studied in bulk experiments with isotopic variants taken into account. As for the other product CH3 radical, the ν2 out-of-plane bending (ν2 = 1 ←2 up to 3 ←4) vibrational distribution was determined by Suzuki and Hirota 22 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. The product speed distributions at different Ecol were measured in crossed molecular beams by a number of groups, 17,18 and the excitation functions were also obtained. The reactions at centerof-mass collision energies in the range of 2.8–3.9 eV were investigated with crossed-beams experiments and with direct dynamics calculations. 18 Both the experiments and calculations provide evidence for previously unobserved reaction pathways which principally lead to O-atom addition and subsequent H-atom elimination. In addition, the expected H-atom abstraction reaction to form OH has been observed, with the modest classical barrier height of 11.76 kcal/mol. Ausfelder et al. 14 used isotope substitution to provide confirmation of the nuclear framework dynamics of this reaction and studied influences on the product rotation with fine-structure state partitioning. The O(3 P) atoms were generated by photolysis of NO2 . The rotational and fine-structure state distributions in the nascent product OH or OD radicals were determined by laser-induced fluorescence. 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. Theoretically, it has received a great attention using different approaches. One of the earlier ab initio studies of the reaction was performed with POL-CI wave functions by Walch and Dunning. 11 The reactant and product electronic asymptotic states are correlated by means of 3 A(×2) (C1 symmetry), 3 A0 +3A00 (Cs symmetry) or 3 E (C3v symmetry) adiabatic sheets, depending on the symmetry considered. The barrier height was obtained to be around 12.0 kcal/mol and the saddle geometry was found to be of C3v symmetry. The


Page 4 of 37

Page 5 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


geometries of stationary points including transition states in this reaction were optimized with the Hartree-Fock method by Gonzalez et al. 23 There were two Cs symmetry transition states in these calculations, and they were almost identical in energy but slightly different in the OHC bending angle (181◦ for 3 A0 and 178◦ for 3A00 ) and the tilt of the methyl group. The calculated barriers on both surfaces are around 15.34 kcal/mol, and the corresponding harmonic vibrational frequencies, performed at UMP2 levels, are 2197i and 2198i in cm−1 , respectively. The ground potential energy surface (PES), constructed by Gonz´ aleza and co-authors 24 using the PMP4 calculations, were employed to obtain the barrier height of 13.3 kcal/mol and the zero-point-inclusive barrier height of 9.4 kcal/mol. The OHC bending angle is 180.0◦ and the imaginary frequency is 2470i cm−1 . Roberto-Neto and co-authors 25 identified the barrier height of 14.2 kcal/mol with high level multireference Møller-Plesset second order perturbation theory calculations, in which a new database was used to assess electronic structure methods. Recently, with the development of theoretical chemistry, high-level ab initio quantum and dynamics calculations have been used in addressing this hydrogen-abstraction reaction. 26–31 The first full dimensional analytical PES was constructed by Corchado and co-workers. 26 In 2012, the high accuracy was achieved by using a large basis set and 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 saddle point, which was taken as a benchmark value. 27 In a recent study, a new full-dimensional analytical PES-2014 for O(3 P) + CH4 was accurately constructed based on CCSD(T)=FULL/aug-cc-pVQZ level by Gonz´ alez-Lavado and co-authors. 32 They showed that the correlation of all electrons (FULL), including core and valence electrons, has a major influence on the description of the barrier. The barrier height obtained from this surface 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) 33–36 and with variational transition state theory calculations on this new surface. 37 The transmission coefficients transition state theory calculations were obtained



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

by the least-action tunneling approximation. In the classical high-temperature limit regime, the RPMD rate theory result agrees with experiment perfectly, although it overestimates the experimental rate coefficients in the low-temperature regime. 33–36 Moreover, calculating accurate reaction rates requires the determination of accurate geometries, energies, and vibrational frequencies for the reactants, products, and transition structure and sometimes also for one or more intermediate complexes. An accurate description of the transition structure is very important for the thermal rate calculations. However, it is somewhat uncertain in the previous theoretical studies, since the approach of O(3 P) along a C–H bond has 3–fold symmetry and leads to a Jahn-Teller conical intersection, 38–40 which corresponds to a 3E state, and then breaking the C3v symmetry splits the degenerate E state into one 3A0 and one 3A00 surface. Therefore, previous calculations on the saddle

3

point at various ab initio levels make this reaction more and more elusive, representing a theoretical challenge. 5,23–27 One of the reasons for the overestimation of the barrier height proposed by Gonzalez et al. 23 was the inadequacy of a single-configuration reference wave function for the post-Hartree-Fock calculations. As the first part of a two-part series, this paper focuses on the title reaction with the 3 00

A transition state. The multireference configuration interaction methods are employed

to an accurate ab initio description of the above two electronic states (3 A0 and 3A00 ), and further kinetic calculations are performed to study the hydrogen-abstraction reaction of CH4 + O(3 P). So far most of the ab initio work concerning this reaction has been based on the so-called frozen core approximation, in which correlation effects involving the electrons in 1s core orbitals of carbon and oxygen are neglected. However, if the goals of a calculation are to obtain chemical accuracy of thermochemical properties, the effects of correlating the electrons in the core orbital generally must be addressed in the calculations. 27,32,41 Not only appropriate active space but also an optimized basis set with additional functions for describing core and core-valence correlation (called CV) effects is employed in our ab initio calculations. The reaction channel is determined. The reaction rate constants are calculated


Page 6 of 37

Page 7 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and compared with the available experimental data and other previous theoretical results. In addition, the kinetic isotope effects are reported. This article is organized as follows. The theoretical methods and computational details describing ab initio calculations, reactionpath analysis and reaction rate constant are presented in Sec.2. In Sec.3, we present the results and discussion. Finally, the conclusions are given in Sec.4.

2

Theoretical Methods

2.1

Ab initio Electronic Structure Calculations

In this work, we present an extensive theoretical study on hydrogen-abstraction reaction of CH4 + O(3 P) based upon the accurate ab initio calculations, which are performed using the MOLPRO 2008.1 program package. 42 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 point), 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. 1, in which the bending angle θ is 6 OCH and φOCHH 0 is the dihedral angle. It is important to generate balanced orbitals for use in subsequent correlation calculations, so state-averaged complete active space SCF (CASSCF) 43,44 calculations including one root in 3 A0 and one in 3A00 symmetries,respectively, are carried out for the orbital optimization using the orbitals obtained with the Unrestricted Hartree-Fock method as the starting guess. Then, utilizing the CASSCF energies as reference values, the energies of each electronic state are computed by the internally contracted multireference configuration interaction method (icMRCI). 45,46 A large amount of correlation energy is descirbed via single and double electron excitations, and further correlation energy due to higher excitations is approximated by the Davidson correction (+Q). 47 There are at least two reaction channels/transition states (3 A0 and 3 A00 ) involved in the title reaction. 11,23 Therefore, two reference 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

states (3 A0 and 3A00 ) are required and used for generating the internally contracted pairs in the icMRCI procedure. The selection of the active space is crucial in the CASSCF and icMRCI+Q calculations. Here various active spaces in CASSCF calculations have been tested (shown in Table S1 and S2) and finally the active space, which consists of twelve orbitals, referred to as (14e,12o), is chosen for the further 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 (CH4 :sp3 hybridization → CH3 :sp2 hybridization) and OH bond-forming accurately. Although, the use of the active space (12e,11o) for CASSCF and icMRCI calculations is limited to relatively small basis sets, many calculations have been performed to check carefully the basis set dependence to ensure reliable theoretical results. Several basis sets, including cc-pCVDZ, aug-cc-pCVTZ, cc-pCVQZ and cc-pCVTZ have been tested to obtain the reasonable reaction endoergicity. (Some results with previous reported calculations are represented in Figure S1.) A comparison with the frozen-core and all-electron results is also provided. The result of all-electron is better than that of frozen-core. The 1s orbitals of carbon and oxygen are not correlated but optimized at SA-CASSCF level. To obtain the CV effects contribution to the energies, especially barriers, 32,48 the cc-pCVDZ basis sets are chosen but modified. The details of this scheme could be found elsewhere, 49–51 and only a brief outline will be given here. The first seven inner s functions of C/O atom are contracted to two functions using the coefficients from the cc-pVDZ basis sets. 52 The outer two s functions are uncontracted as the fourth p functions, while the first seven p functions are contracted to one function. These core basis sets, developed to treat both core and valence correlations, are of the form (9s,4p,1d)/[4s,2p,1d] and designated as optCVDZ, which are smaller than Dunning’s standard cc-pCVDZ (9s,4p,1d)/[4s,3p,1d] and give a good description of core electrons. The geometrical structures of all reactants, transition state, complexes and products are optimized at CASSCF/optCVDZ level. The distance between two fragments of reactants/products is fixed at 30.0 bohr. At the same level, the corresponding frequency


Page 8 of 37

Page 9 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculations of the optimized geometries are also done to prove the characters of minima without imaginary frequencies and the transition state with only one imaginary frequency. 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 with the desired reactant and product. 53 The numerically gradient vector and hessian matrix including rotational partition functions were evaluated at each point along the IRC by harmonic vibrational frequencies calculations. To evaluate relative energy of each point in the IRC, the single point energy is also calculated at subsequent icMRCI+Q calculation on the basis of the CASSCF optimized geometry, which is sufficient to accurately describe the title reaction. 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+Q calculations to account for the core-valence correlations. The rest of the inner electrons are kept frozen and not correlated.

2.2

Kinetic Calculations

Recently, some programs for computing rate coefficients have been developed. For instance, RPMD approach has several desirable features and had been employed in the title reaction. 31,34,37 It gives the correct quantum statistics, which is important for tunneling at lowtemperature, and approaches the classical limit at high-temperature. 33–36 However, it is not suitable for the present work due to the fact that the number of points required to do this is much larger than the number of points available. It will go far beyond our currently reasonable computational resources if we try to get so many points using the present ab initio method. The kinetic calculations in the present work are carried out using the canonical variational transition-state (CVT) theory with general polyatomic rate constants code POLYRATE 2015 program. 54 As for a long-lived potential intermediate species formed along the reaction pathway, there may be more than one kinetic bottleneck beside the one in the central barrier region 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

including a saddle point. The canonical unified statistical (CUS) theory, 55,56 which considers a complex region between two bottlenecks and does not give a bound even in classical mechanics, is a good choice for this kind of kinetic calculation 26,57 in the present work. The general Equation for rate constant calculations in CUS is:

kTCU S = kTCV T RTCU S kTCV T is the rate constant of CVT, which locates the dividing surface between reactants and products at a point s∗,CV T (T ) along the reaction path that minimizes the generalized transition-state theory (TST) rate coefficients for a given temperature T, and it could be given by

CV T kTCV T = [kB T /hΦR T ]qT

where, kB is Boltzmann’s constant, h is Planck’s constant, and ΦR T is the partition function per unit volume (length in the collinear world) for reactants. It could be further expressed as kTCV T = σ

kB T ◦ K exp[−∆G(T, s∗,CV T )/kB T ] h

with σ being the symmetry factor, K ◦ being the reciprocal of the standard-state concentration, and s being the reaction coordinate, measured as the mass-weigthed distance along the reaction path to the saddle point. The CUS theory recrossing factor is

RTCU S

k CV T k CV T = 1 + max + min k k

−1

where kTCV T is calculated at the highest maximum of ∆G, k max at the second highest maximum, and k min at the minimum. The reaction path is calculated by the Euler steepest-descent method from -2.0 bohr on the reactant side to +1.8 bohr on the product side. The rotational partition functions of each


Page 10 of 37

Page 11 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


point on the reaction path are obtained classically by ab initio frequency calculations, and the vibrational modes are treated as quantum mechanical separable harmonic oscillators, with the generalized normal-modes defined in redundant curvilinear coordinates. The title hydrogen-abstraction reaction shows that mH mO , mCH3 , indicating a heavylight-heavy mass combination. Therefore, tunneling at low-temperature plays an important role. 58,59 The following levels of tunneling calculations have been considered: small-curvature tunneling (SCT), large-curvature tunneling (LCT), microcanonical optimized multidimensional tunneling (µOMT) and least-action tunneling (LAT). The LAT method is considered to be the most accurate. 35,37,59 However, LCT/µOMT/LAT could not be achieved in the present work due to the following two reasons. Firstly, instead of the analytical PES, the ab initio saddle point and dozens of points along the IRC are considered to treat explicitly all atoms of the system. Secondly, the potential energy curve 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. So the quantum-tunneling effects are accounted for by the SCT transmission coefficient kTSCT . The rate constants reported are obtained from the CVT/SCT and CUS/SCT calculations.

3 3.1

Results and Discussion Stationary Points and Mechanism

To demonstrate the effective sites and reaction directions, the O(3 P) atom undergoes two series of scans in which the OCH bending angle θ is rotated through 0.0◦ to 180.0◦ and the dihedral angle φOCHH 0 is rotated through 120.0◦ to 180.0◦ . In one series of scan calculations, the fragment CH4 is fixed at the methane equilibrium geometry and RCO is fixed at 6.0 bohr, while in another series of scan calculations the products fragments (CH3 and ROH ) are fixed at the corresponding equilibrium geometries and RCH is fixed at 5.0 bohr. The step sizes are 3.0◦ and 1.0◦ for θ and φ scans, respectively. Since the −CH3 is frozen at C3v symmetry in 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

all scan calculations, the dihedral angle φOCHH 0 changing from 120.0◦ to 180.0◦ is good for a circle, which goes through 3 H atoms of CH4 . The PESs for the O(3 P) + CH4 reaction, which are plotted as a function of θ and φOCHH 0 with RCO fixed at 6.0, are shown in Fig. 2(a) to illustrate the entrance channel. The interaction of the O atom with methane leads to a van der Waals (vdW) potential well H3 CH···O (WellR ,3A00 ) with four geometry configurations of O-H-C(H0 H00 H000 ) (θ=0.0◦ ), O-H0 -C(HH00 H000 ) (θ=109.2◦ ,φOCHH 0 =120.0◦ ), O-H00 -C(HH0 H000 ) (θ=109.2◦ ,φOCHH 0 =240.0◦ ) and O-H000 -C(HH0 H00 ) (θ=109.2◦ ,φOCHH 0 =0.0◦ ). The potential energy curves of the entrance channel for the O(3 P) + CH4 reaction are presented in Fig. 3, in which the red, green and blue lines are corresponding to the PECs of three RCO values of 6.0, 7.5 and 9.0 bohrs, respectively. It can be seen that the system can enter the well region easily when O(3 P) approaches CH4 at all RCO and φOCHH 0 values. When RCO is smaller than 9.0 bohr, the weakly-bound vdW interaction of the ground-state O atom with methane along the O···HCH3 axis is helpful for orienting the system towards WellR with a depth of -0.46 kcal/mol relative to reactants. The geometry parameters of WellR including θ=0.0◦ , ROH =5.976 bohr, RCH =2.100 bohr, are listed in Table 1, in which the other stationarypoint properties from our ab initio calculations and previous literatures are also collected. Czako et al. 27 described this as a very shallow vdW minima with a depth of -0.18 kcal/mol, which could support few vibrational levels. Gonzalez et al. 23 at CCSD(T)=FULL/aug-ccpVQZ//CCSD(T)=FC/cc-pVTZ level found the WellR depth of -0.5 kcal/mol at ROH =5.630 bohr. All these calculations show that WellR does exist in the entrance channel and present a very weak long-range vdW interaction. Under Cs symmetry as the O(3 P) atom approaches the methane along a C-H bond after WellR , a collinear O··H··CH3 transition state (TS) is found with the length of the breaking C-H bond of 2.644 bohr, which is 25.80% (0.542 bohr) longer than the corresponding equilibrium distance in CH4 molecule. The length of the forming O-H bond in the TS geometry is 2.212 bohr, which is 20.68% (0.379 bohr) longer than the corresponding equilibrium distance of OH. Our calculated 6 OCH 0 , 6 OCH 00 and 6 OCH 000 (103.6◦ , 103.6◦ and 104.1◦ ) are close to


Page 12 of 37

Page 13 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


109.2◦ in CH4 , so it is an early barrier with the reactant-like structure. The TS geometric parameters and barrier heights from our ab initio calculations and previous literatures are shown in Fig. 4. As pointed out by Gonzalez et al., 23 for CH4 + O(3 P) reaction, there are two surfaces that must be considered separately: the 3 A0 state with the half-filled p orbital in the HCO plane and the 3A00 state with the half-filled p orbital perpendicular to the plane. In the previous work, the TS geometry was somewhat uncertain and the corresponding calculations at various ab initio levels are not consistent. 5 Czako et al. 27 described the transition state using all-electron CCSDT(Q)/aug-cc-pCVQZ calculations, with the bond lengths being formed ROH and broken RCH at 2.253 bohr and 2.470 bohr, respectively, and the barrier height of 14.08 kcal/mol. They were taken as the the best estimated values to data and served as benchmark values. The PES-2014 for the title reaction obtained at CCSD(T)=FULL/aug-cc-pVQZ level presents the barrier with a height of 14.0 kcal/mol and a lower imaginary frequency for TS of 1900i cm−1 , which is much lower than our calculated value of 2375i cm−1 . Obviously the thinner barrier obtained in our present work facilitates more facile tunnelling. The barrier height of 14.08 kcal/mol from “benchmark” and 14.0 kcal/mol from “PES-2014” are much higher than the experimental value of 8.7– 11.7 kcal/mol 19,60 extracted from Arrhenius plots, though Arrhenius expression does not provide a particularly good fit to the experimental data. Our calculated barrier height is 14.13 kcal/mol at CASSCF(14e,12o)/optCVDZ level, in agreement with the previous reported theoretical values, while the single point calculation at icMRCI+Q(14e,12o)/optCVDZ level gives the barrier height of 12.53 kcal/mol with ROH =2.212 bohr and RCH =2.644 bohr, and including harmonic zero-point vibrational energy reduces the barrier to around 11.07 kcal/mol. An exhaustive geometry optimization for the TS performed at icMRCI level obtained the similar geometry with ROH =2.234 bohr and RCH =2.436 bohr, but this kind of calculations goes far beyond our currently reasonable computational resources, and neither frequency for TS nor geometry optimizations for other stationary points could be achieved. There are several possible O atom attacking directions to CH4 molecule and another



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

two possible vdW wells Well0R and Well00R (reactant-complex-face configuration and reactantcomplex-edge configuration, respectively, called in Reference [ 32,61]) are founded in the entrance channel. The optimized geometries and harmonic vibrational frequencies for these two vdW potential wells are listed in Table 2. Yan et al. 61 described the Well0R and Well00R configurations, with well depths of -0.27 and -0.18 kcal/mol, respectively. Gonz´ alez-Lavado and co-authors 32,37 searched for their configurations using several ab initio methods. They characterized Well0R with a depth well of -0.3 kcal/mol by MP2/6-31G(d,p) calculations. In addition, the well depth reduces to -0.4 kcal/mol on their PES. As for Well00R , they found stable structures at MP2/6-31G(d,p) and CCSD(T)/cc-pVTZ levels, but these structures always present an imaginary frequency. So they concluded that the intermediate complexes in the entrance channel present weak long-range vdW interactions, and that their stability is strongly dependent on the theoretical level. High-level ab initio calculations give an accurate description of the long-range vdW region. In the exit channel, the potential energy surfaces, plotted as a function of θ and φOCHH 0 with RCH fixed at 5.0 bohr and two product fragments at the equilibrium geometries, are shown in Fig. 2(b). It exhibits a shallow vdW well H3 C···OH (WellP ,3A00 ) with ROH =1.858 bohr and RCH 0 /RCH 00 /RCH 000 =2.085 bohr, which is -1.91 kcal/mol lower than the product asymptote and may influence the dynamics of the reaction by randomizing the energy. The geometry parameters of WellP are listed in Table 1. These vdW complexes, in particular WellR , Well0R and Well00R in the entrance channel, may play some role in the reaction dynamics, 62–65 and could be responsible for a number of prominent resonances in the cumulative reaction probability in the tunneling region. For further investigation of the reaction mechanism, it is desirable to find a connecting pathway among these stationary points (reactants, WellR , TS, WellP and products). This pathway is defined as the steepest descent path from the TS to the WellR and WellP , respectively, and is found in mass-weighted cartesian coordinates. The minimum energy path is affirmed by IRC calculation in Fig. 5, and the single point energy of each point in the path is obtained at icMRCI+Q calculation.


Page 14 of 37

Page 15 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2

Thermal 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 idea of the dual-level strategy 26,57,66 is to use two levels of ab initio calculations in the thermal rate constant prediction. The geometries, gradient vector, hessian matrix of all points along IRC are obtained from CASSCF, while the energy is obtained from icMRCI+Q. The spin-orbit (SO) coupling calculations for the O(3 P) atom and CH4 (X1 A1 ) molecular have been performed using spin-orbit pseudopotentials. The obtained energy splitting between two electronic states for OH is 118.2 cm−1 , whereas 140.2 and 259.9 cm−1 for 3 P1 and 3 P0 relative to the 3 P2 of O atom. These are included in the CVT and CUS calculations. According to the geometries, harmonic frequencies, and energy of reactants, products and TSs, the thermal rate constants of hydrogen-abstraction reaction have been calculated at different temperatures ranging from 298K to 2500K, which are plotted in Fig. 6 (a) as a typical Arrhenius behaviour. Normally, all quantum mechanical effects become practically negligible at high-temperature. The present theoretical thermal rate constants calculations are based upon dual-level strategy extensive electronic structure determinations. The discrepancy between the potential energy curve of ab initio single point energies obtained from higher-level calculations and the counterpart from lower-level calculations might not be practically negligible in dynamics calculations at high-temperature. The conventional TST in the present work locates the transition state dividing surface at the saddle points 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 match perfectly the experimental data at high-temperature. 35–37 The calculated rate constants for the title reaction are 9.77×10−18 , 5.93×10−15 , 1.08×10−13 ,



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 37

cm3 molecule−1 s−1 at 300, 600 and 1000 K,respectively, by using CVT/SCT method. The value of the rate constant at 300 K is in better agreement with the experimental result 9.00×10−18 cm3 molecule−1 s−1 than that from Gonzalez’s calculation (1.32×10−17 cm3 molecule−1 s−1 ). 37 Our CUS/SCT results are larger at lower temperatures than theirs with CUS/LAT on PES2014 and smaller at higher temperatures, but both results are in accordance with the available experimental data. The rate coefficient obtained from CUS/CVT is significantly lower than its previous RPMD result at low-temperature below the crossover temperature, 37 suggesting underestimation of the tunneling contribution. 59,67 It may be attributed to the inaccuracies of CUS and CVT. The rate coefficients of the deuterated reaction are also calculated at temperatures ranging from 298K to 1000K, which are plotted in Fig. 6 (b). The ratio of the rate coefficients of the unsubstituted reaction to the deuterated reaction, CH4 /CD4 , is obtained to analyze the kinetics isotope effects (KIEs). (The KIEs at different temperatures are listed in Table S4.) Our calculated KIEs from CVT/SCT and the previous works 29,31,37 show similar trends, but their values do not agree well. In general, the present KIEs values from CVT/SCT are lower than the earlier ones. To the best of our knowledge, there has been no experimental reports on the isotope data for the title reaction available for comparison and we hope that our present work will benefit experiments in the future. Further work on the various accurate dynamical calculations for the hydrogen-abstraction reaction of CH4 + O(3 P) studies, in particular KIEs, will be included.

4

Conclusions

In this work, accurate ab initio and dynamical calculations are performed to investigate the path energy and reaction barrier of the hydrogen-abstraction reaction between O(3 P) atom and methane molecule. These calculations are based upon dual-level strategy extensive electronic structure determinations, including the geometries optimizations, vibrational


Page 17 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


frequencies and IRC calculations at the CASSCF(14e,12o) level, and single point energies calculations at icMRCI+Q(14e,12o) level with a kind of optimized CVDZ basis sets. The core and core-valence correlation effects are included, which are necessary for an accurate quantum chemical description of the transition structure and long-range interactions in the title reaction. Our calculations give a fairly accurate description of the regions around the 3A00 transition state in the O(3 P) attacking a near-collinear H–CH3 direction with a barrier height of 12.53 kcal/mol, which is lower than those reported before. Moreover, several possible vdW wells are obtained in the entrance and exit channels. WellR with a depth of -0.46 kcal/mol relative to reactants presents a very weak long-range vdW interaction, which may play some role in the reaction dynamics. The thermal rate constants for this hydrogen-abstraction reaction are calculated using several dynamical methods with the temperature ranging from 298K to 1000K. 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 value of the rate constant at 300 K is in better agreement with the experimental result 9.00×10−18 cm3 molecule−1 s−1 than that from Gonzalez’s calculation (1.32×10−17 cm3 molecule−1 s−1 ). The kinetics isotope effects of the title reaction are reported. Our calculated KIEs from CVT/SCT show similar trends with the previous works, but they are lower than the earlier ones. The present work provides us a basis for further nonadiabatic dynamical studies, in particular Jahn-Teller conical intersection, and further work for the construction of the global PESs is in progress.

5

Supporting Information

Complete references (42) and (54); Properties of stationary points using CASSCF(10e,10o)/optCVQZ and CASSCF(10e,10o)/optCVTZ on reaction path of O(3 P) + CH4 (X1 A1 ) → OH(X2 Π) +



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

CH3 (X2 A002 ) in comparison with those from other ab initio calculations and from experiments are listed in Table S1 - S2; Thermal rate coefficients for the O(3 P) + CH4 reaction at 300,600 1000 and 2500 K using different dynamical methods are listed in Table S3; H/D kinetic isotope effects for the O(3 P) + CH4 reaction at 300, 600, 1000 and 2500 K using different dynamical methods are listed in Table S4. The comparison of energy difference between EP roducts and EReactants with various ab initio methods are represented in Figure S1.

Acknowledgement This work is supported by the National Natural Science Foundation of China (No. 51604018, 51574016) and the Natural Science Foundation of Beijing (No. 8164060). The authors thank Professor Donald G. Truhlar for providing the POLYRATE 2015 program. Many computations are carried out at Virtual Laboratory of Computational Chemistry, Computer Network Information Center of Chinese Academy of Sciences.


Page 18 of 37

Page 19 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References (1) Bodelier, P. L. E. Sustainability bypassing the methane cycle. Nature 2015, 523, 534– 535. (2) Lawton, T. J.; Rosenzweig, A. C. Methane–make it or break it. Science 2016, 352, 892–893. (3) Jiang, Z.; Jiang, J.; Wang, H.; Du, C.; Luo, J. Experimental study on gas migration law in gob area of fully mechanized caving coal face. Min. Saf. Environ. Prot. 2015, 3, 5–11. (4) Ding, H.; Jiang, Z.; Han, Y. Numerical simulation and application of boreholes along coal seam for methane drainage. J. Univ. Sci. Technol. B. 2008, 30, 1205–1210. (5) Zhang, J.; Liu, K. Imaging the reaction dynamics of O(3 P) + CH4 → OH + CH3 . Chem. Asian J. 2011, 6, 3132–3136. (6) Pan, H.; Liu, K. Communication: imaging the effects of the antisymmetric-stretching excitation in the O(3 P) + CH4 (v3 = 1) reaction. J. Chem. Phys. 2014, 140, 191101. (7) Dianat, A.; Seriani, N.; Ciacchi, L. C.; Bobeth, M.; Cuniberti, G. DFT study of reaction processes of methane combustion on PdO(1 0 0). Chem. Phys. 2014, 443, 53–60. (8) Zhou, C.-W.; Simmie, J. M.; Curran, H. J. Rate constants for hydrogen-abstraction by from n-butanol. Combust. Flam. 2011, 158, 726–731. (9) Jing, F.; Cao, J.; Liu, X.; Hu, Y.; Ma, H.; Bian, W. Theoretical study on mechanism and kinetics of reaction of O(3 P) with propane. Chin. J. Chem. Phys. 2016, 29, 430–436. (10) Martinez, R.; Enriquez, P. A.; Puyuelo, M. P.; Gonzalez, M. Exploring the stereodynamics and microscopic mechanism of the O(3 P) + CH4 , CD4 → OH + CH3 , OD + CD3 combustion reactions. Chem. Phys. 2015, 461, 98–105. 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(11) Walch, S. P.; Dunning, T. H. Calculated barrier to hydrogen atom abstraction from CH4 by O(3 P). J. Chem. Phys. 1980, 72, 3221–3227. (12) Wang, M.-L.; Li, Y.-M.; Zhang, J. Z. H. Application of semirigid vibrating rotor target model to the reaction of O(3 P) + CH4 → CH3 + OH. J. Phys. Chem. A 2001, 105, 2530–2534. (13) Sweeney, G. M.; Watson, A.; McKendrick, K. G. Rotational and spin-orbit effects in the dynamics of O(3 Pj ) + hydrocarbon reactions. I. experimental results. J. Chem. Phys. 1997, 106, 9172–9181. (14) Ausfelder, F.; Kelso, H.; McKendrick, K. G. The dynamics of O(3 P) + deuterated hydrocarbons: influences on product rotation and fine-structure state partitioning. Phys. Chem. Chem. Phys. 2002, 4, 473–481. (15) Troya, D.; Schatz, G. C.; Garton, D. J.; Brunsvold, A. L.; Minton, T. K. Crossed beams and theoretical studies of the O(3 P) + CH4 → H + OCH3 reaction excitation function. J. Chem. Phys. 2004, 120, 731–739. (16) Troya, D.; Garcia-Molina, E. Quasiclassical trajectory study of the O(3 P) + CH4 → OH + CH3 reaction with a specific reaction parameters semiempirical hamiltonian. J. Phys. Chem. A 2005, 109, 3015–3023. (17) Zhang, J.; Lahankar, S. A.; Garton, D. J.; Minton, T. K.; Zhang, W.; Yang, X. Crossedbeams studies of the dynamics of the H-Atom abstraction reaction, O(3 P) + CH4 → OH + CH3 , at hyperthermal collision energies. J. Phys. Chem. A 2011, 115, 10894–10902. (18) Garton, D. J.; Minton, T. K.; Troya, D.; Pascual, R.; Schatz, G. C. Hyperthermal reactions of O(3 P) with alkanes: observations of novel reaction pathways in crossedbeams and theoretical studies. J. Phys. Chem. A 2003, 107, 4583–4587.


Page 20 of 37

Page 21 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(19) Westenberg, A. A.; de Haas, N. Atom-molecule kinetics at high temperature using ESR detection. Technique and Results for O + H2 , O + CH4 , and O + C2 H6 . J. Chem. Phys. 1967, 46, 490–501. (20) 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. et al. Evaluated kinetic data for combustion modelling. J. Phys. Chem. Ref. Data 1992, 21, 411–734. (21) Cohen, N. A reevaluation of low temperature experimental rate data for the reactions of O atoms with methane, ethane, and neopentane. Int. J. Chem. Kinet. 1986, 18, 59–82. (22) Suzuki, T.; Hirota, E. Vibrational distribution of CH3 produced by the reaction of O(1 D) atom with CH4 . J. Chem. Phys. 1993, 98, 2387–2398. (23) Gonzalez, C.; McDouall, J. J. W.; Schlegel, H. B. Ab initio study of the reactions between methane and hydroxyl, hydrogen atom, and triplet oxygen atom. J. Phys. Chem. 1990, 94, 7467–7471. (24) Gonz´ alez, M.; Hernando, J.; Milln, J.; Says, R. Ab initio ground potential energy surface, VTST and QCT study of the O(3 P) + CH4 (X1 A1 ) → OH(X2 Π) + CH3 reaction. J. Chem. Phys. 1999, 110, 7326–7338. (25) Roberto-Neto, O.; Machado, F. B. C.; Truhlar, D. G. Energetic and structural features of the CH4 + O(3 P) → CH3 + OH abstraction reaction: Does perturbation theory from a multiconfiguration reference state (finally) provide a balanced treatment of transition states? J. Chem. Phys. 1999, 111, 10046–10052. (26) Corchado, J. C.; Espinosa-Garcia, J.; Roberto-Neto, O.; Chuang, Y.-Y.; Truhlar, D. G. Dual-level direct dynamics calculations of the reaction rates for a Jahn-Teller reaction: hydrogen abstraction from CH4 or CD4 by O(3 P). J. Phys. Chem. A 1998, 102, 4899– 4910. 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(27) Czako, G.; Bowman, J. M. Dynamics of the O(3 P) + CHD3 (vCH =0,1) reactions on an accurate ab initio potential energy surface. P. Natl. Acad. Sci. USA 2012, 109, 7997–8001. (28) Czako, G. Communication: Direct comparison between theory and experiment for correlated angular and product-state distributions of the ground-state and stretchingexcited O(3 P) + CH4 reactions. J. Chem. Phys. 2014, 140, 231102. (29) Espinosa-Garcia, J.; Garcia-Bernaldez, J. C. Analytical potential energy surface for the CH4 + O(3 P)→CH3 + OH reaction. Thermal rate constants and kinetic isotope effects. Phys. Chem. Chem. Phys. 2000, 2, 2345–2351. (30) Huarte-Larra˜ naga, F.; Manthe, U. Accurate quantum dynamics of a combustion reaction: Thermal rate constants of O(3 P) + CH4 (X1 A1 ) → OH + CH3 . J. Chem. Phys. 2002, 117, 4635–4638. (31) Zhao, H.; Wang, W.; Zhao, Y. Thermal rate constants for the O(3 P) + CH4 → OH + CH3 reaction: the effects of quantum tunneling and potential energy barrier shape. J. Phys. Chem. A 2016, 120, 7589–7597. (32) Gonz´ alez-Lavado, E.; Corchado, J. C.; Espinosa-Garcia, J. The hydrogen abstraction reaction O(3 P) + CH4 : A new analytical potential energy surface based on fit to ab initio calculations. J. Chem. Phys. 2014, 140, 064310. (33) Suleimanov, Y.; Allen, J.; Green, W. RPMDrate: Bimolecular chemical reaction rates from ring polymer molecular dynamics. Comput. Phys. Commun. 2013, 184, 833–840. (34) Li, Y.; Suleimanov, Y. V.; Yang, M.; Green, W. H.; Guo, H. Ring polymer molecular dynamics calculations of thermal rate constants for the O(3 P) + CH4 → OH + CH3 reaction: contributions of quantum effects. J. Phys. Chem. Lett. 2013, 4, 48–52.


Page 22 of 37

Page 23 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(35) Suleimanov, Y. V.; Aoiz, F. J.; Guo, H. Chemical reaction rate coefficients from ring polymer molecular dynamics: theory and practical applications. J. Phys. Chem. A 2016, 120, 8488–8502. (36) Suleimanov, Y. V.; Espinosa-Garcia, J. Recrossing and Tunneling in the Kinetics Study of the OH + CH4 → H2 O + CH3 Reaction. J. Phys. Chem. B 2016, 120, 1418–1428. (37) Gonz´ alez-Lavado, E.; Corchado, J. C.; Suleimanov, Y. V.; Green, W. H.; EspinosaGarcia, J. Theoretical kinetics study of the O(3 P) + CH4 /CD4 hydrogen abstraction reaction: The role of anharmonicity, recrossing effects, and quantum mechanical tunneling. J. Phys. Chem. A 2014, 118, 3243–3252. (38) Cederbaum, L.; Domcke, W.; K¨oppel, H. Jahn-Teller effect induced by non-degenerate vibrational modes in cumulenes. Chem. Phys. 1978, 33, 319–326. (39) Opalka, D.; Segado, M.; Poluyanov, L. V.; Domcke, W. Relativistic Jahn-Teller effect in tetrahedral systems. Phys. Rev. A 2010, 81, 042501. (40) Domcke, W.; Mishra, S.; Poluyanov, L. V. The relativistic E×E Jahn-Teller effect revisited. Chem. Phys. 2006, 322, 405–410. (41) Schirmer, J.; Cederbaum, L.; Domcke, W.; von Niessen, W. Strong correlation effects in inner valence ionization of N2 and CO. Chem. Phys. 1977, 26, 149–153. (42) Werner H.-J., Knowles P. J., Knizia G., Manby F. R., Sch¨ utzM., Celani P., Gy¨orffy W., Kats D., Korona T. and Lindh R. et al., MOLPRO, version 2008.1, a package of ab initio programs, see http://www.molpro.net. (43) Werner, H.-J.; Knowles, P. J. A second order multiconfiguration SCF procedure with optimum convergence. J. Chem. Phys. 1985, 82, 5053–5063. (44) Werner, H.-J.; Knowles, P. J. An efficient second-order MC SCF method for long configuration expansions. Chem. Phys. Lett. 1985, 115, 259–267. 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(45) Werner, H.-J.; Knowles, P. J. An efficient internally contracted multiconfiguration reference configuration-interaction method. J. Chem. Phys. 1988, 89, 5803–5814. (46) Knowles, P. J.; Werner, H.-J. An efficient method for the evaluation of coupling coefficients in configuration interaction calculations. Chem. Phys. Lett. 1988, 145, 514–522. (47) Langhoff, S. R.; Davidson, E. R. Configuration interaction calculations on the nitrogen molecule. Int. J. Quantum Chem 1974, 8, 61–72. (48) Zheng, J.; Zhao, Y.; Truhlar, D. G. The DBH24/08 database and its use to assess electronic structure model chemistries for chemical reaction barrier heights. J. Chem. Theory Comput. 2009, 5, 808–821. (49) Matteoli, E.; Mansoori, G. A. A simple expression for radial distribution functions of pure fluids and mixtures. J. Chem. Phys. 1995, 103, 4672–4677. (50) Peterson, K. A.; Dunning, T. H. Accurate correlation consistent basis sets for molecular corevalence correlation effects: The second row atoms AlAr, and the first row atoms BNe revisited. J. Chem. Phys. 2002, 117, 10548–10560. (51) Partridge, H.; Schwenke, D. W. The determination of an accurate isotope dependent potential energy surface for water from extensive ab initio calculations and experimental data. J. Chem. Phys. 1997, 106, 4618–4639. (52) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (53) Fukui, K. A formulation of reaction coordinate. J. Phys. Chem. 1970, 74, 4161. (54) Zheng J., Zhang S., Lynch B. J., Corchado J. C., Chuang Y. Y. , Fast P. L., Hu W. P., Liu Y. P., Lynch G. C., Nguyen K. A., et al, POLYRATE, version 2015, Computer Program for the Calculation of Chemical Reaction Rates for Polyatomics, see https://comp.chem.umn.edu/polyrate/. 24 ACS Paragon Plus Environment

Page 24 of 37

Page 25 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(55) Hu, W.-P.; Truhlar, D. G. Deuterium kinetic isotope effects and their temperature dependence in the gas-pPhase S2 N reactions X− + CH3 Y → CH3 X + Y− (X, Y = Cl, Br, I). J. Am. Chem. Soc. 1995, 117, 10726–10734. (56) Hu, W.-P.; Truhlar, D. G. Factors affecting competitive ion-molecule reactions: ClO− + C2 H5 Cl and C2 D5 Cl via E2 and S2 N channels. J. Am. Chem. Soc. 1996, 118, 860–869. (57) Ma, H.; Liu, X.; Bian, W.; Meng, L.; Zheng, S. A theoretical study of the mechanism and kinetics of F+N3 reactions. ChemPhysChem 2006, 7, 1786–1794. (58) Garrett, B. C.; Truhlar, D. G. A least-action variational method for calculating multidimensional tunneling probabilities for chemical reactions. J. Chem. Phys. 1983, 79, 4931–4938. (59) Li, Y.; Suleimanov, Y. V.; Green, W. H.; Guo, H. Quantum rate coefficients and kinetic isotope effect for the reaction Cl + CH4 → HCl + CH3 from ring polymer molecular dynamics. J. Phys. Chem. A 2014, 118, 1989–1996. (60) Brabbs, T. A.; Brokaw, R. S. 15th Symp. (Int.) Combust. 1975, 893. (61) Yan, T.; Doubleday, C.; Hase, W. L. A PM3-SRP + analytic function potential energy surface model for O(3 P) reactions with alkanes. Application to O(3 P) + Ethane. J. Phys. Chem. A 2004, 108, 9863–9875. (62) Xie, T.; Wang, D.; Bowman, J. M.; Manolopoulos, D. E. Resonances in the O(3 P) + HCl reaction due to van der Waals minima. J. Chem. Phys. 2002, 116, 7461–7467. (63) Liu, X.; Sobolewski, A. L.; Borrelli, R.; Domcke, W. Computational investigation of the photoinduced homolytic dissociation of water in the pyridine-water complex. Phys. Chem. Chem. Phys. 2013, 15, 5957–5966. (64) Zhang, C.; Fu, M.; Shen, Z.; Ma, H.; Bian, W. Global analytical ab initio ground-state



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

potential energy surface for the C(1 D) + H2 reactive system. J. Chem. Phys. 2014, 140, 234301. (65) Xie, J.; Zhang, J.; Hase, W. L. Is there hydrogen bonding for gas phase SN 2 pre-reaction complexes? Int. J. Mass Spectrom. 2015, 378, 14–19. (66) Liu, J.; Li, Z.; Dai, Z.; Zhang, G.; Sun, C. Dual-level direct dynamics studies for the hydrogen abstraction reaction of 1,1-difluoroethane with O(3 P). Chem. Phys. 2004, 296, 43–51. (67) Zuo, J.; Li, Y.; Guo, H.; Xie, D. Rate Coefficients of the HCl + OH → Cl + H2 O Reaction from Ring Polymer Molecular Dynamics. J. Phys. Chem. A 2016, 120, 3433– 3440.


Page 26 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

WellR b our otherc -0.51 -0.5 (-0.46) -0.18e -0.06f 5.976 5.630 2.100 2.058 2.102 2.056 2.102 2.056 2.102 2.056 0.03 0.0 109.6 109.4 109.4 120.1 119.8 120.1 -

TS(3A00 ) ourb otherc 14.13 14.0 (12.53) 14.08e WellP our otherc 3.22 3.5 (3.85) 3.12e b

our 5.22 (5.76)

b

Products otherc 5.7 5.32e

ROH 30.0 2.212 2.253 1.858 1.843 1.833 1.835 RCH 2.102 2.058 2.062 2.644 2.470 4.833 4.318 30.0 RCH 0 2.102 2.058 2.062 2.094 2.050 2.085 2.041 2.023 2.039 00 Geom RCH 2.102 2.058 2.062 2.094 2.050 2.085 2.041 2.023 2.039 RCH 000 2.102 2.058 2.062 2.094 2.050 2.085 2.041 2.023 2.039 6 OCH 0.00 0.61 0.0 0.00 0.0 0.00 6 OCH 0 109.2 109.5 109.5 103.6 91.0 91.2 00 6 OCH 109.5 109.5 109.5 103.6 89.9 90.0 6 OCH 000 109.5 109.5 109.5 104.1 89.9 90.0 φH 0 CHH 00 120.0 120.0 120.0 120.2 120.0 120.0 φH 00 CHH 000 120.0 120.0 120.0 122.2 119.9 120.0 φH 000 CHH 0 120.0 120.0 120.0 120.2 120.0 120.0 υi 2376i 1900i υ1 1339 1344 1306 60 42 222 404 44 85 599 491 υ2 1339 83 72 373 475 132 123 1393 1430 υ3 1339 83 99 521 582 163 156 1393 υ4 1556 1570 1534 1340 1341 1044 1044 295 293 3019 3126 Freq υ5 1556 1344 1345 1067 1163 324 323 3203 3289 υ6 2939 3035 2917 1344 1345 1154 1210 582 570 3204 υ7 3068 3153 3019 1559 1572 1349 1496 1397 1391 3740 3743 υ8 3068 1559 1572 1425 1518 1411 1418 υ9 3068 2942 3034 2982 3157 3010 3114 υ10 3068 3152 3113 3214 3223 3285 υ11 3068 3152 3135 3288 3487 3305 υ12 3080 3159 3701 3650 a 3 1 The energy of the O( P) + CH4 (X A1 ) asymptote is taken to be zero. b Geometry optimization and vibrational frequency from CASSCF(14e,12o)/optCVDZ level. The values in parentheses correspond to the single point energy calculation from icMRCI+Q(14e,12o)/optCVDZ //CASSCF(14e,12o)/optCVDZ level. c CCSD(T)=FULL/aug-cc-pVQZ//CCSD(T)=FC/cc-pVTZ single point level from Reference [ 32]. d NIST Standard Reference Database Number 69, http://webbook.nist.gov/chemistry/ e Benchmark: all-electron CCSDT(Q)/complete-basis-set quality from Reference [ 27]. f PMP2//UMP2/aug-cc-pVTZ level from Reference [ 61].

Relative Energya

our 0.0

b

Reactants otherc exp.d 0.0 0.0

Table 1: Properties of stationary points on reaction path of O(3 P) + CH4 (X1 A1 ) → OH(X2 Π) + CH3 (X2 A002 ) in comparison with those from other ab initio calculations and from experiments (bond length in bohr, dihedral angle in degrees, energy in kcal/mol and frequencies in cm−1 ).


ACS Paragon Plus Environment

27

3738

3004 3171

606 1402

1.8324 2.039 2.039 2.039 -

exp.d

Page 27 of 37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2: Properties of reactant and product complexes in comparison with those from other ab initio calculations (bond length in bohr, dihedral angle in degrees, energy in kcal/mol and frequencies in cm−1 ). Well0R Well00R b c d our other other ourb otherd Relative Energya -0.56 -0.3 -0.27 -0.40 -0.18 e (-0.50) -0.48 (-0.31) ROH 9.074 6.316 RCO 7.443 RCH 2.102 2.050 RCH 2.102 RCH 0 2.101 RCH 0 2.102 Geom RCH 00 2.101 RCH 00 2.101 RCH 000 2.101 RCH 000 2.101 6 OCH 6 180.0 OCH 120.6 6 OCH 0 6 OCH 0 109.6 120.6 00 6 6 OCH 00 109.4 OCH 54.7 6 OCH 000 6 OCH 000 109.4 54.7 φH 0 CHH 00 120.1 φH 0 COH 0 180.0 00 000 00 φH CHH 120.1 φH COH -89.4 φH 000 CHH 0 120.1 φH 000 COH 179.7 υ1 44 65 υ1 1336 υ2 78 98 υ2 1337 υ3 111 100 υ3 1342 υ4 1395 1401 υ4 1555 υ5 1508 1410 υ5 1558 Freq υ6 1508 1410 υ6 2942 υ7 1766 1777 υ7 3069 υ8 1784 1783 υ8 3072 υ9 3130 3134 υ9 3073 υ10 3279 3278 υ11 3379 3327 υ12 3379 3327 a 3 1 The energy of the O( P) + CH4 (X A1 ) asymptote is taken to be zero. b Geometry optimization and vibrational frequency from CASSCF(14e,12o)/optCVDZ level. The values in parentheses correspond to the single point energy calculation from icMRCI+Q(14e,12o)/optCVDZ //CASSCF(14e,12o)/optCVDZ level. c CCSD(T)=FULL/aug-cc-pVQZ//CCSD(T)=FC/cc-pVTZ single point level from Reference [ 32]. d PMP2//UMP2/aug-cc-pVTZ level from Reference [ 61]. e Benchmark: accurate relative energies obtained at the all-electron CCSDT(Q)/complete-basis-set quality from Reference [ 27].


Page 28 of 37

Page 29 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure Legends

Fig. 1 The coordinate definition geometry of CH4 O.

Fig. 2

The potential energy surface in the entrance channels as a function of θ and

φOCHH 0 with (a): the fragment CH4 is fixed at the methane equilibrium geometry and RCO is fixed at 6.0 bohr; with (b): the fragment CH3 and OH are fixed at the equilibrium geometries and RCH is fixed at 5.0 bohr. The energy of the O(3 P) + CH4 (X1 A1 ) and OH(X2 Π) + CH3 (X2 A002 ) asymptotes are taken to be zero in (a) and (b), respectively. The bending angle θ is 6 OCH and φOCHH 0 is dihedral angle. The coordinates are defined in Fig. 1.

Fig. 3

The bending potential energy curves for O + CH4 reaction system calculated

at the CASSCF(14e,12o)/optCVDZ level. The bending angle θ is 6 OCH and φOCHH 0 is dihedral angle. The red, green and blue lines are corresponding to the PECs of three RCO values of 6.0, 7.5 and 9.0 bohr, respectively. The coordinates are defined in Fig. 1. The energy of the O(3 P) + CH4 (X1 A1 ) asymptote is taken to be zero.

Fig. 4

Structural parameters (distances in bohr, angles in degrees and energies in k-

cal/mol) for the transition state obtained from various ab initio levels. Corchado:geometry at the UMP2/cc-pVTZ level and energy at the PUMP2-SAC//UMP2/cc-pVTZ level from Reference [ 26]; Troya:PMP2/aug-cc-pVDZ from Reference [ 16]; Czako:all-electron CCSD(T)/augcc-pCVQZ from Reference [ 27]; Gonzalez:geometry at the valence-electron CCSD(T)/ccpVTZ level and energy at the all-electron CCSD(T)/cc-pVQZ level from Reference [ 32];

our:geometry at the CASSCF(14e,12o)/ level and energy at the icMRCI+Q(14e,12o)/optCVDZ //CASSCF(



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 5 O(3 P) + CH4 (X1 A1 ) reaction path obtained at the CASSCF(14e,12o)/optCVDZ level. The energy of the O(3 P) + CH4 (X1 A1 ) asymptote is taken to be zero.

Fig. 6

(a) Arrhenius plot for the reaction of CH4 +O(3 P)→ CH3 +OH. (b) Arrhenius

plot for the isotope effects on the rate constants. TST: present study using transition-state theory calculations; CVT/SCT: present study using canonical variational transition-state theory with small-curvature-tunneling (SCT) calculations; CUS/SCT: present study using canonical unified statistical theory with SCT calculations; CUS/LAT(PES-2014): other study using CUS theory with least-action tunneling calculations on PES-2014 from Reference [ 37]; RPMD(PES-2014): other study using ring polymer molecular dynamics approach calculations on PES-2014 from Reference [ 37]; Exp.1: Experimental values from Reference [ 21]; Exp.2: Experimental values from Reference [ 20]. Arrhenius plot for the isotope effects on the rate constants. TSTD : present study on O(3 P) + CD4 reaction using TST calculations; CVTD /SCT: present study O(3 P) + CD4 reaction using CVS theory with SCT calculations; CUSD /SCT: present study O(3 P) + CD4 reaction using CUS theory with SCT calculations.


Page 30 of 37

Page 31 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


O

H θ C H0

H000 φ H00

Figure 1: Peng et al.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

(b)

32



Page 32 of 37

Page 33 of 37

−0.414

−0.416

Relative Energies (kcal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


φ=180o o

φ=180

−0.418

φ=180o

RCO=9.0 bohr

φ=120o −0.42

o

φ=120

RCO=7.5 bohr

φ=150o

−0.422

−0.424 RCO=6.0 bohr −0.426

0

20

φ=120o 40

60

80 100 θOCH (degree)

Figure 3: Peng et al. 33 ACS Paragon Plus Environment

120

140

160

180


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Parameter Corchado:1998 Troya:2005 Czako:2012 Gonzalez:2014 This work Barrier 14.0 14.37 14.08 14.0 12.53

RCH 00 2.033 2.084 2.047 H000 2.050 2.094 H00 RCH 0 2.033 2.084 2.047 2.050 2.094

Page 34 of 37

RCH 000 2.033 2.084 2.047 2.050 2.094 C

H0



RCH 2.362 2.338 2.434 2.470 2.644

∠CHO 180.8 180.0 179.3 180.0 179.0 H

O ROH 2.270 2.328 2.292 2.253 2.212

Page 35 of 37

16

14

TS(3A")

12

Energy/kcal mol −1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10

8

CH3+OH

6

4

WellP 2

0

−2

CH4+O WellR //

−1.8

−1.5

−1.2

−0.9

−0.6

−0.3

0

0.3

0.6

0.9

Intrinsic reaction coordinate S/(aum)1/2 bohr

Figure 5: Peng et al. 35 ACS Paragon Plus Environment

1.2

1.5

//


−7 TST CVT/SCT CUS/SCT CUS/LAT(PES−2014) RPMD(PES−2014) Exp.1 Exp.2

−8 −9

log10(k/cm3molecule−1s−1)

−10 −11 −12 −13 −14 −15 −16 −17 −18

0.5

1

1.5

2 1000/T(K)

2.5

3

(a)

−7 TST CVT/SCT CUS/SCT TSTD

−8 −9

CVTD/SCT −10 log10(k/cm3molecule−1s−1)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 37

CUSD/SCT

−11 −12 −13 −14 −15 −16 −17 −18

0.5

1

1.5

2

2.5

1000/T(K)

(b)

36



3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry 16

14

TS(3A")

12

Energy/kcal mol −1

Page 37 of 37

10

8

CH3+OH

6

4

WellP 2

0

−2

CH4+O WellR //

−1.8

−1.5

−1.2

−0.9

−0.6

−0.3

0

0.3

0.6

0.9

1.2

1.5

//

Intrinsic reaction coordinate S/(aum)1/2 bohr


Mechanism and Kinetics of Methane Combustion, Part I: Thermal Rate

Recommend Documents