Advances in Physical ChemistryVolume 2012 (2012), Article ID 164752, 19 pagesdoi:10.1155/2012/164752

Review Article

Constructing Potential Energy Surfaces for Polyatomic Systems: Recent Progress and New Problems

J. Espinosa-Garcia, M. Monge-Palacios, and J. C. Corchado

Departamento de Química Física, Universidad de Extremadura, 06071 Badajoz, Spain

Received 11 May 2011; Accepted 14 July 2011

Academic Editor: Laimutis Bytautas

Copyright © 2012 J. Espinosa-Garcia et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Different methods of constructing potential energy surfaces in polyatomic systems are reviewed, with the emphasis put on fitting, interpolation, and analytical (defined by functional forms) approaches, based on quantum chemistry electronic structure calculations. The different approaches are reviewed first, followed by a comparison using the benchmark H + CH₄ and the H + NH₃ gas-phase hydrogen abstraction reactions. Different kinetics and dynamics properties are analyzed for these reactions and compared with the available experimental data, which permits one to estimate the advantages and disadvantages of each method. Finally, we analyze different problems with increasing difficulty in the potential energy construction: spin-orbit coupling, molecular size, and more complicated reactions with several maxima and minima, which test the soundness and general applicability of each method. We conclude that, although the field of small systems, typically atom-diatom, is mature, there still remains much work to be done in the field of polyatomic systems.

1. Introduction

At the heart of chemistry lies the knowledge of the reaction mechanism at atomic or molecular levels, that is, the motion of the nuclei in the potential field due to the electrons and nuclei. When the separation of the nuclear and electronic motions is possible, that is, within the Born-Oppenheimer (BO) approximation [1], the potential energy surface (PES) is the electronic energy (which includes electronic kinetic and Coulomb interaction energy) plus nuclear repulsion of a given adiabatic electronic state as a function of geometry. The PES is then the effective potential energy for the nuclear motion when the system is in that electronic state.

The actual knowledge of the chemical reactivity is based on quantum mechanical first principles, and the complete construction of the PES represents a very important challenge in theoretical chemistry. For small reactive systems (three or four atoms), the PES construction is a relatively mature field, although even today new surfaces are still being developed for atom-diatom systems [2–6] using high-level ab initio calculations. For instance, one of the latest atom-diatom surfaces has been developed in the present year by Jiang et al. [6] for the H + HBr hydrogen abstraction reaction. Ab initio calculations were performed based on the multireference configuration interaction (MRCI) method including Davidson’s correction using augmented correlation-consistent-polarized valence X-tuple zeta basis sets, with X up to 5, extrapolating the energies to the complete basis set (CBS) limit, and taking into account the spin-orbit correction. In order to define the complete configuration space for the BrH₂ system, ab initio calculations were performed at more than 12000 geometries. Then, a three-dimensional cubic spline interpolation was employed to yield the ground-state PES for this system. As can be seen, the computational effort (correlation energy + basis set + number of points) is enormous, and doing this for larger systems is still clearly prohibitive.

The extension to larger reactive systems (more than 4 atoms) is an open and promising field, although computationally very costly. Thus, for polyatomic systems of atoms, with internal degrees of freedom, if one needs configurations for each degree of freedom, one would need a total of evaluations of the energy. For instance, taking a typical value of of 10, for the H + CH₄ hydrogen abstraction reaction 10¹² evaluations would be needed. This implies too much computer time to perform ab initio calculations with chemical accuracy (1 kcal mol⁻¹) and becomes totally unaffordable if one aims to achieve spectroscopic accuracy (1 cm⁻¹).

The accuracy of the kinetics and dynamics description of a chemical reaction depends, aside from the quality of the PES, on the dynamics method used. If the motion of the nuclei on the PES is determined by the Hamilton equations in the phase-space configuration, the dynamics methods are classical or quasiclassical trajectory (QCT) approaches [7], and if the system is described by the Schrödinger equation (time dependent or time independent), the dynamics methods are quantum mechanical (QM) [8]. Finally, variational transition-state theory with multidimensional tunneling contribution (VTST/MT) can be derived from a dynamical approach by statistical mechanics and provides an excellent and affordable method to calculate thermal rate and state-selected constants [9].

In sum, electronic structure theory, dynamics methods, and constructed potential energy surfaces are the keystones of the theoretical study of chemical reactivity. Recent years have seen a spectacular development in these strongly related disciplines, paving the way towards chemical accuracy in polyatomic systems and toward spectroscopic accuracy in smaller problems. The scope of this paper is the construction of potential energy surfaces for bimolecular gas-phase polyatomic reactive systems in their electronic ground-state. Section 2 describes different methods of constructing potential energy surfaces, with especial focus on the methodological approach developed by our research group. Section 3 presents a library of PESs available for the scientific community, and Section 4 presents some results for the H + CH₄ benchmark polyatomic system and the H + NH₃ reaction, which permits us to compare results from different methods. Section 5 is devoted to analyzing the treatment of the peculiarities found in some reactions that entail additional challenges to the construction of potential energy surfaces, such as spin-orbit coupling, the molecular size, and the presence of various maxima and minima. Finally, Section 6 summarizes the main conclusions of the paper.

2. Constructing Potential Energy Surfaces

The construction of potential energy surfaces in reactive systems began with the dawn of Quantum Chemistry when different quantum approaches were used, basically empirical or semiempirical. In present days, the data needed for this construction are obtained from high-level ab initio calculations, although density functional theory (DFT) methods are also widely used, especially in the case of larger molecular systems.

As is well known, and this special issue is a clear indication of it, the construction of PESs has a long tradition in theoretical and computational chemistry, and an exhaustive review of the literature is beyond the scope of the present paper (see, for instance, a recent review in [10]). Here we only highlight some important contributions with the aim that the reader will get an overall perspective on this wide field of research.

The most straightforward procedure to describe a reactive system is from electronic structure calculations carried out “on the fly,” sometimes also called “direct dynamics” [11], where every time the dynamics algorithm requires an energy, gradient, or Hessian, it is computed from electronic structure calculations or molecular mechanics methods. The advantage of this approach is the direct application of electronic structure calculations to dynamics problems without intermediaries, but when chemical accuracy is needed, very high-level ab initio calculations are required, and hence the time and computational cost are very high. For polyatomic systems, such an approach is still prohibitive. A recent example is the direct dynamics trajectory study of the formaldehyde cation with molecular hydrogen gas-phase reaction [12]. Even using a modest MP2/6-311G(d,p) ab initio level, which yields a barrier only 0.5 kcal mol⁻¹ below the benchmark value, the ≈9600 trajectories required about 2.5 CPU years on a small Linux cluster. In this regard, one has to keep in mind that QCT calculations require a large number of trajectories to adequately sample the events and obtain results that from a purely statistical point of view are tolerable, that is, with very small errors. Note that the time devoted to the QCT calculation is negligible in this time scale.

The alternative to direct dynamics calculations is to construct a mathematical PES, that is, developing an algorithm that can provide the potential energy for any given geometry of the system without depending on “on-the-fly” electronic structure computations. In this sense, three basic approaches have been considered: fitting, interpolation, or analytical (defined by functional forms). In fact, the fitting and analytical procedures share the idea of functional form, and this analytical function needs to be fitted to theoretical information. This classification is therefore arbitrary, and other possibilities can be found in the literature [10]. We think, however, that it can help the reader get a better general view of this field. The interpolation methods have been widely applied in the construction of potential energy surfaces, and different strategies have been developed [9, 10, 13–15]. Thus, Collins et al. [14, 16, 17] developed and applied a method for generating potential energy surfaces using a modified Shepard interpolation. The PES at any configuration () is represented by a weighted average of Taylor series : where is the number of molecular configurations whose energy and its first (gradient) and second (Hessian) derivatives have been evaluated, is the normalized weighting factor, is a set of internal coordinates, and the Taylor series are expanded around the point : The use of the energy and its first and second derivatives makes the method very robust, but computationally very expensive, especially if very high-level ab initio calculations are needed for a correct description of the reactive system. For instance, very recently a new PES for the H + CH₄ gas-phase reaction has been developed by Collins by interpolation of 30000 data points obtained with high-level ab initio calculations [18]. The ab initio calculations took about 600 days on a workstation with 8 CPU cores in total, while the evaluations of the potential values on all the grids for this reaction took 400 days. Three advantages of this method are the following: first, that the interpolated function agrees precisely with the values of the data points; second, that new electronic structure points may be incorporated in the process to improve the PES; third, that the PES is invariant to the permutation or exchange of indistinguishable nuclei. However, Shepard interpolation often yields slight oscillations in equipotential contours and vibrational frequencies (second derivatives of the energy); that is, the potential energy surface is not smooth enough to provide smooth changes in the vibrational frequencies. This is a limitation when one wants to compute tunneling effects at low energies, since the barrier to tunnel through is unrealistic.

Fitting methods, such as the interpolated ones considered above, have been widely used in the construction of potential energy surfaces [19–21], and they have been recently reviewed [10]. Thus, for instance, a fitting procedure widely used is based on the many-body expansion (MBE) method [22], amplified by Varandas [23] in the so-called double many-body expansion (DMBE) method. The potential is given by a sum of terms corresponding to atoms, diatoms, triatoms, and tetra-atoms and a series of adjustable coefficients. Another fitting procedure has been developed and applied by Bowman et al. [24–27], in which the ab initio data are globally fitted to a permutational symmetry invariant polynomial. The function describing the PES has the form where is an -dimensional vector, depending on internuclear distances, with components , being the distance between nuclei and , and given by and and are polynomials constructed to satisfy the permutation symmetry with respect to the indistinguishable nuclei. For instance, for the CH₅   + polyatomic system [28], the PES was least-squares fitted to 20 639 ab initio energies, obtained at the MP2/cc-pVTZ level. This fit contains 2303 coefficients and an rms fitting error of 51 cm⁻¹. Bowman and coworkers have developed 16 of such PESs for polyatomic systems, and they have been recently revised [27]. In general, the fitting approach is linear least squares, and therefore these surfaces do not exactly reproduce the ab initio data.

The third alternative consists of the analytical surfaces defined by functional forms. In this method the ab initio data are fitted to a valence bond (VB) functional form, augmented with molecular mechanics (MM) terms which give great flexibility to the potential energy surface. These VB/MM functional forms have a long history in the development of potential energy surfaces [29–33], although at first the surfaces were semiempirical in the sense that theoretical and experimental data were used in the fitting procedure. This has been the methodological approach developed by our group [10, 34, 35].

The first surfaces were developed for the H + CH₄ hydrogen abstraction reaction, as a paradigm of polyatomic systems [29–33]. The potential energy for a given geometry, , is given by the sum of three terms: stretching potential, , harmonic bending term, , and anharmonic out-of-plane potential, ,

