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

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     

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.     

Gradient incorporation in one-dimensional applications of interpolating moving least-squares methods for fitting potential energy surfaces

We present several approaches to use gradients in higher degree interpolating moving least squares (IMLS) methods for representing a potential energy surface (PES). General procedures are developed to obtain smooth approximations of the PES and its derivatives from quasi-uniform sets of energy and gradient data points. These methods are illustrated and analyzed for the Morse oscillator and a 1-D slice of the ground-state PES for the HCO radical computed using density functional theory. Variations in the IMLS fits with the number and distribution of points and the degree of the polynomial fitting basis set are examined. We determine the effects of gradient inclusion on the accuracy of the IMLS values of the energy, first and second derivatives for two 1-D test cases. Gradient inclusion reduces the number of data points required by up to 40%.     

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.     

A nested molecule-independent neural network approach for high-quality potential fits

It is shown that neural networks (NNs) are efficient and effective tools for fitting potential energy surfaces. For H2O, a simple NN approach works very well. To fit surfaces for HOOH and H2CO, we develop a nested neural network technique in which we first fit an approximate NN potential and then use another NN to fit the difference of the true potential and the approximate potential. The root-mean-square error (RMSE) of the H2O surface is 1 cm(-1). For the 6-D HOOH and H2CO surfaces, the nested approach does almost as well attaining a RMSE of 2 cm(-1). The quality of the NN surfaces is verified by calculating vibrational spectra. For all three molecules, most of the low-lying levels are within 1 cm(-1) of the exact results. On the basis of these results, we propose that the nested NN approach be considered a method of choice for both simple potentials, for which it is relatively easy to guess a good fitting function, and complicated (e.g., double well) potentials for which it is much harder to deduce an appropriate fitting function. The number of fitting parameters is only moderately larger for the 6-D than for the 3-D potentials, and for all three molecules, decreasing the desired RMSE increases only slightly the number of required fitting parameters (nodes). NN methods, and in particular the nested approach we propose, should be good universal potential fitting tools.     

Formation of hydroxyl radical from the photolysis of frozen hydrogen peroxide

Hydrogen peroxide (HOOH) in ice and snow is an important chemical tracer for the oxidative capacities of past atmospheres. However, photolysis in ice and snow will destroy HOOH and form the hydroxyl radical (*OH), which can react with snowpack trace species. Reactions of *OH in snow and ice will affect the composition of both the overlying atmosphere (e.g., by the release of volatile species such as formaldehyde to the boundary layer) and the snow and ice (e.g., by the *OH-mediated destruction of trace organics). To help understand these impacts, we have measured the quantum yield of *OH from the photolysis of HOOH on ice. Our measured quantum yields (Phi(HOOH --> *OH)) are independent of ionic strength, pH, and wavelength, but are dependent upon temperature. This temperature dependence for both solution and ice data is best described by the relationship ln(Phi(HOOH --> *OH)) = -(684 +/- 17)(1/T) + (2.27 +/- 0.064) (where errors represent 1 standard error). The corresponding activation energy (Ea) for HOOH (5.7 kJ mol(-1)) is much smaller than that for nitrate photolysis, indicating that the photochemistry of HOOH is less affected by changes in temperature. Using our measured quantum yields, we calculate that the photolytic lifetimes of HOOH in surface snow grains under midday, summer solstice sunlight are approximately 140 h at representative sites on the Greenland and Antarctic ice sheets. In addition, our calculations reveal that the majority of *OH radicals formed on polar snow grains are from HOOH photolysis, while nitrate photolysis is only a minor contributor. Similarly, HOOH appears to be much more important than nitrate as a photochemical source of *OH on cirrus ice clouds, where reactions of the photochemically formed hydroxyl radical could lead to the release of oxygenated volatile organic compounds to the upper troposphere.     

Basis set effects on the energy and hardness profiles of the hydrogen fluoride dimer

In earlier work, the present authors have shown that hardness profiles are less dependent on the level of calculation than energy profiles for potential energy surfaces (PESs) having pathological behaviors. At variance with energy profiles, hardness profiles always show the correct number of stationary points. This characteristic has been used to indicate the existence of spurious stationary points on the PESs. In the present work, we apply this methodology to the hydrogen fluoride dimer, a classical difficult case for the density functional theory methods.     

Comparison of hydrated hydroperoxide anion (HOO-)(H2O)n clusters with alkaline hydrogen peroxide (HOOH)(OH-)(H2O)(n-1) clusters, n = 1-8, 20: an ab initio study

