Search

Algorithms for parallel unconstrained minimization of molecular systems are examined. The overall framework of minimization is the same except for the choice of directions for updating the quasi-Newton Hessian. Ideally these directions are chosen so the updated Hessian gives steps that are same as using the Newton method. Three approaches to determine the directions for updating are presented: the straightforward approach of simply cycling through the Cartesian unit vectors (finite difference), a concurrent set of minimizations, and the Lanczos method. We show the importance of using preconditioning and a multiple secant update in these approaches. For the Lanczos algorithm, an initial set of directions is required to start the method, and a number of possibilities are explored. To test the methods we used the standard 50-dimensional analytic Rosenbrock function. Results are also reported for the histidine dipeptide, the isoleucine tripeptide, and cyclic adenosine monophosphate. All of these systems show a significant speed-up with the number of processors up to about eight processors. [ABSTRACT FROM AUTHOR]

Three new algorithms are presented for determining transition state (TS) structures on the reduced potential energy surface, that is, for problems in which a few important degrees of freedom can be isolated. All three methods use constrained optimization to rapidly find the TS without an initial Hessian evaluation. The algorithms highlight how efficiently the TS can be located on a reduced surface, where the rest of the degrees of freedom are minimized. The first method uses a nonpositive definite quasi-Newton update for the reduced degrees of freedom. The second uses Shepard interpolation to fit the Hessian and starts from a set of points that bound the TS. The third directly uses a finite difference scheme to calculate the reduced degrees of freedom of the Hessian of the entire system, and searches for the TS on the full potential energy surface. All three methods are tested on an epoxide hydrolase cluster, and the ring formations of cyclohexane and cyclobutenone. The results indicate that all the methods are able to converge quite rapidly to the correct TS, but that the finite difference approach is the most efficient. [ABSTRACT FROM AUTHOR]

The number of the potential energy calculations required by the quadratic string method (QSM), and the fast marching method (FMM) is significantly reduced by using Shepard interpolation, with a moving least squares to fit the higher-order derivatives of the potential. The derivatives of the potential are fitted up to fifth order. With an error estimate for the interpolated values, this moving least squares enhanced Shepard interpolation scheme drastically reduces the number of potential energy calculations in FMM, often by up 80%. Fitting up through the highest order tested here (fifth order) gave the best results for all grid spacings. For QSM, using enhanced Shepard interpolation gave slightly better results than using the usual second order approximate, damped Broyden-Fletcher-Goldfarb-Shanno updated Hessian to approximate the surface. To test these methods we examined two analytic potentials, the rotational dihedral potential of alanine dipeptide and the SN2 reaction of methyl chloride with fluoride. [ABSTRACT FROM AUTHOR]

To accurately determine the reaction path and its energetics for enzymatic and solution-phase reactions, we present a sequential sampling and optimization approach that greatly enhances the efficiency of the ab initio quantum mechanics/molecular mechanics minimum free-energy path (QM/MM-MFEP) method. In the QM/MM-MFEP method, the thermodynamics of a complex reaction system is described by the potential of mean force (PMF) surface of the quantum mechanical (QM) subsystem with a small number of degrees of freedom, somewhat like describing a reaction process in the gas phase. The main computational cost of the QM/MM-MFEP method comes from the statistical sampling of conformations of the molecular mechanical (MM) subsystem required for the calculation of the QM PMF and its gradient. In our new sequential sampling and optimization approach, we aim to reduce the amount of MM sampling while still retaining the accuracy of the results by first carrying out MM phase-space sampling and then optimizing the QM subsystem in the fixed-size ensemble of MM conformations. The resulting QM optimized structures are then used to obtain more accurate sampling of the MM subsystem. This process of sequential MM sampling and QM optimization is iterated until convergence. The use of a fixed-size, finite MM conformational ensemble enables the precise evaluation of the QM potential of mean force and its gradient within the ensemble, thus circumventing the challenges associated with statistical averaging and significantly speeding up the convergence of the optimization process. To further improve the accuracy of the QM/MM-MFEP method, the reaction path potential method developed by Lu and Yang [Z. Lu and W. Yang, J. Chem. Phys. 121, 89 (2004)] is employed to describe the QM/MM electrostatic interactions in an approximate yet accurate way with a computational cost that is comparable to classical MM simulations. The new method was successfully applied to two example reaction processes, the classical SN2 reaction of Cl-+CH3Cl in solution and the second proton transfer step of the reaction catalyzed by the enzyme 4-oxalocrotonate tautomerase. The activation free energies calculated with this new sequential sampling and optimization approach to the QM/MM-MFEP method agree well with results from other simulation approaches such as the umbrella sampling technique with direct QM/MM dynamics sampling, demonstrating the accuracy of the iterative QM/MM-MFEP method. [ABSTRACT FROM AUTHOR]