The stretching potential is the sum of four London-Eyring-Polanyi (LEP) terms, each one corresponding to a permutation of the four methane hydrogens: where is the distance between the two subscript atoms, Hi stands for one of the four methane hydrogens, and HB is the attacking H atom. Although the functional form of the LEP potential is well known, we will recall it for the sake of completeness:where there are 12 fitting parameters, four for each of the three kinds of bond, , , and . In particular, these are the singlet and triplet dissociation energies, and , the equilibrium bond distance, , and the Morse parameter, . The Morse parameter for the CHi bonds, , however, is not taken as a constant but rather as a function of the distances, with being the average distance, In this way, changes smoothly from its value at methane, , to its value at the methyl radical, , as the reaction evolves. Therefore, 14 parameters are required to describe the stretching potential.

One of the problems with this functional was that the equilibrium C–H distances for the reactants, saddle point, and products are the same, leading to a very rigid surface. Chakraborty et al. [36], for the H + C₂H₆ reaction, included a modification to endow the surface with greater flexibility. The reference C–H bond distance is transformed smoothly from reactant to product using the following equation: where is which is symmetric with respect to all the four hydrogen atoms and goes to zero as one of the hydrogen atoms is abstracted, and is a geometry-dependent switching function, given by where and are adjustable parameters. Therefore, this adds 2 new parameters (total 16 parameters) to describe the stretching potential.

The term is the sum of six harmonic terms, one for each bond angle in methane: where and are force constants and are the reference angles. The force constants are allowed to evolve from their value in methane, , to their value in methyl, , which are two parameters of the fit, by means of switching functions: while is a function of both the and distances: Thus, four adjustable parameters are involved in the definition of .

The reference angles are also allowed to change from their value at methane, , or , to methyl, 120° or radians, by means of switching functions [37]: Finally, the switching functions are given by involving 10 more adjustable parameters. In total, 16 terms need to be fitted for the calibration of the potential.

The potential is a quadratic-quartic term whose aim is to correctly describe the out-of-plane motion of methyl: The force constants, and , have been incorporated into a new switching function which is such that vanishes at the methane limit and which directs the change of the methyl fragment from pyramidal to planar as the reaction evolves: with , , , and being the only parameters of that enter the fitting process. is the angle that measures the deviation from the reference angle: where , , , and are vectors going from the carbon atom to the , , , and l hydrogen atoms, respectively, and are the reference angles defined in (16). The first term to the right of (20) is therefore the angle between the CHi bond and a vector perpendicular to the plane described by the , , and hydrogen atoms and centred at the atom. To correctly calculate , the motion from to has to be clockwise from the point of view of the atom.

In the case of the H + CH₄ reaction the CH₃ radical product is planar. What happens if the product presents a nonplanar geometry? In these cases of CX₃ products, we have modified the reference angle in (13) and (20). Thus, the original expression, where , is replaced by where is the bending angle in the non-planar product, and is related to by the expression, In these cases of non-planar products, CX₃, this correction would add new parameters in the fitting procedure.

The PES, therefore, depends on at least 36 parameters, 16 for the stretching, 16 for the harmonic term, and 4 for the out-of-plane potential. These 36 parameters give great flexibility to the PES, while keeping the VB/MM functional form physically intuitive.

In the case of the H + CH₄ reaction the CH₄ reactant presents symmetry. What happens if the reactant presents a different symmetry? For instance, ammonia presents symmetry , characterized by an inversion mechanism through a planar structure with symmetry . The functional form (5) used for the H + CH₄ reaction cannot be applied without modification to any kind of system. Thus, when this potential is applied to the study of the H + NH₃ → H₂ + NH₂ reaction, a major drawback was observed [38], namely that it wrongly describes the NH₃ inversion reaction, predicting that the planar ammonia ( symmetry), which is a saddle point to the ammonia inversion, is about 9 kcal mol⁻¹ more stable than the pyramidal structure ( symmetry).

In the original expression for the H + CH₄ reaction (5), the term was added to obtain a correct description of the out-of-plane bending in the methyl radical, 580 cm⁻¹. However, Yang and Corchado [38] noted that this term leads to unphysical behaviour along the ammonia inversion path. To avoid this drawback, our laboratory recently developed [39] a new PES where the term is removed: Consequently, each of the 16 parameters is required for the stretching terms and the harmonic bending terms of the PES which were fitted to high-level ab initio CCSD(T)/cc-pVTZ calculations.

In sum, starting from a basic functional form, this form must be adapted to each particular case, looking for the greatest flexibility and suitability for the problem under study. This has been the main aim of the modifications on the original functional form (5), described by (10)–(12), (22), (23), and (24). While this could represent a disadvantage with respect to the fitting or interpolation methods, because new functional forms must be developed in each case, it also represents an advantage, because simple stepwise modifications can give great flexibility to the surface.

Once the functional form is available, the fitting procedure is started. A very popular approach for fitting a function is the least-squares method, which, using some local optimization algorithm, gives values of the parameters that minimize (locally) the function where is the ab initio energy associated with a particular molecular configuration specified by “” and is the energy predicted by the analytical function at the same molecular configuration, which depends on a set of parameters denoted as .

One must note, however, that any fitting procedure has certain limitations. First, due to the large number of parameters, it is very hard to find a global minimum for the fit, which accurately describes the entire surface. Second, due to the nature of the linear least-square method, the result is dependent on the initial parameters; third, since we use a mathematical approach without physical intuition, one usually obtains a number of distinct sets of mathematical parameters, all equally probable and good at reproducing the complete system. To make matters worse, as noted by Banks and Clary [40], while the differences between these parameter sets may be small, the dynamics information obtained from them can vary notably. This is especially true because small changes in the energy derivatives can cause large changes in dynamical properties of the PES. One therefore requires an extremely accurate fit to the topology of the given points to obtain an acceptable set of parameters.

To solve at least partially some of the above problems, we adopt a different approach [41]. Firstly, we will try to obtain the values of the parameters that minimize the function where denotes the gradients (first derivative of the energy), the Hessian elements (second derivatives of the energy), and , , and are weights. In practice, the gradients enter the fitting process only at the stationary points, and the Hessian elements by means of the harmonic vibrational frequencies at selected points. Secondly, the choice of the number of points to be fitted is critical, and this strategy allows this number to be reduced drastically (see more details in the original paper [41]). We begin by choosing the stationary points (reactants, products, saddle point, and intermediate complexes) as data points. The geometry, energy, and frequencies are fitted and, thus, indirectly, the first and second derivatives of the energy. It is necessary to remember that if the information included was only the energy, a mesh of points would be required. However, when the information included is not only the energy at a stationary point but also the first and second derivatives, only one point is needed to reproduce each stationary point configuration. Moreover, additional points on the minimum energy path are included. Thus, with the inclusion of the stationary points and representative points on the reaction path, we search for a good reproduction of the topology of the path connecting reactants to products. Finally, we add the energy of a point not on the reaction path to describe zones of the reaction valley relevant to tunneling dynamics and rotational excitation of the products. The fitting procedure is currently automated on a computer with no user intervention except to analyze the results and, when called for, to change the weights.

This strategy presents certain advantages over other methods. The first is transferability of the functional form and the fitted parameters. Since one can regard VB/MM as some kind of highly specific MM force field, it can be expected that some of the fitting parameters are transferable to a similar system, although obviously the parameter values are system specific. In our group we have used this feature to obtain, for example, the PES for the F + CH₄ reaction [42] using as the starting point the analytical PES for H + CH₄ [43]. Secondly, the PESs are guaranteed to some extent to have the capability of reproducing the energetic interactions of the chemical system in region not included in the fit. For example, high-energy collisions might require knowledge of the PES at very high energies that may not have been sampled in the fitting process. An interpolation method cannot ensure the absence of spurious wells or unphysical repulsive regions when nearby points are absent, while MM force fields can. Thirdly, VB/MM PES parameters can be refitted so as to fine-tune the PES using additional data (e.g., higher ab initio calculations at selected points) at a low computational cost. Furthermore, VB/MM surfaces and their energy derivatives (when they can be analytically calculated) are usually smoother than interpolated surfaces. As was mentioned above, Shepard interpolation sometimes gives discontinuities in the derivatives, which is a highly undesirable feature for kinetics and dynamics studies.

Thus, our research group has developed economical alternatives for constructing analytical PESs of polyatomic systems, which basically are VB/MM-type surfaces. This strategy of seeking an optimal tradeoff of time and computational cost has been successfully used in several gas-phase hydrogen abstraction reactions of five [44, 45], six [42, 43, 46–52], and seven [53] atoms. In general, a good correspondence between kinetics and dynamics theoretical results and experimental measurements is found. In the first phase of our research, theoretical and experimental information was used in the fitting procedure, making the analytical surface semiempirical in nature, which was a serious problem and a limitation in kinetics and dynamics studies. However, in the last few years, our surfaces have been fitted exclusively to very high-level ab initio calculations, avoiding the aforementioned limitations.

3. Library of Potential Energy Surfaces

The FORTRAN codes and the fitted parameters of the polyatomic potential energy surfaces developed by our group are available, free of charge, for the scientific community, and can be downloaded from the POTLIB library [54]: http://comp.chem.umn.edu/potlib/ or http://users.ipfw.edu/DUCHOVIC/POTLIB2001, or obtained from the authors upon request.

4. Applications

