Efficient and precise solvation free energies via alchemical adiabatic molecular dynamics

A new molecular dynamics method for calculating free energies associated with transformations of the thermodynamic state or chemical composition of a system (also known as alchemical transformations) is presented. The new method extends the adiabatic dynamics approach recently introduced by Rosso et al. [J. Chem. Phys. 116, 4389 (2002)] and is based on the use of an additional degree of freedom, lambda, that is used as a switching parameter between the potential energy functions that characterize the two states. In the new method, the coupling parameter lambda is introduced as a fictitious dynamical variable in the Hamiltonian, and a system of switching functions is employed that leads to a barrier in the lambda free energy profile between the relevant thermodynamic end points. The presence of such a barrier, therefore, enhances sampling in the end point (lambda = 0 and lambda = 1) regions which are most important for computing relevant free energy differences. In order to ensure efficient barrier crossing, a high temperature T(lambda) is assigned to lambda and a fictitious mass m(lambda) is introduced as a means of creating an adiabatic separation between lambda and the rest of the system. Under these conditions, it is shown that the lambda free energy profile can be directly computed from the adiabatic probability distribution function of lambda without any postprocessing or unbiasing of the output data. The new method is illustrated on two model problems and in the calculation of the solvation free energy of amino acid side-chain analogs in TIP3P water. Comparisons to previous work using thermodynamic integration and free energy perturbation show that the new lambda adiabatic free energy dynamics method results in very precise free energy calculations using significantly shorter trajectories.     

Self-healing umbrella sampling: a non-equilibrium approach for quantitative free energy calculations

We propose a new approach for the umbrella sampling method in molecular dynamics simulations of complex systems. An accelerated sampling of the slow degrees of freedom is achieved by generating a single self-adaptive trajectory that tends to span uniformly the reaction coordinate using a time dependent bias potential derived from the preceding history of the system. To show the convergent behavior and the efficiency of the method, we present the free energy surface of alanine dipeptide in water as a function of the backbone dihedral angles.     

Essential energy space random walk via energy space metadynamics method to accelerate molecular dynamics simulations

To overcome the possible pseudoergodicity problem, molecular dynamic simulation can be accelerated via the realization of an energy space random walk. To achieve this, a biased free energy function (BFEF) needs to be priori obtained. Although the quality of BFEF is essential for sampling efficiency, its generation is usually tedious and nontrivial. In this work, we present an energy space metadynamics algorithm to efficiently and robustly obtain BFEFs. Moreover, in order to deal with the associated diffusion sampling problem caused by the random walk in the total energy space, the idea in the original umbrella sampling method is generalized to be the random walk in the essential energy space, which only includes the energy terms determining the conformation of a region of interest. This essential energy space generalization allows the realization of efficient localized enhanced sampling and also offers the possibility of further sampling efficiency improvement when high frequency energy terms irrelevant to the target events are free of activation. The energy space metadynamics method and its generalization in the essential energy space for the molecular dynamics acceleration are demonstrated in the simulation of a pentanelike system, the blocked alanine dipeptide model, and the leucine model.     

Minimum free energy pathways and free energy profiles for conformational transitions based on atomistic molecular dynamics simulations

An efficient method for the calculation of minimum free energy pathways and free energy profiles for conformational transitions is presented. Short restricted perturbation-targeted molecular dynamics trajectories are used to generate an approximate free energy surface. Approximate reaction pathways for the conformational change are constructed from one-dimensional line segments on this surface using a Monte Carlo optimization. Accurate free energy profiles are then determined along the pathways by means of one-dimensional adaptive umbrella sampling simulations. The method is illustrated by its application to the alanine "dipeptide." Due to the low computational cost and memory demands, the method is expected to be useful for the treatment of large biomolecular systems.     

Estimating kinetic rates from accelerated molecular dynamics simulations: alanine dipeptide in explicit solvent as a case study

