Pressure Effects on the Relaxation of an Excited Ethane Molecule in High-Pressure Bath Gases.

We use molecular dynamics to calculate the rotational and vibrational energy relaxation of C2H6 in Ar, Kr, and Xe bath gases over a pressure range of 10-400 atm and at temperatures of 300 and 800 K. The C2H6 is instantaneously excited by 80 kcal/mol randomly distributed into both vibrational and rotational modes. The computed relaxation rates show little sensitivity to the identity of the noble gas in the bath. Vibrational relaxation rates show a nonlinear pressure dependence at 300 K. At 800 K the reduced range of bath gas densities covered by the range of pressures does not yet show any nonlinearity in the pressure dependence. Rotational relaxation is characterized with two relaxation rates. The slower rate is comparable to the vibrational relaxation rate. The faster rate has a linear pressure dependence at 300 K but an irregular, nonlinear pressure dependence at 800 K. To understand this, a model was developed based on approximating the periodic box used in the molecular dynamics simulations by an equal-volume collection of cubes where each cube is sized to allow only single occupancy by the noble gas or the molecule. Combinatorial statistics then leads to a pressure- and temperature-dependent analytic distribution of the bath gas species the molecule encounters in a collision. This distribution, the dissociation energy of molecule/bath gas complexes and bath gas clusters, and the computed energy release per collision combine to show that only at 300 K is the energy release sufficient to dissociate likely complexes and clusters. This suggests that persistent and pressure-dependent clusters and complexes at 800 K may be responsible for the nonlinear pressure dependence of rotational relaxation.

Pressure effects on the vibrational and rotational relaxation of vibrationally excited OH (ν, J) in an argon bath.

Quasi-classical molecular dynamics simulations were used to study the energy relaxation of an initially non-rotating, vibrationally excited (ν = 4) hydroxyl radical (OH) in an Ar bath at 300 K and at high pressures from 50 atm to 400 atm. A Morse oscillator potential represented the OH, and two sets of interaction potentials were used based on whether the Ar-H potential was a Buckingham (Exp6) or a Lennard-Jones (LJ) potential. The vibrational and rotational energies were monitored for 25 000-90 000 ps for Exp6 trajectories and 5000 ps for LJ trajectories. Comparisons to measured vibrational relaxation rates show that Exp6 rates are superior. Simulated initial vibrational relaxation rates are linearly proportional to pressure, implying no effect of high-pressure breakdown in the isolated binary collision approximation. The vibrational decay curves upward from single-exponential decay. A model based on transition rates that exponentially depend on the anharmonic energy gap between vibrational levels fits the vibrational decay well at all pressures, suggesting that anharmonicity is a major cause of the curvature. Due to the competition of vibration-to-rotation energy transfer and bath gas relaxation, the rotational energy overshoots and then relaxes to its thermal value. Approximate models with adjustable rates for this competition successfully reproduced the rotational results. These models show that a large fraction of the vibrational energy loss is initially converted to rotational energy but that fraction decreases rapidly as the vibrational energy content of OH decreases. While simulated rates change dramatically between Exp6 and LJ potentials, the mechanisms remain the same.

Normal mode analysis on the relaxation of an excited nitromethane molecule in argon bath.

In our previous work [Rivera-Rivera et al., J. Chem. Phys. 142, 014303 (2015)], classical molecular dynamics simulations followed the relaxation, in a 300 K Ar bath at a pressure of 10-400 atm, of nitromethane (CH3NO2) instantaneously excited by statistically distributing 50 kcal/mol among all its internal degrees of freedom. Both rotational and vibrational energies decayed with nonexponential curves. The present work explores mode-specific mechanisms at work in the decay process. With the separation of rotation and vibration developed by Rhee and Kim [J. Chem. Phys. 107, 1394 (1997)], one can show that the vibrational kinetic energy decomposes only into vibrational normal modes, while the rotational and Coriolis energies decompose into both vibrational and rotational normal modes. The saved CH3NO2 positions and momenta were converted into mode-specific energies whose decay was monitored over 1000 ps. The results identify vibrational and rotational modes that promote/resist energy lost and drive nonexponential behavior.

SIESTA-SIPs: Massively parallel spectrum-slicing eigensolver for an ab initio molecular dynamics package.

