Overcoming free energy barriers using unconstrained molecular dynamics simulations

Association of unconstrained molecular dynamics (MD) and the formalisms of thermodynamic integration and average force [Darve and Pohorille, J. Chem. Phys. 115, 9169 (2001)] have been employed to determine potentials of mean force. When implemented in a general MD code, the additional computational effort, compared to other standard, unconstrained simulations, is marginal. The force acting along a chosen reaction coordinate xi is estimated from the individual forces exerted on the chemical system and accumulated as the simulation progresses. The estimated free energy derivative computed for small intervals of xi is canceled by an adaptive bias to overcome the barriers of the free energy landscape. Evolution of the system along the reaction coordinate is, thus, limited by its sole self-diffusion properties. The illustrative examples of the reversible unfolding of deca-L-alanine, the association of acetate and guanidinium ions in water, the dimerization of methane in water, and its transfer across the water liquid-vapor interface are examined to probe the efficiency of the method.     

Coupling of replica exchange simulations to a non-Boltzmann structure reservoir

Computing converged ensemble properties remains challenging for large biomolecules. Replica exchange molecular dynamics (REMD) can significantly increase the efficiency of conformational sampling by using high temperatures to escape kinetic traps. Several groups, including ours, introduced the idea of coupling replica exchange to a pre-converged, Boltzmann-populated reservoir, usually at a temperature higher than that of the highest temperature replica. This procedure reduces computational cost because the long simulation times needed for extensive sampling are only carried out for a single temperature. However, a weakness of the approach is that the Boltzmann-weighted reservoir can still be difficult to generate. We now present the idea of employing a non-Boltzmann reservoir, whose structures can be generated through more efficient conformational sampling methods. We demonstrate that the approach is rigorous and derive a correct statistical mechanical exchange criterion between the reservoir and the replicas that drives Boltzmann-weighted probabilities for the replicas. We test this approach on the trpzip2 peptide and demonstrate that the resulting thermal stability profile is essentially indistinguishable from that obtained using very long (>100 ns) standard REMD simulations. The convergence of this reservoir-aided REMD is significantly faster than for regular REMD. Furthermore, we demonstrate that modification of the exchange criterion is essential; REMD simulations using a standard exchange function with the non-Boltzmann reservoir produced incorrect results.     

Efficiency of alchemical free energy simulations. II. Improvements for thermodynamic integration

We attempt to optimize the efficiency of thermodynamic integration, as defined by the minimal number of unphysical intermediate states required for the computation of accurate and precise free energy differences. The suitability of various numerical quadrature methods is tested. In particular, we compare the trapezoidal rule, Simpson's rule, Gauss-Legendre, Gauss-Kronrod-Patterson, and Clenshaw-Curtis integration, as well as integration based on a cubic spline approximation of the integrand. We find that Simpson's rule and spline integration are already significantly more efficient that the trapezoidal rule, i.e., correct free energy differences can be obtained using fewer λ-states. We demonstrate that Simpson's rule can be used advantageously with nonequidistant values of the abscissa, which increases the flexibility of the method. Efficiency is enhanced even further if higher order methods, such as Gauss-Legendre, Gauss-Kronrod-Patterson, or Clenshaw-Curtis integration, are used; no more than seven λ-states, which in the case of Clenshaw-Curtis integration include the physical end states, were required for accurate results in all test problems studied. Thus, the performance of thermodynamic integration can equal that of Bennett's acceptance ratio method. We also show, however, that the high efficiency found here relies on the particular functional form of the soft-core potential used; overall, thermodynamic integration is more susceptible to the details of the hybrid Hamiltonian used than Bennett's acceptance ratio method. Therefore, we recommend Bennett's acceptance ratio method as the most robust method to compute alchemical free energy differences; nevertheless, scenarios when thermodynamic integration may be preferable are discussed.     

Overcoming sustainability barriers within the chemical industry

1. Download : Download high-res image (109KB)
2. Download : Download full-size image

     