Molecular dynamics (MD) simulation is the standard computational technique used to obtain information on the time evolution of the conformations of proteins and many other molecular systems. However, for most biological systems of interest, the time scale for slow conformational transitions is still inaccessible to standard MD simulations. Several sampling methods have been proposed to address this issue, including the accelerated molecular dynamics method. In this work, we study the extent of sampling of the phi/psi space of alanine dipeptide in explicit water using accelerated molecular dynamics and present a framework to recover the correct kinetic rate constant for the helix to beta-strand transition. We show that the accelerated MD can drastically enhance the sampling of the phi/psi conformational phase space when compared to normal MD. In addition, the free energy density plots of the phi/psi space show that all minima regions are accurately sampled and the canonical distribution is recovered. Moreover, the kinetic rate constant for the helix to beta-strand transition is accurately estimated from these simulations by relating the diffusion coefficient to the local energetic roughness of the energy landscape. Surprisingly, even for such a low barrier transition, it is difficult to obtain enough transitions to accurately estimate the rate constant when one uses normal MD.     

Metadynamics in essential coordinates: free energy simulation of conformational changes

We propose an approach that combines an extraction of collective motions of a molecular system with a sampling of its free energy surface. A recently introduced method of metadynamics allows exploration of the free energy surface of a molecular system by means of coarse-grained dynamics with flooding of free energy minima. This free energy surface is defined as a function of a set of collective variables (e.g., interatomic distances, angles, torsions, and others). In this study, essential coordinates determined by essential dynamics (principle component analysis) were used as collective variables in metadynamics. First, dynamics of the model system (explicitly solvated alanine dipeptide, Ace-Ala-Nme) was simulated by a classical molecular dynamics simulation. The trajectory (1 ns) was then analyzed by essential dynamics to obtain essential coordinates. The free energy surface as a function of the first and second essential coordinates was then explored by metadynamics. The resulting free energy surface is in agreement with other studies of this system. We propose that a combination of these two methods (metadynamics and essential dynamics) has great potential in studies of conformational changes in peptides and proteins.     

Using Selectively Applied Accelerated Molecular Dynamics to Enhance Free Energy Calculations

Accelerated molecular dynamics (aMD) has been shown to enhance conformational space sampling relative to classical molecular dynamics; however, the exponential reweighting of aMD trajectories, which is necessary for the calculation of free energies relating to the classical system, is oftentimes problematic, especially for systems larger than small poly peptides. Here, we propose a method of accelerating only the degrees of freedom most pertinent to sampling, thereby reducing the total acceleration added to the system and improving the convergence of calculated ensemble averages, which we term selective aMD. Its application is highlighted in two biomolecular cases. First, the model system alanine dipeptide is simulated with classical MD, all-dihedral aMD, and selective aMD, and these results are compared to the infinite sampling limit as calculated with metadynamics. We show that both forms of aMD enhance the convergence of the underlying free energy landscape by 5-fold relative to classical MD; however, selective aMD can produce improved statistics over all-dihedral aMD due to the improved reweighting. Then we focus on the pharmaceutically relevant case of computing the free energy of the decoupling of oseltamivir in the active site of neuraminidase. Results show that selective aMD greatly reduces the cost of this alchemical free energy transformation, whereas all-dihedral aMD produces unreliable free energy estimates.     

An algorithm to find minimum free-energy paths using umbrella integration

The calculation of free-energy barriers by umbrella sampling and many other methods is hampered by the necessity for an a priori choice of the reaction coordinate along which to sample. We avoid this problem by providing a method to search for saddle points on the free-energy surface in many coordinates. The necessary gradients and Hessians of the free energy are obtained by multidimensional umbrella integration. We construct the minimum free-energy path by following the gradient down to minima on the free-energy surface. The change of free energy along the path is obtained by integrating out all coordinates orthogonal to the path. While we expect the method to be applicable to large systems, we test it on the alanine dipeptide in vacuum. The minima, transition states, and free-energy barriers agree well with those obtained previously with other methods.     

