#### 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_{4} and the H + NH_{3} 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_{2} system, *ab initio* calculations were performed at more than 12000 geometries. Then, a three-dimensional cubic spline interpolation was employed to yield the ground-state PES for this system. As can be seen, the computational effort (correlation energy + basis set + number of points) is enormous, and doing this for larger systems is still clearly prohibitive.

The extension to larger reactive systems (more than 4 atoms) is an open and promising field, although computationally very costly. Thus, for polyatomic systems of atoms, with internal degrees of freedom, if one needs configurations for each degree of freedom, one would need a total of evaluations of the energy. For instance, taking a typical value of of 10, for the H + CH_{4} hydrogen abstraction reaction 10^{12} 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 + CH_{4} benchmark polyatomic system and the H + NH_{3} 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 :
where is the number of molecular configurations whose energy and its first (gradient) and second (Hessian) derivatives have been evaluated, is the normalized weighting factor, is a set of internal coordinates, and the Taylor series are expanded around the point :
The use of the energy and its first and second derivatives makes the method very robust, but computationally very expensive, especially if very high-level *ab initio* calculations are needed for a correct description of the reactive system. For instance, very recently a new PES for the H + CH_{4} gas-phase reaction has been developed by Collins by interpolation of 30000 data points obtained with high-level *ab initio* calculations [18]. The *ab initio* calculations took about 600 days on a workstation with 8βCPU cores in total, while the evaluations of the potential values on all the grids for this reaction took 400 days. Three advantages of this method are the following: first, that the interpolated function agrees precisely with the values of the data points; second, that new electronic structure points may be incorporated in the process to improve the PES; third, that the PES is invariant to the permutation or exchange of indistinguishable nuclei. However, Shepard interpolation often yields slight oscillations in equipotential contours and vibrational frequencies (second derivatives of the energy); that is, the potential energy surface is not smooth enough to provide smooth changes in the vibrational frequencies. This is a limitation when one wants to compute tunneling effects at low energies, since the barrier to tunnel through is unrealistic.

Fitting methods, such as the interpolated ones considered above, have been widely used in the construction of potential energy surfaces [19β21], and they have been recently reviewed [10]. Thus, for instance, a fitting procedure widely used is based on the many-body expansion (MBE) method [22], amplified by Varandas [23] in the so-called double many-body expansion (DMBE) method. The potential is given by a sum of terms corresponding to atoms, diatoms, triatoms, and tetra-atoms and a series of adjustable coefficients. Another fitting procedure has been developed and applied by Bowman et al. [24β27], in which the *ab initio* data are globally fitted to a permutational symmetry invariant polynomial. The function describing the PES has the form
where is an -dimensional vector, depending on internuclear distances, with components , being the distance between nuclei and , and given by
and and are polynomials constructed to satisfy the permutation symmetry with respect to the indistinguishable nuclei. For instance, for the CH_{5}^{βββ+} 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 + CH_{4} hydrogen abstraction reaction, as a paradigm of polyatomic systems [29β33]. The potential energy for a given geometry, , is given by the sum of three terms: stretching potential, , harmonic bending term, , and anharmonic out-of-plane potential, ,

The stretching potential is the sum of four London-Eyring-Polanyi (LEP) terms, each one corresponding to a permutation of the four methane hydrogens:
where is the distance between the two subscript atoms, 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:where there are 12 fitting parameters, four for each of the three kinds of bond, , , and . In particular, these are the singlet and triplet dissociation energies, and , the equilibrium bond distance, , and the Morse parameter, . The Morse parameter for the CH_{i} bonds, , however, is not taken as a constant but rather as a function of the distances,
with being the average distance,
In this way, changes smoothly from its value at methane, , to its value at the methyl radical, , as the reaction evolves. Therefore, 14 parameters are required to describe the stretching potential.

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

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

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

The potential is a quadratic-quartic term whose aim is to correctly describe the out-of-plane motion of methyl:
The force constants, and , have been incorporated into a new switching function which is such that vanishes at the methane limit and which directs the change of the methyl fragment from pyramidal to planar as the reaction evolves:
with , , , and being the only parameters of that enter the fitting process. is the angle that measures the deviation from the reference angle:
where , , , and are vectors going from the carbon atom to the , , , and l hydrogen atoms, respectively, and are the reference angles defined in (16). The first term to the right of (20) is therefore the angle between the 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_{4} reaction the CH_{3} radical product is planar. What happens if the product presents a nonplanar geometry? In these cases of CX_{3} products, we have modified the reference angle in (13) and (20). Thus, the original expression,
where , is replaced by
where is the bending angle in the non-planar product, and is related to by the expression,
In these cases of non-planar products, CX_{3}, 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_{4} reaction the CH_{4} reactant presents symmetry. What happens if the reactant presents a different symmetry? For instance, ammonia presents symmetry , characterized by an inversion mechanism through a planar structure with symmetry . The functional form (5) used for the H + CH_{4} reaction cannot be applied without modification to any kind of system. Thus, when this potential is applied to the study of the H + NH_{3} β H_{2} + NH_{2} reaction, a major drawback was observed [38], namely that it wrongly describes the NH_{3} inversion reaction, predicting that the planar ammonia ( symmetry), which is a saddle point to the ammonia inversion, is about 9βkcalβmol^{β1} more stable than the pyramidal structure ( symmetry).

