Interpolating moving least-squares methods for fitting potential energy surfaces: computing high-density potential energy surface data from low-density ab initio data points

A highly accurate and efficient method for molecular global potential energy surface (PES) construction and fitting is demonstrated. An interpolating-moving-least-squares (IMLS)-based method is developed using low-density ab initio Hessian values to compute high-density PES parameters suitable for accurate and efficient PES representation. The method is automated and flexible so that a PES can be optimally generated for classical trajectories, spectroscopy, or other applications. Two important bottlenecks for fitting PESs are addressed. First, high accuracy is obtained using a minimal density of ab initio points, thus overcoming the bottleneck of ab initio point generation faced in applications of modified-Shepard-based methods. Second, high efficiency is also possible (suitable when a huge number of potential energy and gradient evaluations are required during a trajectory calculation). This overcomes the bottleneck in high-order IMLS-based methods, i.e., the high cost/accuracy ratio for potential energy evaluations. The result is a set of hybrid IMLS methods in which high-order IMLS is used with low-density ab initio Hessian data to compute a dense grid of points at which the energy, Hessian, or even high-order IMLS fitting parameters are stored. A series of hybrid methods is then possible as these data can be used for neural network fitting, modified-Shepard interpolation, or approximate IMLS. Results that are indicative of the accuracy, efficiency, and scalability are presented for one-dimensional model potentials as well as for three-dimensional (HCN) and six-dimensional (HOOH) molecular PESs.     

Interpolating moving least-squares methods for fitting potential energy surfaces: Analysis of an application to a six-dimensional system

The basic formal and numerical aspects of different degree interpolated moving least-squares (IMLS) methods are applied to a six-dimensional potential energy surface (PES) of the HOOH molecule, for which an analytic ("exact") potential is available in the literature. The results of systematic investigations of the effects of weight function parameters, the degree and partial degree of IMLS, the number of data points allowed, and the optimal automatic point selection of data points up to full third-degree IMLS fits are reported. With partial reduction of cross terms and automatic point selection the full six-dimensional HOOH PES can be fit over a range of 100 kcal/mol to an accuracy of less than 1 kcal/mol with approximately 1350 ab initio points.     

Interpolating moving least-squares methods for fitting potential energy surfaces: a strategy for efficient automatic data point placement in high dimensions

An accurate and efficient method for automated molecular global potential energy surface (PES) construction and fitting is demonstrated. An interpolating moving least-squares (IMLS) method is developed with the flexibility to fit various ab initio data: (1) energies, (2) energies and gradients, or (3) energies, gradients, and Hessian data. The method is automated and flexible so that a PES can be optimally generated for trajectories, spectroscopy, or other applications. High efficiency is achieved by employing local IMLS in which fitting coefficients are stored at a limited number of expansion points, thus eliminating the need to perform weighted least-squares fits each time the potential is evaluated. An automatic point selection scheme based on the difference in two successive orders of IMLS fits is used to determine where new ab initio data need to be calculated for the most efficient fitting of the PES. A simple scan of the coordinate is shown to work well to identify these maxima in one dimension, but this search strategy scales poorly with dimension. We demonstrate the efficacy of using conjugate gradient minimizations on the difference surface to locate optimal data point placement in high dimensions. Results that are indicative of the accuracy, efficiency, and scalability are presented for a one-dimensional model potential (Morse) as well as for three-dimensional (HCN), six-dimensional (HOOH), and nine-dimensional (CH4) molecular PESs.     

Interpolating moving least-squares methods for fitting potential-energy surfaces: further improvement of efficiency via cutoff strategies

In standard applications of interpolating moving least squares (IMLS) for fitting a potential-energy surface (PES), all available ab initio points are used. Because remote ab initio points negligibly influence IMLS accuracy and increase IMLS time-to-solution, we present two methods to locally restrict the number of points included in a particular fit. The fixed radius cutoff (FRC) method includes ab initio points within a hypersphere of fixed radius. The density adaptive cutoff (DAC) method includes points within a hypersphere of variable radius depending on the point density. We test these methods by fitting a six-dimensional analytical PES for hydrogen peroxide. Both methods reduce the IMLS time-to-solution by about an order of magnitude relative to that when no cutoff method is used. The DAC method is more robust and efficient than the FRC method.     