Calibration of force-field dependency in free energy landscapes of peptide conformations by quantum chemical calculations

The free energy landscapes of peptide conformations were calibrated by ab initio quantum chemical calculations, after the enhanced conformational diversity search using the multicanonical molecular dynamics simulations. Three different potentials of mean force for an isolated dipeptide were individually obtained by the multicanonical molecular dynamics simulations using the conventional force fields, AMBER parm94, AMBER parm96, and CHARMm22. Each potential of mean force was then calibrated based upon the umbrella sampling algorithm from the adiabatic energy map that was calculated separately by the ab initio molecular orbital method, and all of the calibrated potentials of mean force coincided well. The calibration method was also applied to the simulations of a peptide dimer in explicit water models, and it was shown that the calibrated free energy landscapes did not depend on the force field used in the classical simulations, as far as the conformational space was sampled well. The current calibration method fuses the classical free energy calculation with the quantum chemical calculation, and it should generally make simulations for biomolecular systems much more reliable when combining with enhanced conformational sampling.     

Theoretical modeling of postionization fragmentation of rare-gas trimer cations

The dynamics of ionic rare-gas trimers (Ar(3) (+), Kr(3) (+), and Xe(3) (+)) produced by a sudden ionization of neutral precursors is investigated theoretically with a hybrid classical-quantum method for solving the equations of motion governed by a Hamiltonian obtained from a previously tested diatomics-in-molecules model. Initial conditions are selected with Monte Carlo sampling. Two possibilities for generating the initial electronic state are considered: diabatic (local) and adiabatic (delocalized). The dynamics generally leads to fragmentation, producing either monomer ions or dimer ions in a relatively short time; however, a large number of long-lived metastable trimer ions are also seen in some cases. We have analyzed the dynamics with respect to the fraction of monomer ions produced, the distribution of the kinetic energy of the products, and the distribution of fragmentation times of the trimers. Initial diabatic ionization is associated with much faster fragmentation than adiabatic ionization. Spin-orbit coupling plays an important role in the fragmentation dynamics.     

Bennett's acceptance ratio and histogram analysis methods enhanced by umbrella sampling along a reaction coordinate in configurational space

Free energy perturbation, a method for computing the free energy difference between two states, is often combined with non-Boltzmann biased sampling techniques in order to accelerate the convergence of free energy calculations. Here we present a new extension of the Bennett acceptance ratio (BAR) method by combining it with umbrella sampling (US) along a reaction coordinate in configurational space. In this approach, which we call Bennett acceptance ratio with umbrella sampling (BAR-US), the conditional histogram of energy difference (a mapping of the 3N-dimensional configurational space via a reaction coordinate onto 1D energy difference space) is weighted for marginalization with the associated population density along a reaction coordinate computed by US. This procedure produces marginal histograms of energy difference, from forward and backward simulations, with higher overlap in energy difference space, rendering free energy difference estimations using BAR statistically more reliable. In addition to BAR-US, two histogram analysis methods, termed Bennett overlapping histograms with US (BOH-US) and Bennett-Hummer (linear) least square with US (BHLS-US), are employed as consistency and convergence checks for free energy difference estimation by BAR-US. The proposed methods (BAR-US, BOH-US, and BHLS-US) are applied to a 1-dimensional asymmetric model potential, as has been used previously to test free energy calculations from non-equilibrium processes. We then consider the more stringent test of a 1-dimensional strongly (but linearly) shifted harmonic oscillator, which exhibits no overlap between two states when sampled using unbiased Brownian dynamics. We find that the efficiency of the proposed methods is enhanced over the original Bennett's methods (BAR, BOH, and BHLS) through fast uniform sampling of energy difference space via US in configurational space. We apply the proposed methods to the calculation of the electrostatic contribution to the absolute solvation free energy (excess chemical potential) of water. We then address the controversial issue of ion selectivity in the K(+) ion channel, KcsA. We have calculated the relative binding affinity of K(+) over Na(+) within a binding site of the KcsA channel for which different, though adjacent, K(+) and Na(+) configurations exist, ideally suited to these US-enhanced methods. Our studies demonstrate that the significant improvements in free energy calculations obtained using the proposed methods can have serious consequences for elucidating biological mechanisms and for the interpretation of experimental data.     

