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

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 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 + 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 𝑇_𝑖(𝑅):𝑉=^𝑁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) 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 𝑉=𝑝(𝑥)+_𝑖<𝑗𝑞_𝑖,𝑗(𝑥)⋅𝑦_𝑖,𝑗,(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 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, 𝑉_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=⁴_𝑖=1𝑉₃𝑅_CHi,𝑅_CHB,𝑅_HiHB,(6) where 𝑅 is the distance between the two subscript atoms, H_i stands for one of the four methane hydrogens, and H_B 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:^{𝑉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 CH_i bonds, 𝛼_CH, however, is not taken as a constant but rather as a function of the CH distances,𝛼_CH=𝑎_CH+𝑏_CH⎛⎜⎜⎝tanh𝑐_CH𝑅−𝑅^𝑒_CH+12⎞⎟⎟⎠,(8) with 𝑅 being the average 𝑅_CHi distance, 𝑅=14⁴_𝑖=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 + 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:𝑅^𝑜_CH=𝑃₁𝑅^𝑜_CH,𝑅+1−𝑃₁𝑅^𝑜_CH,𝑃,(10) where 𝑃₁ is𝑃₁=⁴_𝑖=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 𝑇₁ is a geometry-dependent switching function, given by𝑇₁𝑅_CHi=1−tanh𝑤₁𝑅_CHi−𝑤₂,(12) where 𝑤₁ and 𝑤₂ 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=12³_4𝑖=1_𝑗=𝑖+1𝑘^0𝑖𝑗𝑘_𝑖𝑘_𝑗𝜃_𝑖𝑗−𝜃^0𝑖𝑗²,(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𝑆¹𝑅_CHi𝑆₁𝑅^CH_j−1_CH4−𝑘^CH3𝑆₂𝑅_CHi𝑆₂𝑅^CH_j−1,(14) while 𝑘_𝑖 is a function of both the 𝑅_CHi and 𝑅_HiHB distances:𝑘_𝑖𝐴=₁𝐴₁exp−𝐴₂𝑅_CHi−𝑅^𝑒_CH²,𝐴₂=1−exp−𝑎𝑎₁𝑅_HiHB²,=𝑎𝑎₂+𝑎𝑎₃exp−𝑎𝑎₄𝑅^H_i^H_B−𝑅^𝑒_HiHB².(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𝑖𝑗=acos−1+acos3−1+acos−13−𝜋2𝑆^𝜑𝑅_CHi𝑆_𝜑𝑅^CH_j−13−23𝜋𝑆^𝜃𝑅_CHk𝑆_𝜃𝑅_CHl−1.(16) Finally, the switching functions are given by𝑆𝑆¹₂𝑅_CHi=1−tanh𝛼^𝑠1𝑅_CHi−𝑅^𝑒_CH𝑅_CHi−𝛽^𝑠1⁸,𝑆_𝜑𝑅_CHi=1−tanh𝛼^𝑠2𝑅_CHi−𝑅^𝑒_CH𝑅_CHi−𝛽^𝑠2⁶,𝑆𝑅_CHi=1−tanh𝐴_𝜑𝑅_CHi−𝑅^𝑒_CHexp𝐵_𝜑𝑅_CHi−𝐶_𝜑³_𝜃,𝑅_CHi=1−tanh𝐴_𝜃𝑅_CHi−𝑅^𝑒_CHexp𝐵_𝜃𝑅_CHi−𝐶_𝜃³,(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=⁴_𝑖=1𝑓_Δ𝑖⁴_{𝑗=1𝑗≠𝑖}Δ_𝑖𝑗²+⁴_𝑖=1ℎ_Δ𝑖⁴_{𝑗=1𝑗≠𝑖}Δ_𝑖𝑗⁴.(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−𝑆₃𝑅_CHi𝑓_Δ,𝑆₃=1−𝑆₃𝑅_CHiℎ_Δ,𝑅_CHi=1−tanh𝛼^𝑠3𝑅_CHi−𝑅^𝑒_CH𝑅_CHi−𝛽^𝑠3²,(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 CH_i 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,𝜃^𝑜𝑖𝑗=𝜏+𝜏−𝜋+𝜏−2𝜋2𝑆^𝜑𝑅_CX𝑖⋅𝑆_𝜑𝑅^CX_𝑗−13𝑆^𝜗𝑅_CX𝑘⋅𝑆_𝜗𝑅_CX𝑙−1,(21) where 𝜏=109.47^∘, is replaced by𝜃^𝑜𝑖𝑗=+𝜏+𝜏−𝜏𝜏−₂𝜏₁𝑆_𝜑𝑅_CX𝑖⋅𝑆_𝜑𝑅^CX_𝑗−1𝑆_𝜗𝑅_CX𝑘⋅𝑆_𝜗𝑅_CX𝑙−1,(22) where 𝜏₂ is the bending angle in the non-planar product, and 𝜏₁ is related to 𝜏₂ by the expression, 𝜏₁=𝜋−arcsinsin𝜏²/2sin(𝜋/3).(23) 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 𝐶_3𝑣, characterized by an inversion mechanism through a planar structure with symmetry 𝐷_3ℎ. 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 (𝐷_3h symmetry), which is a saddle point to the ammonia inversion, is about 9 kcal mol⁻¹ more stable than the pyramidal structure (𝐶_3v symmetry).

In the original expression for the H + CH₄ reaction (5), the 𝑉_op 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 𝑉_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𝑅=_𝑥||𝐸(𝑥)−𝐹(𝑥,𝑝)||²,(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𝑅=_𝑥𝑤^𝑒𝑥||𝐸(𝑥)−𝐹(𝑥,𝑝)||²+_𝑥𝑤^𝑔𝑥|||𝑔(𝑥)−𝜕𝐹(𝑥,𝑝)𝜕𝑥|||²+_𝑥𝑤^𝐻𝑥||||𝐻(𝑥)−𝜕²𝐹(𝑥,𝑝)𝜕𝑥²||||²,(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 + 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₄ Reaction
The 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 𝑣^=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 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

Table 2 

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

Figure 1 

Arrhenius plots of ln𝑘 (cm3 molecule−1 s−1) for the forward thermal rate coefficients of the H + CH4 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 

CD3 product angular distribution (with respect to the incident CD4) for the H + CD4 → HD + CD3 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 (𝑎2𝑜) versus collision energy (eV) for the H + CH4 → H2 + CH3 reaction on the CBE (black line), ZBB3 (blue line), and ZFWCZ (red line) surfaces.

(b) The H + NH₃ 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 + 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 (𝐶_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⁻¹, 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 (𝐷_3ℎ symmetry) is about 9 kcal mol⁻¹ more stable than the pyramidal structure (𝐶_3𝑣 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 (𝑉_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⁻¹ 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, 𝑘_tunn, which is in turn obtained from the imaginary action integral, 𝜃(𝐸), using the WKB approximation [79]: Δ𝐸₀=𝑘^tunn𝑘2𝑐_tunn,=2𝑐𝜈²𝜋exp[−𝜃(𝐸)],(27) 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 kcal⋅mol⁻¹ 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 Δ𝐸=0.37 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

Table 4 

Eigenvalues of ammonia inversion (in cm⁻¹).

(a)

(b)

Figure 4 

Potential energy (b) and vibrational ground-state energy (a) along the reaction path of the H + NH3 → H2 + NH2 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 ln𝑘 (cm3 molecule−1 s−1) for the forward thermal rate coefficients of the H + NH3 → H2 + NH2 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 (𝑎20) versus the collision energy (kcal mol−1) for the H + NH3 → H2 + NH2 reaction computed using the analytical EC-2009 (solid line) and interpolated MC (dashed line) surfaces.

Figure 7 

H2 product angular distribution (with respect to the incident H) for the H + NH3 → H2 + NH2 reaction at 25 kcal·mol−1 (solid line) and 40 kcal mol−1 (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, X(²P), with molecular systems, R–H: X²P+H-R⟶XH+R(X=F,Cl,Br).(28) The halogen atom presents two spin-orbit electronic states, ²P_3/2 and ²P_1/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(²P_3/2, ²P_1/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(²P_3/2, ²P_1/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, ²P_1/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(²P_3/2) + H₂) will dominate the adiabatically forbidden reaction (Cl(²P_1/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𝑄𝑒=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⁻¹, 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 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(²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 (𝑛>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].

(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, ²P_1/2 and ²P_3/2, with a separation of 𝜀=882 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.