A new method, referred to as the sequential quadratic programming method, is presented for determining minimum energy paths. The 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 can efficiently determine the reaction mechanism, 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 Müller-Brown potential, amide hydrolysis, and an 89 atom cluster taken from the active site of 4-oxalocrotonate tautomerase for the reaction which catalyzes 2-oxo-4-hexenedioate to the intermediate 2-hydroxy-2,4-hexadienedioate. [ABSTRACT FROM AUTHOR]

The diagonally implicit Runge-Kutta framework 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 nonembedded, diagonally implicit Runge-Kutta methods available from the numerical analysis literature and these are reviewed for orders two, three, and four. 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: the Müller-Brown [Theor. Chem. Acta 53, 75 (1979)] potential and two gas phase chemical reactions. The results show that many of the methods are capable of integrating efficiently while reliably keeping the error bound within the desired tolerance. This allows the reaction path to be determined through automatic integration by only specifying the desired accuracy and transition state. [ABSTRACT FROM AUTHOR]

We present the use of an optimal combined explicit-implicit method for following the reaction path to high accuracy. This is in contrast to most purely implicit reaction path integration algorithms, which are only efficient on stiff ordinary differential equations. The defining equation for the reaction path is considered to be stiff, however, we show here that the reaction path is not uniformly stiff and instead is only stiff near stationary points. The optimal algorithm developed in this work is a combination of explicit and implicit methods with a simple criterion to switch between the two. Using three different chemical reactions, we combine and compare three different integration methods: the implicit trapezoidal method [C. Gonzalez and H. Schlegel, J. Chem. Phys. 90, 2154 (1989)], an explicit stabilized third order algorithm [A. A. Medovikov, BIT 38, 372 (1998)] implemented in the code DUMKA3 and the traditional explicit fourth order Runge-Kutta method written in the code RKSUITE. The results for high accuracy show that when the implicit trapezoidal method is combined with either explicit method the number of energy and gradient calculations can potentially be reduced by almost a half compared with integrating either method alone. Finally, to explain the improvements of the combined method we expand on the concepts of stability and stiffness and relate them to the efficiency of integration methods. [ABSTRACT FROM AUTHOR]

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 Müller-Brown potential, a seven-atom Lennard-Jones cluster, and the enolation of acetaldehyde to vinyl alcohol. [ABSTRACT FROM AUTHOR]

Wepresent a quantum mechanics/molecular mechanics (QM/MM) methodto evaluate the partial charges of amino acid residues for use inMM potentials based on their protein environment. For each residueof interest, the nearby residues are included in the QM system whilethe rest of the protein is treated at the MM level of theory. Aftera short structural optimization, the partial charges of the centralresidue are fit to the electrostatic potential using the restrainedelectrostatic potential (RESP) method. The resulting charges and electrostaticpotential account for the individual environment of the residue, althoughthey lack the transferable nature of library partial charges. To evaluatethe quality of the QM/MM RESP charges, thermodynamic integration isused to measure the pKashift of the asparticacid residues in three different proteins, turkey egg lysozyme, beta-cryptogein,and Thioredoxin. Compared to the AMBER ff99SB library values, theQM/MM RESP charges show better agreement between the calculated andexperimental pKavalues for almost allof the residues considered. [ABSTRACT FROM AUTHOR]