(a) The H + CH₄ ReactionThe gas-phase H + CH₄ → H₂ + CH₃ hydrogen abstraction reaction, as well as its deuterated isotopomers, is the prototype polyatomic reactive system and has been widely studied both theoretically and experimentally [10]. The construction of its PES has a long history, and it is one of the few reactive systems for which different approaches to the construction have been developed, which permits a direct comparison.We will focus attention, except otherwise stated, on the three most recent and accurate surfaces for this reactive system, which were constructed with different strategies. Chronologically, in 2006 Zhang et al. [55, 56] developed the family of ZBBi surfaces, using the invariant polynomial method, based on the fitting to more than 20000 ab initio energies at the RCCSD(T)/aug-cc-pVTZ level. In 2009 we developed [41] an analytical PES, CBE surface, which is a VB/MM functional form, and the 36 parameters are fitted using exclusively high-level electronic structure calculations at the CCSD(T)/cc-pVTZ level. Very recently, in 2011, Zhou et al. [18] constructed a full-dimensional surface using the modified Shepard interpolation scheme based on 30000 data points at the CCSD(T)/aug-cc-pVTZ level—the ZFWCZ surface. The energy, geometry and vibrational frequencies of the saddle point are summarized in Table 1 for the three surfaces. Comparing the wide range of properties, we conclude that the three surfaces present excellent agreement, with small differences in the CBE with respect to the other two surfaces in the <H–C–H′ bending angle, ≈3°, and in the ZBB3 compared to the other two surfaces in the imaginary frequency, ≈100 cm⁻¹. The three surfaces present a colinear saddle point, 180.0°, with a similar barrier height, with small differences, ±0.2 kcal mol⁻¹, that is, within the chemical accuracy, and very close to the best estimation predicted at the CCSD(T)/aug-cc-pVQZ level, 14.87 kcal mol⁻¹ [57].The best macroscopic measure of the accuracy of a potential energy surface is probably the rate constant, at least in the thermal bottleneck region. When we compare different dynamics methods using the same surface, the dynamics approach is tested, and when we compare theoretical and experimental results, both the dynamics method and the surface are tested. Figure 1 plots the thermal rate coefficients computed with accurate fulldimensional quantum dynamics approaches, on the analytical CBE surface; a very accurate interpolation surface developed by Wu et al., WWM surface [57–59], which used the modified Shepard interpolated method developed by Collins et al., based on CCSD(T)/cc-pVQZ or CCSD(T)/aug-cc-pVQZ ab initio calculations; and a fitted surface, named ZBB2, an earlier version of the ZBB3 surface developed by Bowman et al. In the same figure, there also appear the results obtained with the variational transition-state theory with multidimensional tunneling effect [41] and the experimental data [60] for comparison. First, with quantum dynamics approaches, the rate coefficients obtained with the three surfaces agree almost perfectly in the common temperature range. Second, when the rate coefficients are obtained using different dynamics approaches, VTST/MT and MCDTH (multiconfigurational time-dependent Hartree approach [58, 59]) on the same CBE surface, excellent agreement is found, and both reproduce the experimental information in the common temperature range. These results indicate, first, that VTST/MT is a powerful and computationally economic tool for the kinetics study of polyatomic systems, with results comparable to those obtained with computationally more expensive quantum dynamical methods. Second, the CBE surface is accurate at least in the region of low energies, which is the most relevant region for thermal rate coefficient calculations and, therefore, provides a satisfactory description of the reaction path and the transition-state region.Next, we analyze some dynamics properties where different surfaces have been used and compared with the sparse experimental data. We focus on the H + CD₄ → HD + CD₃ reaction because there is more experimental information available for comparison. The product energy partitioning has been experimentally obtained by Valentini’s group. [61] only for the vibration and rotation of the HD product, and this same group found that more than 95% of the HD product is formed in the and vibrational states. These results appear in Table 2, together with the QCT results on the analytical CBE and the fitted ZBB1 [55] surfaces. The QCT calculations on these two PESs show excellent agreement for the HD product energy partitioning as well as for the HD vibrational distribution. In this latter case, the QCT results reproduce the experimental evidence, but, in the first case, they strongly contrast with the experimental measurements, which measure an internal excitation of the HD product of 7% and 9% for vibration and rotation, respectively. The agreement between the results from the two surfaces leads us to think that the discrepancies with experiment are mainly due to the dynamical method, that is, to limitations of the QCT approach. Experimental problems, however, cannot be totally ruled out, as was recently observed by Hu et al. [62], who suggested that the conclusions from Valentini et al.’s CARS experimental study might need to be reinterpreted.The product angular distribution is, doubtless, one of the most sensitive dynamics features with which to test the quality of the potential energy surface, but experimentally it is very difficult to measure in some cases. When the Photoloc technique is used, the laboratory speed depends on both the scattering angle and the speed of the CD₃ product, which is influenced by the HD coproduct internal energy distribution. Uncertainties in this quantity could produce errors in the scattering angle.Camden et al. [63] reported the first study of the state-to-state dynamics differential cross-section at high energies (1.95 eV) for the H + CD₄ gas-phase reaction using the Photoloc technique. They found that the CD₃ products are sideways/forward scattered with respect to the incident CD₄, suggesting a stripping mechanism (note that in the original papers [63–65] the CD₃ product is measured with respect to the incident H). Later, this same laboratory [64, 65] reported new experimental studies, also at high energy (1.2 eV), finding the same experimental behaviour. Experimentally, state-to-state dynamics studies are difficult to perform at low energies for the title reaction, because the H atoms, which are produced in a photolysis process, are hot. Only very recently have Zhang et al. [66] reported crossed molecular beam experiments for the H + CD₄ reaction at lower collision energies, ranging from 0.72 to 1.99 eV. Note that the lower value, 0.72 eV, is close to the barrier height, and consequently its dynamics will be influenced mainly by the transition-state region. Figure 2(a) plots these experimental results. A backward angle is clearly observed at low energies (rebound mechanism), changing towards sideways when the energy increases (stripping mechanism).This is an excellent opportunity to test the quality of the PES and the dynamics methods (Figures 2(b)–2(d)). Only two surfaces have been used to study this problem theoretically: an analytical surface from our group (versions 2002 [43, 67] and 2008, labeled here as CBE [41, 68]), using both QCT and QM calculations, and an “on-the-fly” B3LYP/631G(d,p) density functional theory surface [65], using QCT calculations. At low energies, ≈0.7 eV, our analytical surfaces, PES-2002 and CBE, using QCT and QM methods (Figures 2(b) and 2(c)), show backward scattering, associated with a rebound mechanism, reproducing the recent experimental data [66]. The B3LYP “on-the-fly” surface using QCT calculations (Figure 2(d)) yields a more sideways scattering, with large uncertainties, in contrast with experiment. These differences could be due to the poor statistics on the B3LYP surface [65] and to the severe underestimation of the barrier heights, about 5 kcal mol⁻¹ lower than the best ab initio calculations. This low barrier artificially permits reactive trajectories with larger impact parameters, favouring the sideways scattering region. When this error in the barrier height is corrected, our PES-2002 and, more noticeably in the better CBE surface, the low impact parameters are favoured, and this explains the rebound mechanism. Note that, interestingly, this experimental behaviour was already predicted by our group in 2006 [67] before the experimental data were available.At higher collision energies, as they increase from 1.06 to 1.99 eV, Zhang et al. [66] found a shift of the product angular distribution from backwards to sideways. This behaviour agrees qualitatively with the previous observations of Camden et al. [63–65], although these latter workers found a clear sideways distribution at 1.21 and 1.95 eV, with practical extinction of the backward signal. This may have been due to the use of the Photoloc technique, which neglects the internal energy contribution of the HD co-product, and, as the authors themselves recognized in 2006, “clearly a more detailed picture of the differential cross section is desirable but it will have to await more experimental work.” The QCT results using the analytical CBE surface (Figure 2(b)) reproduce the new experimental evidence. QM calculations on the same CBE surface [68] (Figure 2(c)) give more sideways scattering than the QCT calculations and experiment, although it is not clear whether this is due to significant quantum effects or simply are an artifact of the reduced dimensionality approach in the QM calculations.Note that QCT calculations on the old PES-2002 surface also predicted the subsequently experimentally observed behaviour [66]—backwards-sideways scattering. However, they contradicted the experimental information available at the time the study was carried out—in 2006. In addition, since direct-dynamics QCT calculations at the B3LYP/631G(d,p) level showed sideways scattering, reproducing the experimental measurements [63–65], the validity of the analytical PES-2002 was questioned. The recent experiments of Zhang et al. [66], however, changed the whole picture and now the direct-dynamics QCT calculations at low collision energies are questioned because of, as noted above, the presence of reactive trajectories with erroneously large impact parameters due to the underestimation of the barrier height. In the case of higher collision energies, however, this error in the barrier should be of less concern than at lower energies.Very recently, Zhou et al. [18] have performed an exhaustive analysis of the total reaction probabilities and integral cross-section using quantum dynamics calculations on the three most recent surfaces: fitted ZBB3, analytical CBE, and interpolated ZFWCZ surfaces, with collision energies up to 1.7 eV. Figure 3 plots the integral cross-section results. At collision energies up to 1.0 eV, the three surfaces show satisfactory agreement, while at higher energies of collision, the CBE surface overestimates this dynamics property, while the remaining two surfaces show good agreement. This could be attributed to deficiencies of the CBE surface at high energies. This is not surprising since the calibration of the CBE surface was done thinking of thermal behaviour. The information used during the fit focused on the reaction path and reaction valley, and higher energy areas were neither sampled nor weighted sufficiently. Therefore, as the collision energy increases, the accuracy of the CBE surface diminishes.Finally, another severe test of the quality of the PES is the study of the effect of the vibrational excitation on the dynamics. In fact, the dynamics of a vibrationally excited polyatomic reaction presents a challenge both theoretically and experimentally. Camden et al. [69] carried out the first experimental study on the effect of the C–H stretch excitation on the gas-phase H + CH₄ hydrogen abstraction reaction. They found that the excitation of the asymmetric C–H stretch mode enhances the reaction cross-section by a factor of 3.0 ± 1.5 with respect to the ground-state methane, and this enhancement is practically independent of the collision energies for the three cases analyzed—1.52, 1.85, and 2.20 eV. In the following year, 2006, two theoretical papers on this issue were published: one by Xie et al. [56] using QCT calculations on the ZBB3 surface, and another from our laboratory [70] also using QCT calculations on the older analytical PES-2002 surface. At 1.52 eV, the theoretical results are close, with computed enhancement factors of 2.3 and 1.9, respectively, both of them within the experimental uncertainties. Note that our more recent CBE surface also predicts an enhancement factor of 1.9 (unpublished results).

Table 1: CH₅ saddle point properties.a

Fitting (ZBB3)b
Analytical (CBE)c
Interpolated (ZFWCZ)d

Barrier height
14.78
15.01
15.03

     Geometry

R(C–H′)
1.409
1.390
1.399

R(H′–HB)
0.901
0.973
0.895

<H-C-H′
102.8
106.8
103.1

<C–H′–HB
180.0
180.0
180.0

     Vibrational frequencies

3289
3173
3223

3289
3173
3223

3104
3036
3073

1726
1833
1784

1434
1443
1442

1434
1443
1442

1078
1173
1115

1078
1173
1115

1062
1085
1062

440
542
522

440
542
522

1320i
1488i
1467i

Energy in kcal mol⁻¹, geometry in and degrees, and vibrational frequency in cm⁻¹; [55]; [41]; [18].

Table 2: HD product energy partitioning and vibrational distribution (percentages) for the H + CD₄ reaction at 1.52 eV.

Surface
(HD)
(HD)
HD ()
HD ()

Analyticala
18
18
84
15

Fittedb
22
18
78
22

Exp.c
7
9
≥95

CBE surface, [41] ZBB1 surface, [55] Valentini et al. [61].

Figure 1: Arrhenius plots of (cm³ molecule⁻¹ s⁻¹) for the forward thermal rate coefficients of the H + CH₄ reaction against the reciprocal of temperature (K), in the range 250–1000 K. Black line: MCDTH quantum calculations on the CBE surface; red line: MCDTH quantum calculations on the WWM surface; blue line: quantum calculations on the ZBB2 surface; black dashed line: VTST/MT calculations on the CBE surface; crosses: experimental values [60].

