On coordinate transformations in steepest descent path and stationary point locations

Two aspects of the problems of calculating steepest descent paths and locating stationary points on surfaces E( X ), which are sources of some confusion in the literature, are addressed. These include writing proper expressions for the gradient and Hessian, and their transformation properties relative to coordinate transformations, based on the invariance of the surface E( X ). The appropriate transformation is derived, based on a constrained energy minimization condition, to achieve what we call the Hessian eigenvalue representation. This not only allows decoupling of the variables, but also points to the minimization direction and preserves the eigenvalues of the Hessian. These results allow one to use the steepest descent path and stationary point location algorithms in any coordinate system and obtain invariant results. The validity of these considerations are also confirmed through numerical examples. The stationary condition with constrained kinematic path length is also shown to yield a Hessian eigenvalue representation for the normal modes for small vibrations. Lastly, we have constructed a mathematically consistent definition of mass-weighted Cartesians where the intrinsic reaction path of Fukui is a steepest descent path. © 1992 John Wiley & Sons, Inc.     

Search for stationary points of arbitrary index by augmented Hessian method

Convergence properties of the augmented Hessian (AH ) method when searching for stationary points of an arbitrary fixed index are investigated. It is shown that the displacement vector of this method is proportional to one of the Hessian eigenvectors if the current point is far from a stationary one of the required index. A simple and reliable criterion for nearness of the current point to a stationary one of the desired index is proposed. The efficiency of a new one-dimensional optimization scheme that uses this criterion is studied. The case of coincidence of Hessian eigenvalues, which is a bottleneck of the standard AH method, is analyzed. A relation of the AH method to those by Poppinger and Wales is outlined. The correctness of the results obtained is illustrated on an example of a model surface. © 1995 John Wiley & Sons, Inc.     

Are most of the stationary points in a molecular association minima? Application of Fraga's potential to benzene–benzene

The importance of characterizing the stationary points of the intermolecular potential by means of Hessian eigenvalues is illustrated for the calculation of the benzene–benzene interaction using an atom-to-atom pair potential proposed by Fraga (FAAP). Two models, the standard one-center-per atom and another using three-centers-per atom due to Hunter and Sanders, are used to evaluate the electrostatic contributions and the results are compared. It is found in both cases that although using low-gradient thresholds allows optimization procedures to avoid many stationary points that are not true minima computing time considerations makes the usual procedure of using high-gradient thresholds [say, 10^?2 kj/(mol Å)] as the most efficient. Moreover, this later procedure can be recommended because the actual minima can be characterized by means of Hessian eigenvalues even if these high-gradient thresholds are used, and further decreasing of the convergence criterion does not imply significant modifications in the geometric parameters of the minima. The possible advantages of using the three-centers-per-atom model for the calculation of molecular associations between delocalized systems are also discussed on the basis of the agreement of the benzene–benzene results with experimental and theoretical data taken from the literature. © 1993 John Wiley & Sons, Inc.     

The search for stationary points on a quantum mechanical/molecular mechanical potential-energy surface

 We propose a methodology to locate stationary points on a quantum mechanical/molecular mechanical potential-energy surface. This algorithm is based on a suitable approximation of an initial full Hessian matrix, either a modified Broyden–Fletcher–Goldfarg–Shanno or a Powell update formula for the location of, respectively, a minimum or a transition state, and the so-called rational function optimization. The latter avoids the Hessian matrix inversion required by a quasi-Newton–Raphson method. Some examples are presented and analyzed. Received: 16 July 2001 / Accepted: 9 October 2001 / Published online: 9 January 2002     

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     

Analysis of the Valley-Ridge inflection points through the partitioning technique of the Hessian eigenvalue equation

The Valley-Ridge inflection (VRI) points are related to the branching of a reaction valley or reaction channel. These points are a special class of points of the potential energy surface (PES). They are also special points of the Valley-Ridge borderline of the PES. The nature of the VRI points and their differences with respect to the other points of the Valley-Ridge borderline is analyzed using the Löwdin’s partitioning technique applied to the eigenvalue equation of the Hessian matrix. Eigenvalues and eigenvectors of the Hessian are better imaginable than the former used adjoint matrix.     