On searching neighbors in computer simulations of macromolecular systems

Various algorithms for evaluating nonbonded interactions in molecular dynamics (MD ) simulations of macromolecular systems are considered, and a combination of two techniques using a space grid for finding neighbor atoms is proposed. The application of grid search techniques to nonrectangular periodic systems is discussed. Finally, the computing time required by different algorithms is compared on a Cray-1 vector processing computer as well as on a CDC Cyber 170/760 sequential computer. Neighbor list techniques turn out to be faster than the grid search techniques for the systems considered here. However, storage requirements may exclude the use of neighbor list techniques for large systems. Finally, MD of complex macromolecular systems turned out to be about 4–14 times faster on a Cray-1 than on a Cyber 170/760.     

Gradient SHAKE: An improved method for constrained energy minimization in macromolecular simulations

Current macromolecular energy minimization algorithms become inefficient and prone to failure when bond length constraints are imposed. They are required to relieve steric stresses in biomolecules prior to a molecular dynamics simulation. Unfortunately, the latter often require constraints, leading to difficulties in initiating trajectories from unconstrained energy minima. This difficulty was overcome by requiring that the components of the energy gradient vanish along the constrained bonds. The modified energy minimization algorithm converges to a lower energy in a fewer number of iterations and is more robust than current implementations. The method has been successfully applied to the Dickerson DNA dodecamer, CGCGAATTCGCG. © 1995 John Wiley & Sons, Inc.     

Grain-boundary free energy via thermodynamic integration

In a previous publication by Lusk and Beale [Phys. Rev. E 69, 026117 (2004)], fluctuating cell (FC) theory was used to estimate the free energy of symmetric tilt grain boundaries in an assembly of nearly hard disks. The FC method is much faster than the more traditional thermodynamic integration, but the accuracy of the algorithm has not been assessed in association with persistent defect structures. This motivated the present work wherein the FC free energies are compared directly with the data obtained via thermodynamic integration from an Einstein crystal to an assembly of hard disks. This comparison is made over the range of possible misorientations for symmetric tilt boundaries and indicates that the FC method gives quantitatively accurate estimates for grain-boundary free energy. We also demsonstrate that the FC approximation is quantitatively accurate at determining the free-energy contribution of each particle whether in the bulk or the grain boundary. The FC calculation is about two orders of magnitude faster than a full thermodynamic integration. This approach may offer a numerically efficient means of estimating the free energy of persistent defect structures to greater accuracy than is afforded by the quasiharmonic and local harmonic approximations.     

Soft-core potentials in thermodynamic integration: comparing one- and two-step transformations

Molecular dynamics-based free energy calculations allow the determination of a variety of thermodynamic quantities from computer simulations of small molecules. Thermodynamic integration (TI) calculations can suffer from instabilities during the creation or annihilation of particles. This "singularity" problem can be addressed with "soft-core" potential functions which keep pairwise interaction energies finite for all configurations and provide smooth free energy curves. "One-step" transformations, in which electrostatic and van der Waals forces are simultaneously modified, can be simpler and less expensive than "two-step" transformations in which these properties are changed in separate calculations. Here, we study solvation free energies for molecules of different hydrophobicity using both models. We provide recommended values for the two parameters α(LJ) and β(C) controlling the behavior of the soft-core Lennard-Jones and Coulomb potentials and compare one- and two-step transformations with regard to their suitability for numerical integration. For many types of transformations, the one-step procedure offers a convenient and accurate approach to free energy estimates.     

Converting solution macromolecular thermodynamic properties into gas-phase mass spectrometry observations

Solution dissociation constants and changes in free energies associated with ligand binding to proteins have been measured in the gas phase using mass spectrometry.     

The self-referential method combined with thermodynamic integration

The self-referential method [M. B. Sweatman, Phys. Rev. E 72, 016711 (2005)] for calculating the free energy of crystalline solids via molecular simulation is combined with thermodynamic integration to produce a technique that is convenient and efficient. Results are presented for the chemical potential of hard sphere and Lennard-Jones face centered cubic crystals that agree well with this previous work. For the small system sizes studied, this technique is about 100 times more efficient than the parameter hopping technique used previously.     