Integration of Shift-and-Invert Parallel Spectral Transformation (SIPs) eigensolver (as implemented in the SLEPc library) into an ab initio molecular dynamics package, SIESTA, is described. The effectiveness of the code is demonstrated on applications to polyethylene chains, boron nitride sheets, and bulk water clusters. For problems with the same number of orbitals, the performance of the SLEPc eigensolver depends on the sparsity of the matrices involved, favoring reduced dimensional systems such as polyethylene or boron nitride sheets in comparison to bulk systems like water clusters. For all problems investigated, performance of SIESTA-SIPs exceeds the performance of SIESTA with default solver (ScaLAPACK) at the larger number of cores and the larger number of orbitals. A method that improves the load-balance with each iteration in the self-consistency cycle by exploiting the emerging knowledge of the eigenvalue spectrum is demonstrated. © 2018 Wiley Periodicals, Inc.

Anharmonic Rovibrational Partition Functions for Fluxional Species at High Temperatures via Monte Carlo Phase Space Integrals.

Monte Carlo phase space integration (MCPSI) is used to compute full dimensional and fully anharmonic, but classical, rovibrational partition functions for 22 small- and medium-sized molecules and radicals. Several of the species considered here feature multiple minima and low-frequency nonlocal motions, and efficiently sampling these systems is facilitated using curvilinear (stretch, bend, and torsion) coordinates. The curvilinear coordinate MCPSI method is demonstrated to be applicable to the treatment of fluxional species with complex rovibrational structures and as many as 21 fully coupled rovibrational degrees of freedom. Trends in the computed anharmonicity corrections are discussed. For many systems, rovibrational anharmonicities at elevated temperatures are shown to vary consistently with the number of degrees of freedom and with temperature once rovibrational coupling and torsional anharmonicity are accounted for. Larger corrections are found for systems with complex vibrational structures, such as systems with multiple large-amplitude modes and/or multiple minima.

Shift-and-invert parallel spectral transformation eigensolver: Massively parallel performance for density-functional based tight-binding.

The Shift-and-invert parallel spectral transformations (SIPs), a computational approach to solve sparse eigenvalue problems, is developed for massively parallel architectures with exceptional parallel scalability and robustness. The capabilities of SIPs are demonstrated by diagonalization of density-functional based tight-binding (DFTB) Hamiltonian and overlap matrices for single-wall metallic carbon nanotubes, diamond nanowires, and bulk diamond crystals. The largest (smallest) example studied is a 128,000 (2000) atom nanotube for which ∼330,000 (∼5600) eigenvalues and eigenfunctions are obtained in ∼190 (∼5) seconds when parallelized over 266,144 (16,384) Blue Gene/Q cores. Weak scaling and strong scaling of SIPs are analyzed and the performance of SIPs is compared with other novel methods. Different matrix ordering methods are investigated to reduce the cost of the factorization step, which dominates the time-to-solution at the strong scaling limit. A parallel implementation of assembling the density matrix from the distributed eigenvectors is demonstrated.

Pressure effects on the relaxation of an excited nitromethane molecule in an argon bath.

Classical molecular dynamics simulations were performed to study the relaxation of nitromethane in an Ar bath (of 1000 atoms) at 300 K and pressures 10, 50, 75, 100, 125, 150, 300, and 400 atm. The molecule was instantaneously excited by statistically distributing 50 kcal/mol among the internal degrees of freedom. At each pressure, 1000 trajectories were integrated for 1000 ps, except for 10 atm, for which the integration time was 5000 ps. The computed ensemble-averaged rotational energy decay is ∼100 times faster than the vibrational energy decay. Both rotational and vibrational decay curves can be satisfactorily fit with the Lendvay-Schatz function, which involves two parameters: one for the initial rate and one for the curvature of the decay curve. The decay curves for all pressures exhibit positive curvature implying the rate slows as the molecule loses energy. The initial rotational relaxation rate is directly proportional to density over the interval of simulated densities, but the initial vibrational relaxation rate decreases with increasing density relative to the extrapolation of the limiting low-pressure proportionality to density. The initial vibrational relaxation rate and curvature are fit as functions of density. For the initial vibrational relaxation rate, the functional form of the fit arises from a combinatorial model for the frequency of nitromethane "simultaneously" colliding with multiple Ar atoms. Roll-off of the initial rate from its low-density extrapolation occurs because the cross section for collision events with L Ar atoms increases with L more slowly than L times the cross section for collision events with one Ar atom. The resulting density-dependent functions of the initial rate and curvature represent, reasonably well, all the vibrational decay curves except at the lowest density for which the functions overestimate the rate of decay. The decay over all gas phase densities is predicted by extrapolating the fits to condensed-phase densities.

