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     

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     

Finding transition states using reduced potential-energy surfaces

 We propose a method to locate saddle points that is based on the interplay between the driving coordinate and the restricted quasi-Newton algorithm. The method locates the transition state using a reduced potential-energy surface. The reduced potential-energy surface is characterized by the set of driving coordinates. The proposed algorithm starts at a point on the surface that is slightly perturbed from either reactant or product and, in principle, converges to the transition state. Finally we give a special type of update Hessian matrix formula that should be applied in optimizations carried out on reduced potential-energy surfaces. Received: 29 September 2000 / Accepted: 3 January 2001 / Published online: 3 April 2001     

On the determination of the minimum on the crossing seam of two potential energy surfaces

An improved gradient-based algorithm is presented for the determination of the minimum energy point on the crossing seam hypersurface between two arbitrary potential energy hypersurfaces. The Hessian matrix is updated employing the gradient information. The method is demonstrated in a study of some representative cases including charge-transfer states of a typical molecular-device molecule (a rigid spiro π – σ – π molecular cation) with, as well as without, an external electric field.     

An analytical computation of Christoffel symbols for reaction coordinate and trajectory treatments under internal coordinates

We examine the Hessian matrix of the potential energy under internal coordinates. We report all Christoffel symbols which exist for molecules if we use the known coordinates such as bond distances, bond angles, torsion angles, and out-of-plane angles. We use as an example triatomic HCN in an extended geometry.     

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     

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.     

Technique for incorporating the density functional Hessian into the geometry optimization of biomolecules, solvated molecules, and large floppy molecules

Traditional geometry optimization methods require the gradient of the potential surface, together with a Hessian which is often approximated. Approximation of the Hessian causes difficulties for large, floppy molecules, increasing the number of steps required to reach the minimum. In this article, the costly evaluation of the exact Hessian is avoided by expanding the density functional to second order in both the nuclear and electronic variables, and then searching for the minimum of the quadratic functional. The quadratic search involves the simultaneous determination of both the geometry step and the associated change in the electron density matrix. Trial calculations on Taxol indicate that the cost of the quadratic search is comparable to the cost of the density functional energy plus gradient. While this procedure circumvents the bottleneck coupled-perturbed step in the evaluation of the full Hessian, the second derivatives of the electron-repulsion integrals are still required for atomic-orbital-based calculations, and they are presently more expensive than the energy plus gradient. Hence, we anticipate that the quadratic optimizer will initially find application in fields in which existing optimizers breakdown or are inefficient, particularly biochemistry and solvation chemistry.     

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     

Quadratic string method for determining the minimum-energy path based on multiobjective optimization

Based on a multiobjective optimization framework, we develop a new quadratic string method for finding the minimum-energy path. In the method, each point on the minimum-energy path is minimized by integration in the descent direction perpendicular to path. Each local integration is done on a quadratic surface approximated by a damped Broyden-Fletcher-Goldfarb-Shanno updated Hessian, allowing the algorithm to take many steps between energy and gradient calls. The integration is performed with an adaptive step-size solver, which is restricted in length to the trust radius of the approximate Hessian. The full algorithm is shown to be capable of practical superlinear convergence, in contrast to the linear convergence of other methods. The method also eliminates the need for predetermining such parameters as step size and spring constants, and is applicable to reactions with multiple barriers. The effectiveness of this method is demonstrated for the Muller-Brown potential, a seven-atom Lennard-Jones cluster, and the enolation of acetaldehyde to vinyl alcohol.     

Improved RGF method to find saddle points

The predictor-corrector method for following a reduced gradient (RGF) to determine saddle points [Quapp, W. et al., J Comput Chem 1998, 19, 1087] is further accelerated by a modification allowing an implied corrector step per predictor but almost without additional costs. The stability and robustness of the RGF method are improved, and the new version in addition reduces the number of gradient and Hessian calculations.     

Extremely localized molecular orbitals (ELMO): a non-orthogonal Hartree-Fock method

A new optimization method for extremely localized molecular orbitals (ELMO) is derived in a non-orthogonal formalism. The method is based on a quasi Newton-Raphson algorithm in which an approximate diagonal-blocked Hessian matrix is calculated through the Fock matrix. The Hessian matrix inverse is updated at each iteration by a variable metric updating procedure to account for the intrinsically small coupling between the orbitals. The updated orbitals are obtained with approximately n ² operations. No n ³ processes such as matrix diagonalization, matrix multiplication or orbital orthogonalization are employed. The use of localized orbitals allows for the creation of high-quality initial “guess” orbitals from optimized molecular orbitals of small systems and thus reduces the number of iterations to converge. The delocalization effects are included by a Jacobi correction (JC) which allows the accurate calculation of the total energy with a limited number of operations. This extension, referred to as ELMO(JC), is a variational method that reproduces the Hartree-Fock (HF) energy with an error of less than 2 kcal/mol for a reduced total cost compared to standard HF methods. The small number of variables, even for a very large system, and the limited number of operations potentially makes ELMO a method of choice to study large systems. Received: 30 December 1996 / Accepted: 5 June 1997     

