On the convergence improvement in the metadynamics simulations: a Wang-Landau recursion approach

As a popular tool in exploring free energy landscapes, the metadynamics method has been widely applied to elucidate various chemical or biochemical processes. As deeply discussed by Laio et al. [J. Phys. Chem. B 109, 6714 (2005)], the size of the updating Gaussian function is pivotal to the free energy convergence toward the target free energy surface. For instance, a greater Gaussian height can facilitate the quick visit of a conformation region of interest; however, it may lead to a larger error of the calculated free energy surface. In contrast, a lower Gaussian height can guarantee a better resolution of the calculated free energy surface; however, it will take longer time for such a simulation to navigate through the defined conformational region. In order to reconcile such confliction, the authors present a method by implementing the Wang-Landau recursion scheme in the metadynamics simulations to adaptively update the height of the unit Gaussian function. As demonstrated in their model studies on both a toy system, and a realistic molecular system treated with the hybrid quantum mechanical and molecular mechanical (QMMM) potential, the present approach can quickly result in more decently converged free energy surfaces, compared with the classical metadynamics simulations employing the fixed Gaussian heights.     

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.     

A scalable parallel Monte Carlo method for free energy simulations of molecular systems

We present a method of parallelizing flat histogram Monte Carlo simulations, which give the free energy of a molecular system as an output. In the serial version, a constant probability distribution, as a function of any system parameter, is calculated by updating an external potential that is added to the system Hamiltonian. This external potential is related to the free energy. In the parallel implementation, the simulation is distributed on to different processors. With regular intervals the modifying potential is summed over all processors and distributed back to every processor, thus spreading the information of which parts of parameter space have been explored. This implementation is shown to decrease the execution time linearly with added number of processors.     

A highly efficient hybrid method for calculating the hydration free energy of a protein

We develop a new method for calculating the hydration free energy (HFE) of a protein with any net charge. The polar part of the energetic component in the HFE is expressed as a linear combination of four geometric measures (GMs) of the protein structure and the generalized Born (GB) energy plus a constant. The other constituents in the HFE are expressed as linear combinations of the four GMs. The coefficients (including the constant) in the linear combinations are determined using the three‐dimensional reference interaction site model (3D‐RISM) theory applied to sufficiently many protein structures. Once the coefficients are determined, the HFE and its constituents of any other protein structure are obtained simply by calculating the four GMs and GB energy. Our method and the 3D‐RISM theory give perfectly correlated results. Nevertheless, the computation time required in our method is over four orders of magnitude shorter.     

Achieving fast convergence of ab initio free energy perturbation calculations with the adaptive force-matching method

This paper studies the possibility of improving the convergence of ab initio free energy perturbation (FEP) calculations by developing customized force fields with the adaptive force-matching (AFM) method. The ab initio FEP method relies on a molecular mechanics (MM) potential to sample configuration space. If the Boltzmann weight of the MM sampling is close to that of the ab initio method, the efficiency of ab initio FEP will be optimal. The difference in the Boltzmann weights can be quantified by the relative energy difference distribution (REDD). The force field developed through AFM significantly improves the REDD when compared with standard MM models, thus improving the convergence of the ab initio FEP calculation. The static dielectric constant ε_s of ice-Ih was studied with PW-91 through ab initio FEP. With a customized force field developed through AFM, we were able to converge ε_s to 80 ± 4 with 3,600 configurations. A similar ab initio FEP calculation with the TIP4P model would require 220 times more configurations to achieve the same accuracy. Our study indicates that the PW-91 functional underestimates ice-Ih ε_s by about 20%.     

Simulated scaling method for localized enhanced sampling and simultaneous "alchemical" free energy simulations: a general method for molecular mechanical, quantum mechanical, and quantum mechanical/molecular mechanical simulations

A potential scaling version of simulated tempering is presented to efficiently sample configuration space in a localized region. The present "simulated scaling" method is developed with a Wang-Landau type of updating scheme in order to quickly flatten the distributions in the scaling parameter lambdam space. This proposal is meaningful for a broad range of biophysical problems, in which localized sampling is required. Besides its superior capability and robustness in localized conformational sampling, this simulated scaling method can also naturally lead to efficient "alchemical" free energy predictions when dual-topology alchemical hybrid potential is applied; thereby simultaneously, both of the chemically and conformationally distinct portions of two end point chemical states can be efficiently sampled. As demonstrated in this work, the present method is also feasible for the quantum mechanical and quantum mechanical/molecular mechanical simulations.     

A hybrid approach combining energy density analysis with the interaction energy decomposition method

We propose a new analysis technique for characterizing molecular interactions that combines an energy decomposition scheme, such as the Kitaura-Morokuma decomposition method, with energy density analysis, which partitions the total energy of the system into atomic contributions. The combined scheme, termed Interaction-EDA, enables us to estimate the local contribution of interaction energy components, such as electrostatic, exchange, polarization, and charge transfer. The evaluation of the local interaction energy is rather important in large systems. For a numerical assessment, the Interaction-EDA method is applied to the process of CO adsorption on Si(100) - (2 x 1) surface.     