In the original expression for the H + CH_{4} reaction (5), the 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 term is removed:
Consequently, each of the 16 parameters is required for the stretching terms and the harmonic bending terms of the PES which were fitted to high-level *ab initio* CCSD(T)/cc-pVTZ calculations.

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

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

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

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

This strategy presents certain advantages over other methods. The first is transferability of the functional form and the fitted parameters. Since one can regard VB/MM as some kind of highly specific MM force field, it can be expected that some of the fitting parameters are transferable to a similar system, although obviously the parameter values are system specific. In our group we have used this feature to obtain, for example, the PES for the F + CH_{4} reaction [42] using as the starting point the analytical PES for H + CH_{4} [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 _{4} Reaction*

The gas-phase H + CH

_{4}β H

_{2}+ CH

_{3}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 + CD

_{4}β HD + CD

_{3}reaction because there is more experimental information available for comparison. The product energy partitioning has been experimentally obtained by Valentiniβs group. [61] only for the vibration and rotation of the HD product, and this same group found that more than 95% of the HD product is formed in the and vibrational states. These results appear in Table 2, together with the QCT results on the analytical CBE and the fitted ZBB1 [55] surfaces. The QCT calculations on these two PESs show excellent agreement for the HD product energy partitioning as well as for the HD vibrational distribution. In this latter case, the QCT results reproduce the experimental evidence, but, in the first case, they strongly contrast with the experimental measurements, which measure an internal excitation of the HD product of 7% and 9% for vibration and rotation, respectively. The agreement between the results from the two surfaces leads us to think that the discrepancies with experiment are mainly due to the dynamical method, that is, to limitations of the QCT approach. Experimental problems, however, cannot be totally ruled out, as was recently observed by Hu et al. [62], who suggested that the conclusions from Valentini et al.βs CARS experimental study might need to be reinterpreted.

The product angular distribution is, doubtless, one of the most sensitive dynamics features with which to test the quality of the potential energy surface, but experimentally it is very difficult to measure in some cases. When the Photoloc technique is used, the laboratory speed depends on both the scattering angle and the speed of the CD

_{3}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

_{4}gas-phase reaction using the Photoloc technique. They found that the CD

_{3}products are sideways/forward scattered with respect to the incident CD

_{4}, suggesting a stripping mechanism (note that in the original papers [63β65] the CD

_{3}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

_{4}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 + CH

_{4}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).

**(a)**

**(b)**

**(c)**

**(d)**

*(b) The H + NH _{3} 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

_{3}β H

_{2}+ NH

_{2}hydrogen abstraction reaction is similar to the H + CH

_{4}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 + CH

_{4}analogue. Also, the inversion of ammonia between two pyramidal structures ( symmetry) passing through a planar structure ( symmetry, which is a saddle point) is an example of intramolecular dynamics [71] with a barrier height in the range 5.20β5.94βkcalβmol

^{β1}, depending on the level of calculation.

Only three potential energy surfaces have been developed for the H + NH

_{3}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

_{3}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

_{3}inversion motion, predicting incorrectly that the planar ammonia ( symmetry) is about 9βkcalβmol

^{β1}more stable than the pyramidal structure ( symmetry). To correct this behaviour of the NH

_{3}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 () and the ground-state vibrationally adiabatic potential curve (), which is defined as the sum of the potential and vibrational zero-point energies along the reaction path. Note that is the reaction coordinate, being zero at the saddle point, positive in the product channel, and negative in the reactant channel. Firstly, the curve is smooth and shows no oscillations in either surface. However, shows significant oscillations when computed for the interpolated surface. The reason is that the potential energy surface is not smooth enough to provide smooth changes in the frequencies (second derivatives of the energy). The situation becomes worse as one moves away from the saddle point and approaches areas where little

*ab initio*information for interpolation is available. Thus, the bump in at around βa.u. is due to wild oscillations of the frequencies that lead to unphysical values of for . This is a limitation when one wants to compute the tunneling effect at low energies (see below), since the barrier to tunnel through is totally unrealistic.

Figure 5 plots the thermal rate coefficients computed with the VTST/MT approach on both surfaces, where tunneling was estimated using the least-action tunneling (LAT) method [75, 76], together with experimental values [77] for comparison. In the common temperature range, 490β1780βK, both surfaces reproduce the experimental data, which is a test of both the surface and the dynamical method. The differences between the rates for the two surfaces are relatively small, 40% at 600βK, and diminishing as temperature increases. At low temperatures, the differences between the two surfaces increase, with the rates computed on the MC surface being 92% larger than the ones with the EC-2009 surface at 200βK. In this low temperature regime, where tunneling is important, the analytical EC-2009 surface is more accurate due to the unrealistic oscillations in the adiabatic reaction path on the MC surface, which negatively influence the tunneling.

The integral cross-sections in the range 10β30βkcalβmol

^{β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 H

_{2}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, , which is in turn obtained from the imaginary action integral, , using the WKB approximation [79]: with being the speed of light and the eigenvalue associated with the inversion mode of ammonia, 1113βcm

^{β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 kcalmol

^{β1}give rise to a factor of four in the computed splitting. With respect to the effect of isotopic substitution, for the ND

_{3}case we obtain a value β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.

**(a)**

**(b)**

#### 5. When the Problems Increase

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

##### 5.1. The Spin-Orbit Problem and Multisurface Dynamics

This is a typical problem, for instance, in reactions involving halogen atoms, , with molecular systems, RβH:
The halogen atom presents two spin-orbit electronic states, ^{2}P_{3/2} and ^{2}P_{1/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(^{2}P_{3/2}, ^{2}P_{1/2}) + H_{2} 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_{2} 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(^{2}P_{3/2}, ^{2}P_{1/2}) + H_{2} 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, ^{2}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_{2} reaction; that is, the excited chlorine atom (Cl*) is more reactive to H_{2} 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_{2} 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(^{2}P_{3/2}) + H_{2}) will dominate the adiabatically forbidden reaction (Cl(^{2}P_{1/2}) + H_{2}). Hence, these results are in direct contrast with the experiment of Liu et al. [85β87], suggesting that this experiment claiming high reactivity of Cl* needs to be reexamined.

In the case of polyatomic systems, this level of sophistication has not been achieved and would still be computationally prohibitive. Thus, some approaches have been considered to take into account, indirectly, the s-o effect in the construction of the PES and the dynamics study. First, for thermochemical or rate coefficient calculations, the s-o effect on the multiple electronic states is taken into account by the electronic partition function of the reactants in the usual expression
where is the s-o splitting of the halogen atom. Second, there is an additional effect on the barrier height (Figure 10). In fact, if we assume that along the entire reaction path the states are fully quenched, considering the s-o effect would lower the energy of the s-o ground state of the halogen atom by 1/3β below its nonrelativistic energy, increasing the barrier height by this amount. This represents 0.38, 0.83, and 1.17βkcalβmol^{β1}, respectively, for F, Cl, and Br. Our group considered these approaches in the construction of the surface for the F(^{2}P) + CH_{4} [42], Cl(^{2}P) + CH_{4} [49], and Cl(^{2}P) + NH_{3} [94, 95] reactions. Note, however, that this approach is to some extent incoherent in the sense that as we move away from the saddle point towards reactants the energy asymptotically tends to zero, the energy of the lower state of the halogen atom. However, as we move from the saddle point towards reactants, it should tend asymptotically to until it reaches the point where the two states interact. From that point towards the reactants, the surface has to tend to zero. There is therefore a gap of up to between the exact result and our approximation results. However, we can assume that this gap is sufficiently small and is located so far from the dynamically important regions of the surface that we can safely neglect it. Moreover, in our analytical surfaces, the changes in the energy can be fitted in order to change the slope of the reaction path so that this asymptotic behaviour can be corrected. In addition, in the latter case of the Cl(^{2}P) + NH_{3} 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 + CH_{4} benchmark reaction involves five light hydrogen atoms with only eleven electrons, when third-row atoms are considered, for instance, H + SiH_{4}, nineteen electrons are involved, or when larger systems are considered, for instance, H + CCl_{4}, 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_{3} [39, 44], 11 electrons, F + NH_{3} [45], 19 electrons, and Cl(^{2}P) + NH_{3} [94, 95], 27 electrons; six-body systems, H + SiH_{4} [46], 19 electrons, H + GeH_{4} [96], 37 electrons, Cl(^{2}P) + CH_{4} [49], 27 electrons, Br(^{2}P) + CH_{4} [51], 45 electrons, H + CCl_{4} [52, 97], 75 electrons; one of seven bodies, OH + CH_{4} [53], 19 electrons. In addition, we have also studied reactions with asymmetrically substituted methane, H + CH_{3}Cl and Cl + CHClF_{2} [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_{4} [41], H + CCl_{4} [97], H + NH_{3} [39], and Cl(^{2}P) + NH_{3} [95].

It is noteworthy that the functional form has remained almost unchanged, being that of the H + CH_{4} 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_{4}, and for the seven-body system, OH + CH_{4}, 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_{3} surface by allowing the equilibrium NβH distance to vary along the reaction path [95].

##### 5.3. Reactions with More Complicated Topology

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

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

**(a)**

**(b)**

**(c)**

**(d)**

We finish this section by describing the most complicated surface analyzed by our group [95]: the analytical potential energy surface for the Cl(^{2}P) + NH_{3} 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, ^{2}P_{1/2} and ^{2}P_{3/2}, with a separation of β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 + NH_{3} 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_{4} 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.