Interpolating moving least-squares methods for fitting potential energy surfaces: Improving efficiency via local approximants

The local interpolating moving least-squares (IMLS) method for constructing potential energy surfaces is investigated. The method retains the advantageous features of the IMLS approach in that the ab initio derivatives are not required and high degree polynomials can be used to provide accurate fits, while at the same time it is much more efficient than the standard IMLS approach because the least-squares solutions need to be calculated only once at the data points. Issues related to the implementation of the local IMLS method are investigated and the accuracy is assessed using HOOH as a test case. It is shown that the local IMLS method is at the same level of accuracy as the standard IMLS method. In addition, the scaling of the method is found to be a power law as a function of number of data points N, N(-q). The results suggest that when fitting only to the energy values for a d-dimensional system by using a Qth degree polynomial the power law exponent q approximately Qd when the energy range fitted is large (e.g., E<100 kcalmol for HOOH), and q>Qd when the energy range fitted is smaller (E<30 kcalmol) and the density of data points is higher. This study demonstrates that the local IMLS method provides an efficient and accurate means for constructing potential energy surfaces.     

Interpolating moving least-squares methods for fitting potential energy surfaces: an application to the H2CN unimolecular reaction

Classical trajectories have been used to compute rates for the unimolecular reaction H2CN-->H+HCN on a fitted ab initio potential energy surface (PES). The ab initio energies were obtained from CCSD(T)/aug-cc-pvtz electronic structure calculations. The ab initio energies were fitted by the interpolating moving least-squares (IMLS) method. This work continues the development of the IMLS method for producing ab initio PESs for use in molecular dynamics simulations of many-atom systems. A dual-level scheme was used in which the preliminary selection of data points was done using a low-level theory and the points used for fitting the final PES were obtained at the desired higher level of theory. Classical trajectories were used on various low-level IMLS fits to tune the fit to the unimolecular reaction under study. Procedures for efficiently picking data points, selecting basis functions, and defining cutoff limits to exclude distant points were investigated. The accuracy of the fitted PES was assessed by comparing interpolated values of quantities to the corresponding ab initio values. With as little as 330 ab initio points classical trajectory rate constants were converged to 5%-10% and the rms error over the six-dimensional region sampled by the trajectories was a few tenths of a kcal/mol.     

Interpolating moving least-squares methods for fitting potential energy surfaces: applications to classical dynamics calculations

As a continuation of our efforts to develop efficient and accurate interpolating moving least-squares (IMLS) methods for generating potential energy surfaces, we carry out classical trajectories and compute kinetics properties on higher degree IMLS surfaces. In this study, we have investigated the choice of coordinate system, the range of points (i.e., the cutoff radius) used in fitting, and strategies for selections of data points and basis elements. We illustrate and test the method by applying it to hydrogen peroxide (HOOH). In particular, reaction rates for the O-O bond breaking in HOOH are calculated on fitted surfaces using the classical trajectory approach to test the accuracy of the IMLS method for providing potentials for dynamics calculations.     

A local interpolation scheme using no derivatives in potential sampling: application to O(1D) + H2 system

We recently proposed a local interpolation scheme, in which interpolant moving least squares (IMLS) and Shepard interpolation are employed to describe potential energy surfaces. This IMLS/Shepard scheme is used to interpolate quantum chemical potential energy surfaces for which analytical derivatives are not available. In this study, we apply the scheme to the highly exothermic O((1)D) + H(2) --> H + OH reaction and compare it with results based on Shepard interpolation using second-order Taylor expansions. An analytical surface is used to define the potential function so that errors in the interpolation function may accurately be determined. We find that the present scheme reproduces the correct reactive cross-sections more accurately than the Shepard scheme, and with rms errors for energy and gradients that are significantly smaller than those from Shepard interpolation. This occurs even though the present scheme does not utilize derivative and Hessian information, whereas the Shepard interpolation does. The Bayesian approach proposed by Bettens and Collins does not improve the IMLS/Shepard results significantly, although it does the Shepard-only approach. The accuracy of the IMLS/Shepard scheme is surprising, but can be explained by the more global nature of the interpolation.     