Measurement of nonequilibrium entropy from space-time thermodynamic integration

The entropy of a system transiently driven out of equilibrium by a time-inhomogeneous stochastic dynamics is first expressed as a transient response function generalizing the nonlinear Kawasaki-Crooks response. This function is then reformulated into three statistical averages defined over ensembles of nonequilibrium trajectories. The first average corresponds to a space-time thermodynamic perturbation relation, while the two following ones correspond to space-time thermodynamic integration relations. Provided that trajectories are initiated starting from a distribution of states that is analytically known, the ensemble averages are computationally amenable to Markov chain Monte Carlo methods. The relevance of importance sampling in path ensembles is confirmed in practice by computing the nonequilibrium entropy of a driven toy system. We finally study a situation where the dynamics produces entropy. In this case, we observe that space-time thermodynamic integration still yields converged estimates, while space-time thermodynamic perturbation turns out to converge very slowly.     

On the enhanced sampling over energy barriers in molecular dynamics simulations

We present here calculations of free energies of multidimensional systems using an efficient sampling method. The method uses a transformed potential energy surface, which allows an efficient sampling of both low and high energy spaces and accelerates transitions over barriers. It allows efficient sampling of the configuration space over and only over the desired energy range(s). It does not require predetermined or selected reaction coordinate(s). We apply this method to study the dynamics of slow barrier crossing processes in a disaccharide and a dipeptide system.     

Pseudorotation barriers of biological oxyphosphoranes: a challenge for simulations of ribozyme catalysis

Pseudorotation reactions of biologically relevant oxyphosphoranes were studied by using density functional and continuum solvation methods. A series of 16 pseudorotation reactions involving acyclic and cyclic oxyphosphoranes in neutral and monoanionic (singly deprotonated) forms were studied, in addition to pseudorotation of PF5. The effect of solvent was treated by using three different solvation models for comparison. The barriers to pseudorotation ranged from 1.5 to 8.1 kcal mol(-1) and were influenced systematically by charge state, apicophilicity of ligands, intramolecular hydrogen bonding, cyclic structure and solvation. Barriers to pseudorotation for monoanionic phosphoranes occur with the anionic oxo ligand as the pivotal atom, and are generally lower than for neutral phosphoranes. The OCH3 groups were observed to be more apicophilic than OH groups, and hence pseudorotations that involve axial OCH3/equatorial OH exchange had higher reaction and activation free energy values. Solvent generally lowered barriers relative to the gas-phase reactions. These results, together with isotope 18O exchange experiments, support the assertion that dianionic phosphoranes are not sufficiently long-lived to undergo pseudorotation. Comparison of the density functional results with those from several semiempirical quantum models highlight a challenge for new-generation hybrid quantum mechanical/molecular mechanical potentials for non-enzymatic and enzymatic phosphoryl transfer reactions: the reliable modeling of pseudorotation processes.     

镜像配基筛选:克服“酶屏障”的多肽、寡核苷酸药物研发技术

生物体内存在的酶屏障是功能性多肽和寡核苷酸成药的重要瓶颈.以镜像噬菌体展示技术(mirror-image phage display)和镜像适配体筛选技术(Spiegelmer technology)为代表的镜像配基筛选技术能有效地克服这一屏障,获得对靶点具有特异识别能力且在生理条件下高度稳定的D型多肽和L型寡核苷酸配基,有望成为多肽或寡核苷酸类特异性药物研发的有利工具.本文主要综述镜像配基筛选技术的原理及其在药学领域的研究进展.     

Isolated pentagon rule in buckybowls: a computational study on thermodynamic stabilities and bowl-to-bowl inversion barriers

