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 + CH4 and the H + NH3 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 BrH2 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 3𝑛−6 internal degrees of freedom, if one needs 𝑚 configurations for each degree of freedom, one would need a total of 𝑚(3𝑛−6) evaluations of the energy. For instance, taking a typical value of 𝑚 of 10, for the H + CH4 hydrogen abstraction reaction 1012 evaluations would be needed. This implies too much computer time to perform ab initio calculations with chemical accuracy (1 kcal mol−1) and becomes totally unaffordable if one aims to achieve spectroscopic accuracy (1 cm−1).

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 + CH4 benchmark polyatomic system and the H + NH3 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−1 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 𝑇𝑖(𝑅):𝑉=𝑁data𝑖=1𝑤𝑖(𝑍)𝑇𝑖(𝑅),(1) where 𝑁data 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 𝑖:𝑇𝑖(𝑅)=𝑉(𝑅)+3𝑛−6𝑘=1𝜕𝑉𝜕𝑅𝑘𝑅𝑘−𝑅𝑘+1(𝑖)2!3𝑛−6𝑘=13𝑛−6𝑗=1𝑅𝑘−𝑅𝑘𝜕(𝑖)2𝑉𝜕𝑅𝑘𝜕𝑅𝑗𝑅𝑗−𝑅𝑗.(𝑖)(2) 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 + CH4 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 𝑉=𝑝(𝑥)+𝑖<𝑗𝑞𝑖,𝑗(𝑥)⋅𝑦𝑖,𝑗,(3) where 𝑥 is an 𝑛-dimensional vector, depending on 𝑛 internuclear distances, with components 𝑥𝑖,𝑗=exp(−𝑟𝑖𝑗), 𝑟𝑖𝑗 being the distance between nuclei 𝑖 and 𝑗, and 𝑦𝑖,𝑗 given by𝑦𝑖,𝑗=𝑒−𝑟𝑖,𝑗𝑟𝑖,𝑗,(4) and 𝑝(𝑥) and 𝑞(𝑥) are polynomials constructed to satisfy the permutation symmetry with respect to the indistinguishable nuclei. For instance, for the CH5   + 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−1. 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 + CH4 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, 𝑉stretch, harmonic bending term, 𝑉harm, and anharmonic out-of-plane potential, 𝑉op, 𝑉=𝑉stretch+𝑉harm+𝑉op.(5)

The stretching potential is the sum of four London-Eyring-Polanyi (LEP) terms, each one corresponding to a permutation of the four methane hydrogens:𝑉stretch=4𝑖=1𝑉3𝑅CHi,𝑅CHB,𝑅HiHB,(6) 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 𝑉3 LEP potential is well known, we will recall it for the sake of completeness:𝑉3𝑅CHi,𝑅CHB,𝑅HiHB𝑅=𝑄CHi𝑅+𝑄CHB𝑅+𝑄HiHB−12𝐽𝑅CHi𝑅−𝐽HiHB2+12𝐽𝑅HiHB𝑅−𝐽CHB2+12𝐽𝑅CHB𝑅−𝐽CHi2,𝑄𝑅𝑋𝑌=𝐸1𝑅𝑋𝑌+𝐸3𝑅𝑋𝑌2,𝐽𝑅𝑋𝑌=𝐸1𝑅𝑋𝑌−𝐸3𝑅𝑋𝑌2,𝐸1𝑅𝑋𝑌=𝐷1𝑋𝑌exp−2𝛼𝑋𝑌𝑅𝑋𝑌−𝑅𝑒𝑋𝑌−2exp−𝛼𝑋𝑌𝑅𝑋𝑌−𝑅𝑒𝑋𝑌,𝐸3𝑅𝑋𝑌=𝐷3𝑋𝑌exp−2𝛼𝑋𝑌𝑅𝑋𝑌−𝑅𝑒𝑋𝑌+2exp−𝛼𝑋𝑌𝑅𝑋𝑌−𝑅𝑒𝑋𝑌,(7)where there are 12 fitting parameters, four for each of the three kinds of bond, 𝑅CHi, 𝑅CHB, and 𝑅HiHB. In particular, these are the singlet and triplet dissociation energies, 𝐷1𝑋𝑌 and 𝐷3𝑋𝑌, the equilibrium bond distance, 𝑅𝑒𝑋𝑌, and the Morse parameter, 𝛼𝑋𝑌. The Morse parameter for the CHi bonds, 𝛼CH, however, is not taken as a constant but rather as a function of the CH distances,𝛼CH=𝑎CH+𝑏CH⎛⎜⎜⎝𝑐tanhCH𝑅−𝑅𝑒CH+12⎞⎟⎟⎠,(8) with 𝑅 being the average 𝑅CHi distance, 1𝑅=44𝑖=1𝑅CHi.(9) In this way, 𝛼CH changes smoothly from its value at methane, 𝑎CH+(𝑏CH/2), to its value at the methyl radical, 𝑎CH+𝑏CH, 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 + C2H6 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:𝑅𝑜CH=𝑃1𝑅𝑜CH,𝑅+1−𝑃1𝑅𝑜CH,𝑃,(10) where 𝑃1 is𝑃1=4𝑖=1𝑇1𝑅CHi,(11) which is symmetric with respect to all the four hydrogen atoms and goes to zero as one of the hydrogen atoms is abstracted, and 𝑇1 is a geometry-dependent switching function, given by𝑇1𝑅CHi𝑤=1−tanh1𝑅CHi−𝑤2,(12) where 𝑤1 and 𝑤2 are adjustable parameters. Therefore, this adds 2 new parameters (total 16 parameters) to describe the stretching potential.