Using swarm intelligence for finding transition states and reaction paths

We describe an algorithm that explores potential energy surfaces (PES) and finds approximate reaction paths and transition states. A few (≈6) evolving atomic configurations ("climbers") start near a local minimum M1 of the PES. The climbers seek a shallow ascent, low energy, path toward a saddle point S12, cross over to another valley of the PES, and climb down to a new minimum M2 that was not known beforehand. Climbers use both energy and energy derivatives to make individual decisions, and they use relative fitness to make team-based decisions. In sufficiently long runs, they keep exploring and may go through a sequence M1-S12-M2-S23-M3 ... of minima and saddle points without revisiting any of the critical points. We report results on eight small test systems that highlight advantages and disadvantages of the method. We also investigated the PES of Li(8), Al(7)(+), Ag(7), and Ag(2)NH(3) to illustrate potential applications of this new method.     

An approach to reaction path branching using valley–ridge inflection points of potential-energy surfaces

Valley–ridge inflection (VRI) points of a potential-energy surface (PES) may have a strong relation to the occurrence of bifurcations along reaction pathways of molecular rearrangements. We discuss two different definitions of VRI points in the literature. The calculation of symmetric VRI points has already been reported [W. Quapp et al. (1998) Theor. Chem. Acc. 100: 285–299]. Here, we in addition calculate special asymmetric VRI points which are placed on gradient extremals (GE). Following a GE opens the possibility to find the VRI point on it. An application is presented to search for asymmetric VRI points near the isomerization valley of the PES of the HCN molecule. A new method for GE-following is based on a mathematical connection between the following of a reduced gradient and the calculation of GEs. The tangent search method to follow a GE to the smallest eigenvalue [W. Quapp et al. (2000) Theor. Chem. Acc. 105: 145–155] is extended to follow also GEs to higher eigenvalues in order to find a VRI point. The new method needs gradient and second derivatives of the PES only.     

Searching for saddle points of potential energy surfaces by following a reduced gradient

The old coordinate driving procedure to find transition structures in chemical systems is revisited. The well-known gradient criterion, ∇E( x )= 0 , which defines the stationary points of the potential energy surface (PES), is reduced by one equation corresponding to one search direction. In this manner, abstract curves can be defined connecting stationary points of the PES. Starting at a given minimum, one follows a well-selected coordinate to reach the saddle of interest. Usually, but not necessarily, this coordinate will be related to the reaction progress. The method, called reduced gradient following (RGF), locally has an explicit analytical definition. We present a predictor–corrector method for tracing such curves. RGF uses the gradient and the Hessian matrix or updates of the latter at every curve point. For the purpose of testing a whole surface, the six-dimensional PES of formaldehyde, H₂CO, was explored by RGF using the restricted Hartree–Fock (RHF) method and the STO-3G basis set. Forty-nine minima and saddle points of different indices were found. At least seven stationary points representing bonded structures were detected in addition to those located using another search algorithm on the same level of theory. Further examples are the localization of the saddle for the HCN⇌CNH isomerization (used for steplength tests) and for the ring closure of azidoazomethine to 1H-tetrazole. The results show that following the reduced gradient may represent a serious alternative to other methods used to locate saddle points in quantum chemistry. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1087–1100, 1998     

Efficiency considerations in the construction of interpolated potential energy surfaces for the calculation of quantum observables by diffusion Monte Carlo