A method to estimate the Gibbs free energy of non-equilibrium alloys by thermal analysis

A method has been developed to estimate the Gibbs free energy $ \left( {G_{\text{S}}^{\text{NE}} } \right) $ of the non-equilibrium solid alloys with multicomponents based on differential scanning calorimetry (DSC) analysis. In this method, the DSC curves of the non-equilibrium and equilibrium alloys during heating up to fully melting and those of the alloys during solidifying were measured. Then the thermal effects of the solid phase transformations from non-equilibrium to equilibrium states and the equilibrium solidification could be calculated. By evolving the traditional equal-G curve principle to equal-G point, the Gibbs free energy of the equilibrium solid alloy with multicomponents could be obtained on condition that the free energy of the liquid alloy was known. Considering the thermal effects of the solid phase transformations from non-equilibrium to equilibrium states, the Gibbs free energy value of the non-equilibrium alloys with a given composition could be achieved although the phase constitution of the equilibrium solid alloys and the Gibbs free energy of each phase were not known, and the calculation errors could be reduced by dividing the alloys into many infinitesimal virtual pure metals. The Gibbs free energy of the non-equilibrium Al?CSi?CMn alloys was calculated by using this method, confirming the validity of this method.     

A simple mathematical method to determine the efficiencies of log-conical detectors

A new type of photon detector, log-conical, is proposed. The average path length traveled by an incident photon of arbitrary energy as well as the geometrical solid angle are calculated in a mathematical expression to determine the efficiencies of this detector for an arbitrarily positioned isotropic radiating point source. The off-axis effect of the source position was analyzed to demonstrate the powerful capability of the proposed method. The results are compared with those obtained using a standard 3″×3″ cylindrical detector of the same volume in order to show the enhanced efficiency of the log-conical detector.     

From A to B in free energy space

The authors present a new method for searching low free energy paths in complex molecular systems at finite temperature. They introduce two variables that are able to describe the position of a point in configurational space relative to a preassigned path. With the help of these two variables the authors combine features of approaches such as metadynamics or umbrella sampling with those of path based methods. This allows global searches in the space of paths to be performed and a new variational principle for the determination of low free energy paths to be established. Contrary to metadynamics or umbrella sampling the path can be described by an arbitrary large number of variables, still the energy profile along the path can be calculated. The authors exemplify the method numerically by studying the conformational changes of alanine dipeptide.     

Calculation of the entropy and free energy of peptides by molecular dynamics simulations using the hypothetical scanning molecular dynamics method

Hypothetical scanning (HS) is a method for calculating the absolute entropy S and free energy F from a sample generated by any simulation technique. With this approach each sample configuration is reconstructed with the help of transition probabilities (TPs) and their product leads to the configuration's probability, hence to the entropy. Recently a new way for calculating the TPs by Monte Carlo (MC) simulations has been suggested, where all system interactions are taken into account. Therefore, this method--called HSMC--is in principle exact where the only approximation is due to insufficient sampling. HSMC has been applied very successfully to liquid argon, TIP3P water, self-avoiding walks on a lattice, and peptides. Because molecular dynamics (MD) is considered to be significantly more efficient than MC for a compact polymer chain, in this paper HSMC is extended to MD simulations as applied to peptides. Like before, we study decaglycine in vacuum but for the first time also a peptide with side chains, (Val)(2)(Gly)(6)(Val)(2). The transition from MC to MD requires implementing essential changes in the reconstruction process of HSMD. Results are calculated for three microstates, helix, extended, and hairpin. HSMD leads to very stable differences in entropy TDeltaS between these microstates with small errors of 0.1-0.2 kcal/mol (T=100 K) for a wide range of calculation parameters with extremely high efficiency. Various aspects of HSMD and plans for future work are discussed.     

A hybrid potential reaction path and free energy study of the chorismate mutase reaction.

We present a combination of two techniques--QM/MM statistical simulation methods and QM/MM internal energy minimizations--to get a deeper insight into the reaction catalyzed by the enzyme chorismate mutase. Structures, internal energies and free energies, taken from the paths of the reaction in solution and in the enzyme have been analyzed in order to estimate the relative importance of the reorganization and preorganization effects. The results we obtain for this reaction are in good agreement with experiment and show that chorismate mutase achieves its catalytic efficiency in two ways; first, it preferentially binds the active conformer of the substrate and, second, it reduces the free energy of activation for the reaction relative to that in solution by providing an environment which stabilizes the transition state.     

Free energy simulations: use of reverse cumulative averaging to determine the equilibrated region and the time required for convergence