Theoretical/experimental comparison of deep tunneling decay of quasi-bound H(D)OCO to H(D) + CO2.

The measured H(D)OCO survival fractions of the photoelectron-photofragment coincidence experiments by the Continetti group are qualitatively reproduced by tunneling calculations to H(D) + CO2 on several recent ab initio potential energy surfaces for the HOCO system. The tunneling calculations involve effective one-dimensional barriers based on steepest descent paths computed on each potential energy surface. The resulting tunneling probabilities are converted into H(D)OCO survival fractions using a model developed by the Continetti group in which every oscillation of the H(D)-OCO stretch provides an opportunity to tunnel. Four different potential energy surfaces are examined with the best qualitative agreement with experiment occurring for the PIP-NN surface based on UCCSD(T)-F12a/AVTZ electronic structure calculations and also a partial surface constructed for this study based on CASPT2/AVDZ electronic structure calculations. These two surfaces differ in barrier height by 1.6 kcal/mol but when matched at the saddle point have an almost identical shape along their reaction paths. The PIP surface is a less accurate fit to a smaller ab initio data set than that used for PIP-NN and its computed survival fractions are somewhat inferior to PIP-NN. The LTSH potential energy surface is the oldest surface examined and is qualitatively incompatible with experiment. This surface also has a small discontinuity that is easily repaired. On each surface, four different approximate tunneling methods are compared but only the small curvature tunneling method and the improved semiclassical transition state method produce useful results on all four surfaces. The results of these two methods are generally comparable and in qualitative agreement with experiment on the PIP-NN and CASPT2 surfaces. The original semiclassical transition state theory method produces qualitatively incorrect tunneling probabilities on all surfaces except the PIP. The Eckart tunneling method uses the least amount of information about the reaction path and produces too high a tunneling probability on PIP-NN surface, leading to survival fractions that peak at half their measured values.

Improved multidimensional semiclassical tunneling theory.

We show that the analytic multidimensional semiclassical tunneling formula of Miller et al. [Miller, W. H.; Hernandez, R.; Handy, N. C.; Jayatilaka, D.; Willets, A. Chem. Phys. Lett. 1990, 172, 62] is qualitatively incorrect for deep tunneling at energies well below the top of the barrier. The origin of this deficiency is that the formula uses an effective barrier weakly related to the true energetics but correctly adjusted to reproduce the harmonic description and anharmonic corrections of the reaction path at the saddle point as determined by second order vibrational perturbation theory. We present an analytic improved semiclassical formula that correctly includes energetic information and allows a qualitatively correct representation of deep tunneling. This is done by constructing a three segment composite Eckart potential that is continuous everywhere in both value and derivative. This composite potential has an analytic barrier penetration integral from which the semiclassical action can be derived and then used to define the semiclassical tunneling probability. The middle segment of the composite potential by itself is superior to the original formula of Miller et al. because it incorporates the asymmetry of the reaction barrier produced by the known reaction exoergicity. Comparison of the semiclassical and exact quantum tunneling probability for the pure Eckart potential suggests a simple threshold multiplicative factor to the improved formula to account for quantum effects very near threshold not represented by semiclassical theory. The deep tunneling limitations of the original formula are echoed in semiclassical high-energy descriptions of bound vibrational states perpendicular to the reaction path at the saddle point. However, typically ab initio energetic information is not available to correct it. The Supporting Information contains a Fortran code, test input, and test output that implements the improved semiclassical tunneling formula.

A classical trajectory study of the dissociation and isomerization of C2H5.

Motivated by photodissociation experiments in which non-RRKM nanosecond lifetimes of the ethyl radical were reported, we have performed a classical trajectory study of the dissociation and isomerization of C2H5 over the energy range 100-150 kcal/mol. We used a customized version of the AIREBO semiempirical potential (Stuart, S. J.; et al. J. Chem. Phys. 2000, 112, 6472-6486) to more accurately describe the gas-phase decomposition of C2H5. This study constitutes one of the first gas-phase applications of this potential form. At each energy, 10,000 trajectories were run and all underwent dissociation in less than 100 ps. The calculated dissociation rate constants are consistent with RRKM models; no evidence was found for nanosecond lifetimes. An analytic kinetics model of isomerization/dissociation competition was developed that incorporated incomplete mode mixing through a postulated divided phase space. The fits of the model to the trajectory data are good and represent the trajectory results in detail through repeated isomerizations at all energies. The model correctly displays single exponential decay at lower energies, but at higher energies, multiexponential decay due to incomplete mode mixing becomes more apparent. At both ends of the energy range, we carried out similar trajectory studies on CD2CH3 to examine isotopic scrambling. The results largely support the assumption that a H or a D atom is equally likely to dissociate from the mixed-isotope methyl end of the molecule. The calculated fraction of products that have the D atom dissociation is ∼20%, twice the experimental value available at one energy within our range. The calculated degree of isotopic scrambling is non-monotonic with respect to energy due to a non-monotonic ratio of the isomerization to dissociation rate constants.