(a)

(b)

(c)

(d)

Figure 2: CD₃ product angular distribution (with respect to the incident CD₄) for the H + CD₄ → HD + CD₃ reaction at different collision energies. (a) Experimental results from [64]; (b) QCT angular distribution on the CBE analytical surface [66]; (c) QM angular distribution on the CBE analytical surface [66]; (d) QCT angular distribution on the B3LYP “on-the-fly” surface [65].

Figure 3: Quantum mechanical integral cross-section () versus collision energy (eV) for the H + CH₄ → H₂ + CH₃ reaction on the CBE (black line), ZBB3 (blue line), and ZFWCZ (red line) surfaces.

(b) The H + NH₃ ReactionThe reaction of hydrogen atom with ammonia is a typical five-body reactive system, and presents a rare opportunity to study both intermolecular and intramolecular dynamics. For the intermolecular case, the gas-phase H + NH₃ → H₂ + NH₂ hydrogen abstraction reaction is similar to the H + CH₄ reaction. It presents a barrier height of 14.5 kcal mol⁻¹ and a reaction exoergicity of 5.0 kcal mol⁻¹, as compared to 14.87 and an endoergicity of 2.88 kcal mol⁻¹, respectively, for the H + CH₄ analogue. Also, the inversion of ammonia between two pyramidal structures ( symmetry) passing through a planar structure ( symmetry, which is a saddle point) is an example of intramolecular dynamics [71] with a barrier height in the range 5.20–5.94 kcal mol⁻¹, depending on the level of calculation.Only three potential energy surfaces have been developed for the H + NH₃ system. In 2005, Moyano and Collins [72] developed an interpolation potential energy surface for the ammonia inversion, and the hydrogen abstraction and exchange reactions of H + NH₃ using a modified Shepard interpolated scheme based on 2000 data points calculated at the CCSD(T)/aug-cc-pVDZ level (PES1 version) or as single-point calculations at the CCSD(T)/aug-cc-pVTZ (PES2 version) level. The first and second derivatives of the energy were calculated by finite differences in the energy. Previously, in 1997, our group constructed the first surface for the hydrogen abstraction reaction exclusively, CE-1997 [44], which was fitted to a combination of experimental and theoretical information that is, it was semiempirical in nature, which represents a limitation for dynamics studies. Moreover, as was previously noted, recently Yang and Corchado [38] reported a major drawback of the CE-1997, namely that it describes incorrectly the NH₃ inversion motion, predicting incorrectly that the planar ammonia ( symmetry) is about 9 kcal mol⁻¹ more stable than the pyramidal structure ( symmetry). To correct this behaviour of the NH₃ inversion, together with its semiempirical character, a new analytical potential energy surface, named EC-2009, was recently developed by our group [39], describing simultaneously the hydrogen abstraction and ammonia inversion reactions. This EC-2009 surface is basically a valence bond-molecular mechanics (VB/MM) surface, given by (24), and was fitted exclusively to very high level ab initio calculations at the CCSD(T)/cc-pVTZ level. Note that the first derivatives of this surface are analytical, which implies a significant reduction in the computer time required for dynamical calculations as well as more accurate derivatives than is possible with numerical methods.We begin by analyzing the hydrogen abstraction reaction. Table 3 lists the energy, geometry, and vibrational frequencies of the saddle point, with the interpolated MC (PES2 version) and the analytical EC-2009 surfaces. Using as target the CCSD(T)/cc-pVTZ ab initio level, the two surfaces present similarities, although important differences must be noted. First, while the interpolated MC surface reproduces the ab initio N–H′–H bend angle, the analytical surface yields a collinear approach. However, in previous papers [73, 74], we demonstrated that this was not a serious problem in the kinetics and dynamics description of the system. Second, the imaginary frequency obtained with the MC surface differs by about 300 cm⁻¹ from the ab initio value.For the two surfaces, Figure 4 shows the energy along the minimum energy path () and the ground-state vibrationally adiabatic potential curve (), which is defined as the sum of the potential and vibrational zero-point energies along the reaction path. Note that is the reaction coordinate, being zero at the saddle point, positive in the product channel, and negative in the reactant channel. Firstly, the curve is smooth and shows no oscillations in either surface. However, shows significant oscillations when computed for the interpolated surface. The reason is that the potential energy surface is not smooth enough to provide smooth changes in the frequencies (second derivatives of the energy). The situation becomes worse as one moves away from the saddle point and approaches areas where little ab initio information for interpolation is available. Thus, the bump in at around  a.u. is due to wild oscillations of the frequencies that lead to unphysical values of for . This is a limitation when one wants to compute the tunneling effect at low energies (see below), since the barrier to tunnel through is totally unrealistic.Figure 5 plots the thermal rate coefficients computed with the VTST/MT approach on both surfaces, where tunneling was estimated using the least-action tunneling (LAT) method [75, 76], together with experimental values [77] for comparison. In the common temperature range, 490–1780 K, both surfaces reproduce the experimental data, which is a test of both the surface and the dynamical method. The differences between the rates for the two surfaces are relatively small, 40% at 600 K, and diminishing as temperature increases. At low temperatures, the differences between the two surfaces increase, with the rates computed on the MC surface being 92% larger than the ones with the EC-2009 surface at 200 K. In this low temperature regime, where tunneling is important, the analytical EC-2009 surface is more accurate due to the unrealistic oscillations in the adiabatic reaction path on the MC surface, which negatively influence the tunneling.The integral cross-sections in the range 10–30 kcal mol⁻¹ have been evaluated using quasi-classical trajectory (QCT) calculations on both surfaces [72, 78] and are plotted in Figure 6. As apparent from this figure, both surfaces agree reasonably, showing typical threshold behaviour, starting from 10 kcal mol⁻¹ and increasing with the translational energy. Unfortunately, there is no experimental data for comparison, and we think that these theoretical results might stimulate experimental work on this little studied system.Finally, the angular distributions of the H₂ product with respect to the incident H atom has only been determined using the EC-2009 surface [78]. Figure 7 plots this property for collision energies of 25 and 40 kcal mol⁻¹. At 25 kcal mol⁻¹, the scattering distribution is in the sideways-backward hemisphere, associated with a rebound mechanism and low impact parameters. When the collision energy increases, 40 kcal mol⁻¹, the scattering is shifted slightly towards the sideways hemisphere, due to larger impact parameters.As was mentioned above, the EC-2009 surface describes, in addition to the aforementioned hydrogen abstraction reaction, the ammonia inversion, an example of interesting intramolecular dynamics. Figure 8 plots the equipotential contours in the two significant coordinates, r(N–H) bond length and <H–N–H angle, for the analytical EC-2009 PES and the CCSD(T)/cc-pVTZ ab initio level, for comparison. As seen, the contours are smooth and the analytical PES reproduces the fitted ab initio data points [39].A very stringent test of the quality of this surface is the ammonia splitting, which demands spectroscopic accuracy. The inversion motion is represented by a symmetric double-well potential. As a consequence of the perturbation originating this double well, a splitting of each degenerate vibrational level into two levels appears, , due to quantum mechanical tunneling [71], and the splitting increases rapidly with the vibrational number. The splitting can be computed from the tunneling rate of inversion, , which is in turn obtained from the imaginary action integral, , using the WKB approximation [79]: with being the speed of light and the eigenvalue associated with the inversion mode of ammonia, 1113 cm⁻¹. Figure 9 plots this path for ammonia where the first two pairs of split eigenvalues are superimposed. The computed splitting is listed in Table 4 together with experimental values [80] for comparison. The EC-2009 results overestimate the experimental values. The reason for the discrepancy mainly lies in the shape of the PES, although other factors such as the tunneling calculation cannot be discarded. Indeed, we found [39] that the overestimation we obtain in tunneling splitting is due to our barrier to inversion being slightly lower and thinner than those of other studies. Unfortunately, tunneling splitting is so sensitive to the shape of the PES that a few tenths of kcalmol⁻¹ give rise to a factor of four in the computed splitting. With respect to the effect of isotopic substitution, for the ND₃ case we obtain a value  cm⁻¹. Although this is greater than the experimentally reported value, 0.05 cm⁻¹ [81], it correctly predicts about one order of magnitude reduction upon isotopic substitution. In sum, despite the enormous effort required for the construction of the potential energy surface, it does not suffice to obtain spectroscopic accuracy, and our surface can only give a qualitative description of the splitting in ammonia.

Table 3: NH₄ hydrogen abstraction saddle point properties.a

Interpolated (MC)
Analytical (EC-2009)
Ab initio

Barrier height
14.64
14.48
14.73

     Geometry

R(N–H′)
1.323
1.279
1.308

R(H′–HB)
0.900
0.868
0.890

<C–H′–HB
160.8
180.0
158.4

    Vibrational frequencies

3450
3444
3478

3350
3373
3384

2057
1861
1888

1527
1623
1566

1174
1497
1280

1047
1080
1063

650
622
677

629
581
506

1371i
1602i
1662i

Energy in kcal mol⁻¹, geometry in and degrees, and vibrational frequency in cm⁻¹; [72]; [39]; [39], at the CCSD(T)/cc-pVTZ level.

Table 4: Eigenvalues of ammonia inversion (in cm⁻¹).

ν ₂
EC-2009a
Exp.b

0⁺
0.00
0.00

0⁻
3.64
0.79

1⁺
861.39
932.43

1⁻
962.90
968.12

2⁺
1533.36
1598.47

2⁻
1942.29
1822.18

a[39]; b[80].

(a)

(b)

Figure 4: Potential energy (b) and vibrational ground-state energy (a) along the reaction path of the H + NH₃ → H₂ + NH₂ reaction computed using the EC-2009 (solid lines) and MC (dashed lines) surfaces. The zero of energy is set to the equilibrium potential energy of the reactants.

Figure 5: Arrhenius plots of (cm³ molecule⁻¹ s⁻¹) for the forward thermal rate coefficients of the H + NH₃ → H₂ + NH₂ reaction against the reciprocal of temperature (K) in the range 200–2000 K. Solid black line: analytical EC-2009; dashed blue line: interpolated MC; dotted red line: experimental values from [77].

Figure 6: QCT reaction cross-section () versus the collision energy (kcal mol⁻¹) for the H + NH₃ → H₂ + NH₂ reaction computed using the analytical EC-2009 (solid line) and interpolated MC (dashed line) surfaces.

Figure 7: H₂ product angular distribution (with respect to the incident H) for the H + NH₃ → H₂ + NH₂ reaction at 25 kcal·mol⁻¹ (solid line) and 40 kcal mol⁻¹ (dashed line), computed using QCT calculations on the EC-2009 surface.

Figure 8: Ammonia inversion reaction. Contour plots of the analytical EC-2009 (upper panel) and CCSD(T)/ccpVTZ ab initio surface (lower panel).