Harmonic Fourier beads method for studying rare events on rugged energy surfaces

We present a robust, distributable method for computing minimum free energy paths of large molecular systems with rugged energy landscapes. The method, which we call harmonic Fourier beads (HFB), exploits the Fourier representation of a path in an appropriate coordinate space and proceeds iteratively by evolving a discrete set of harmonically restrained path points-beads-to generate positions for the next path. The HFB method does not require explicit knowledge of the free energy to locate the path. To compute the free energy profile along the final path we employ an umbrella sampling method in two generalized dimensions. The proposed HFB method is anticipated to aid the study of rare events in biomolecular systems. Its utility is demonstrated with an application to conformational isomerization of the alanine dipeptide in gas phase.     

Generalized essential energy space random walks to more effectively accelerate solute sampling in aqueous environment

Molecular dynamics sampling can be enhanced via the promoting of potential energy fluctuations, for instance, based on a Hamiltonian modified with the addition of a potential-energy-dependent biasing term. To overcome the diffusion sampling issue, which reveals the fact that enlargement of event-irrelevant energy fluctuations may abolish sampling efficiency, the essential energy space random walk (EESRW) approach was proposed earlier. To more effectively accelerate the sampling of solute conformations in aqueous environment, in the current work, we generalized the EESRW method to a two-dimension-EESRW (2D-EESRW) strategy. Specifically, the essential internal energy component of a focused region and the essential interaction energy component between the focused region and the environmental region are employed to define the two-dimensional essential energy space. This proposal is motivated by the general observation that in different conformational events, the two essential energy components have distinctive interplays. Model studies on the alanine dipeptide and the aspartate-arginine peptide demonstrate sampling improvement over the original one-dimension-EESRW strategy; with the same biasing level, the present generalization allows more effective acceleration of the sampling of conformational transitions in aqueous solution. The 2D-EESRW generalization is readily extended to higher dimension schemes and employed in more advanced enhanced-sampling schemes, such as the recent orthogonal space random walk method.     

An application of coupled reference interaction site model/molecular dynamics to the conformational analysis of the alanine dipeptide

We present an application of our recently proposed coupled reference interaction site model (RISM) molecular dynamics (MD) solvation free energy methodology [Freedman and Truong, Chem. Phys. Lett. 381, 362 (2003); J. Chem. Phys. 121, 2187 (2004)] to study the conformational stability of alanine dipeptide in aqueous solution. In this methodology, radial distribution functions obtained from a single MD simulation are substituted into a RISM expression for solvation free energy. Consequently, iterative solution of the RISM equation is not needed. The relative solvation free energies of seven different conformations of the alanine dipeptide in aqueous solution are calculated. Results from the coupled RISM/MD methodology are in good agreement with those from earlier simulations using the accurate free energy perturbation approach, showing that the alphaR conformation is most stabilized by solution. This study establishes a framework for applying this coupled RISM/MD method to larger biological systems.     

Sensitivity field for nonautonomous molecular rotors

We propose a numerical approach to quantify the control of a nonautonomous molecular rotor motion. Unlike straightforward molecular dynamics simulations in an explicitly time-dependent framework, our method is based on the theory of geometric phases. This theory allows us to define a sensitivity field (SF) in control parameter space that characterizes average motion of a molecule induced by a cyclic perturbation. We show that the SF can be obtained using only equilibrium free energy sampling techniques. A density plot of the SF quantifies response of a molecule to an arbitrary cyclic adiabatic evolution of parameters. For demonstration, we numerically find the SFs for two surface mounted molecular rotor molecules that can be driven, in practice, by strong time-dependent electric fields of a STM tip.     

