Escaping freeenergy minima
See allHide authors and affiliations

Communicated by David Chandler, University of California, Berkeley, CA (received for review May 8, 2002)
Abstract
We introduce a powerful method for exploring the properties of the multidimensional free energy surfaces (FESs) of complex manybody systems by means of coarsegrained nonMarkovian dynamics in the space defined by a few collective coordinates. A characteristic feature of these dynamics is the presence of a historydependent potential term that, in time, fills the minima in the FES, allowing the efficient exploration and accurate determination of the FES as a function of the collective coordinates. We demonstrate the usefulness of this approach in the case of the dissociation of a NaCl molecule in water and in the study of the conformational changes of a dialanine in solution.
Molecular dynamics (MD) and the Monte Carlo simulation method have had a very deep influence on the most diverse fields, from materials science to biology and from astrophysics to pharmacology. Yet, despite their success, these simulation methods suffer from limitations that reduce the scope of their applications. A severe constraint is the limited time scale that presentday computer technology and sampling algorithms explore. In particular, there are many circumstances where the free energy surface (FES) has several local minima separated by large barriers. Examples of these situations include conformational changes in solution, protein folding, firstorder phase transitions, and chemical reactions. In such circumstances a simulation started in one minimum will be able to move spontaneously to the next minimum only under very favorable circumstances. A host of methods have been suggested to lift this restriction and explore the FES (1–13) or to characterize the transition state (14, 15). Here we propose a solution to this problem by combining the ideas of coarsegrained dynamics (16, 17) on the FES (10, 12) with those of adaptive bias potential methods (2, 11), obtaining a procedure that allows the system to escape from local minima in the FES and, at the same time, permits a quantitative determination of the FES as a byproduct of the integrated process.
Methodology
We shall assume here that there exists a finite number of relevant collective coordinates s_{i}, i = 1,n where n is a small number, and we consider the dependence of the free energy ℱ(s) on these parameters. Practical examples of appropriate choices of these variables will be given below. The exploration of the FES is guided by the forces F = −∂ℱ/∂s. To estimate these forces efficiently, we introduce an ensemble of P replicas of the system, each obeying the constraint that the collective coordinates have a preassigned value s_{i} = s, and each evolved independently at the same temperature T. Since the P replicas are statistically independent, the estimate of thermodynamic observables (e.g., the forces on the constraints) is improved with respect to an evaluation on a single replica, and it can be parallelized in a straightforward manner. The constraints are imposed on each replica via the standard methods of constrained molecular dynamics (18) by adding to the Lagrangean a term Σ_{i}_{=1,}_{n}λ_{i}(s_{i} −s) where λ_{i} are Lagrange multipliers. Averaging over the time and over the replicas we can evaluate the derivative of the free energy relative to s as F = 〈λ_{i}〉 (1, 4). For the sake of simplicity we neglect here the small kinematic correction term discussed in ref. 4, although this can be added easily. Taking a cue from the work of Kevrikidis et al. (16), we use these forces, determined by microscopic dynamics, to define coarsegrained dynamics in the space of the s_{i}s. In our case the definition of these dynamics is rather arbitrary and designed only to explore the FES efficiently. The dynamics is defined from the discretized evolution equation: In Eq. 1 we have introduced the scaled variables σ = s/Δs_{i} and the scaled forces φ = FΔs_{i}, where Δs_{i} is the estimated size of the FES well in the direction s_{i}, φ^{t} is the modulus of the nth dimensional vector φ and δσ is a dimensionless stepping parameter. It should be stressed that Eq. 1 has the form of a steepest descent equation in the direction given by the forces φ_{i} and does not imply any real dynamical evolution. The index t is only used to label the states. After the collective coordinates are updated by using Eq. 1, a new ensemble of replicas of the system with values σ is prepared, and new forces F are calculated for the next iteration. At the same time the driving forces are evaluated from the microscopic Hamiltonian in short standard microscopic MD runs. Since there is no dynamical continuity between the P replicas at different iterations, one can use large values of δσ and move very efficiently in the space of the collective coordinates.
Clearly Eq. 1 alone cannot guarantee an efficient exploration of the FES, nor is it useful to determine the FES. This task can be achieved if we replace the forces in Eq. 1 with a historydependent term (2, 11): where the height and the width of the Gaussian W and δσ are chosen to provide a reasonable compromise between accuracy and efficiency in exploring the FES, as we will show below. The component of the forces coming from the Gaussian will discourage the system from revisiting the same spot and encourage an efficient exploration of the FES. As the system diffuses although the FES, the Gaussian potentials accumulate and fill the FES well, allowing the system to migrate from well to well. After a while the sum of the Gaussian terms will almost exactly compensate the underlying FES well. A typical example of this behavior is shown in Fig. 1, in which the dynamics in Eq. 1 is used to explore a onedimensional potential energy surface (PES) V(σ) with three wells. The dynamics start from a local minimum that is filled by the Gaussians in ≈20 steps. Then the dynamics escapes from the well from the lowest energy saddle point, filling the second well in ≈80 steps. The second highest saddle point is reached in ≈160 steps, and the full PES is filled in a total of ≈320 steps. Hence, in the case of this example, since the form of the potential is known, it can be verified that, for large t and if the width of the Gaussians is sufficiently small with respect to the length of a typical variation of V, modulus an additive constant. The same property can be verified for any multidimensional potential if the variables σ_{i} are allowed to vary in a finite region. Hence, we assume that Eq. 3 holds also if the method is used to explore the FES, and that the free energy can be estimated from Eq. 3 for large t. This procedure resembles the algorithm recently proposed by Wang and Landau (11), in which a timedependent bias is introduced to modify the density of states to produce locally flat histograms.
We observe an interesting isomorphism between our dynamics and the ordered deposition of the beads of a polymer chain on the FES. In fact, we can regard σ as the position of a monomer in the ndimensional parameter space. Each monomer is held at a distance δσ from the previous one and monomers repel each other with the Gaussian term of Eq. 2. Neglecting terms of order δσ^{2} each monomer is placed by Eq. 1 in the position of free energy minimum compatible with the requirement of minimum overlap with previous beads. Hence, the polymer chain generated in this manner fills the wells in the FES. This isomorphism helps to clarify how our approach works and why it can be made more precise by reducing δσ.
The efficiency of the method in filling a well in the PES (or in the FES) can be estimated by the number of Gaussians that are needed to fill the well. This number is proportional to (1/δσ)^{n}, where n is the dimensionality of the problem. Hence, the efficiency of the method scales exponentially with the number of dimensions involved. If n is large, the only way to obtain reasonable efficiencies is to use Gaussians with a size comparable to that of the well. On the other hand, a sum of Gaussians can only reproduce features of the FES on a scale larger than ≈δσ. A judicious choice of Δs_{i}, W, and δσ will ensure the right compromise between accuracy and sampling efficiency, and the optimal height and width of the Gaussians are determined by the typical variations of the FES.
If prior information on the nature of the freeenergy well is not available, the scaling parameters Δs_{i} are chosen by performing short coarsegrained dynamic runs without bias potential (see Eq. 1). In such a case the system moves around the minimum. The scaling parameters are chosen so that the elongations have approximately the same value in all directions. This amounts to an empirical form of preconditioning that makes the FES minimum nearly spherical in n dimensions and easy to fill with ndimensional Gaussians.
The metastep length δσ determines the efficiency of the method in a manner that is exponential in n and it should be made as large as possible. However, the dynamics in Eq. 1 is able to reconstruct details of the FES only on the length scale of δσ. Moreover, an overlarge value of δσ would cause the system to jump irregularly from one basin to the other. The value of δσ is chosen requiring the metatrajectory to remain localized in the FES minimum if the bias potential is not applied. With this choice of δσ, a single step of the dynamics in Eq. 1 cannot lead a system from a FES minimum to another, and the initial state at each new iteration can be generated from the last MD step of the previous iteration without creating major overlaps. Since the shake algorithm is used to impose a new set of σ values, large forces on the constraints are generated at the beginning of each microscopic dynamics, and the initial part of the trajectory has to be discarded for the calculation of F.
The height of the Gaussians W can be estimated as follows. We notice that we want to reach a situation in which the modified FES is flat. In such a case the forces coming from the FES and that coming from the Gaussian approximately balance each other. If we require the maximum value of the Gaussian forces to be smaller than the typical FES force, we arrive at the relation W/δσe^{−1/2} = α〈F〉^{1/2} with α < 1. In all cases so far studied a value of α close to 0.5 allows a fast escape from the local minima in the FES without significant loss of the underlying structure.
In the future adaptive ways of determining all of these parameters should be considered for adapting the procedure in an optimal manner to the local shape of the FES.
Results
Dissociation of NaCl in Water.
We first discuss an application of the method to the dissociation of a NaCl molecule in water. The most stable state is the dissociated one, whereas the undissociated contact ion pair corresponds to a metastable local minimum. Here we shall study how the system escapes from this metastable state. We model the system by using the AMBER95 force field (20) and solvate the NaCl complex in 215 TIP3P water molecules with periodic boundary conditions. The electrostatic interaction between classical atoms is taken into account by the P3M method (19). As collective coordinates, we use the distance r_{Na}_{}_{Cl} between Na and Cl, and, to take into account the dynamics of the solvation shell during the dissociation, we also use the electric fields V_{Na} and V_{Cl} on the Na and on the Cl due to the water molecules within ≈6.5 Å of the ions. V_{Na} and V_{Cl} depend significantly on the hydration pattern around the Na and on the Cl ion. If, for example, a hydrogen bond to one of the two ions is formed or broken, the field on the ion changes by several kcal/mol. A dynamic of the form in Eq. 1 was performed on the system, with δσ = 0.25 and W = 0.3 kcal/mol. The scaling parameters Δs_{i} were 0.53 Å, 32 and 32 kcal/mol for r_{Na}_{}_{Cl}, V_{Na}, and V_{Cl}, respectively. The forces in Eq. 1 were evaluated by short MD runs with a time step of 0.7 fs performed on six replicas of the system. The replicas were equilibrated for 200 MD steps at 300 K, and the forces were evaluated by averaging over the following 500 MD steps. Starting from a distance of 2.6 Å, corresponding to a contact pair, the ions dissociate after 120 coarsegrained iterations. The dynamics was then continued for another 350 steps, imposing a maximum separation of the ion pair of 5 Å, to study the structure of the FES around the transition state (see Fig. 2). Seven recrossings for the coarsegrained dynamics were observed. The transition state is located r_{Na}_{}_{Cl} ≈4.02 Å, V_{Na} ≈−120 kcal/mol, and V_{Cl} ≈81.5 kcal/mol (see Fig. 2). The overall topology of the FES confirms the importance of the solvent degrees of freedom for describing the reaction, since the transition region is in transverse orientation relative to the axis of the coordinates (21), and hence neither the distance nor the fields alone can provide an exhaustive description of the dissociation. The transition barrier (S in Fig. 2), estimated with thermodynamic integration (4) along r_{Na}_{}_{Cl}, is of 3.3 kcal/mol. The onedimensional free energy as a function of r_{Na}_{}_{Cl} can also be obtained by integrating the threedimensional FES (F_{3}) with respect to V_{Na} and V_{Cl} as F_{1} = 1/βlog[∫dV_{Na}dV_{Cl} exp(−βF_{3})]. This leads to a barrier of 3.4 kcal/mol. Hence, the method is able to reproduce the onedimensional free energy profile estimated with thermodynamic integration.
We have verified that these collective coordinates identify the transition state and the reaction coordinate also in the dynamical sense (14). To this effect we prepared a set of replicas of the system with the critical values of the collective coordinates. We then followed the procedure described in ref. 21, assigning to the particles values of the velocity chosen from a Maxwellian distribution at 300 K, and following the trajectory backward and forward in time. By this procedure, we found that, within a time scale of 200 fs, the overwhelming majority of the replicas fell into one of the two basins of attraction without any recrossing. Moreover, the timereversed dynamics leads systematically to the opposite basin of attraction (14). Finally, if we prepare the systems with the shake algorithm (18) for a given value of the collective coordinates, and if the particle velocities are not reassigned before letting the system evolve freely, the time needed to go in either of the attraction basins is of the order of a few ps, and several recrossings are observed. This is due to the fact that the shake algorithm imposes zero velocities in the direction of the constraint and it is only after thermalization, which takes place in a time scale of picoseconds, that the collective coordinates acquire sufficient velocity to evolve toward the FES minima. This is a strong confirmation that we have identified the correct transition state.
To characterize the solvent structure during the dissociation, we have computed the Na–water pair correlation function at the transition state S. The coordination number, obtained by integrating the pair correlation function up to a distance of 3 Å, is 5.28. This value is consistent with the picture depicted in ref. 21, in which the sodium at the transition state is shown to be typically only 5fold coordinated.
Isomerization of Alanine Dipeptide in Water.
As a second application, we used our method to explore the FES of an alanine dipeptide as a function of the backbone dihedral angles Φ and Ψ (ref. 22 and references therein). The system is described by the allatom AMBER95 force field (20) and solvated in 287 TIP3P water molecules. The parameters of the coarsegrained dynamics of Eq. 1 are δσ = 0.25 and W = 0.25 kcal/mol, while the scaling parameters Δs_{i} are 60° for both Φ and Ψ. The other settings are the same as for the simulation of the NaCl complex. The FES of alanine dipeptide in water has been studied in detail by other authors [refs. 22 (and references therein) and 23], who identified several minima. Our dynamics (1) is able to capture the overall features of the FES extremely quickly. Starting from a configuration close to the state C_{7}_{eq} [refs. 22 (and references therein) and 23], the dynamics fills this attraction basin in approximately 45 steps. The minimum is localized at Φ ≈ −80° and Ψ ≈ 150°. The estimated free energy at the saddle point is 2 kcal/mol higher than the free energy in C_{7}_{eq}. This compares to a value of 2.1 kcal/mol obtained in ref. 22 by thermodynamic integration. Our estimate is obtained with 45 steps of coarsegrained dynamics on six replicas, corresponding to a total of only ≈120 ps of microscopic dynamics. After another 150 steps, the dynamics fills the attraction basin of the state α_{r}, which corresponds to the equilibrium configuration of the system, and recrosses to C_{7}_{eq}. Finally, the dynamics visits the Φ > 0 region of the configuration space after ≈300 steps. In Fig. 3 we plot the free energy as a function of Φ and Ψ after 40, 100, and 200 steps of coarsegrained dynamics. Also small details of the free energy landscape, like the location of the saddle points and the presence of secondary minima [e.g., the state C_{5} (23)], can be detected graphically after a few steps of coarsegrained dynamics. The onedimensional free energy profile as a function of Ψ obtained by analytic integration of the FES reproduces the profile obtained in ref. 22 within 0.2 kcal/mol for all of the values of Ψ. Concerning the choice of the collective coordinates for this system, it should be remarked that the dihedral angles Φ and Ψ do not provide a complete description of the dialanine isomerization reaction (ref. 22 and references therein), and the real reaction coordinate also should take into account the solvent degrees of freedom. Still, our results show that the method is able to reproduce the overall feature of the FES even if relevant degrees of freedom are not explicitly included in the collective coordinate space. This shows that the neglected degrees of freedom, although relevant for determining the reaction coordinate, are associated with relatively small barriers and are sampled efficiently during the microscopic dynamics for each value of Φ and Ψ, despite the relatively short runs. This is due to the use of the replicas and to the ability of the dynamics to retrace the same region of parameter space during the dynamics.
Discussion
The method is based on the combination of the ideas of coarsegrained dynamics (16) in the space defined by a few collective coordinates and the introduction of a historydependent bias (2, 11). Constructing dynamics on a FES that depends on a few collective coordinates, we simplify enormously the complexity of the problem, which depends exponentially on the number of degrees of freedom. The FES is much smoother than the underlying PES and it is topologically simpler, with a greatly reduced number of local minima. The historydependent bias prevents the system from visiting regions it has already explored. This term is crucial for the efficiency of the method: for example, Langevin dynamics in the space of the collective coordinates would provide scant information on the underlying FES and it would be slower in locating global minima. This has been checked by performing a Langevin dynamics for the dissociation of a NaCl molecule in water, with the same collective variables introduced in Results and adding to Eq. 1 a Gaussian noise at a temperature T that is gradually reduced (24). The number of evaluations of the forces required to converge to the minimum of the FES (i.e., the dissociated state) is well above 1,000, if the temperature is decreased with a logarithic schedule, as required to ensure quasiergodicity in the collective coordinate space (25). A historydependent bias potential as defined in Eq. 2 but applied in a regular MD simulation without coarse graining the dynamics would be efficient in finding escapes from the local minima, as has been shown (2), but it would not provide quantitative information about the FES. In particular, Eq. 3 holds only because the dynamic is performed in the collective coordinates space using the derivatives of the free energy. Hence, both the coarsegrained dynamics on the FES and the historydependent bias are essential ingredients of the method, and only their combination allows an efficient and accurate determination of the FES.
Another advantage of our search algorithm is that it can be easily parallelized over the P replicas, thus optimally exploiting presentday computer architectures and, moreover, does not require a very accurate determination of the forces. In fact, in model calculations on an analytic PES, we have found that the algorithm is effective even for noise levels as large as 20–30% of the forces, since the height and the width of the Gaussians is such that the coarsegrained dynamics explores a region of space of the size of a Gaussian several times. and the system is only discouraged (and not prevented) from revisiting the same spot of configuration space. Hence, the errors in the forces tend to even out during the evolution. The dynamics (1) is “selfhealing,” i.e., it is in principle capable to compensate the effect of completely wrongly located Gaussians or of a wrong additional bias potential. We can therefore use quite short runs and a small number of replicas to evaluate the forces. This is at variance with thermodynamic integration (18) and the blue moon ensemble method (1, 4), where a very precise determination of the forces is needed and where treating multidimensional reaction coordinates is a daunting task. Other techniques, like the umbrella sampling or bias potential methods (5, 6, 9, 18), allow an accurate estimate of the FES in many dimensions. The drawback of these methods is that they are efficient only if a good approximation for the FES is known a priori, whereas this is not required in our method.
A last advantage of the method is its ability to provide qualitative information on the free energy of a system in a very short time. For example, the overall topology of a FES can be determined by very few coarsegrained dynamics steps by using large Gaussians. Subsequently, our qualitative knowledge of the FES can be improved by using smaller Gaussians, eventually reducing the dimensionality of the problem by exploiting the topological information obtained with the large Gaussians.
Acknowledgments
M.P. thanks Yannis Kevrikidis for a very illuminating discussion on the coarsegrained method and David Chandler for having patiently shared his understanding of the statistical mechanics of transition paths. A.L. thanks Joost VandeVondele for many precious suggestions.
Abbreviations

FES, free energy surface

PES, potential energy surface

MD, molecular dynamics
 Received May 8, 2002.
 Accepted July 18, 2002.
 Copyright © 2002, The National Academy of Sciences
References
 ↵
 ↵
 Huber T.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Bolhuis P. G.
 ↵
 ↵
Kevrekidis, I. G., Theodoropoulos, C. & Gear, C. W. (2002) Comput. Chem. Eng., in press.
 ↵
Weinan, E. & Engquist, B. (2002) arXiv ePrint archive, http://arxiv.org/abs/physics/0205048.
 ↵
 Allen M. P.
 ↵
 Hünenberger P.
 ↵
 ↵
 Geissler P. L.
 ↵
 Bolhuis P. G.
 ↵
 ↵
 Kirkpatrick S.
 ↵