The 𝑉harm term is the sum of six harmonic terms, one for each bond angle in methane:𝑉harm=1234𝑖=1𝑗=𝑖+1𝑘0𝑖𝑗𝑘𝑖𝑘𝑗𝜃𝑖𝑗−𝜃0𝑖𝑗2,(13) where 𝑘0𝑖𝑗 and 𝑘𝑖 are force constants and 𝜃0𝑖𝑗 are the reference angles. The 𝑘0𝑖𝑗 force constants are allowed to evolve from their value in methane, 𝑘CH4, to their value in methyl, 𝑘CH3, which are two parameters of the fit, by means of switching functions: 𝑘0𝑖𝑗=𝑘CH4+𝑘CH4𝑆1𝑅CHi𝑆1𝑅CHj+𝑘−1CH4−𝑘CH3𝑆2𝑅CHi𝑆2𝑅CHj,−1(14) while 𝑘𝑖 is a function of both the 𝑅CHi and 𝑅HiHB distances:𝑘𝑖=𝐴1exp−𝐴2𝑅CHi−𝑅𝑒CH2,𝐴1=1−exp−𝑎𝑎1𝑅HiHB2,𝐴2=𝑎𝑎2+𝑎𝑎3exp−𝑎𝑎4𝑅HiHB−𝑅𝑒HiHB2.(15) Thus, four adjustable parameters are involved in the definition of 𝑘𝑖.

The reference angles are also allowed to change from their value at methane, 𝜏=109.47∘, or arccosine(−1/3), to methyl, 120° or 2𝜋/3 radians, by means of switching functions [37]: 𝜃0𝑖𝑗−1=acos3+−1acos3−𝜋2𝑆𝜑𝑅CHi𝑆𝜑𝑅CHj+−1−1acos3−2𝜋3𝑆𝜃𝑅CHk𝑆𝜃𝑅CHl.−1(16) Finally, the switching functions are given by𝑆1𝑅CHi𝛼=1−tanh𝑠1𝑅CHi−𝑅𝑒CH𝑅CHi−𝛽𝑠18,𝑆2𝑅CHi𝛼=1−tanh𝑠2𝑅CHi−𝑅𝑒CH𝑅CHi−𝛽𝑠26,𝑆𝜑𝑅CHi𝐴=1−tanh𝜑𝑅CHi−𝑅𝑒CH𝐵exp𝜑𝑅CHi−𝐶𝜑3,𝑆𝜃𝑅CHi𝐴=1−tanh𝜃𝑅CHi−𝑅𝑒CH𝐵exp𝜃𝑅CHi−𝐶𝜃3,(17) involving 10 more adjustable parameters. In total, 16 terms need to be fitted for the calibration of the 𝑉harm potential.