A modified Shepard interpolation scheme is used to construct global potential energy surfaces (PES) in order to calculate quantum observables--vibrationally averaged internal coordinates, fully anharmonic zero-point energies and nuclear radial distribution functions--for a prototypical loosely bound molecular system, the water dimer. The efficiency of PES construction is examined with respect to (a) the method used to sample configurational space, (b) the method used to choose which points to add to the PES data set, and (c) the use of either a one- or two-part weight function. The most efficient method for constructing the PES is found to require a quantum sampling regime, a combination of both h-weight and rms methods for choosing data points and use of the two-part weight function in the interpolation. Using this regime, the quantum diffusion Monte Carlo zero-point energy converges to the exact result within addition of 50 data points. The vibrationally averaged O-O distance and O-O radial distribution function, however, converge more slowly and require addition of over 500 data points. The methods presented here are expected to be applicable to both other loosely bound complexes as well as tightly bound molecular species. When combined with high quality ab initio calculations, these methods should be able to accurately characterize the PES of such species.     

Haptic quantum chemistry

We present an implementation designed to physically experience quantum mechanical forces between reactants in chemical reactions. This allows one to screen the profile of potential energy surfaces for the study of reaction mechanisms. For this, we have developed a interface between the user and a virtual laboratory by means of a force‐feedback haptic device. Potential energy surfaces of chemical reactions can be explored efficiently by rendering in the haptic device the gradients calculated with first‐principles methods. The underlying potential energy surface is accurately fitted on the fly by the interpolating moving least‐squares (IMLS) scheme to a grid of quantum chemical electronic energies (and geometric gradients). In addition, we introduce a new IMLS‐based method to locate minimum‐energy paths between two points on a potential energy surface. © 2009 Wiley Periodicals, Inc. J Comput Chem 2009     

Geometrical properties of the potential energy of the soft-sphere binary mixture

We report a detailed study of the stationary points (zero-force points) of the potential energy surface (PES) of a model structural glassformer. We compare stationary points found with two different algorithms (eigenvector following and square gradient minimization), and show that the mapping between instantaneous configuration and stationary points defined by those algorithms is as different as to strongly influence the instability index K versus temperature plot, which relevance in analyzing the liquid dynamics is thus questioned. On the other hand, the plot of K versus energy is much less sensitive to the algorithm employed, showing that the energy is the good variable to discuss geometric properties of the PES. We find new evidence of a geometric transition between a minima-dominated phase and a saddle-point-dominated one. We analyze the distances between instantaneous configurations and stationary points, and find that above the glass transition, the system is closer to saddle points than to minima.     

Six-dimensional potential energy surface of the dissociative chemisorption of HCl on Au(111) using neural networks

We constructed a six-dimensional potential energy surface(PES)for the dissociative chemisorption of HCl on Au(111)using the neural networks method based on roughly 70000 energies obtained from extensive density functional theory(DFT)calculations.The resulting PES is accurate and smooth,based on the small fitting errors and good agreement between the fitted PES and the direct DFT calculations.Time-dependent wave packet calculations show that the potential energy surface is very well converged with respect to the number of DFT data points,as well as to the fitting process.The dissociation probabilities of HCl initially in the ground rovibrational state from six-dimensional quantum dynamical calculations are quite diferent from the four-dimensional fixed-site calculations,indicating it is essential to perform full-dimensional quantum dynamical studies for the title molecule-surface interaction system.     

Leading coordinate analysis of reaction pathways in proton chain transfer: Application to a two-proton transfer model for the green fluorescent protein

The ‘leading coordinate’ approach to computing an approximate reaction pathway, with subsequent determination of the true minimum energy profile, is applied to a two-proton chain transfer model based on the chromophore and its surrounding moieties within the green fluorescent protein (GFP). Using an ab initio quantum chemical method, a number of different relaxed energy profiles are found for several plausible guesses at leading coordinates. The results obtained for different trial leading coordinates are rationalized through the calculation of a two-dimensional relaxed potential energy surface (PES) for the system. Analysis of the 2-D relaxed PES reveals that two of the trial pathways are entirely spurious, while two others contain useful information and can be used to furnish starting points for successful saddle-point searches. Implications for selection of trial leading coordinates in this class of proton chain transfer reactions are discussed, and a simple diagnostic function is proposed for revealing whether or not a relaxed pathway based on a trial leading coordinate is likely to furnish useful information.     

Ab initio calculation and quasi‐classical dynamics study of the two lowest potential energy surfaces of the O(1D)+HBr system*