Figure 9: Classical potential for the ammonia inversion path obtained from the EC-2009 surface. The first two pairs of eigenvalues are shown.

5. When the Problems Increase

The benchmark H + CH₄ hydrogen abstraction reaction, with five light atoms and a single heavy atom, which allows a large number of very high-level ab initio calculations to be performed, gives the impression of the process being a “piece of cake” regarding polyatomic bimolecular reactions. However, it still presents kinetics and dynamics differences depending on the potential energy used and discrepancies with the experimental measures. Then, what will the case be when more complicated systems are studied? In this section we will analyze some problems which can appear alone or in combination in the study of polyatomic systems and strongly complicate the construction of potential energy surfaces.

5.1. The Spin-Orbit Problem and Multisurface Dynamics

This is a typical problem, for instance, in reactions involving halogen atoms, , with molecular systems, R–H: The halogen atom presents two spin-orbit electronic states, ²P3/2 and ²P1/2, with a splitting of 404 cm⁻¹ (1.1 kcal mol⁻¹), 882 cm⁻¹ (2.5 kcal mol⁻¹), and 3685 cm⁻¹ (3.5 kcal mol⁻¹) for F, Cl, and Br, respectively. A priori, the smaller the separation, the greater the possibility of the reaction coming from the two states, which complicates the PES construction and the dynamics study.

This problem affects all the previously described theoretical methods to develop surfaces, because it is a problem intrinsic to the initial information required: the quantum mechanical calculations. For this problem, relativistic calculations would be needed, which would immensely increase the computational cost and would make these calculations impractical in polyatomic systems. In addition, new functional forms for the analytical functions need to be developed to include the coupling between states and its dependence on coordinates [82, 83], in order to make it possible to include nonadiabatic effects and hopping between surfaces in the dynamics study of these multistate systems.

In the case of atom-diatom reactions, some results have shed light on the spin-orbit problem, although some theory/experiment controversies still persist. For instance, for the well-studied F(²P3/2, ²P1/2) + H₂ reaction, Alexander et al. [82, 84] found that the reactivity of the excited s-o state of F is small, 10–25% of the reactivity of the ground s-o state, and concluded that the overall dynamics of the F + H₂ reaction could be well described by calculations on a single, electronically adiabatic PES, although for a direct comparison with experiment the coupling between the ground and excited s-o surfaces must be considered.

For the analogue Cl(²P3/2, ²P1/2) + H₂ reaction, because of the larger energy separation, one would expect that the reaction could evolve on the ground-state adiabatic surface, with the contribution of the chlorine excited state, ²P1/2, being practically negligible according to the Born-Oppenheimer (BO) approach. However, recently a theory/experiment controversy has arisen on this issue. Thus, Lee and Liu [85–87] demonstrated experimentally the contrary for the Cl + H₂ reaction; that is, the excited chlorine atom (Cl*) is more reactive to H₂ than the ground-state chlorine (Cl) by a factor of at least ≈6. In any case, even taking into account experimental error bars, the authors were confident about the reactivity relationship Cl* > Cl, which could have a significant effect on the temperature dependence of the thermal rate coefficients. The authors interpreted this surprising result by postulating a nonadiabatic transition (breakdown of the BO approximation) from the excited Cl* to the ground-state Cl surface by either electrostatic or spin-orbit coupling in the entrance channel. This experimental study initiated a major theoretical and experimental debate on the influence of the excited Cl* in the reactivity of the Cl + H₂ reaction. Different laboratories [88–93] performed theoretical/experimental studies of the validity of the BO approximation in this reaction, concluding that the adiabatically allowed reaction (Cl(²P3/2) + H₂) will dominate the adiabatically forbidden reaction (Cl(²P1/2) + H₂). Hence, these results are in direct contrast with the experiment of Liu et al. [85–87], suggesting that this experiment claiming high reactivity of Cl* needs to be reexamined.

In the case of polyatomic systems, this level of sophistication has not been achieved and would still be computationally prohibitive. Thus, some approaches have been considered to take into account, indirectly, the s-o effect in the construction of the PES and the dynamics study. First, for thermochemical or rate coefficient calculations, the s-o effect on the multiple electronic states is taken into account by the electronic partition function of the reactants in the usual expression where is the s-o splitting of the halogen atom. Second, there is an additional effect on the barrier height (Figure 10). In fact, if we assume that along the entire reaction path the states are fully quenched, considering the s-o effect would lower the energy of the s-o ground state of the halogen atom by 1/3  below its nonrelativistic energy, increasing the barrier height by this amount. This represents 0.38, 0.83, and 1.17 kcal mol⁻¹, respectively, for F, Cl, and Br. Our group considered these approaches in the construction of the surface for the F(²P) + CH₄ [42], Cl(²P) + CH₄ [49], and Cl(²P) + NH₃ [94, 95] reactions. Note, however, that this approach is to some extent incoherent in the sense that as we move away from the saddle point towards reactants the energy asymptotically tends to zero, the energy of the lower state of the halogen atom. However, as we move from the saddle point towards reactants, it should tend asymptotically to until it reaches the point where the two states interact. From that point towards the reactants, the surface has to tend to zero. There is therefore a gap of up to between the exact result and our approximation results. However, we can assume that this gap is sufficiently small and is located so far from the dynamically important regions of the surface that we can safely neglect it. Moreover, in our analytical surfaces, the changes in the energy can be fitted in order to change the slope of the reaction path so that this asymptotic behaviour can be corrected. In addition, in the latter case of the Cl(²P) + NH₃ reaction, because of the presence of wells in the reactant channel (see below), this gap occurs before the system reaches the well in the region connecting the well with reactants, and its influence on the kinetics and dynamics is entirely negligible.

Figure 10: Schematic representation of the potential energy along the reaction path for a reaction with spin/orbit effects on the reactants. Red line: nonrelativistic calculations.

5.2. The Molecular “Size” Problem

Obviously, the cost of calculating the quantum chemical information needed to build the PES increases exponentially with the number of electrons involved, and it is still a prohibitive task for large molecules and heavy atoms. Thus, for instance, while the H + CH₄ benchmark reaction involves five light hydrogen atoms with only eleven electrons, when third-row atoms are considered, for instance, H + SiH₄, nineteen electrons are involved, or when larger systems are considered, for instance, H + CCl₄, four heavy chlorine atoms and 75 electrons must be included in the calculations. This represents an enormous computational effort, and that is prohibitive if high-level ab initio calculations are used to obtain chemical accuracy. Obviously, fitting or interpolation approaches, based on grids of 20000–30000 data points or direct dynamics calculations are still unaffordable, although more economical alternatives could be used to calculate the quantum mechanical data, such as the “dual level” technique, where the geometries and vibrational frequencies are calculated at a lower ab initio or DFT level and the energies are calculated as single points on these geometries at a higher quantum mechanical level. Even so, the computational effort would be enormous.

In these complicated cases with a large number of electrons, our strategy to build the PES, based on a smaller number of ab initio calculations, represents an interesting and practical alternative. Thus, our laboratory has constructed surfaces for several five-body systems, H + NH₃ [39, 44], 11 electrons, F + NH₃ [45], 19 electrons, and Cl(²P) + NH₃ [94, 95], 27 electrons; six-body systems, H + SiH₄ [46], 19 electrons, H + GeH₄ [96], 37 electrons, Cl(²P) + CH₄ [49], 27 electrons, Br(²P) + CH₄ [51], 45 electrons, H + CCl₄ [52, 97], 75 electrons; one of seven bodies, OH + CH₄ [53], 19 electrons. In addition, we have also studied reactions with asymmetrically substituted methane, H + CH₃Cl and Cl + CHClF₂ [98, 99], which represent another challenge in the PES construction for polyatomic systems, because in addition to the aforementioned problems, the possibility of several reaction channels needs to be taken into account. Although the older surfaces were semiempirical, that is, they combined theoretical and experimental information in the fitting procedure due to computational limitations at that time, the newer surfaces, those developed from 2007 onward, are based exclusively on quantum mechanical information: H + CH₄ [41], H + CCl₄ [97], H + NH₃ [39], and Cl(²P) + NH₃ [95].

It is noteworthy that the functional form has remained almost unchanged, being that of the H + CH₄ reaction, adding different modifications following the requirements of the systems under study. Thus, for the five-body systems we removed the dependency on one of the hydrogen atoms of CH₄, and for the seven-body system, OH + CH₄, we added additional Morse and harmonic terms to describe the OH bond and <HOH angle, while for the asymmetric reactions we allowed the four atoms bonded to the carbon atom to vary independently. In addition, we improved the analytical form when building the Cl + NH₃ surface by allowing the equilibrium N–H distance to vary along the reaction path [95].

5.3. Reactions with More Complicated Topology

The last problem analyzed in this section is that associated with the presence of several maxima and minima in the polyatomic reaction. As was recently noted by Clary [100], “reactions with several maxima and minima in the potential energy surface present the severest challenge to calculating and fitting potential energy surfaces and carrying out quantum dynamics calculations.”

Basically, taking into account the topology of the potential energy surface and independently of whether the reaction is exothermic, endothermic, or thermoneutral, the bimolecular reactions can be classified into three broad categories (Figure 11): reactions with a single barrier (Figure 11(a)); barrierless reactions (Figure 11(b)), and reactions with more complicated potentials (e.g., Figures 11(c) and 11(d)). We focus on this last case with polyatomic systems (). Obviously, the existence of wells in the entry and exit channels is associated with the presence of heavy atoms, especially halogen atoms, which favour the formation of intermediate complexes, either van der Waals or hydrogen-bonded complexes. Thus, bimolecular reactions with a complicated topology, as shown in Figures 11(c) or 11(d), also tend to pose several of the problems considered above, that is, the increase of the number of electrons involved and the spin-orbit problem. In addition, the presence of several maxima and minima requires a fine evaluation of the energies, gradients, and Hessians needed for a correct description of the PES. This, obviously, represents a complication in the construction of the PES, and only very few global surfaces have as yet been developed [26, 95].

(a)

(b)

(c)

(d)

Figure 11: Schematic representation of the potential energy along the reaction path for different types of bimolecular reactions. R, P and C denote, respectively, reactants, products, and intermediate complexes. In all the cases, the vertical axis represents the energy and the horizontal axis the reaction coordinate.

