MINIMUM ENERGY AND STEEPEST DESCENT PATH ALGORITHMS FOR QM/MM APPLICATONS

Abstract

A number of new methods are presented to determine the reaction path both for chemical systems where the transition state(TS) is known and for the more complicated case when only the endpoints are available. To determine the minimum energy path(MEP) two algorithms were developed.The first MEP method is a quadratic string method (QSM) which is based on a multiobjective optimization framework. In the method, each point on the MEP is integrated in the descent direction perpendicular to path. Each local integration is done on an approximate quadratic surface with an 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 method is demonstrated for the Muller Brown potential, a 7-atom Lennard-Jones cluster and the enolation of acetaldehyde to vinyl alcohol.The second MEP method is referred to as the Sequential Quadratic Programming Method (SQPM). This method is based on minimizing the points representing the path in the subspace perpendicular to the tangent of the path while using a penalty term to prevent kinks from forming. Rather than taking one full step, the minimization is divided into a number of sequential steps on an approximate quadratic surface. The resulting method is shown to be capable of super-linear convergence. However, the emphasis of the algorithm is on its robustness and its ability to determine the reaction mechanism efficiently, from which transition state can be easily identified and refined with other methods. To improve the resolution of the path close to the transition state, points are clustered close to this region with a reparametrization scheme. The usefulness of the algorithm is demonstrated for the Mu$ller Brown potential, amide hydrolysis and an 89 atom cluster taken from the active site of 4-Oxalocrotonate tautomerase (4-OT) for the reaction which catalyzes 2-oxo-4-hexenedioate to the intermediate 2-hydroxy-2,4-hexadienedioate.When the TS is known we present two methods for integrating the steepest descent path (SDP). Also the concepts of stability and stiffness are elaborated upon. The first SDP method is an optimally combined explicit-implicit method for following the reaction path to high accuracy. Although the SDP is generally considered to be a stiff ODE, it is shown that the reaction path is not uniformly stiff and instead is only stiff near stationary points. The optimal algorithm is developed by combining the explicit and implicit methods with a simple criterion, based on the stiffness, to switch between the two. Using two different methods an algorithm is developed to efficiently integrate the SDP. This method is tested on a number of small molecules.The final method given is based on the diagonally implicit Runge-Kutta framework, which is shown to be a general form for constructing stable, efficient steepest descent reaction path integrators, of any order. With this framework tolerance driven, adaptive step-size methods can be constructed by embedding methods to obtain error estimates of each step without additional computational cost. There are many embedded and non-embedded, diagonally implicit Runge-Kutta methods available from the numerical analysis literature and these are reviewed for orders 2,3 and 4. New embedded methods are also developed which are tailored to the application of reaction path following. All integrators are summarized and compared for three systems.