The lowest singlet 1¹A′ and 1¹A″ potential energy surfaces (PES) of the O(¹D)+HBr system have been ab initio computed. The complete active space self‐consistent field (CASSCF) method was used in most of the calculations, considering all the valence orbitals as active. The calculations were complemented with both analytical gradient calculations to characterize the stationary points and multireference configuration interaction (MRCI) calculations at selected nuclear geometries to improve the determination of the barrier heights and of the energetics. Electronic energy values for both PESs were then independently fitted by polynomial expansions in bond order coordinates. On the fitted surfaces quasi‐classical trajectories were separately run. Single‐surface calculations behave qualitatively different for the ground and the excited PES at low collision energies. A satisfactory agreement with existing experimental data was obtained by using the ground PES while calculations performed on the excited 1¹A″ PES worsened the agreement. However, when collision energy is increased, detailed experimental distributions are less well reproduced by calculations on the ground PES. This may imply the participation via nonadiabatic transitions of the 2¹A′ PES at higher energies while the adiabatic ground singlet PES well describes reactive scattering at low collision energy. © 2001 John Wiley & Sons, Inc. Int J Quantum Chem, 2001     

Following the streambed reaction on potential-energy surfaces: a new robust method

A simple procedure with low computational efforts is proposed to follow the reaction path of the potential-energy hypersurface (PES) starting from minima or saddle points. The method uses a modification of the so-called “following the reduced gradient” [Quapp W, Hirsch M, Imig O, Heidrich D (1998) J Comput Chem 19:1087]. The original method connects points where the gradient has a constant direction. In the present article the procedure is replaced by taking iterative varying directions of the gradient controlled by the last tangent of the searched curve. The resulting minimum energy path is that valley floor gradient extremal (GE) which belongs to the smallest (absolute) eigenvalue of the Hessian and, hence, that GE which usually leads along the streambed of a chemical reaction. The new method avoids third derivatives of the PES and obtains the GE of least ascent by second-order calculations only. Nevertheless, we are able to follow the streambed GE uphill or downhill. We can connect a minimum with its saddles if the streambed leads up to a saddle, or we find a turning point or a bifurcation point. The effectiveness and the characteristic properties of the new algorithm are demonstrated by using polynomial test surfaces, an ab initio PES of H₂O, and the analytic potentials of Lennard-Jones (LJ) clusters. By tracing the streambeds we located previously identified saddle points for LJ _N with N=3, 7, 8, and 55. Saddles for LJ _N with N=15, 20, and 30 as presented here are new results. Received: 8 March 2000 / Accepted: 17 July 2000 / Published online: 24 October 2000     

Newton trajectories for finding stationary points on molecular potential energy surfaces

We present a new algorithm for computing Newton trajectories based on the Quadratic String Method (QSM) and explain how this can be used to find key stationary points on the molecular potential energy surface (PES). This method starts by using the intersections of Newton trajectories to locate stationary points on the PES. These points could then be used to determine the minimum energy path. The new method, called QSM-NT, is shown to be efficient and reliable for both analytical potential energy surfaces and potential energy surfaces computed from quantum chemistry calculations. The advantages and pitfalls of this method for exploring PES are discussed. In particular, the problem of discontinuous Newton trajectories is elucidated.     

Interpolated potential energy surfaces: How accurate do the second derivatives have to be?

A global potential energy surface for the water dimer is constructed using the modified Shepard interpolation scheme of Collins et al. According to this interpolation scheme, the energy at an arbitrary geometry is expressed as a weighted sum of Taylor series expansions from neighboring data points, where the energy and derivative data required are obtained from ab initio calculations. For some ab initio methods, errors are introduced into the second derivative matrix, either by numerical differencing of ab initio energies or numerical integration during the ab initio calculation. Therefore, we test the accuracy required of the second derivative data by truncation of the exact second derivatives to a series of approximate second derivatives, and assess the effect on the results of a quantum diffusion Monte Carlo (QDMC) simulation. Our results show that the calculated zero-point energy and wave function histograms converge to within the numerical uncertainty of the QDMC simulation by inclusion of either three significant figures or three decimal places in the second derivatives.