Partial hessian fitting for determining force constant parameters in molecular mechanics

We present a new protocol for deriving force constant parameters that are used in molecular mechanics (MM) force fields to describe the bond‐stretching, angle‐bending, and dihedral terms. A 3 × 3 partial matrix is chosen from the MM Hessian matrix in Cartesian coordinates according to a simple rule and made as close as possible to the corresponding partial Hessian matrix computed using quantum mechanics (QM). This partial Hessian fitting (PHF) is done analytically and thus rapidly in a least‐squares sense, yielding force constant parameters as the output. We herein apply this approach to derive force constant parameters for the AMBER‐type energy expression. Test calculations on several different molecules show good performance of the PHF parameter sets in terms of how well they can reproduce QM‐calculated frequencies. When soft bonds are involved in the target molecule as in the case of secondary building units of metal‐organic frameworks, the MM‐optimized geometry sometimes deviates significantly from the QM‐optimized one. We show that this problem is rectified effectively by use of a simple procedure called Katachi that modifies the equilibrium bond distances and angles in bond‐stretching and angle‐bending terms. © 2016 Wiley Periodicals, Inc.     

A conceptual realization of stationary points using the Hessian matrix

In this study we use surface-fitting equations to generate the energy and its first derivative in terms of two torsional angles in the methanediol model, using GAUSSIAN 88 at the 4-31G level of approximation. The Hessian matrix is further used to locate the stationary points, and a Gaussian fit of the absolute values of the sum of the eigenvalues of the Hessian is used in order to generate a surface in which all the stationary points are identified.     

A survey of optimization procedures for stable structures and transition states

We examine a variety of methods for obtaining the stable geometry of molecules and the transition states of simple systems and summarize some of our findings. We find the most efficient methods for optimizing structure to be those based on calculated gradients and estimated second derivative (Hessian) matrices, the later obtained either from the Broyden–Fletcher–Goldfarb–Shanno (BFGS ) quasi-Newton update method or from approximations to the coupled perturbed Hartree–Fock method. For uncovering transition states we find particularly useful a variety of the augmented Hessian theory used to uncover regions of the potential energy hypersurface with one and only one negative eigenvalue of the Hessian matrix characterizing the catchment region of the transition state. Once this region is found we minimize the norm of the gradient vector to catch the nearest extreme point of the surface. Examples of these procedures are given.     

Multipath variational transition state theory: rate constant of the 1,4-hydrogen shift isomerization of the 2-cyclohexylethyl radical

We propose a new formulation of variational transition state theory called multipath variational transition state theory (MP-VTST). We employ this new formulation to calculate the forward and reverse thermal rate constant of the 1,4-hydrogen shift isomerization of the 2-cyclohexylethyl radical in the gas phase. First, we find and optimize all the local-minimum-energy structures of the reaction, product, and transition state. Then, for the lowest-energy transition state structures, we calculate the reaction path by using multiconfiguration Shepard interpolation (MSCI) method to represent the potential energy surface, and, from this representation, we also calculate the ground-state vibrationally adiabatic potential energy curve, the reaction-path curvature vector, and the generalized free energy of activation profile. With this information, the path-averaged generalized transmission coefficients <γ> are evaluated. Then, thermal rate constant containing the multiple-structure anharmonicity and torsional anharmonicity effects is calculated using multistructural transition state theory (MS-TST). The final MP-VTST thermal rate constant is obtained by multiplying k(MS-T)(MS-TST) by <γ>. In these calculations, the M06 density functional is utilized to compute the energy, gradient, and Hessian at the Shepard points, and the M06-2X density functional is used to obtain the structures (conformers) of the reactant, product, and the saddle point for computing the multistructural anharmonicity factors.     

Comment “On the quadratic reaction path evaluated in a reduced potential energy surface model and the problem to locate transition states” [by J. M. Anglada,E. Besalú, J. M. Bofill,and R. Crehuet,J Comput Chem 2001, 22, 4, 387–406]

The difference is explained between steepest ascent and following a reduced gradient (distinguished coordinate method) for the location of saddle points. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 537–540, 2001     

Analysis of the updated Hessian matrices for locating transition structures

We present an analysis of the behavior of different updating Hessian formulas when they are used for the location and optimization of transition structures. The analysis is based on the number of iterations, the minimum of the weighted Euclidean matrix norm, and first-order perturbation theory applied to each type of Hessian correction. Finally, we give a derivation of a family of updated Hessians from the variational method proposed by Greenstadt. We conclude that the proposed family of updated Hessians is useful for the optimization of transition structures. © 1995 John Wiley & Sons, Inc.