A classical trajectory study of the intramolecular dynamics, isomerization, and unimolecular dissociation of HO2.

The classical dynamics and rates of isomerization and dissociation of HO2 have been studied using two potential energy surfaces (PESs) based on interpolative fittings of ab initio data: An interpolative moving least-squares (IMLS) surface [A. Li, D. Xie, R. Dawes, A. W. Jasper, J. Ma, and H. Guo, J. Chem. Phys. 133, 144306 (2010)] and the cubic-spline-fitted PES reported by Xu, Xie, Zhang, Lin, and Guo (XXZLG) [J. Chem. Phys. 127, 024304 (2007)]. Both PESs are based on similar, though not identical, internally contracted multi-reference configuration interaction with Davidson correction (icMRCI+Q) electronic structure calculations; the IMLS PES includes complete basis set (CBS) extrapolation. The coordinate range of the IMLS PES is limited to non-reactive processes. Surfaces-of-section show similar generally regular phase space structures for the IMLS and XXZLG PESs with increasing energy. The intramolecular vibrational energy redistribution (IVR) at energies above and below the threshold of isomerization is slow, especially for O-O stretch excitations, consistent with the regularity in the surfaces-of-section. The slow IVR rates lead to mode-specific effects that are prominent for isomerization (on both the IMLS and XXZLG) and modest for unimolecular dissociation to H + O2 (accessible only on the XXZLG PES). Even with statistical distributions of initial energy, slow IVR rates result in double exponential decay for isomerization, with the slower rate correlated with slow IVR rates for O-O vibrational excitation. The IVR and isomerization rates computed for the IMLS and XXZLG PESs are quantitatively, but not qualitatively, different from one another with the largest differences ascribed to the ~2 kcal/mol difference in the isomerization barrier heights. The IMLS and XXZLG results are compared with those obtained using the global, semi-empirical double-many-body expansion DMBE-IV PES [M. R. Pastrana, L. A. M. Quintales, J. Brandão, and A. J. C. Varandas, J. Chem. Phys. 94, 8073 (1990)], for which the surfaces-of-section display more irregular phase space structure, much faster IVR rates, and significantly less mode-specific effects in isomerization and unimolecular dissociation. The calculated IVR results for all three PESs are reasonably well represented by an analytic, coupled three-mode energy transfer model.

Ab initio wavenumber accurate spectroscopy: 1CH2 and HCN vibrational levels on automatically generated IMLS potential energy surfaces.

We report here calculated J = 0 vibrational frequencies for (1)CH(2) and HCN with root-mean-square error relative to available measurements of 2.0 cm(-1) and 3.2 cm(-1), respectively. These results are obtained with DVR calculations with a dense grid on ab initio potential energy surfaces (PESs). The ab initio electronic structure calculations employed are Davidson-corrected MRCI calculations with double-, triple-, and quadruple-zeta basis sets extrapolated to the complete basis set (CBS) limit. In the (1)CH(2) case, Full CI tests of the Davidson correction at small basis set levels lead to a scaling of the correction with the bend angle that can be profitably applied at the CBS limit. Core-valence corrections are added derived from CCSD(T) calculations with and without frozen cores. Relativistic and non-Born-Oppenheimer corrections are available for HCN and were applied. CBS limit CCSD(T) and CASPT2 calculations with the same basis sets were also tried for HCN. The CCSD(T) results are noticeably less accurate than the MRCI results while the CASPT2 results are much poorer. The PESs were generated automatically using the local interpolative moving least-squares method (L-IMLS). A general triatomic code is described where the L-IMLS method is interfaced with several common electronic structure packages. All PESs were computed with this code running in parallel on eight processors. The L-IMLS method provides global and local fitting error measures important in automatically growing the PES from initial ab initio seed points. The reliability of this approach was tested for (1)CH(2) by comparing DVR-calculated vibrational levels on an L-IMLS ab initio surface with levels generated by an explicit ab initio calculation at each DVR grid point. For all levels ( approximately 200) below 20 000 cm(-1), the mean unsigned difference between the levels of these two calculations was 0.1 cm(-1), consistent with the L-IMLS estimated mean unsigned fitting error of 0.3 cm(-1). All L-IMLS PESs used in this work have comparable mean unsigned fitting errors, implying that fitting errors have a negligible role in the final errors of the computed vibrational levels with experiment. Less than 500 ab initio calculations of the energy and gradients are required to achieve this level of accuracy.

