The hydrogen abstraction reaction O(3P) + CH4: A new analytical potential energy surface based on fit to ab initio calculations Eloisa González-Lavado, Jose C. Corchado, and Joaquin Espinosa-Garcia Citation: The Journal of Chemical Physics 140, 064310 (2014); doi: 10.1063/1.4864358 View online: http://dx.doi.org/10.1063/1.4864358 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/6?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Accurate ab initio potential energy surface, thermochemistry, and dynamics of the Br(2P, 2P3/2) + CH4 HBr + CH3 reaction J. Chem. Phys. 138, 134301 (2013); 10.1063/1.4797467 A new ab initio intermolecular potential energy surface and predicted rotational spectra of the KrH2O complex J. Chem. Phys. 137, 224314 (2012); 10.1063/1.4770263 Anharmonic force field and vibrational dynamics of CH2F2 up to 5000 cm1 studied by Fourier transform infrared spectroscopy and state-of-the-art ab initio calculations J. Chem. Phys. 136, 214302 (2012); 10.1063/1.4720502 Accurate ab initio potential energy surface, thermochemistry, and dynamics of the Cl(2P, 2P3/2) + CH4 HCl + CH3 and H + CH3Cl reactions J. Chem. Phys. 136, 044307 (2012); 10.1063/1.3679014 The hydrogen abstraction reaction H + CH 4 . I. New analytical potential energy surface based on fitting to ab initio calculations J. Chem. Phys. 130, 184314 (2009); 10.1063/1.3132223
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 194.57.123.86 On: Mon, 07 Apr 2014 07:34:55
THE JOURNAL OF CHEMICAL PHYSICS 140, 064310 (2014)
The hydrogen abstraction reaction O(3 P) + CH4 : A new analytical potential energy surface based on fit to ab initio calculations Eloisa González-Lavado, Jose C. Corchado, and Joaquin Espinosa-Garciaa) Departamento de Química Física, Universidad de Extremadura, 06071 Badajoz, Spain
(Received 22 December 2013; accepted 25 January 2014; published online 12 February 2014) Based exclusively on high-level ab initio calculations, a new full-dimensional analytical potential energy surface (PES-2014) for the gas-phase reaction of hydrogen abstraction from methane by an oxygen atom is developed. The ab initio information employed in the fit includes properties (equilibrium geometries, relative energies, and vibrational frequencies) of the reactants, products, saddle point, points on the reaction path, and points on the reaction swath, taking especial caution respecting the location and characterization of the intermediate complexes in the entrance and exit channels. By comparing with the reference results we show that the resulting PES-2014 reproduces reasonably well the whole set of ab initio data used in the fitting, obtained at the CCSD(T) = FULL/aug-ccpVQZ//CCSD(T) = FC/cc-pVTZ single point level, which represents a severe test of the new surface. As a first application, on this analytical surface we perform an extensive dynamics study using quasi-classical trajectory calculations, comparing the results with recent experimental and theoretical data. The excitation function increases with energy (concave-up) reproducing experimental and theoretical information, although our values are somewhat larger. The OH rotovibrational distribution is cold in agreement with experiment. Finally, our results reproduce experimental backward scattering distribution, associated to a rebound mechanism. These results lend confidence to the accuracy of the new surface, which substantially improves the results obtained with our previous surface (PES-2000) for the same system. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4864358] I. INTRODUCTION
The construction of potential energy surfaces (PES) represents a process in continuous evolution, closely related to the development of high-level theoretical methods, functional forms more adequate to represent the nuclear motion, and the appearance of new experimental data. Hence, the PESs play a central role in the analysis of the kinetics and dynamics of reactive systems; and in the case of polyatomic systems this is no trivial task, and it is time-consuming. In 2000 our group1 reported an analytical surface for the title reaction (PES-2000), which was symmetric with respect to any permutation of the four methane hydrogens, and yielded rate constants in good agreement with the experimental data. Since this surface provides not only the energy of the reactive system but also the analytical gradients (i.e., the analytical first energy derivatives), in recent years it has been used as a testing bench for different dynamics methods. For instance, variational transition state theory with multidimensional tunneling effect (VTST/MT),1 multiconfiguration time-dependent Hartree (MCTDH),2 quasi-classical trajectory (QCT),3 reduced-dimensional quantum dynamics (QD),4 or more recently ring polymer molecular dynamics (RPMD).5 However, since we constructed this surface, we have observed some drawbacks in the symmetric PES-2000. First, the barrier height (13.0 kcal mol−1 ) is too low with respect to high-level ab initio calculations6 (14.1 kcal mol−1 ) and the top of the barrier is too broad (imaginary frequency of a) [email protected]
0021-9606/2014/140(6)/064310/12/$30.00
1549 cm−1 ) with respect to high-level ab initio calculations1 (imaginary frequency around 2000 cm−1 ). Second, experimentally this reaction yields a rotationally cold OH product,7 while PES-2000 gives a rotationally hotter OH product.8 Third, in the calibration process we used the LCT3 method9 to describe the tunneling effect of large curvature, but later, Truhlar et al. found that this method overestimates this effect, proposing new methods to correct this deficiency, LCT4,10 and LAT.11, 12 Finally, the presence of intermediate complexes in the entrance and exit channels were not considered. In addition, the most serious drawback is that this PES-2000 surface is semiempirical in nature, i.e., it was fitted to reproduce theoretical and experimental data. The present study of this reaction has been divided into two parts. In the first part, we report the construction of a new analytical surface for the title reaction that corrects the major deficiencies of PES-2000, which is based exclusively on highlevel ab initio calculations. It is named PES-2014 and is also symmetric with respect to the permutations of the H atoms. In the second part, as a test of quality, we report a dynamics study based on this surface. The O(3 P) + CH4 → OH + CH3 reaction is an interesting system from an experimental point of view, since it is of substantial importance in the chemistry of hydrocarbon combustion (high temperatures), and also from a theoretical point of view, because it presents a heavy-light-heavy mass combination which is a good candidate for large tunnelling effects (low temperatures). The thermal rate constants have been measured experimentally with a wide variety of methods.13–32 The most recent recommended expression for the rate
140, 064310-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 194.57.123.86 On: Mon, 07 Apr 2014 07:34:55
064310-2
González-Lavado, Corchado, and Espinosa-Garcia
constant for the temperature range 300–2500 K is:32 k(T) = 1.15 × 10−15 T1.56 exp(−4270/T) cm3 molecule−1 s−1 , although below 400 K, the value of k is considerably less reliable due to uncertainties in the reaction stoichiometry.30 The dynamics of this reaction (nascent vibrational distributions of the CH3 product, nascent rotational distributions of the OH product, and differential cross section) has also been widely studied with different techniques.7, 33–38 Due to the great experimental interest of this reaction, theoretically, it has received a great attention using different approaches: semiempirical direct dynamics,3, 8 reduced-dimensional quantum dynamics,39–43 and QCT calculations.3, 44–46 In parallel with the develop of dynamics methods, great efforts have been made in the construction of its potential energy surface,1, 6, 45, 47 highlighting the recent full-dimensional ab initio PES reported by Czakó and Bowman (CB) in 20126 based on a permutationally invariant fit of 17 212 accurate energy points obtained by an efficient composite method. In addition, this reaction represents a theoretical challenge because the approach of O(3 P) along a CH bond has 3-fold symmetry and leads to a Jahn-Teller conical intersection rather than a saddle point. The conical intersection corresponds to a 3 E state, and breaking the C3v symmetry splits this into two surfaces, 3 A and 3 A . Walch and Dunning48 and Schlegel et al.49 using ab initio molecular orbital calculations, obtained optimized geometries and frequencies at the stationary points, and the barrier height and heat of reaction. The saddle point geometry is found to be of Cs symmetry, but close to C3v and the predicted barrier heights for the two surfaces (3 A and 3 A ) present a very small difference (0.2 kcal mol−1 ). Therefore, we can assume that the PES2014 surface can describe either of the two surfaces. This approach has been employed in previous studies.1 In summary, in the present paper to correct the deficiencies of PES-2000, we report the construction of a new analytical surface for the title reaction based exclusively on highlevel ab initio calculations, named PES-2014. The paper is structured as follows. In Sec. II, high-level electronic structure calculations are outlined, with especial caution taken with respect to the barrier height, the complexes in the entrance and exit channels, and the topology of the reaction path from reactants to products. In Sec. III, a detailed description of the PES for the O(3 P) + CH4 hydrogen abstraction reaction and the fitting approach is presented. The results of the application of the fitting method are presented in Sec. IV and tested against theoretical electronic structure data. The dynamics results using quasi-classical trajectory calculations are presented in Secs. V and VI and compared with the available experimental evidence. Finally the conclusions are presented in Sec. VII. II. ELECTRONIC STRUCTURE CALCULATIONS
In this work we develop the PES-2014 to describe the hydrogen abstraction reaction in the title. For the fitting procedure the input information is exclusively based on high-level ab initio calculations describing the topology of the reaction from reactants to products, with special care taken in the description of the intermediate complexes.
J. Chem. Phys. 140, 064310 (2014)
We took the coupled-cluster singles and doubles approach, including a quasi-perturbative estimate of the connected triple excitations, correlating only valence electrons, CCSD(T), using the correlation consistent polarized triplezeta basis set, cc-pVTZ,50 abbreviated CCSD(T) = FC/ccpVTZ level, where FC means frozen core. We began optimizing and characterizing all the stationary point geometries (reactants, products, complexes in the entrance and exit channels, and saddle point). Starting from the saddle point geometry and going downhill to both the asymptotic reactant and product channels in mass-weighted Cartesian coordinates, we constructed the minimum energy path (MEP) at this level. In addition, we computed the gradients and Hessians for 60 points along the MEP. Since at this level the vibrational frequency calculations are numerical, to compute a Cartesian coordinate Hessian approximately 324 energies (324 is 18 squared, and 18 is 3 coordinates per atom times 6 atoms) are necessary. As 60 points were considered along the reaction path, 19 940 (324 times 60) energy points were calculated. By adding optimization and vibrational frequency calculations at the stationary points, we calculated about 20 000 energies at the CCSD(T)/cc-pVTZ level of theory, which were then used in the fitting procedure. Finally, to improve the energy description of the reactive system, we performed single-point calculations at a higher level, CCSD(T) = FULL/aug-cc-pVQZ (see below). The calculations were performed using the GAUSSRATE code,51 which serves as an interface between the GAUSSIAN09 systems of programs52 and POLYRATE-2010.53 Barrier height is undoubtedly one of the most difficult energy properties of a reaction to estimate accurately, and as Zhang and Liu38 noted, “the saddle-point geometry is somewhat uncertain. The predictions of the relative order of the lengths of the breaking C–H bond and the forming O– H bond obtained at various ab initio levels are not consistent.” We characterized the saddle point properties by using ab initio calculations at the CCSD(T)/cc-pVTZ level, and they are listed in Table I while Fig. 1 shows the optimized geometry. It is well known that the classic barrier height is strongly dependent on the basis set. Recently, Truhlar et al.54 used a new database to assess electronic structure methods for barrier heights. They showed that the costly CCSD(T) = FC/cc-pVTZ level presents a mean unsigned error (MUE) of 2.86 kcal mol−1 , and that the correlation of all electrons (FULL) and the use of a larger basis set have a major influence on the description of the barrier. Thus, using the single point calculation at the CCSD(T) = FULL/aug-cc-pVQZ level, we obtain a barrier height of 14.0 kcal mol−1 , i.e., 2.7 kcal mol−1 lower than that obtained with the CCSD(T) = FC/cc-pVTZ method. This lack of convergence when changing the level of accuracy, both method and basis, may be due to the existence of a conical intersection giving rise to the saddle point. This effect was also found in other related systems as F + OH → O(3 P) + HF, where the three states correlating to the O(3 P) intersect several times along the reaction path.55 These results illustrate the dramatic influence of electronic correlation and basis set on the correct description of the barrier. These will be the values used to fit the analytical PES (Sec. III). The single-point calculations with the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 194.57.123.86 On: Mon, 07 Apr 2014 07:34:55
064310-3
González-Lavado, Corchado, and Espinosa-Garcia
J. Chem. Phys. 140, 064310 (2014)
TABLE I. Properties of the reactants, products, and saddle point.a Geometry Parameter
Frequency
PES
Ab initio
PES
R(C–H)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 194.57.123.86 On: Mon, 07 Apr 2014 07:34:55
THE JOURNAL OF CHEMICAL PHYSICS 140, 064310 (2014)
The hydrogen abstraction reaction O(3 P) + CH4 : A new analytical potential energy surface based on fit to ab initio calculations Eloisa González-Lavado, Jose C. Corchado, and Joaquin Espinosa-Garciaa) Departamento de Química Física, Universidad de Extremadura, 06071 Badajoz, Spain
(Received 22 December 2013; accepted 25 January 2014; published online 12 February 2014) Based exclusively on high-level ab initio calculations, a new full-dimensional analytical potential energy surface (PES-2014) for the gas-phase reaction of hydrogen abstraction from methane by an oxygen atom is developed. The ab initio information employed in the fit includes properties (equilibrium geometries, relative energies, and vibrational frequencies) of the reactants, products, saddle point, points on the reaction path, and points on the reaction swath, taking especial caution respecting the location and characterization of the intermediate complexes in the entrance and exit channels. By comparing with the reference results we show that the resulting PES-2014 reproduces reasonably well the whole set of ab initio data used in the fitting, obtained at the CCSD(T) = FULL/aug-ccpVQZ//CCSD(T) = FC/cc-pVTZ single point level, which represents a severe test of the new surface. As a first application, on this analytical surface we perform an extensive dynamics study using quasi-classical trajectory calculations, comparing the results with recent experimental and theoretical data. The excitation function increases with energy (concave-up) reproducing experimental and theoretical information, although our values are somewhat larger. The OH rotovibrational distribution is cold in agreement with experiment. Finally, our results reproduce experimental backward scattering distribution, associated to a rebound mechanism. These results lend confidence to the accuracy of the new surface, which substantially improves the results obtained with our previous surface (PES-2000) for the same system. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4864358] I. INTRODUCTION
The construction of potential energy surfaces (PES) represents a process in continuous evolution, closely related to the development of high-level theoretical methods, functional forms more adequate to represent the nuclear motion, and the appearance of new experimental data. Hence, the PESs play a central role in the analysis of the kinetics and dynamics of reactive systems; and in the case of polyatomic systems this is no trivial task, and it is time-consuming. In 2000 our group1 reported an analytical surface for the title reaction (PES-2000), which was symmetric with respect to any permutation of the four methane hydrogens, and yielded rate constants in good agreement with the experimental data. Since this surface provides not only the energy of the reactive system but also the analytical gradients (i.e., the analytical first energy derivatives), in recent years it has been used as a testing bench for different dynamics methods. For instance, variational transition state theory with multidimensional tunneling effect (VTST/MT),1 multiconfiguration time-dependent Hartree (MCTDH),2 quasi-classical trajectory (QCT),3 reduced-dimensional quantum dynamics (QD),4 or more recently ring polymer molecular dynamics (RPMD).5 However, since we constructed this surface, we have observed some drawbacks in the symmetric PES-2000. First, the barrier height (13.0 kcal mol−1 ) is too low with respect to high-level ab initio calculations6 (14.1 kcal mol−1 ) and the top of the barrier is too broad (imaginary frequency of a) [email protected]
0021-9606/2014/140(6)/064310/12/$30.00
1549 cm−1 ) with respect to high-level ab initio calculations1 (imaginary frequency around 2000 cm−1 ). Second, experimentally this reaction yields a rotationally cold OH product,7 while PES-2000 gives a rotationally hotter OH product.8 Third, in the calibration process we used the LCT3 method9 to describe the tunneling effect of large curvature, but later, Truhlar et al. found that this method overestimates this effect, proposing new methods to correct this deficiency, LCT4,10 and LAT.11, 12 Finally, the presence of intermediate complexes in the entrance and exit channels were not considered. In addition, the most serious drawback is that this PES-2000 surface is semiempirical in nature, i.e., it was fitted to reproduce theoretical and experimental data. The present study of this reaction has been divided into two parts. In the first part, we report the construction of a new analytical surface for the title reaction that corrects the major deficiencies of PES-2000, which is based exclusively on highlevel ab initio calculations. It is named PES-2014 and is also symmetric with respect to the permutations of the H atoms. In the second part, as a test of quality, we report a dynamics study based on this surface. The O(3 P) + CH4 → OH + CH3 reaction is an interesting system from an experimental point of view, since it is of substantial importance in the chemistry of hydrocarbon combustion (high temperatures), and also from a theoretical point of view, because it presents a heavy-light-heavy mass combination which is a good candidate for large tunnelling effects (low temperatures). The thermal rate constants have been measured experimentally with a wide variety of methods.13–32 The most recent recommended expression for the rate
140, 064310-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 194.57.123.86 On: Mon, 07 Apr 2014 07:34:55
064310-2
González-Lavado, Corchado, and Espinosa-Garcia
constant for the temperature range 300–2500 K is:32 k(T) = 1.15 × 10−15 T1.56 exp(−4270/T) cm3 molecule−1 s−1 , although below 400 K, the value of k is considerably less reliable due to uncertainties in the reaction stoichiometry.30 The dynamics of this reaction (nascent vibrational distributions of the CH3 product, nascent rotational distributions of the OH product, and differential cross section) has also been widely studied with different techniques.7, 33–38 Due to the great experimental interest of this reaction, theoretically, it has received a great attention using different approaches: semiempirical direct dynamics,3, 8 reduced-dimensional quantum dynamics,39–43 and QCT calculations.3, 44–46 In parallel with the develop of dynamics methods, great efforts have been made in the construction of its potential energy surface,1, 6, 45, 47 highlighting the recent full-dimensional ab initio PES reported by Czakó and Bowman (CB) in 20126 based on a permutationally invariant fit of 17 212 accurate energy points obtained by an efficient composite method. In addition, this reaction represents a theoretical challenge because the approach of O(3 P) along a CH bond has 3-fold symmetry and leads to a Jahn-Teller conical intersection rather than a saddle point. The conical intersection corresponds to a 3 E state, and breaking the C3v symmetry splits this into two surfaces, 3 A and 3 A . Walch and Dunning48 and Schlegel et al.49 using ab initio molecular orbital calculations, obtained optimized geometries and frequencies at the stationary points, and the barrier height and heat of reaction. The saddle point geometry is found to be of Cs symmetry, but close to C3v and the predicted barrier heights for the two surfaces (3 A and 3 A ) present a very small difference (0.2 kcal mol−1 ). Therefore, we can assume that the PES2014 surface can describe either of the two surfaces. This approach has been employed in previous studies.1 In summary, in the present paper to correct the deficiencies of PES-2000, we report the construction of a new analytical surface for the title reaction based exclusively on highlevel ab initio calculations, named PES-2014. The paper is structured as follows. In Sec. II, high-level electronic structure calculations are outlined, with especial caution taken with respect to the barrier height, the complexes in the entrance and exit channels, and the topology of the reaction path from reactants to products. In Sec. III, a detailed description of the PES for the O(3 P) + CH4 hydrogen abstraction reaction and the fitting approach is presented. The results of the application of the fitting method are presented in Sec. IV and tested against theoretical electronic structure data. The dynamics results using quasi-classical trajectory calculations are presented in Secs. V and VI and compared with the available experimental evidence. Finally the conclusions are presented in Sec. VII. II. ELECTRONIC STRUCTURE CALCULATIONS
In this work we develop the PES-2014 to describe the hydrogen abstraction reaction in the title. For the fitting procedure the input information is exclusively based on high-level ab initio calculations describing the topology of the reaction from reactants to products, with special care taken in the description of the intermediate complexes.
J. Chem. Phys. 140, 064310 (2014)
We took the coupled-cluster singles and doubles approach, including a quasi-perturbative estimate of the connected triple excitations, correlating only valence electrons, CCSD(T), using the correlation consistent polarized triplezeta basis set, cc-pVTZ,50 abbreviated CCSD(T) = FC/ccpVTZ level, where FC means frozen core. We began optimizing and characterizing all the stationary point geometries (reactants, products, complexes in the entrance and exit channels, and saddle point). Starting from the saddle point geometry and going downhill to both the asymptotic reactant and product channels in mass-weighted Cartesian coordinates, we constructed the minimum energy path (MEP) at this level. In addition, we computed the gradients and Hessians for 60 points along the MEP. Since at this level the vibrational frequency calculations are numerical, to compute a Cartesian coordinate Hessian approximately 324 energies (324 is 18 squared, and 18 is 3 coordinates per atom times 6 atoms) are necessary. As 60 points were considered along the reaction path, 19 940 (324 times 60) energy points were calculated. By adding optimization and vibrational frequency calculations at the stationary points, we calculated about 20 000 energies at the CCSD(T)/cc-pVTZ level of theory, which were then used in the fitting procedure. Finally, to improve the energy description of the reactive system, we performed single-point calculations at a higher level, CCSD(T) = FULL/aug-cc-pVQZ (see below). The calculations were performed using the GAUSSRATE code,51 which serves as an interface between the GAUSSIAN09 systems of programs52 and POLYRATE-2010.53 Barrier height is undoubtedly one of the most difficult energy properties of a reaction to estimate accurately, and as Zhang and Liu38 noted, “the saddle-point geometry is somewhat uncertain. The predictions of the relative order of the lengths of the breaking C–H bond and the forming O– H bond obtained at various ab initio levels are not consistent.” We characterized the saddle point properties by using ab initio calculations at the CCSD(T)/cc-pVTZ level, and they are listed in Table I while Fig. 1 shows the optimized geometry. It is well known that the classic barrier height is strongly dependent on the basis set. Recently, Truhlar et al.54 used a new database to assess electronic structure methods for barrier heights. They showed that the costly CCSD(T) = FC/cc-pVTZ level presents a mean unsigned error (MUE) of 2.86 kcal mol−1 , and that the correlation of all electrons (FULL) and the use of a larger basis set have a major influence on the description of the barrier. Thus, using the single point calculation at the CCSD(T) = FULL/aug-cc-pVQZ level, we obtain a barrier height of 14.0 kcal mol−1 , i.e., 2.7 kcal mol−1 lower than that obtained with the CCSD(T) = FC/cc-pVTZ method. This lack of convergence when changing the level of accuracy, both method and basis, may be due to the existence of a conical intersection giving rise to the saddle point. This effect was also found in other related systems as F + OH → O(3 P) + HF, where the three states correlating to the O(3 P) intersect several times along the reaction path.55 These results illustrate the dramatic influence of electronic correlation and basis set on the correct description of the barrier. These will be the values used to fit the analytical PES (Sec. III). The single-point calculations with the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 194.57.123.86 On: Mon, 07 Apr 2014 07:34:55
064310-3
González-Lavado, Corchado, and Espinosa-Garcia
J. Chem. Phys. 140, 064310 (2014)
TABLE I. Properties of the reactants, products, and saddle point.a Geometry Parameter
Frequency
PES
Ab initio
PES
R(C–H)