Computation of stationary points via a homotopy method

A homotopy method is presented that locates both minimizers and saddle points of energy functions in an efficient manner. In contrast to other methods, it makes possible the exploration of large parts of potential energy surfaces. Along a homotopy path stationary points of odd and even order occur alternately. A path tracing procedure requiring only gradients and at most one evaluation of the Hessian matrix is given. Test results on a model potential and three MINDO/3 potentials are reported. Received: 6 May 1996 / Accepted: 2 April 1998 / Published online: 23 June 1998     

How good is a Broyden–Fletcher–Goldfarb–Shanno-like update Hessian formula to locate transition structures? Specific reformulation of Broyden–Fletcher–Goldfarb–Shanno for optimizing saddle points

Based on a study of the Broyden–Fletcher–Goldfarb–Shanno (BFGS) update Hessian formula to locate minima on a hypersurface potential energy, we present an updated Hessian formula to locate and optimize saddle points of any order that in some sense preserves the initial structure of the BFGS update formula. The performance and efficiency of this new formula is compared with the previous updated Hessian formulae such as the Powell and MSP ones. We conclude that the proposed update is quite competitive but no more efficient than the normal updates normally used in any optimization of saddle points. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 349–362, 1998     

Remarks on the updated Hessian matrix methods

Optimizing a function with respect to a set of variables using the quasi‐Newton–Raphson method implies updating the Hessian matrix at each iteration. The Broyden–Fletcher–Goldfarb–Shanno update formula is used for minimization and the Murtagh–Sargent–Powell update formula for optimization of first‐order saddle points. Two new formulae are proposed to update the Hessian matrix. One of these formulae is derived using exponential weights and should be used to locate first‐order saddle points. The second formula is a modification of the TS–Broyden–Fletcher–Goldfarb–Shanno update and could used for both minimum and first‐order saddle point optimizations. These two update Hessian matrix formulae present a performance that is the same and in many cases better that the Broyden–Fletcher–Goldfarb–Shanno and Murtagh–Sargent–Powell formulae. © 2003 Wiley Periodicals, Inc. Int J Quantum Chem 94: 324–332, 2003     

Steepest descent reaction path integration using a first-order predictor-corrector method

The theoretical treatment of chemical reactions inevitably includes the integration of reaction pathways. After reactant, transition structure, and product stationary points on the potential energy surface are located, steepest descent reaction path following provides a means for verifying reaction mechanisms. Accurately integrated paths are also needed when evaluating reaction rates using variational transition state theory or reaction path Hamiltonian models. In this work an Euler-based predictor-corrector integrator is presented and tested using one analytic model surface and five chemical reactions. The use of Hessian updating, as a means for reducing the overall computational cost of the reaction path calculation, is also discussed.     

On the degeneracy of the Hessian matrix

A widespread notion in the computational chemistry literature about the Hessian matrix has been revisited, namely, that the Hessian matrix over Cartesian space is sixfold degenerate due to the three translational and three rotational degrees of freedom. It has been shown that this is true only at critical points on the potential energy hypersurface, otherwise the Hessian matrix is only threefold degenerate. The rotational degrees of freedom generally do not cause degeneracy in the Hessian matrix away from critical points.On leave until January 1993 from the Department of General and Analytical Chemistry, Technical University Budapest, Szt. Gellért 4, H-1111 Budapest, Hungary.     

Comparison of methods for finding saddle points without knowledge of the final states