Interpolating moving least-squares methods for fitting potential energy surfaces: using classical trajectories to explore configuration space.

We develop two approaches for growing a fitted potential energy surface (PES) by the interpolating moving least-squares (IMLS) technique using classical trajectories. We illustrate both approaches by calculating nitrous acid (HONO) cis-->trans isomerization trajectories under the control of ab initio forces from low-level HF/cc-pVDZ electronic structure calculations. In this illustrative example, as few as 300 ab initio energy/gradient calculations are required to converge the isomerization rate constant at a fixed energy to approximately 10%. Neither approach requires any preliminary electronic structure calculations or initial approximate representation of the PES (beyond information required for trajectory initial conditions). Hessians are not required. Both approaches rely on the fitting error estimation properties of IMLS fits. The first approach, called IMLS-accelerated direct dynamics, propagates individual trajectories directly with no preliminary exploratory trajectories. The PES is grown "on the fly" with the computation of new ab initio data only when a fitting error estimate exceeds a prescribed tight tolerance. The second approach, called dynamics-driven IMLS fitting, uses relatively inexpensive exploratory trajectories to both determine and fit the dynamically accessible configuration space. Once exploratory trajectories no longer find configurations with fitting error estimates higher than the designated accuracy, the IMLS fit is considered to be complete and usable in classical trajectory calculations or other applications.

Benchmark Calculations of Energetic Properties of Groups 4 and 6 Transition Metal Oxide Nanoclusters Including Comparison to Density Functional Theory.

The heats of formation and the normalized clustering energies (NCEs) for the group 4 and group 6 transition metal oxide (TMO) trimers and tetramers have been calculated by the Feller-Peterson-Dixon (FPD) method. The heats of formation predicted by the FPD method do not differ much from those previously derived from the NCEs at the CCSD(T)/aT level except for the CrO3 nanoclusters. New and improved heats of formation for Cr3O9 and Cr4O12 were obtained using PW91 orbitals instead of Hartree-Fock (HF) orbitals. Diffuse functions are necessary to predict accurate heats of formation. The fluoride affinities (FAs) are calculated with the CCSD(T) method. The relative energies (REs) of different isomers, NCEs, electron affinities (EAs), and FAs of (MO2)n (M = Ti, Zr, Hf, n = 1-4) and (MO3)n (M = Cr, Mo, W, n = 1-3) clusters have been benchmarked with 55 exchange-correlation density functional theory (DFT) functionals including both pure and hybrid types. The absolute errors of the DFT results are mostly less than ±10 kcal/mol for the NCEs and the EAs and less than ±15 kcal/mol for the FAs. Hybrid functionals usually perform better than the pure functionals for the REs and NCEs. The performance of the two types of functionals in predicting EAs and FAs is comparable. The B1B95 and PBE1PBE functionals provide reliable energetic properties for most isomers. Long range corrected pure functionals usually give poor FAs. The standard deviation of the absolute error is always close to the mean errors, and the probability distributions of the DFT errors are often not Gaussian (normal). The breadth of the distribution of errors and the maximum probability are dependent on the energy property and the isomer.

Vectorization of a classical trajectory code on a floating point systems, Inc. Model 164 attached processor.

A triatomic classical trajectory code has been modified by extensive vectorization of the algorithms to achieve much improved performance on an FPS 164 attached processor. Extensive timings on both the FPS 164 and a VAX 11/780 with floating point accelerator are presented as a function of the number of trajectories simultaneously run. The timing tests involve a potential energy surface of the LEPS variety and trajectories with 1000 time steps. The results indicate that vectorization results in timing improvements on both the VAX and the FPS. For larger numbers of trajectories run simultaneously, up to a factor of 25 improvement in speed occurs between VAX and FPS vectorized code.

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: 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: 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.