A fast, empirical method, Mut-pKa, is presented for predicting the p Ka values of ionizable residues in proteins based on mutation. The method compares the effect of mutating each residue that may act as a hydrogen bond donor or acceptor for the ionizable residue. The energetic effect of each type of mutation, along with a desolvation measure and the overall background charge, is fit against p Ka data for histidine and carboxyl residues. A total of 214 residues from 35 different proteins were used in the dataset. Using 11 parameters for each type of ionizable residue, a root mean squared error (RMSE) of 0.78 and 1.12 pH units were obtained for carboxyl and histidines residues, respectively, using leave one out cross validation (LOOCV). The results were particularly promising for buried residues, which had RMSE values of 0.99 and 1.13 for carboxyl and histidine residues, respectively. A number of desolvation measures were tested. The simplest measure, the number of atoms surrounding the residue, was found to work best. The effect of using dynamics was also studied using short molecular dynamics runs, followed by minimization of the structures. Mut-pKa has significantly fewer parameters than, but similar performance to, other empirical methods. Because of this and the LOOCV results, we believe the model is robust and that overfitting is not a problem. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011 [ABSTRACT FROM AUTHOR]

Recognizing the limits of trying to achieve chemical accuracy for p K calculations with a purely electrostatic model, we include empirical corrections into the Poisson-Boltzmann solver macroscopic electrostatics with atomic detail (Bashford, Biochemistry 1990;29:10219-10225), to improve the reliability and accuracy of the model. The total number of parameters is kept to a minimum to maximize the robustness of the model for compounds outside of the fitting dataset. The parameters are based on: (a) the electrostatic interaction between functional groups close to the titratable site, (b) the electrostatic work required to desolvate the residue, and (c) the site-to-site interactions. These interactions are straightforward to calculate once the electrostatic field has been solved for each residue using the linearized Poisson-Boltzmann equation and are assumed to be linearly related to the intrinsic p K. Two hundred and eighty-six residues from 30 proteins are used to determine the empirical parameters, which result in a root mean square error (RMSE) of 0.70 for the entire set. Eight proteins with 46 experimentally known values were excluded from the parameterization to test the model. This test set had a RMSE of 1.08. We show that the parameterized model improves the results over other models, although like other models the error is strongly correlated with the degree to which a residue is buried. The parameters themselves indicate that local effects are most important for determining the p K, whereas site-to-site interactions are found to be less significant. Proteins 2011; © 2011 Wiley-Liss, Inc. [ABSTRACT FROM AUTHOR]

A fast, empirical method, Mut-pKa, is presented for predicting the pKa values of ionizable residues in proteins based on mutation. The method compares the effect of mutating each residue that may act as a hydrogen bond donor or acceptor for the ionizable residue. The energetic effect of each type of mutation, along with a desolvation measure and the overall background charge, is fit against pKa data for histidine and carboxyl residues. A total of 214 residues from 35 different proteins were used in the dataset. Using 11 parameters for each type of ionizable residue, a root mean squared error (RMSE) of 0.78 and 1.12 pH units were obtained for carboxyl and histidines residues, respectively, using leave one out cross validation (LOOCV). The results were particularly promising for buried residues, which had RMSE values of 0.99 and 1.13 for carboxyl and histidine residues, respectively. A number of desolvation measures were tested. The simplest measure, the number of atoms surrounding the residue, was found to work best. The effect of using dynamics was also studied using short molecular dynamics runs, followed by minimization of the structures. Mut-pKa has significantly fewer parameters than, but similar performance to, other empirical methods. Because of this and the LOOCV results, we believe the model is robust and that overfitting is not a problem., (Copyright © 2011 Wiley Periodicals, Inc.)

We investigate the effect of systematically applying molecular dynamics (MD) and quantum mechanics/molecular mechanics (QM/MM) to docked poses in an attempt to improve the correspondence between theoretical prediction and experimental observation. The proposed scheme involves running a short time scale MD simulation on a docked ligand pose (and any known structurally important crystal structure waters in the active site), followed by QM/MM minimization. Both of these steps are relatively fast for moderately sized ligands; longer time scale MD involving the protein is not found to improve the results. The final binding energy is given in terms of the QM/MM total energy, a van der Waals correction, and a term to account for desolvation effects. This methodology is first tested with a trypsin inhibitor, for which we establish the importance of running MD before reoptimizing with QM/MM. The method is then applied to cytochrome c peroxidase using a set of binders and decoys. In this example, the proposed methodology affords much better discrimination between binders and decoys than the traditional docking approach used. For both systems presented, application of this protocol results in a significantly better energetic ranking and a smaller root mean squared deviation from known crystallographic ligand poses. This work highlights the importance of including polarization effects through QM/MM and of sampling with MD to refine a set of initial docked poses.