The adaptive biasing force method: everything you always wanted to know but were afraid to ask

Abstract

In the host of numerical schemes devised to calculate free energy differences by way of geometric transformations, the adaptive biasing force algorithm has emerged as a promising route to map complex free-energy landscapes. It relies upon the simple concept that as a simulation progresses, a continuously updated biasing force is added to the equations of motion, such that in the long-time limit it yields a Hamiltonian devoid of an average force acting along the transition coordinate of interest. This means that sampling proceeds uniformly on a flat free-energy surface, thus providing reliable free-energy estimates. Much of the appeal of the algorithm to the practitioner is in its physically intuitive underlying ideas and the absence of any requirements for prior knowledge about free-energy landscapes. Since its inception in 2001, the adaptive biasing force scheme has been the subject of considerable attention, from in-depth mathematical analysis of convergence properties to novel developments and extensions. The method has also been successfully applied to many challenging problems in chemistry and biology. In this contribution, the method is presented in a comprehensive, self-contained fashion, discussing with a critical eye its properties, applicability, and inherent limitations, as well as introducing novel extensions. Through free-energy calculations of prototypical molecular systems, many methodological aspects are examined, from stratification strategies to overcoming the so-called hidden barriers in orthogonal space, relevant not only to the adaptive biasing force algorithm but also to other importance-sampling schemes. On the basis of the discussions in this paper, a number of good practices for improving the efficiency and reliability of the computed free-energy differences are proposed.

Figure 1

Upper row: (A) original two-dimensional, double-well potential displayed as level sets; (C) time trajectory of the first coordinate x in the stochastic process, showing oscillations between the two metastable wells. Lower row: (B) level sets of the same potential biased by the free energy associated with the transition coordinate ξ(x,y) = x; (D) time trajectory of the transition coordinate x in the adaptively biased dynamics, showing no metastability.

Figure 2

Upper row: (A) the original two-dimensional potential is zero inside the hourglass shape, and +∞ outside; (C) time trajectory of the first coordinate x in the stochastic process, showing oscillations between the two metastable wells. Lower row: (B) free energy along the transition coordinate ξ(x,y) = x, featuring a purely entropic barrier; (D) time trajectory of the transition coordinate x in the free-energy-biased dynamics, showing no metastability.

Figure 3

Stratification strategies for the translationally invariant toy model of a tagged water molecule diffusing in a bulk aqueous medium. The transition path spans 20 Å and is handled in a single window (A), in two 10 Å windows (B), in four 5 Å windows (C), and in eight 2.5 Å windows (D). For each stratification strategy, a potential of mean force calculation is carried out until the root-mean-square deviation with respect to the accurate zero free-energy profile is less than 0.1 kcal/mol.

Figure 4

Coordinate dependence of sampling and the system force distribution in reversible folding of deca-alanine. (A) Samples in each bin for a 10 ns adaptive biasing force calculation on the domain ξ ∈ [4, 32] Å (black curve) and a 100 ns calculation on the domain ξ ∈ [4, 32] Å (red curve). Note the logarithmic scale on the vertical axis. (B) Distribution of the ξ-component of the instantaneous system force for different ranges of ξ. (C) Standard deviation of the ξ-component of the instantaneous system force as a function of ξ.

Figure 5

Coordinate dependence of correlations in the system force in reversible folding of deca-alanine. (A) Autocorrelation functions of the random component of the instantaneous system force for different ranges of ξ. The functions are normalized by the variance ⟨ΔF_ξ²⟩_i so that correlation at t = 0 is unity. (B) Correlation time for ξ-ranges in panel A.

Figure 6

Propagation of error in the mean force to the error of free-energy differences. (A) Mean system force on ξ for a 10 ns adaptive biasing force calculation, with error bars determined according to eq 22. (B) Error of free-energy differences A(ξ) – A(z_a) with the reference position z_a = 14 Å. This position is denoted by a violet dashed line.

Figure 7

N_step × Var(t) as a function of the number of molecular dynamics steps, N_step. This quantity was calculated for the permeation of K⁺ through a transmembrane hexametric channel of a peptaibol, trichotoxin in a window along the normal to the membrane spanning the range between z equal to −15 and −9 Å. z = 0 is located in the middle of the membrane. Note that for N_step between 0.6 × 10⁶ and 2.3 × 10⁶, N_step × Var(t) differs from its average value in this range by no more than 10%. The inset: the number of force samples along z in the window as a function of N_step.

Figure 8

Common free-energy landscape featuring parallel valleys, collinear to the transition coordinate, ξ. These valleys are separated by substantial free-energy barriers in the direction ζ, orthogonal to ξ. ζ can be interpreted as a slow degree of freedom coupled to ξ and hampering progression along the latter direction. Excessively stratified reaction pathways preclude spontaneous crossing of high barriers, typically ΔA₁. Wider windows should allow diffusion toward values of ξ, where the barrier separating valleys in the direction ζ is smaller, typically ΔA₂.

Figure 9

Improving the rate of convergence with multiple-walker strategies. (A) Potentials of mean force obtained for different multiple-walker strategies and total simulation times of 8 ns. For reference, the black curve shows the result for a total simulation time of 128 ns. “Single” refers to a single long 8 ns simulation, and “independent” denotes the result of combining force samples for 16 walkers after they ran independently for 0.5 ns. For the curves marked “shared” and “selection”, force samples were synchronized among walkers every 20 ps. “Selection” also included the application of walker selection, as described in the text, on the same interval. (B) Potentials of mean force obtained for different multiple-walker strategies and total simulation times of 128 ns.

Figure 10

Three-dimensional adaptive biasing force calculation. (A) Three reference structures of deca-alanine. The root-mean-square deviation of selected atoms from their positions in each of the three reference structures defines each of the three transition coordinates. (B) Free-energy isosurfaces as a function of the three transition coordinates. The violet, pink, and gray surfaces contain all points for which the free energy is less than 5, 10, and 20 kcal/mol, respectively. In all cases the tick marks represent a root-mean-square deviation of 1 Å. The minimum value of the free energy, which occurs at (RMSD_α RMSD₃₁₀ RMSD_ω) = (0.3 2.6 2.7) Å, is defined to be zero. (C) Another view of the surface shown in panel B, with representative structures shown for the two energetically favorable regions.

Figure 11

Combining thermodynamics and kinetics for reversible unfolding of deca-alanine. (A) Comparison of the system force as determined by adaptive biasing force and that calculated from the Bayesian scheme on the same trajectory. (B) Diffusivity as a function of the end-to-end distance of deca-alanine for different observation time intervals Δt.