A method is proposed for improving the accuracy and efficiency of free energy simulations. The essential idea is that the convergence of the relevant measure (e.g., the free energy derivative in thermodynamic integration) is monitored in the reverse direction starting from the last frame of the trajectory, instead of the usual approach, which begins with the first frame and goes in the forward direction. This simple change in the use of the simulation data makes it straightforward to eliminate the contamination of the averages by contributions from the equilibrating region. A statistical criterion is introduced for distinguishing the equilibrated (production) region from the equilibrating region. The proposed method, called reverse cumulative averaging, is illustrated by its application to the well-studied case of the alchemical free energy simulation of ethane to methanol.     

A diffusional alternative to the Marcus free energy gap law

The results of the precise kinetic fitting of the highly exergonic electron transfer from excited perylene to tetracyanoethylene in acetonitrile were used to estimate the Stern-Volmer constant of perylene quenching by double channel electron transfer (to the ground and excited states of the radical ion pair). It appears that the Stern-Volmer constant is exactly the same as the diffusional height of the Rehm-Weller plateau, substituting the exergonic wing of the Marcus free energy gap law. Even the single channel transfer is shown to be fast enough that the quenching should be under true diffusional control at the highest available exergonicity.     

Toward efficient chemical potential calculations by expanded ensemble simulations; to make the free energy pathway fairly level

A scheme is suggested of how to construct good bias potentials ("balancing factors") to be used in expanded ensemble (EE) calculations of chemical potentials of solutions. A combination of two strategies are used: (i) to use a pathway for particle insertions that avoids large variations in free energy and (ii) to use calculated free energy derivatives to construct a bias potential that makes the pathway fairly level. Only a few very short simulations are needed to accomplish the latter, and then, a full EE simulation is done to obtain the chemical potential. By practical calculations of the chemical potential of benzene, cyclohexane, and benzylamine in water, it is shown that this method is at least equally efficient to the recent adaptive EE (AEE) method by Aberg et al. (J. Chem. Phys. 2004, 120, 3370). Furthermore, the new method provides an alternative strategy that complements existing EE methods.     

Efficiency of alchemical free energy simulations. I. A practical comparison of the exponential formula, thermodynamic integration, and Bennett's acceptance ratio method

We investigate the relative efficiency of thermodynamic integration, three variants of the exponential formula, also referred to as thermodynamic perturbation, and Bennett's acceptance ratio method to compute relative and absolute solvation free energy differences. Our primary goal is the development of efficient protocols that are robust in practice. We focus on minimizing the number of unphysical intermediate states (λ-states) required for the computation of accurate and precise free energy differences. Several indicators are presented which help decide when additional λ-states are necessary. In all tests Bennett's acceptance ratio method required the least number of λ-states, closely followed by the "double-wide" variant of the exponential formula. Use of the exponential formula in only strict "forward" or "backward" mode was not found to be competitive. Similarly, the performance of thermodynamic integration in terms of efficiency was rather poor. We show that this is caused by the use of the trapezoidal rule as method of numerical quadrature. A systematic study focusing on the optimization of thermodynamic integration is presented in a companion paper.     

A density functional theory based estimation of the anharmonic contributions to the free energy of a polypeptide helix

We have employed density functional theory to determine the temperature dependence of the intrinsic stability of an infinite poly-L-alanine helix. The most relevant helix types, i.e., the α- and the 3(10)-helix, and several unfolded conformations, which serve as reference for the stability analysis, have been included. For the calculation of the free energies for the various chain conformations we have explicitly included both, harmonic and anharmonic contributions. The latter have been calculated by means of a thermodynamic integration approach employing stochastic Langevin molecular dynamics, which is shown to provide a dramatic increase in the computational efficiency as compared to commonly employed deterministic molecular dynamics schemes. Employing this approach we demonstrate that the anharmonic part of the free energy amounts to the order of 0.1-0.4 kcal/mol per peptide unit for all analysed conformations. Although small, the anharmonic contribution stabilizes the helical conformations with respect to the fully extended structure.     

Calculating distribution coefficients based on multi-scale free energy simulations: an evaluation of MM and QM/MM explicit solvent simulations of water-cyclohexane transfer in the SAMPL5 challenge

One of the central aspects of biomolecular recognition is the hydrophobic effect, which is experimentally evaluated by measuring the distribution coefficients of compounds between polar and apolar phases. We use our predictions of the distribution coefficients between water and cyclohexane from the SAMPL5 challenge to estimate the hydrophobicity of different explicit solvent simulation techniques. Based on molecular dynamics trajectories with the CHARMM General Force Field, we compare pure molecular mechanics (MM) with quantum-mechanical (QM) calculations based on QM/MM schemes that treat the solvent at the MM level. We perform QM/MM with both density functional theory (BLYP) and semi-empirical methods (OM1, OM2, OM3, PM3). The calculations also serve to test the sensitivity of partition coefficients to solute polarizability as well as the interplay of the quantum-mechanical region with the fixed-charge molecular mechanics environment. Our results indicate that QM/MM with both BLYP and OM2 outperforms pure MM. However, this observation is limited to a subset of cases where convergence of the free energy can be achieved.