We finish this section by describing the most complicated surface analyzed by our group [95]: the analytical potential energy surface for the Cl(²P) + NH₃ polyatomic reaction, which presents several wells in the entry and exit channels, with a topology showed in Figure 11(d). Different intermediate complexes were found in the entry and exit channels at the CCSD(T)/cc-pVTZ level, and the intrinsic reaction path was calculated (energies, gradients, and Hessians) starting from the saddle point. With respect to the barrier height, one of the most difficult energy properties of a PES to estimate, different laboratories report electronic structure calculations using different levels (correlation energy and basis sets) [94, 95, 101–103]. Gao et al. [102] constructed the reaction path using the MPWB1K density functional (DFT) method [104], finding a barrier of 5.2 kcal mol⁻¹, and, using different correlation energy levels and basis sets, despite, they reported values in the range 4.8–6.2 kcal mol⁻¹. Xu and Lin [103] performed a computational study of the mechanisms and kinetics of the reaction. The geometries of the stationary points were optimized using the B3LYP DFT method, and their energies were refined with the modified Gaussian-2 (G2M) theory. They obtained a value of 7.2 kcal mol⁻¹, in contrast with the preceding values. Finally, we [94] obtained a barrier height of 7.6 kcal mol⁻¹ at the CCSD(T)/cc-pVTZ level and of 5.8 kcal mol⁻¹ [95] at a higher level, CCSD(T) = FULL/aug-cc-pVTZ, close to Gao et al.’s result. These results illustrate the dramatic influence of the electronic correlation and basis set on the correct description of the barrier, and it must be borne in mind that small differences in the saddle point produce large deviations in the kinetics and dynamics analysis. To further complicate the study of this reaction, the spin-orbit coupling must be considered, since the chlorine atom has two low-lying fine structure electronic states, ²P1/2 and ²P3/2, with a separation of  cm⁻¹  ≈2.5 kcal mol⁻¹. As discussed above, the spin-orbit coupling was taken into account in our nonrelativistic calculations in two ways: first, in the electronic partition function of the reactant (29), and second, by adding one-third of the split between the two states, that is, 0.8 kcal mol⁻¹, to the barrier height, where we have assumed that the s-o coupling is essentially fully quenched at the saddle point. With this correction, our best estimate of the barrier height was 6.6 kcal mol⁻¹.

With all this information, an analytical PES was constructed and kinetics information was obtained using VTST/MT methods [95]. The results agree well with experimental values of the rate coefficients and equilibrium constants, showing that the wells have little influence on the kinetics of the reaction. This is mainly due to the fact that VTST/MT for this reaction, whose reaction path is lower in energy than the products, assumes that tunneling has no effect. Therefore, kinetics is controlled by the properties of a transition state, which is located near the saddle point. In addition, calculations of the recrossing of the transition state that the presence of wells could cause showed that it is very small. Therefore, only the saddle point region determines the reaction probabilities.

An exhaustive dynamics study using QCT and QM methods is currently in progress in our group and should be published soon. From the results it seems that the reaction cross-sections show significant values at very low collision energies using both methodologies. In order to analyze to what extent the presence of wells (especially the well on the reactant side) is responsible for this behaviour, further studies will be carried out with a model surface similar to the Cl + NH₃ from which the wells have been removed. Note that the latter study is possible because of the availability of an analytical surface that can be refitted to remove the wells without significantly modifying other regions of the PES, which is an added value to this kind of analytical PES. These kinds of study can help us to understand the dynamics of such a complicated system and whether transition-state theory, which (as noted above) assumes that only the saddle point region is significant, needs to be challenged.

6. Final Remarks

The construction of potential energy surfaces in polyatomic reactive systems represents a major theoretical challenge, with a very high computational and human time cost. Based on high-level electronic structure calculations, fitting, interpolation, and analytical (defined by functional forms) approaches have been developed and applied in the kinetics and dynamics (classical, quasi-classical, and quantum mechanical) study of these reactions. However, in spite of the enormous progress in the last 20–30 years in theoretical algorithms and computational power, the construction of potential energy surfaces in polyatomic systems has still not reached the level of accuracy achieved for the triatomic systems.

The quality of these surfaces is still an open and debatable question, and even for the benchmark H + CH₄ polyatomic reaction, which involves only five light hydrogen atoms and eleven electrons, small differences are found depending on the PES construction and the dynamics method. Unfortunately, the problems will increase for other important chemical systems, where some effects are present such as spin-orbit coupling, increase in the number of electrons and molecular size, or potentials with more complicated topology. In these cases, the kinetics and dynamics results will be even more strongly dependent on the quality of the PES. In the last few decades, there has been much theoretical effort on the part of various laboratories, but there is still much to do in this research field.

Throughout this paper, the emphasis has been put on the strategy developed by our group which has constructed about 15 surfaces for polyatomic systems of five, six and seven bodies. These are freely available for download from the POTLIB websites, http://comp.chem.umn.edu/POTLIB/ or http://users.ipfw.edu/DUCHOVIC/POTLIB2001/, or can be requested from the authors.

The analytical surfaces developed in our group are based on a basic functional form with slight modifications to suit it to each of the systems studied. In this sense, they are very specialized force fields, with mathematical functions that allow further improvements to be introduced. For example, adding terms to describe anharmonicity, mode-mode coupling, or an improved dependence of the constants of the fit on the reaction coordinate are possibilities, which can improve the functional form.

The inclusion of additional reaction channels, such as the exchange reactions in methane or ammonia, is another pending matter in work on our surfaces. This would open up the possibility of analyzing competitive channels and could lead to a better understanding of the behaviour of complex reactions.

An additional advantage of this approach is the negligible computational cost of the evaluation of the potential energy surface and its derivatives, which is a very desirable feature when one wants to apply expensive QM methods for the study of these reactions. Undoubtedly, with the evolution of computer resources, eventually a point will be reached when direct dynamics calculations using high ab initio levels will be feasible. However, until then, the possibility of having an analytical function that can describe the potential energy with reasonable accuracy is the fastest option. In recent years ad hoc surfaces have been constructed which only reproduce a single kinetics or dynamics property but do not provide a complete description of the reaction system. However, the aim of our work is to build a PES that reproduces at least qualitatively all the kinetics and dynamics features of the reactive system. Obviously a quantitative description of all kinetics and dynamics properties will be the final goal, and with this target in mind in our group we have adopted a flexible functional form that can be easily improved and adapted to other similar systems.

Acknowledgments

This work has been partially supported in recent years by the Junta de Extremadura, Spain, and Fondo Social Europeo (Projects nos. 2PR04A001, PRI07A009, and IB10001). One of the authors (M.Monge-Palacios) thanks Junta de Extremadura (Spain) for a scholarship.

References

M. Born and R. Oppenheimer, “Zur quantentheorie der molekeln,” Annals of Physics, vol. 84, pp. 457–484, 1927.

C. N. Ramachandran, D. De Fazio, S. Cavalli, F. Tarantelli, and V. Aquilanti, “Revisiting the potential energy surface for the He + → HeH⁺ + H reaction at the full configuration interaction level,” Chemical Physics Letters, vol. 469, no. 1–3, pp. 26–30, 2009.

F. George and S. Kumar, “Ab initio ground and the first excited adiabatic and quasidiabatic potential energy surfaces of H⁺ + CO system,” Chemical Physics, vol. 373, no. 3, pp. 211–218, 2010.

T. Tanaka and T. Takayanagi, “Quantum reactive scattering calculations of H + F₂ and Mu + F₂ reactions on a new ab initio potential energy surface,” Chemical Physics Letters, vol. 496, no. 4–6, pp. 248–253, 2010.

W. Skomorowski, F. Pawlowski, T. Korona, R. Moszynski, P. S. Zuchowski, and J. M. Hutson, “Interaction between LiH molecule and Li atom from the state-of-the-art electronic structure calculations,” Journal of Chemical Physics, vol. 134, p. 114109, 2011.

B. Jiang, Ch. Xie, and D. Xie, “New ab iniito potential energy surface for BrH₂ and rate constants for the H + HBr → H₂ + Br abstraction reaction,” Journal of Chemical Physics, vol. 134, p. 114301, 2011.

M. Karplus, R. N. Porter, and R. D. Sharma, “Exchange reactions with activation energy. I. Simple barrier potential for (H, H₂),” Journal of Chemical Physics, vol. 43, no. 9, pp. 3259–3287, 1965.

D.J. Kouri, Y. Sun, et al., “New time-dependent and time-independent computational methods for molecular collisions,” in Mathematical Frontiers in Computational Chemical Physics, D. G. Truhlar, Ed., pp. 207–244, Springer, New York, NY, USA, 1988.

A. Fernandez-Ramos, J. A. Miller, S. J. Klippenstein, and D. G. Truhlar, “Modeling the kinetics of bimolecular reactions,” Chemical Reviews, vol. 106, pp. 4518–4584, 2006.

T. V. Albu, J. Espinosa-García, and D. G. Truhlar, “Computational chemistry of polyatomic reaction kinetics and dynamics: the quest for an accurate CH₅ potential energy surface,” Chemical Reviews, vol. 107, no. 11, pp. 5101–5132, 2007.

W. L. Hase, K. Song, and M. S. Gordon, “Direct dynamics simulations,” Computing in Science and Engineering, vol. 5, no. 4, pp. 36–44, 2003.

J. Liu, K. Song, W. L. Hase, and S. L. Anderson, “Direct dynamics trajectory study of the reaction of formaldehyde cation with D₂: vibrational and zero-point energy effects on quasiclassical trajectories,” Journal of Physical Chemistry A, vol. 109, no. 50, pp. 11376–11384, 2005.

T.-S. Ho and H. Rabitz, “A general method for constructing multidimensional molecular potential energy surfaces from ab initio calculations,” Journal of Chemical Physics, vol. 104, no. 7, pp. 2584–2597, 1996.

M. A. Collins, “Molecular potential-energy surfaces for chemical reaction dynamics,” Theoretical Chemistry Accounts, vol. 108, no. 6, pp. 313–324, 2002.

T.-S. Ho and H. Rabitz, “Reproducing kernel Hilbert space interpolation methods as a paradigm of high dimensional model representations: application to multidimensional potential energy surface construction,” Journal of Chemical Physics, vol. 119, no. 13, pp. 6433–6442, 2003.

J. Ischtwan and M. A. Collins, “Molecular potential energy surfaces by interpolation,” Journal of Chemical Physics, vol. 100, no. 11, pp. 8080–8088, 1994.

M. A. Addicoat and M. A. Collins, “Potential energy surfaces: the forces of chemistry,” in Tutorials in Molecular Reaction Dynamics, M. Brouard and C. Vallance, Eds., p. 28, RSC Publishing, Cambridge, UK, 2010.

Y. Zhou, B. Fu, C. Wang, M. A. Collins, and D. H. Zhang, “Ab initio potential energy surface and quantum dynamics for the H + CH₄ → H₂ + CH₃,” Journal of Chemical Physics, vol. 134, pp. 64323–64330, 2011.

J. N. Murrell, S. Carter, S. C. Farantos, P. Huxley, and A. J. C. Varandas, Molecular Potential Energy Functions, Wiley, Chistester, UK, 1984.