Within the harmonic approximation to transition state theory, the biggest challenge involved in finding the mechanism or rate of transitions is the location of the relevant saddle points on the multidimensional potential energy surface. The saddle point search is particularly challenging when the final state of the transition is not specified. In this article we report on a comparison of several methods for locating saddle points under these conditions and compare, in particular, the well-established rational function optimization (RFO) methods using either exact or approximate Hessians with the more recently proposed minimum mode following methods where only the minimum eigenvalue mode is found, either by the dimer or the Lanczos method. A test problem involving transitions in a seven-atom Pt island on a Pt(111) surface using a simple Morse pairwise potential function is used and the number of degrees of freedom varied by varying the number of movable atoms. In the full system, 175 atoms can move so 525 degrees of freedom need to be optimized to find the saddle points. For testing purposes, we have also restricted the number of movable atoms to 7 and 1. Our results indicate that if attempting to make a map of all relevant saddle points for a large system (as would be necessary when simulating the long time scale evolution of a thermal system) the minimum mode following methods are preferred. The minimum mode following methods are also more efficient when searching for the lowest saddle points in a large system, and if the force can be obtained cheaply. However, if only the lowest saddle points are sought and the calculation of the force is expensive but a good approximation for the Hessian at the starting position of the search can be obtained at low cost, then the RFO approaches employing an approximate Hessian represent the preferred choice. For small and medium sized systems where the force is expensive to calculate, the RFO approaches employing an approximate Hessian is also the more efficient, but when the force and Hessian can be obtained cheaply and only the lowest saddle points are sought the RFO approach using an exact Hessian is the better choice. These conclusions have been reached based on a comparison of the total computational effort needed to find the saddle points and the number of saddle points found for each of the methods. The RFO methods do not perform very well with respect to the latter aspect, but starting the searches further away from the initial minimum or using the hybrid RFO version presented here improves this behavior considerably in most cases.     

Two-surface Monte Carlo with basin hopping: quantum mechanical trajectory and multiple stationary points of water cluster

The efficiency of the two-surface monte carlo (TSMC) method depends on the closeness of the actual potential and the biasing potential used to propagate the system of interest. In this work, it is shown that by combining the basin hopping method with TSMC, the efficiency of the method can be increased by several folds. TSMC with basin hopping is used to generate quantum mechanical trajectory and large number of stationary points of water clusters.     

On the restricted step method coupled with the augmented Hessian for the search of stationary points of any continuous function

In any optimization using the augmented Hessian technique, the step is not restricted to any length. Since the restriction of the step at each iteration is very important in order to achieve good convergence, we present a coupled method such that the augmented Hessian automatically gives both the adequate length of the step and the correct Hessian structure. The method is showed for the minima and saddle points of any order. © 1997 John Wiley & Sons, Inc.     

A study of coronene—coronene association using atom—atom pair potentials

A study of the coronene—coronene association using different interaction potentials based on an atom-atom pair potential proposed by Fraga has been performed. The interaction potentials employed differ in the way the electrostatic and/or dispersion contributions are computed. The influence of both contributions on the geometries predicted for the coronene dimer is discussed in order to analyze the effectiveness of the different interaction potentials. The stationary points found in each interaction energy hypersurface are characterized by calculating the Hessian eigenvalues. Results are discussed in the light of those previously reported for the benzene dimer. Stacked-displaced structures are suggested to be the preferred conformations for the coronene—coronene association. © 1996 John Wiley & Sons, Inc.     

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.     

Calculating rate constants with updated Hessians using variational transition state theory with multidimensional tunneling

Variational transition state theory with multidimensional tunneling (VTST/MT) has been used for calculating the rate constants of reactions. The updated Hessians have been used to reduce the computational costs for both geometry optimization and trajectory following procedures. In this paper, updated Hessians are used to reduce the computational costs while calculating the rate constants applying VTST/MT. Although we found that directly applying the updated Hessians will not generate good vibrational frequencies along the minimum energy path (MEP), however, we can either re-compute the full Hessian matrices at fixed intervals or calculate the Block Hessians, which is constructed by numerical one-side difference for the Hessian elements in the "critical" region and Bofill updating scheme for the rest of the Hessian elements. Due to the numerical instability of the Bofill update method near the saddle point region, we have suggested a simple strategy in which we follow the MEP until certain percentage of the classical barrier height from the barrier top with full Hessians computed and then performing rate constant calculation with the extended MEP using Block Hessians. This strategy results a mean unsigned percentage deviation (MUPD) around 10% with full Hessians computed till the point with 80% classical barrier height for four studied reactions. This proposed strategy is attractive not only it can be implemented as an automatic procedure but also speeds up the VTST/MT calculation via embarrassingly parallelization to a personal computer cluster.     