Computations on a series of sumanene and pinakene isomers reveal that the isolated pentagon species are found to be more stable than non-isolated pentagon structures. However, the present study indicates that there is no correlation between the bowl-to-bowl inversion barrier and the relative stability. Unexpectedly, the least stable isomer has the lowest bowl-to-bowl inversion barrier in pinakene isomers. Thus, curvature cannot be taken as a measure to evaluate the stability of buckybowls. The relative energy ordering and HOMO-LUMO gap correlate well in the series of pinakene isomers considered.     

Overcoming barriers in trajectory space: mechanism and kinetics of rare events via Wang-Landau enhanced transition path sampling

Within the framework of transition path sampling (TPS), activation energies can be computed as path ensemble averages without a priori information about the reaction mechanism [C. Dellago and P. G. Bolhuis, Mol. Simul. 30, 795 (2004)]. Activation energies computed for different conditions can then be used to determine by numerical integration the rate constant for a system of interest from the rate constant known for a reference system. However, in systems with complex potential energy surfaces, multiple reaction pathways may exist making ergodic sampling of trajectory space difficult. Here, we present a combination of TPS with the Wang-Landau (WL) flat-histogram algorithm for an efficient sampling of the transition path ensemble. This method, denoted by WL-TPS, has the advantage that from one single simulation, activation energies at different temperatures can be determined even for systems with multiple reaction mechanisms. The proposed methodology for rate constant calculations does not require the knowledge of the reaction coordinate and is generally applicable to Arrhenius and non-Arrhenius processes. We illustrate the applicability of this technique by studying a two-dimensional toy system consisting of a triatomic molecule immersed in a fluid of repulsive soft disks. We also provide an expression for the calculation of activation volumes from path averages such that the pressure dependence of the rate constant can be obtained by numerical integration.     

Free energy difference calculations by thermodynamic integration: Difficulties in obtaining a precise value

Free energy difference calculations have been performed by the “slow growth” method of thermodynamic integration of the AMBER 3.0 molecular dynamics program for the mutation of a conformationally restricted threonine dipeptide, N-acetyl threonyl-N-methylamide, to the corresponding alanyl dipeptide. By varying the total simulation length, it has been determined that precise free energy values are obtained only for simulations of greater than 100 ps total simulation time length. By varying the starting configurations for simulations of the same length, it has been determined that averaging the free energies obtained from shorter simulations may not give precise answers. Possible reasons for this behavior are discussed.     

Direct dynamics simulations using Hessian-based predictor-corrector integration algorithms

In previous research [J. Chem. Phys. 111, 3800 (1999)] a Hessian-based integration algorithm was derived for performing direct dynamics simulations. In the work presented here, improvements to this algorithm are described. The algorithm has a predictor step based on a local second-order Taylor expansion of the potential in Cartesian coordinates, within a trust radius, and a fifth-order correction to this predicted trajectory. The current algorithm determines the predicted trajectory in Cartesian coordinates, instead of the instantaneous normal mode coordinates used previously, to ensure angular momentum conservation. For the previous algorithm the corrected step was evaluated in rotated Cartesian coordinates. Since the local potential expanded in Cartesian coordinates is not invariant to rotation, the constants of motion are not necessarily conserved during the corrector step. An approximate correction to this shortcoming was made by projecting translation and rotation out of the rotated coordinates. For the current algorithm unrotated Cartesian coordinates are used for the corrected step to assure the constants of motion are conserved. An algorithm is proposed for updating the trust radius to enhance the accuracy and efficiency of the numerical integration. This modified Hessian-based integration algorithm, with its new components, has been implemented into the VENUS/NWChem software package and compared with the velocity-Verlet algorithm for the H(2)CO-->H(2)+CO, O(3)+C(3)H(6), and F(-)+CH(3)OOH chemical reactions.     

Molecular dynamics and thermodynamic simulations of segregation phenomena in binary metal nanoparticles

Journal of Thermal Analysis and Calorimetry - An approach combining atomistic molecular dynamics (MD) and thermodynamic simulations has been applied to predict the distribution of components in...