G. C. Schatz, “The analytical representation of electronic potential-energy surfaces,” Reviews of Modern Physics, vol. 61, no. 3, pp. 669–688, 1989.

J. M. Bowman and G. C. Schatz, “Theoretical studies of polyatomic bimolecular reaction dynamics,” Annual Review of Physical Chemistry, vol. 46, no. 1, pp. 169–195, 1995.

K. Sorbie and J. N. Murrell, “Theoretical study of the O(¹D)+H₂(¹) reactive quenching process,” Molecular Physics, vol. 31, pp. 905–920, 1976.

A. J. C. Varandas, “Four-atom bimolecular reactions with relevance in environmental chemistry,” International Reviews in Physical Chemistry, vol. 19, no. 2, pp. 199–245, 2000.

A. Brown, B. J. Braams, K. Christoffel, Z. Jin, and J. M. Bowman, “Classical and quasiclassical spectral analysis of using an ab initio potential energy surface,” Journal of Chemical Physics, vol. 119, no. 17, pp. 8790–8793, 2003.

Z. Jin, B. J. Braams, and J. M. Bowman, “An ab initio based global potential energy surface describing → + H₅,” Journal of Physical Chemistry A, vol. 110, no. 4, pp. 1569–1574, 2006.

B. J. Braams and J. M. Bowman, “Permutationally invariant potential energy surfaces in high dimensionality,” International Reviews in Physical Chemistry, vol. 28, no. 4, pp. 577–606, 2009.

J. M. Bowman, G. Czako, and B. Fu, “High-dimensional ab initio potential energy surfaces for reaction dynamics calculations,” Physical Chemistry Chemical Physics, vol. 13, no. 18, pp. 8094–8111, 2011.

A. B. McCoy, B. J. Braams, A. Brown, X. Huang, Z. Jin, and J. M. Bowman, “Ab initio diffusion monte carlo calculations of the quantum behavior of in full dimensionality,” Journal of Physical Chemistry A, vol. 108, no. 23, pp. 4991–4994, 2004.

T. Valencich and D. L. Bunker, “Energy-dependent cross sections for the tritium-methane hot atom reactions,” Chemical Physics Letters, vol. 20, no. 1, pp. 50–52, 1973.

L. M. Raff, “Theoretical investigations of the reaction dynamics of polyatomic systems: chemistry of the hot atom (T* + CH₄) and (T* + CD₄) systems,” Journal of Chemical Physics, pp. 2220–2244, 1970.

R. Steckler, K. J. Dykema, F. B. Brown, G. C. Hancock, D. G. Truhlar, and T. Valencich, “A comparative study of potential energy surfaces for CH₃ + H₂↔ CH₄ + H,” Journal of Chemical Physics, vol. 87, no. 12, pp. 7024–7035, 1987.

T. Joseph, R. Steckler, and D. G. Truhlar, “A new potential energy surface for the CH₃ + H₂↔ CH₄ + H reaction: calibration and calculations of rate constants and kinetic isotope effects by variational transition state theory and semiclassical tunneling calculations,” Journal of Chemical Physics, vol. 87, no. 12, pp. 7036–7049, 1987.

M. J. T. Jordan and R. G. Gilbert, “Classical trajectory studies of the reaction CH₄ +H → CH₃ + H₂,” Journal of Chemical Physics, vol. 102, no. 14, pp. 5669–5682, 1995.

J. Espinosa-García and J. C. Corchado, “Global surfaces and reaction-path potentials. Applications of variational transition-state theory,” Recent Research Development in Physical Chemistry, vol. 1, pp. 165–192, 1997.

J. Espinosa-Garcia, “Analysis of economical strategies for the construction of potential energy surfaces,” Trends in Physical Chemistry, vol. 8, pp. 49–73, 2000.

A. Chakraborty, Y. Zhao, H. Lin, and D. G. Truhlar, “Combined valence bond-molecular mechanics potential-energy surface and direct dynamics study of rate constants and kinetic isotope effects for the H + C₂H₆ reaction,” Journal of Chemical Physics, vol. 124, no. 4, pp. 44315–44328, 2006.

R. J. Duchovic, W. L. Hase, and H. B. Schlegel, “Analytic function for the H + CH₃⇄ CH₄ potential energy surface,” Journal of Physical Chemistry, vol. 88, no. 7, pp. 1339–1347, 1984.

M. Yang and J. C. Corchado, “Seven-dimensional quantum dynamics study of the H + NH₃ → H₂ + NH₂ reaction,” Journal of Chemical Physics, vol. 126, no. 21, pp. 214312–214321, 2007.

J. Espinosa-Garcia and J. C. Corchado, “Analytical potential energy surface and kinetics of the NH₃ + H → NH2 + H₂ hydrogen abstraction and the ammonia inversion reactions,” Journal of Physical Chemistry A, vol. 114, no. 12, pp. 4455–4463, 2010.

S. T. Banks and D. C. Clary, “Reduced dimensionality quantum dynamics of Cl + CH₄ → HCl + CH₃ on an ab initio potential,” Physical Chemistry Chemical Physics, vol. 9, no. 8, pp. 933–943, 2007.

J. C. Corchado, J. L. Bravo, and J. Espinosa-Garcia, “The hydrogen abstraction reaction H + CH₄. I. New analytical potential energy surface based on fitting to ab initio calculations,” Journal of Chemical Physics, vol. 130, no. 18, pp. 184314–184323, 2009.

C. Rángel, M. Navarrete, and J. Espinosa-García, “Potential energy surface for the F(, ) + CH₄ hydrogen abstraction reaction. Kinetics and dynamics study,” Journal of Physical Chemistry A, vol. 109, no. 7, pp. 1441–1448, 2005.

J. Espinosa-García, “New analytical potential energy surface for the CH₄ + H hydrogen abstraction reaction: thermal rate constants and kinetic isotope effects,” Journal of Chemical Physics, vol. 116, no. 24, pp. 10664–10673, 2002.

J. C. Corchado and J. Espinosa-García, “Analytical potential energy surface for the NH₃ + H ↔ NH₂ + H₂ reaction: application of variational transition-state theory and analysis of the equilibrium constants and kinetic isotope effects using curvilinear and rectilinear coordinates,” Journal of Chemical Physics, vol. 106, no. 10, pp. 4013–4021, 1997.

J. Espinosa-García and J. C. Corchado, “Analytical surface for the reaction with no saddle-point NH₃ + F → NH₂ + FH. Application of variational transition state theory,” Journal of Physical Chemistry A, vol. 101, no. 40, pp. 7336–7344, 1997.

J. Espinosa-García, J. Sansón, and J. C. Corchado, “The SiH₄ + H → SiH₃ + H₂ reaction: potential energy surface, rate constants, and kinetic isotope effects,” Journal of Chemical Physics, vol. 109, no. 2, pp. 466–473, 1998.

J. Espinosa-García and J. C. García-Bernáldez, “Analytical potential energy surface for the CH₄ + O(³P) → CH₃ + OH reaction. Thermal rate constants and kinetic isotope effects,” Physical Chemistry Chemical Physics, vol. 2, no. 10, pp. 2345–2351, 2000.

J. Espinosa-García, J. L. Bravo, and C. Rangel, “New analytical potential energy surface for the F(²P) + CH₄ hydrogen abstraction reaction: kinetics and dynamics,” Journal of Physical Chemistry A, vol. 111, no. 14, pp. 2761–2771, 2007.

J. C. Corchado, D. G. Truhlar, and J. Espinosa-García, “Potential energy surface, thermal, and state-selected rate coefficients, and kinetic isotope effects for Cl + CH₄ → HCl + CH₃,” Journal of Chemical Physics, vol. 112, no. 21, pp. 9375–9389, 2000.

C. Rangel, M. Navarrete, J. C. Corchado, and J. Espinosa-García, “Potential energy surface, kinetics, and dynamics study of the Cl + CH₄ → HCl + CH₃ reaction,” Journal of Chemical Physics, vol. 124, no. 12, pp. 124306–124324, 2006.

J. Espinosa-García, “Potential energy surface for the CH₃ + HBr → CH₄ + Br hydrogen abstraction reaction: thermal and state-selected rate constants, and kinetic isotope effects,” Journal of Chemical Physics, vol. 117, no. 5, pp. 2076–2086, 2002.

C. Rangel and J. Espinosa-García, “Potential energy surface for the CCl₄ + H → CCl₃ + ClH reaction: kinetics and dynamics study,” Journal of Chemical Physics, vol. 122, no. 13, pp. 134315–134319, 2005.

J. Espinosa-García and J. C. Corchado, “Potential energy surface for a seven-atom reaction. Thermal rate constants and kinetic isotope effects for CH₄ + OH,” Journal of Chemical Physics, vol. 112, no. 13, pp. 5731–5739, 2000.

R.J. Duchovic, Y.L. Volobuev, G.C. Lynch, et al., “POLTLIB,” 2001.

X. Zhang, B. J. Braams, and J. M. Bowman, “An ab initio potential surface describing abstraction and exchange for H + CH₄,” Journal of Chemical Physics, vol. 124, no. 2, pp. 021104–021108, 2006.

Z. Xie, J. M. Bowman, and X. Zhang, “Quasiclassical trajectory study of the reaction H + CH₄ ( =0,1) → CH₃ + H₂ using a new ab initio potential energy surface,” Journal of Chemical Physics, vol. 125, no. 13, pp. 133120–133127, 2006.

T. Wu, H. J. Werner, and U. Manthe, “Accurate potential energy surface and quantum reaction rate calculations for the H + CH₄ → H₂ + CH₃ reaction,” Journal of Chemical Physics, vol. 124, no. 16, pp. 164307–164318, 2006.

F. Huarte-Larrañaga and U. Manthe, “Quantum dynamics of the CH₄ + H → CH₃ + H₂ reaction: full-dimensional and reduced dimensionality rate constant calculations,” Journal of Physical Chemistry A, vol. 105, no. 12, pp. 2522–2529, 2001.

G. Schiffel, U. Manthe, and G. Nyman, “Full-dimensional quantum reaction rate calculations for H + CH₄ → H₂ + CH₃ on a recent potential energy surface,” Journal of Physical Chemistry A, vol. 114, no. 36, pp. 9617–9622, 2010.

J. W. Sutherland, M. C. Su, and J. V. Michael, “Rate constants for H + CH₄, CH₃ + H₂, and CH₄ dissociation at high temperature,” International Journal of Chemical Kinetics, vol. 33, no. 11, pp. 669–684, 2001.

G. J. Germann, Y. D. Huh, and J. J. Valentini, “State-to-state dynamics of atom + polyatom abstraction reactions. II. The H + CD₄ → HD(v',J') + CD₃ reaction,” Journal of Chemical Physics, vol. 96, no. 8, pp. 1957–1966, 1992.