Hydroperoxide anion (HOO(-)), the conjugate base of hydrogen peroxide (HOOH), has been relatively little studied despite the importance of HOOH in commercial processes, atmospheric science, and biology. The anion has been shown to exist as a stable species in alkaline water. This project explored the structure of gas phase (HOO(-))(H(2)O)(n) clusters and identified the lowest energy configurations for n ≤ 8 at the B3LYP/6-311++G** level of theory and for n ≤ 6 at the MP2/aug-cc-pVTZ level of theory. As a start toward understanding equilibration between HOO(-) and HOOH in an alkaline environment, (HOOH)(OH(-))(H(2)O)(n-1) clusters were likewise examined, and the lowest energy configurations were determined for n ≤ 8 (B3LYP/6-311++G**) and n ≤ 6 (MP2/aug-cc-pVTZ). Some studies were also done for n = 20. The two species have very different solvation behaviors. In low energy (HOOH)(OH(-))(H(2)O)(n-1) clusters, HOOH sits on the surface of the cluster, is 4-coordinated (each O is donor once and acceptor once), and donates to the hydroxide ion. In contrast, in low energy (HOO(-))(H(2)O)(n) clusters, (HOO(-)) takes a position in the cluster center surrounded on all sides by water molecules, and its optimum coordination number appears to be 7 (one O is donor-acceptor-acceptor while the other is a 4-fold acceptor). For n ≤ 6 the lowest (HOOH)(OH(-))(H(2)O)(n-1) cluster lies 1.0-2.1 kcal/mol below the lowest (HOO(-))(H(2)O)(n) cluster, but the lowest clusters found for n = 20 favor (HOO(-))(H(2)O)(20). The results suggest that ambient water could act as a substantial kinetic brake that slows equilibration between (HOOH)(OH(-)) and (HOO(-))(H(2)O) because extensive rearrangement of solvation shells is necessary to restabilize either species after proton transfer.     

Iterative solutions with energy selected bases for highly excited vibrations of tetra-atomic molecules

The use of energy selected bases (ESB) with iterative diagonalization of the Hamiltonian matrix is described for vibrations of tetra-atomic systems. The performance of the method is tested by computing vibrational states of HOOH below 10,000 cm(-1) (1296 A+ symmetry states) and H(2)CO below 13,500 cm(-1) (729 A(1) symmetry states). For iterative solutions, we tested both the implicitly restarted Lanczos method (IRLM) and the standard (nonreorthogonalizing) Lanczos approach. Comparison with other contracted basis approach as well as direct product grid representation shows superior performance of the ESB/IRLM approach. Of the two systems, H(2)CO is found to be more challenging than HOOH since it has much stronger couplings among vibrational modes, which leads to a drastically larger primitive basis set. For H(2)CO we also discuss some interesting behavior of the molecule in the high internal energy regime.     

Testing the performance of density functionals for the calculation of energetic properties of complex-forming radical-molecule reactions

The performance of some popular and some more recent density functional methods for the calculation of energies of stationary points on the potential surfaces of radical-molecule reactions was examined. The functionals studied are B3-LYP, BH&H, BH&H-LYP, MPW1K, MPWB1K, TPSS, TPSSh, BB1K, M05 and M05-2X, in conjunction with nine different AO basis sets. The reaction energies, barrier heights and the relative energies of the pre-and post-reaction complexes were compared with those obtained at the CCSD(T)/CBS limit for the reactions of OH radicals with HOOH and CH₃OOH. Very poor barrier heights are provided by the B3-LYP, TPSS and TPSSh functionals. The best overall performance was obtained with the BB1K, MPW1K and MPWB1K functionals. In these reactions all of the studied functionals provide converged results only if they are used with large basis sets like aug-cc-pVTZ and def2-TZVP. The data show that before relying on a functional for a specific reaction, it is desirable to make some test calculations on the performance. The same functional can predict some relative energies very well and some others very poorly even in systems including chemically similar reactants.     

Reaction dynamics simulations of the identity S(N)2 reaction H(2)O + HOOH(2)(+)--> H(2)OOH(+)+ H(2)O. Requirements for reaction and competition with proton transfer

Despite the fact that the transition structure of the gas phase S(N)2 reaction H(2)O + HOOH(2)(+)--> HOOH(2)(+)+ H(2)O is well below the reactants in potential energy, the reaction has not yet been observed by experiment. Variational transition state RRKM theory reveals a strong preference for the competing proton transfer reaction H(2)O + HOOH(2)(+)--> H(3)O(+)+ HOOH due to entropy factors. Born-Oppenheimer reaction dynamics simulations confirm these results. However, by increasing the collision energy to around 7.5 eV the probability for nucleophilic substitution increases relative to proton transfer. These observations are explained by the presence of the key common intermediate HOO(H)[dot dot dot]H-OH(2)(+) which leads to effective proton transfer, but can be avoided with increasing collision energy. However, the S(N)2 probability remains below 0.2 since successful passage through the TS requires optimum initial orientation of the reactants, excitation of the relative translational motion and good phase correlation between the O-O vibration and the motion of the incoming water.     