The 𝑉op potential is a quadratic-quartic term whose aim is to correctly describe the out-of-plane motion of methyl:𝑉op=4𝑖=1𝑓Δ𝑖4𝑗=1𝑗≠𝑖Δ𝑖𝑗2+4𝑖=1ℎΔ𝑖4𝑗=1𝑗≠𝑖Δ𝑖𝑗4.(18) 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:𝑓Δ𝑖=1−𝑆3𝑅CHi𝑓Δ,ℎΔ𝑖=1−𝑆3𝑅CHiℎΔ,𝑆3𝑅CHi𝛼=1−tanh𝑠3𝑅CHi−𝑅𝑒CH𝑅CHi−𝛽𝑠32,(19) with 𝑓Δ, ℎΔ, 𝛼𝑠3, and 𝛽𝑠3 being the only parameters of 𝑉op that enter the fitting process. Δ𝑖𝑗 is the angle that measures the deviation from the reference angle:Δ𝑖𝑗=acos⃗𝑟𝑘−⃗𝑟𝑗×⃗𝑟𝑙−⃗𝑟𝑗‖‖⃗𝑟𝑘−⃗𝑟𝑗×⃗𝑟𝑙−⃗𝑟𝑗‖‖⃗𝑟𝑖‖‖⃗𝑟𝑖‖‖−𝜃0𝑖𝑗,(20) where ⃗𝑟𝑖, ⃗𝑟𝑗, ⃗𝑟𝑘, and ⃗𝑟𝑙 are vectors going from the carbon atom to the 𝑖, 𝑗, 𝑘, and l hydrogen atoms, respectively, and 𝜃0𝑖𝑗 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 + CH4 reaction the CH3 radical product is planar. What happens if the product presents a nonplanar geometry? In these cases of CX3 products, we have modified the reference angle 𝜃𝑜𝑖𝑗 in (13) and (20). Thus, the original expression,𝜃𝑜𝑖𝑗𝜋=𝜏+𝜏−2𝑆𝜑𝑅CX𝑖⋅𝑆𝜑𝑅CX𝑗+−1𝜏−2𝜋3𝑆𝜗𝑅CX𝑘⋅𝑆𝜗𝑅CX𝑙,−1(21) where 𝜏=109.47∘, is replaced by𝜃𝑜𝑖𝑗=𝜏+𝜏−𝜏1𝑆𝜑𝑅CX𝑖⋅𝑆𝜑𝑅CX𝑗+−1𝜏−𝜏2𝑆𝜗𝑅CX𝑘⋅𝑆𝜗𝑅CX𝑙,−1(22) where 𝜏2 is the bending angle in the non-planar product, and 𝜏1 is related to 𝜏2 by the expression, 𝜏1𝜏=𝜋−arcsinsin2/2sin(𝜋/3).(23) In these cases of non-planar products, CX3, 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 + CH4 reaction the CH4 reactant presents 𝑇𝑑 symmetry. What happens if the reactant presents a different symmetry? For instance, ammonia presents symmetry 𝐶3𝑣, characterized by an inversion mechanism through a planar structure with symmetry 𝐷3ℎ. The functional form (5) used for the H + CH4 reaction cannot be applied without modification to any kind of system. Thus, when this potential is applied to the study of the H + NH3 → H2 + NH2 reaction, a major drawback was observed [38], namely that it wrongly describes the NH3 inversion reaction, predicting that the planar ammonia (𝐷3h symmetry), which is a saddle point to the ammonia inversion, is about 9 kcal mol−1 more stable than the pyramidal structure (𝐶3v symmetry).