W. Hu, G. Lendvay, D. Troya et al., “H + CD4 abstraction reaction dynamics: product energy partitioning,” Journal of Physical Chemistry A, vol. 110, no. 9, pp. 3017–3027, 2006.

J. P. Camden, H. A. Bechtel, and R. N. Zare, “Dynamics of the simplest reaction of a carbon atom in a tetrahedral environment,” Angewandte Chemie, vol. 42, no. 42, pp. 5227–5230, 2003.

J. P. Camden, H. A. Bechtel, D. J. A. Brown et al., “A reinterpretation of the mechanism of the simplest reaction at an sp 3-hybridized carbon atom: H + CD₄ → CD₃ + HD,” Journal of the American Chemical Society, vol. 127, no. 34, pp. 11898–11899, 2005.

J. P. Camden, W. Hu, H. A. Bechtel et al., “H + CD₄ abstraction reaction dynamics: excitation function and angular distributions,” Journal of Physical Chemistry A, vol. 110, no. 2, pp. 677–686, 2006.

W. Zhang, Y. Zhou, G. Wu et al., “Depression of reactivity by the collision energy in the single barrier H + CD₄ → HD + CD₃ reaction,” Proceedings of the National Academy of Sciences of the United States of America, vol. 107, no. 29, pp. 12782–12785, 2010.

C. Rangel, J. Sansón, J. C. Corchado, J. Espinosa-Garcia, and G. Nyman, “Product angular distribution for the H + CD₄ → HD + CD₃ reaction,” Journal of Physical Chemistry A, vol. 110, no. 37, pp. 10715–10719, 2006.

J. Espinosa-García, G. Nyman, and J. C. Corchado, “The hydrogen abstraction reaction H+ CH₄. II. Theoretical investigation of the kinetics and dynamics,” Journal of Chemical Physics, vol. 130, no. 18, pp. 184315–184323, 2009.

J. P. Camden, H. A. Bechtel, D. J. A. Brown, and R. N. Zare, “Effects of C-H stretch excitation on the H + CH₄ reaction,” Journal of Chemical Physics, vol. 123, no. 13, pp. 134301–134309, 2005.

C. Rangel, J. C. Corchado, and J. Espinosa-García, “Quasi-classical trajectory calculations analyzing the reactivity and dynamics of asymmetric stretch mode excitations of methane in the H + CH₄ reaction,” Journal of Physical Chemistry A, vol. 110, no. 35, pp. 10375–10383, 2006.

C. H. Townes and A. L. Schawlow, Microwave Spectroscopy, McGraw-Hill, New York, NY, USA, 1955.

G. E. Moyano and M. A. Collins, “Interpolated potential energy surface for abstraction and exchange reactions of NH₃ + H and deuterated analogues,” Theoretical Chemistry Accounts, vol. 113, no. 4, pp. 225–232, 2005.

J. Espinosa-García, “Capability of LEPS surfaces to describe the kinetics and dynamics of non-collinear reactions,” Journal of Physical Chemistry A, vol. 105, no. 1, pp. 134–139, 2001.

J. Espinosa-García, “Capability of LEP-type surfaces to describe noncollinear reactions. 2. Polyatomic systems,” Journal of Physical Chemistry A, vol. 105, no. 38, pp. 8748–8755, 2001.

B. C. Garrett and D. G. Truhlar, “A least-action variational method for calculating multidimensional tunneling probabilities for chemical reactions,” Journal of Chemical Physics, vol. 79, no. 10, pp. 4931–4938, 1983.

R. Meana-Pañeda, D. G. Truhlar, and A. Fernández-Ramos, “Least-action tunnelling transmission coefficient for polyatomic reactions,” Journal of Chemical Theory and Computation, vol. 6, no. 1, pp. 6–17, 2010.

T. Ko, P. Marshall, and A. Fontijn, “Rate coefficients for the H + NH₃ reaction over a wide temperature range,” Journal of Physical Chemistry, vol. 94, no. 4, pp. 1401–1404, 1990.

J. Espinosa-García and J. C. Corchado, “Quasi-classical trajectory calculations of the hydrogen abstraction reaction H + NH₃,” Journal of Physical Chemistry A, vol. 114, no. 21, pp. 6194–6200, 2010.

A. M. Hapern, B. R. Ramachandran, and E. D. Glendening, “The inversion potential of ammonia: an intrinsic reaction coordinate calculation for student investigation,” Journal of Chemical Education, vol. 84, no. 6, pp. 1067–1072, 2007.

V. Špirko, “Vibrational anharmonicity and the inversion potential function of NH₃,” Journal of Molecular Spectroscopy, vol. 101, no. 1, pp. 30–47, 1983.

V. M. Devi, P. P. Das, K. N. Rao, S. Urban, D. Papouseck, and V. Spirko, “The diode laser spectrum of the band of ¹⁴ND₃ and ¹⁵ND₃,” Journal of Molecular Spectroscopy, vol. 88, pp. 293–299, 1981.

M. H. Alexander, D. E. Manolopoulos, and H. J. Werner, “An investigation of the F + H₂ reaction based on a full ab initio description of the open-shell character of the F(²P) atom,” Journal of Chemical Physics, vol. 113, no. 24, pp. 11084–11100, 2000.

G. Capecchi and H. J. Werner, “Ab initio calculations of coupled potential energy surfaces for the Cl(,) + H₂ reaction,” Physical Chemistry Chemical Physics, vol. 6, no. 21, pp. 4975–4983, 2004.

M. H. Alexander, H. J. Werner, and D. E. Manolopoulos, “Spin-orbit effects in the reaction of F(²P) with H₂,” Journal of Chemical Physics, vol. 109, no. 14, pp. 5710–5713, 1998.

S. H. Lee, “State-specific excitation function for Cl(²P) + H₂ (v = 0,j): effects of spin-orbit and rotational states,” Journal of Chemical Physics, vol. 110, no. 17, pp. 8229–8232, 1999.

S.-H. Lee and K. Liu, “Exploring the spin-orbit reactivity in the simplest chlorine atom reaction,” Journal of Chemical Physics, vol. 111, no. 14, pp. 6253–6259, 1999.

F. Dong, S. H. Lee, and K. Liu, “Direct determination of the spin-orbit reactivity in CI(, ) + H₂/D₂/HD reactions,” Journal of Chemical Physics, vol. 115, no. 3, pp. 1197–1204, 2001.

M. H. Alexander, G. Capecchi, and H. J. Werner, “Theoretical study of the validity of the Born-Oppenheimer approximation in the CL + H₂ → HCL + H reaction,” Science, vol. 296, no. 5568, pp. 715–718, 2002.

N. Balucani, D. Skouteris, L. Cartechini et al., “Differential cross sections from quantum calculations on coupled ab initio potential energy surfaces and scattering experiments for Cl(²P) + H₂ reactions,” Physical Review Letters, vol. 91, no. 1, pp. 013201–013204, 2003.

M. H. Alexander, G. Capecchi, and H.-J. Werner, “Details and consequences of the nonadiabatic coupling in the Cl(²P) + H² reaction,” Faraday Discussions, vol. 127, pp. 59–72, 2004.

E. Garand, J. Zhou, D. E. Manolopoulos, M. H. Alexander, and D. M. Neumark, “Nonadiabatic interactions in the Cl + H₂ reaction probed by and photoelectron imaging,” Science, vol. 319, no. 5859, pp. 72–75, 2008.

D. M. Neumark, “Slow electron velocity-map imaging of negative ions: applications to spectroscopy and dynamics,” Journal of Physical Chemistry A, vol. 112, no. 51, pp. 13287–13301, 2008.

X. Wang, W. Dong, C. Xiao et al., “The extent of non-born-oppenheimer coupling in the reaction of Cl(²P) with para-H₂,” Science, vol. 322, no. 5901, pp. 573–576, 2008.

M. Monge-Palacios and J. Espinosa-Garcia, “Reaction-path dynamics calculations of the Cl + NH₃ hydrogen abstraction reaction: the role of the intermediate complexes,” Journal of Physical Chemistry A, vol. 114, no. 12, pp. 4418–4426, 2010.

M. Monge-Palacios, C. Rangel, J. C. Corchado, and J. Espinosa-Garcia, “Analytical potential energy surface for the reaction with intermediate complexes NH₃ + Cl → NH₂ + HCl: application to the kinetics study,” International Journal of Quantum Chemistry. In press

J. Espinosa-García, “Analytical potential energy surface for the GeH₄ + H → GeH₃ + H₂ reaction: thermal and vibrational-state selected rate constants and kinetic isotope effects,” Journal of Chemical Physics, vol. 111, no. 20, pp. 9330–9336, 1999.

J. Espinosa-García, C. Rangel, M. Monge-Palacios, and J. C. Corchado, “Kinetics and dynamics study of the H + CCl₄ → HCl(v', j')+ CCl₃ reaction,” Theoretical Chemistry Accounts, vol. 128, pp. 743–755, 2010.

C. Rangel and J. Espinosa-García, “Analytical potential energy surface describing abstraction reactions in asymmetrically substituted polyatomic systems of type CX₃Y + A → products,” Journal of Physical Chemistry A, vol. 110, no. 2, pp. 537–547, 2006.

C. Rangel and J. Espinosa-García, “Potential energy surface for asymmetrically substituted reactions of type CWXYZ + A. kinetics study,” Journal of Physical Chemistry A, vol. 111, no. 23, pp. 5057–5062, 2007.

D. C. Clary, “Theoretical studies on bimolecular reaction dynamics,” Proceedings of the National Academy of Sciences of the United States of America, vol. 105, no. 35, pp. 12649–12653, 2008.

S. Kondo, K. Tokuhashi, A. Takahashi, and M. Kaise, “Ab initio study of reactions between halogen atoms and various fuel molecules by Gaussian-2 theory,” Journal of Hazardous Materials, vol. 79, no. 1-2, pp. 77–86, 2000.

Y. Gao, I. M. Alecu, P.-C. Hsieh, B. P. Morgan, P. Marshall, and L. N. Krasnoperov, “Thermochemistry is not a lower bound to the activation energy of endothermic reactions: a kinetic study of the gas-phase reaction of atomi-chlorine with ammonia,” Journal of Physical Chemistry A, vol. 110, no. 21, pp. 6844–6850, 2006.

Z. F. Xu and M. C. Lin, “Computational studies on the kinetics and mechanisms for NH₃ reactions with ClOx (x = 0-4) radicals,” Journal of Physical Chemistry A, vol. 111, no. 4, pp. 584–590, 2007.

Y. Zhao and D. G. Truhlar, “Hybrid meta density functional theory methods for thermochemistry, thermochemical kinetics, and noncovalent interactions: the MPW1B95 and MPWB1K models and comparative assessments for hydrogen bonding and van der Waals interactions,” Journal of Physical Chemistry A, vol. 108, no. 33, pp. 6908–6918, 2004.