Intramolecular hydroarylation of alkynes catalyzed by platinum or gold: mechanism and endo selectivity

The cyclization of differently substituted aryl alkynes with PtII or AuI catalysts proceeds by endo-dig pathways. When AgI was used to generate reactive cationic AuI catalysts, 2H-chromenes dimerize to form cyclobutane derivatives by a AgI-catalyzed process. A DFT study on the cyclization mechanism shows a kinetic and thermodynamic preference for 6-endo-dig versus 5-exo-dig cyclizations in PtII-catalyzed processes. Calculations indicate that although Friedel-Crafts and the cyclopropanation processes via metal cyclopropyl carbenes show very similar activation energies, platinum cyclopropyl carbenes are the stationary points with the lowest energy.     

TheRate: Program for ab initio direct dynamics calculations of thermal and vibrational-state-selected rate constants

We introduce TheRate (THEoretical RATEs), a complete application program with a graphical user interface (GUI) for calculating rate constants from first principles. It is based on canonical variational transition-state theory (CVT) augmented by multidimensional semiclassical zero and small curvature tunneling approximations. Conventional transition-state theory (TST) with one-dimensional Wigner or Eckart tunneling corrections is also available. Potential energy information needed for the rate calculations are obtained from ab initio molecular orbital and/or density functional electronic structure theory. Vibrational-state-selected rate constants may be calculated using a diabetic model. TheRate also introduces several technical advancements, namely the focusing technique and energy interpolation procedure. The focusing technique minimizes the number of Hessian calculations required by distributing more Hessian grid points in regions that are critical to the CVT and tunneling calculations and fewer Hessian grid points elsewhere. The energy interpolation procedure allows the use of a computationally less demanding electronic structure theory such as DFT to calculate the Hessians and geometries, while the energetics can be improved by performing a small number of single-point energy calculations along the MEP at a more accurate level of theory. The CH₄+H↔CH₃+H₂ reaction is used as a model to demonstrate usage of the program, and the convergence of the rate constants with respect to the number of electronic structure calculations. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1039–1052, 1998     

Computation of the amide I band of polypeptides and proteins using a partial Hessian approach

A partial Hessian approximation for the computation of the amide I band of polypeptides and proteins is introduced. This approximation exploits the nature of the amide I band, which is largely localized on the carbonyl groups of the backbone amide residues. For a set of model peptides, harmonic frequencies computed from the Hessian comprising only derivatives of the energy with respect to the displacement of the carbon, oxygen, and nitrogen atoms of the backbone amide groups introduce mean absolute errors of 15 and 10 cm(-1) from the full Hessian values at the Hartree-Fock/STO-3G and density functional theory EDF16-31G(*) levels of theory, respectively. Limiting the partial Hessian to include only derivatives with respect to the displacement of the backbone carbon and oxygen atoms yields corresponding errors of 24 and 22 cm(-1). Both approximations reproduce the full Hessian band profiles well with only a small shift to lower wave number. Computationally, the partial Hessian approximation is used in the solution of the coupled perturbed Hartree-Fock/Kohn-Sham equations and the evaluation of the second derivatives of the electron repulsion integrals. The resulting computational savings are substantial and grow with the size of the polypeptide. At the HF/STO-3G level, the partial Hessian calculation for a polypeptide comprising five tryptophan residues takes approximately 10%-15% of the time for the full Hessian calculation. Using the partial Hessian method, the amide I bands of the constituent secondary structure elements of the protein agitoxin 2 (PDB code 1AGT) are calculated, and the amide I band of the full protein estimated.