In the original expression for the H + CH4 reaction (5), the 𝑉op term was added to obtain a correct description of the out-of-plane bending in the methyl radical, 580 cm−1. 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 𝑉op term is removed: 𝑉=𝑉stretch+𝑉harm.(24) 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𝑅=𝑥||||𝐸(𝑥)−𝐹(𝑥,𝑝)2,(25) 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𝑅=𝑥𝑤𝑒𝑥||||𝐸(𝑥)−𝐹(𝑥,𝑝)2+𝑥𝑤𝑔𝑥|||𝑔(𝑥)−𝜕𝐹(𝑥,𝑝)|||𝜕𝑥2+𝑥𝑤𝐻𝑥||||𝐻𝜕(𝑥)−2𝐹(𝑥,𝑝)𝜕𝑥2||||2,(26) 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 𝑚3𝑛−6 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 + CH4 reaction [42] using as the starting point the analytical PES for H + CH4 [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 + CH4 Reaction
The gas-phase H + CH4 → H2 + CH3 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−1. The three surfaces present a colinear saddle point, 180.0°, with a similar barrier height, with small differences, ±0.2 kcal mol−1, 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−1 [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 + CD4 → HD + CD3 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 𝑣=0 and 𝑣=1 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 CD3 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 + CD4 gas-phase reaction using the Photoloc technique. They found that the CD3 products are sideways/forward scattered with respect to the incident CD4, suggesting a stripping mechanism (note that in the original papers [63–65] the CD3 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 + CD4 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−1 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 + CH4 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).

(b) The H + NH3 Reaction
The 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 + NH3 → H2 + NH2 hydrogen abstraction reaction is similar to the H + CH4 reaction. It presents a barrier height of 14.5 kcal mol−1 and a reaction exoergicity of 5.0 kcal mol−1, as compared to 14.87 and an endoergicity of 2.88 kcal mol−1, respectively, for the H + CH4 analogue. Also, the inversion of ammonia between two pyramidal structures (𝐶3𝑣 symmetry) passing through a planar structure (𝐷3ℎ 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−1, depending on the level of calculation.
Only three potential energy surfaces have been developed for the H + NH3 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 + NH3 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 NH3 inversion motion, predicting incorrectly that the planar ammonia (𝐷3ℎ symmetry) is about 9 kcal mol−1 more stable than the pyramidal structure (𝐶3𝑣 symmetry). To correct this behaviour of the NH3 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−1 from the ab initio value.
For the two surfaces, Figure 4 shows the energy along the minimum energy path (𝑉MEP) 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 𝑉MEP 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 𝑠=−1 a.u. is due to wild oscillations of the frequencies that lead to unphysical values of 𝑉𝐺𝑎 for |𝑠|>1. 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−1 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−1 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 H2 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−1. At 25 kcal mol−1, 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−1, 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, 𝑘tunn, which is in turn obtained from the imaginary action integral, 𝜃(𝐸), using the WKB approximation [79]: Δ𝐸0=𝑘tunn,𝑘2𝑐tunn=2𝑐𝜈2𝜋[],exp−𝜃(𝐸)(27) with 𝑐 being the speed of light and 𝑣2 the eigenvalue associated with the inversion mode of ammonia, 1113 cm−1. 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 kcal⋅mol−1 give rise to a factor of four in the computed splitting. With respect to the effect of isotopic substitution, for the ND3 case we obtain a value Δ𝐸=0.37 cm−1. Although this is greater than the experimentally reported value, 0.05 cm−1 [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.

5. When the Problems Increase

The benchmark H + CH4 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, X(2P), with molecular systems, R–H: X2P+H-R⟶XH+R(X=F,Cl,Br).(28) The halogen atom presents two spin-orbit electronic states, 2P3/2 and 2P1/2, with a splitting of 404 cm−1 (1.1 kcal mol−1), 882 cm−1 (2.5 kcal mol−1), and 3685 cm−1 (3.5 kcal mol−1) 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(2P3/2, 2P1/2) + H2 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 + H2 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(2P3/2, 2P1/2) + H2 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, 2P1/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 + H2 reaction; that is, the excited chlorine atom (Cl*) is more reactive to H2 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 + H2 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(2P3/2) + H2) will dominate the adiabatically forbidden reaction (Cl(2P1/2) + H2). 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−𝜀𝑄𝑒=4+2exp𝑘𝐵𝑇,(29) 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−1, respectively, for F, Cl, and Br. Our group considered these approaches in the construction of the surface for the F(2P) + CH4 [42], Cl(2P) + CH4 [49], and Cl(2P) + NH3 [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 1/3𝜀 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 1/3𝜀 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(2P) + NH3 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.

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 + CH4 benchmark reaction involves five light hydrogen atoms with only eleven electrons, when third-row atoms are considered, for instance, H + SiH4, nineteen electrons are involved, or when larger systems are considered, for instance, H + CCl4, 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 + NH3 [39, 44], 11 electrons, F + NH3 [45], 19 electrons, and Cl(2P) + NH3 [94, 95], 27 electrons; six-body systems, H + SiH4 [46], 19 electrons, H + GeH4 [96], 37 electrons, Cl(2P) + CH4 [49], 27 electrons, Br(2P) + CH4 [51], 45 electrons, H + CCl4 [52, 97], 75 electrons; one of seven bodies, OH + CH4 [53], 19 electrons. In addition, we have also studied reactions with asymmetrically substituted methane, H + CH3Cl and Cl + CHClF2 [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 + CH4 [41], H + CCl4 [97], H + NH3 [39], and Cl(2P) + NH3 [95].

It is noteworthy that the functional form has remained almost unchanged, being that of the H + CH4 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 CH4, and for the seven-body system, OH + CH4, 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 + NH3 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 (𝑛>4). 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].

We finish this section by describing the most complicated surface analyzed by our group [95]: the analytical potential energy surface for the Cl(2P) + NH3 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−1, and, using different correlation energy levels and basis sets, despite, they reported values in the range 4.8–6.2 kcal mol−1. 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−1, in contrast with the preceding values. Finally, we [94] obtained a barrier height of 7.6 kcal mol−1 at the CCSD(T)/cc-pVTZ level and of 5.8 kcal mol−1 [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, 2P1/2 and 2P3/2, with a separation of 𝜀=882 cm−1  ≈2.5 kcal mol−1. 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−1, 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−1.

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 + NH3 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 + CH4 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.


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.