Computer simulation of electron transfer processes across the electrode|electrolyte interface: a treatment of solvent and electrode polarizability

A polarizable molecular dynamics model for adiabatic electron transfer across the electrode|electrolyte interface is presented. The electronic polarizability of the water and of the metal electrode is accounted for by a dynamical fluctuating charge algorithm, image charges, and the Ewald summation adapted for a conducting interface. The effects of the solvent electronic polarizability are studied by computing the diabatic and adiabatic free energy curves for both polarizable and non-polarizable water models. This represents the first effort to compute the adiabatic free energy curves from simulation for a fully polarizable electrochemical system.     

Free energy derivatives: A new method for probing the convergence problem in free energy calculations

The convergence behavior of free energy calculations has been explored in more detail than in any previously reported work, using a model system of two neon atoms in a periodic box of water. We find that for thermodynamic integration-type free energy calculations as much as a nanosecond or more molecular dynamics sampling is required to obtain a fully converged value for a single λ point of the integrand. The concept of “free energy derivatives” with respect to the individual parameters of the force field is introduced. This formalism allows the total convergence of the simulation to be deconvoluted into components. A determination of the statistical “sampling ratio” from these simulations indicates that for window-type free energy calculations carried out in a periodic waterbox of typical size at least 0.6 ps of sampling should be performed at each window (0.7 ps if constraint contributions to the free energy are being determined). General methods to estimate and reduce the error in thermodynamic integration and free energy perturbation calculations are discussed. We show that the difficulty in applying such methods is determining a reliable estimate of the correlation length from a short series of data. © 1994 by John Wiley & Sons, Inc.     

BayesWHAM: A Bayesian approach for free energy estimation,reweighting, and uncertainty quantification in the weighted histogram analysis method

The weighted histogram analysis method (WHAM) is a powerful approach to estimate molecular free energy surfaces (FES) from biased simulation data. Bayesian reformulations of WHAM are valuable in proving statistically optimal use of the data and providing a transparent means to incorporate regularizing priors and estimate statistical uncertainties. In this work, we develop a fully Bayesian treatment of WHAM to generate statistically optimal FES estimates in any number of biasing dimensions under arbitrary choices of the Bayes prior. Rigorous uncertainty estimates are generated by Metropolis‐Hastings sampling from the Bayes posterior. We also report a means to project the FES and its uncertainties into arbitrary auxiliary order parameters beyond those in which biased sampling was conducted. We demonstrate the approaches in applications of alanine dipeptide and the unthreading of a synthetic mimic of the astexin‐3 lasso peptide. Open‐source MATLAB and Python implementations of our codes are available for free public download. © 2017 Wiley Periodicals, Inc.     

An Anderson impurity model for efficient sampling of adiabatic potential energy surfaces of transition metal complexes

We present a model intended for rapid sampling of ground and excited state potential energy surfaces for first-row transition metal active sites. The method is computationally inexpensive and is suited for dynamics simulations where (1) adiabatic states are required "on-the-fly" and (2) the primary source of the electronic coupling between the diabatic states is the perturbative spin-orbit interaction among the 3d electrons. The model Hamiltonian we develop is a variant of the Anderson impurity model and achieves efficiency through a physically motivated basis set reduction based on the large value of the d-d Coulomb interaction U(d) and a Lanczos matrix diagonalization routine to solve for eigenvalues. The model parameters are constrained by fits to the partial density of states obtained from ab initio density functional theory calculations. For a particular application of our model we focus on electron transfer occurring between cobalt ions solvated by ammonium, incorporating configuration interaction between multiplet states for both metal ions. We demonstrate the capability of the method to efficiently calculate adiabatic potential energy surfaces and the electronic coupling factor we have calculated compares well to previous calculations and experiment. (