Multidentate Anchors for Surface Functionalization

The bottom‐up functionalization of solid surfaces shows increasing importance for a wide range of interdisciplinary applications. Multidentate anchors with more than two contact points can bind to solid surfaces with strong chemisorption, well‐defined upright configuration, and tailored functionality. The surface functionalization using multidentate anchors with three (tripodal), four (quadripodal), or more binding points is summarized herein, with a focus on those beyond classical tripodal anchors. In particular, the molecular design on how to achieve multisite interaction between anchor and substrate and the introduction of functional groups to thin films are discussed.     

Quantum-classical dynamics of scattering processes in adiabatic and diabatic representations

We demonstrate the workability of a TDDVR based [J. Chem. Phys. 118, 5302 (2003)], novel quantum-classical approach, for simulating scattering processes on a quasi-Jahn-Teller model [J. Chem. Phys. 105, 9141 (1996)] surface. The formulation introduces a set of DVR grid points defined by the Hermite part of the basis set in each dimension and allows the movement of grid points around the central trajectory. With enough trajectories (grid points), the method converges to the exact quantum formulation whereas with only one grid point, we recover the conventional molecular dynamics approach. The time-dependent Schrodinger equation and classical equations of motion are solved self-consistently and electronic transitions are allowed anywhere in the configuration space among any number of coupled states. Quantum-classical calculations are performed on diabatic surfaces (two and three) to reveal the effects of symmetry on inelastic and reactive state-to-state transition probabilities, along with calculations on an adiabatic surface with ordinary Born-Oppenheimer approximation. Excellent agreement between TDDVR and DVR results is obtained in both the representations.     

Quantum mechanical/molecular mechanical study of mechanisms of heme degradation by the enzyme heme oxygenase: the strategic function of the water cluster

Heme degradation by heme oxygenase (HO) enzymes is important in maintaining iron homeostasis and prevention of oxidative stress, etc. In response to mechanistic uncertainties, we performed quantum mechanical/molecular mechanical investigations of the heme hydroxylation by HO, in the native route and with the oxygen surrogate donor H2O2. It is demonstrated that H2O2 cannot be deprotonated to yield Fe(III)OOH, and hence the surrogate reaction starts from the FeHOOH complex. The calculations show that, when starting from either Fe(III)OOH or Fe(III)HOOH, the fully concerted mechanism involving O-O bond breakage and O-C(meso) bond formation is highly disfavored. The low-energy mechanism involves a nonsynchronous, effectively concerted pathway, in which the active species undergoes first O-O bond homolysis followed by a barrier-free (small with Fe(III)HOOH) hydroxyl radical attack on the meso position of the porphyrin. During the reaction of Fe(III)HOOH, formation of the Por+*FeIV=O species, compound I, competes with heme hydroxylation, thereby reducing the efficiency of the surrogate route. All these conclusions are in accord with experimental findings (Chu, G. C.; Katakura, K.; Zhang, X.; Yoshida, T.; Ikeda-Saito, M. J. Biol. Chem. 1999, 274, 21319). The study highlights the role of the water cluster in the distal pocket in creating "function" for the enzyme; this cluster affects the O-O cleavage and the O-Cmeso formation, but more so it is responsible for the orientation of the hydroxyl radical and for the observed alpha-meso regioselectivity of hydroxylation (Ortiz de Montellano, P. R. Acc. Chem. Res. 1998, 31, 543). Differences/similarities with P450 and HRP are discussed.     

Biproton transfer in associates of the formic acid molecule with H2O...YH systems (YH=H2O,HOF and HOOH)

An ab initio (STO-3G basis set) method is used to calculate the pathways and activation barriers of cooperative biproton migrations in the associates of the formic acid molecule with water and YH (YH=H₂O, HOF, and HOOH). The associate with YH=H₂O exhibits a cooperative triproton transfer with an activation barrier of 4.46 (C_s symmetry) or 3.08 kcal/mole (C₂ symmetry); in the systems with Y=HOF and HOOH, the biproton transfer is accompanied by synchronous shifts of two hydroxyl groups and has the activation barriers of 75.03 (C₂, Y=HOF), 75.74 (C_s, Y=HOF), and 89.62 kcal/mole (Y=HOOH). The hypothesis of stereochemical similarity for low-barrier cooperative processes in molecular associates is proposed. Scientific Research Institute of Physical and Organic Chemistry, Rostov State University. Translated fromZhurnal Strukturnoi Khimii, Vol. 37, No. 5, pp. 834–844, September–October, 1996. Translated by I. Izvekova