Exploration vs Convergence Speed in Adaptive-Bias Enhanced Sampling.

In adaptive-bias enhanced sampling methods, a bias potential is added to the system to drive transitions between metastable states. The bias potential is a function of a few collective variables and is gradually modified according to the underlying free energy surface. We show that when the collective variables are suboptimal, there is an exploration-convergence tradeoff, and one must choose between a quickly converging bias that will lead to fewer transitions or a slower to converge bias that can explore the phase space more efficiently but might require a much longer time to produce an accurate free energy estimate. The recently proposed on-the-fly probability enhanced sampling (OPES) method focuses on fast convergence, but there are cases where fast exploration is preferred instead. For this reason, we introduce a new variant of the OPES method that focuses on quickly escaping metastable states at the expense of convergence speed. We illustrate the benefits of this approach in prototypical systems and show that it outperforms the popular metadynamics method.

Introduction

Molecular dynamics has become a valuable tool in the study of a variety of phenomena in physics, chemistry, biology, and materials science. One of the long-standing challenges in this important field is the sampling of rare events, such as chemical reactions or conformational changes in biomolecules. To simulate effectively such systems, many enhanced sampling methods have been developed. An important class of such methods is based on an adaptive-bias approach and includes adaptive umbrella sampling,[1] metadynamics (MetaD),[2,3] and the recently developed on-the-fly probability enhanced sampling (OPES).[4−6] Adaptive-bias methods operate by adding to the system’s energy U(R) an external bias potential V = V(s), which is a function of a set of collective variables (CVs), s. The CVs, s = s(R), depend on the atomic coordinates R and are meant to describe the slow modes associated with the rare event under study. They also define a free energy surface (FES), F(s) = −1/β log P(s), where β = (kBT)−1 is the inverse Boltzmann factor and P(s) the marginal s distribution, P(s) ∝ ∫e–βδ[s – s(R)]dR. The bias is periodically updated until it converges to a chosen form. A popular choice is to have it exactly offset the underlying FES, V(s) = −F(s), so that the resulting s distribution is uniform.

The main limitation of adaptive-bias methods is that finding good collective variables is sometimes difficult and a bad choice of CVs might not promote the desired transitions in an affordable computer time. In practical applications, one generally has to live with suboptimal CVs[7] that still can drive transitions but do not include some of the slow modes. In this case, applying a static bias cannot speed up the slow modes that are not accounted for, and thus transitions remain quite infrequent. It is sometimes possible to achieve a faster transition rate using a rapidly changing bias, which can push the system out of a metastable state through a high free energy pathway, different from the energetically favored one. However, unless one wishes to deal explicitly with out-of-equilibrium statistics,[8−10] it is not possible to obtain reliable information about the system while the bias changes in a nonadiabatic fashion. To estimate the FES and other observables, one must let the adaptive-bias method approach convergence, and as the bias becomes quasi-static, transitions inevitably become less frequent.

We refer to this situation as an exploration–convergence tradeoff that every adaptive-bias enhanced sampling method has to deal with when suboptimal CVs are used. Some methods, like OPES, focus more on quickly converging to a quasi-static bias potential and thus obtaining an efficiently reweighted FES, while others, like metadynamics, focus more on escaping metastable states and exploring the phase space. We will demonstrate this qualitative difference in some prototypical systems. For simplicity, in the paper, we only consider the well-tempered variant of metadynamics,[3] but in the Supporting Information, we provide examples that use the original nontempered MetaD[2] and other popular variants, such as parallel-bias MetaD.[11]

We propose here a variant of OPES, named OPES-explore, which focuses on rapid exploration, rather than on fast convergence. It shares many features with the original OPES, and is designed to be an easy-to-use tool requiring few input parameters. To this end, we also introduce an adaptive bandwidth algorithm that can be used in both OPES variants and further reduces the number of input parameters that need to be specified. The detailed description of the adaptive bandwidth algorithm is given in the Supporting Information. All OPES simulations presented make use of this algorithm.

OPES Method

The enhanced sampling method OPES works by adding an adaptive-bias potential to the energy of the system so as to modify the Boltzmann probability distribution into the desired target one. Most adaptive-bias methods aim at sampling uniformly the CV space, but it has been shown that choosing a different target distribution could be advantageous.[12,13] There are two different classes of target distributions that can be sampled with OPES: metadynamics-like and replica-exchange-like. We will consider here only the former type, introduced in ref (4), but the interested reader can find information about OPES for replica-exchange-like sampling in ref (5).

To define a metadynamics-like target distribution, one has to choose a set of collective variables, s = s(R). As stated in the introduction, the unbiased marginal probability along such CVs is P(s) ∝ ∫e–β δ[s – s(R)]dR, where U(R) is the potential energy. The target distribution is then defined by requiring a specific marginal probability distribution over the CVs, ptg(s). Consequently, the desired bias potential is written asso that ∫e–β[  δ[s – s(R)]dR ∝ ptg(s). A typical choice for ptg(s) is the well-tempered distribution[3]where the bias factor γ > 1 controls how much the original distribution is smoothed out. In the limit of γ = ∞, one targets a uniform distribution.

The core idea of OPES is to update self-consistently the estimate of the probability distributions and of the bias potential in an on-the-fly fashion similar to self-healing umbrella sampling.[14] The estimate of the unbiased probability is obtained via a weighted kernel density estimation so that at step n, one haswhere the weights w are given by w = eβ, and the Gaussian kernels G(s, s′) = h exp[−1/2(s – s′)∑–1(s – s′)] have a diagonal covariance matrix ∑ = σ2 δ and fixed height . The number of kernels to represent P(s) would grow linearly with simulation time, but this is avoided thanks to an on-the-fly kernel compression algorithm,[15] as described in detail in the supporting information of ref (4). The compression algorithm also allows for the bandwidth of the kernels to shrink over time, as the effective sample size Neff( = (∑w)2/∑w2 grows. The idea is to start with a coarse estimate of P(s) and then refine it as more data are available. The kernel bandwidth of the i-th CV at step n iswhere d is the total number of CVs.

The instantaneous bias is based on the probability estimate P(s); following eq and using the approximation pWT(s) ∝ [P(s)]1/γ, one haswhere ϵ is a regularization term that limits the maximum possible absolute value of the bias potential and Z can be understood as normalization of P(s) over the CV space thus far explored, Ω

The integral is calculated approximately as a sum of P over the compressed kernels, as explained in the supporting information of ref (4). The intuitive idea is that new kernels are added to the compressed representation only when a new region of CV space is sampled (otherwise they are merged with the existing ones), thus the explored CV-space volume |Ω| is approximately proportional to the total number of compressed kernels.

The introduction of the Z term is one of the key innovations of OPES. In similar methods, once a new metastable state is found, one often sees a dramatic increase in the exit time, compared to the first one[16] (see Supporting Information, Figure S8). This exit time problem is present also when the CVs are optimal and should not be confused with the exploration–convergence tradeoff that is the primary concern of this paper. Other convergence-focused methods introduce extra parameters to tackle this problem, for example in transition-tempered metadynamics[17] prior knowledge of the position of all metastable states is required. Instead, OPES avoids the exit time problem by taking into account the expansion of the CV space via the Z term, which allows the bias to adjust more quickly when a new CV-space region is sampled.[4]

At the start of an OPES simulation, only a handful of parameters need to be chosen, namely the initial kernel bandwidth, the pace at which the bias is updated, and the approximate FES barrier height that needs to be overcome. From this last information, a prescription is given to automatically set the values γ and ϵ. The number of parameters can be reduced even further if one uses, as we shall do here, the adaptive bandwidth algorithm discussed in the Supporting Information.

Exploratory OPES Variant

We present now a new OPES variant called OPES-explore which, compared to the original OPES formulation, leads to a faster exploration of the phase space at the cost of slower convergence. We have recalled that the Z term allows OPES to quickly adapt the bias when a metastable state is found in a previously unexplored region of CV space. However, if the CVs used are suboptimal, it may happen that a new metastable state is found in an already explored s region.[17,18] In such a case, the Z term remains constant and is therefore ineffective in accelerating the exit time. Instead, to encourage a rapid exit, one would need a method that allows the bias to significantly change shape again. Fortunately, it is possible to achieve this exploratory behavior simply by making a minimal change to the OPES protocol, which gives rise to the OPES-explore variant.

In formulating OPES-explore, we restrict ourselves to the case of using as a target the well-tempered distribution, ptg(s) = pWT(s), eq . In OPES, the bias is expressed as a function of P(s), the on-the-fly estimate of the unknown equilibrium distribution P(s). At the beginning of the simulation, this estimate is not reliable but it improves over time and converges in a self-consistent way. In OPES-explore instead, one builds the bias starting from the on-the-fly estimate of the distribution that is being sampled in the biased simulationwhere s is the CV value sampled at step k. As the simulation converges, pWT(s) approaches the target well-tempered distribution pWT(s). Thus, analogously to Section , we use the approximation P(s) ∝ [pWT(s)]γ and write the bias according to eq where ϵ and Z have been added for the same reasons as in eq . We notice that the expressions in eqs and 7, which define the probability estimates used in the two OPES schemes, converge, respectively, to P(s) and pWT(s) only within the self-consistent scheme where the simulation runs with a bias that is updated on-the-fly according to eqs and 8, respectively. Both OPES variants are applications of the more general eq , but OPES estimates on-the-fly P(s) and uses it to calculate the bias, while OPES-explore does the same but with pWT(s) ∝ [P(s)]1/γ.

The free energy surface as a function of the CVs can be estimated in two distinct ways, either directly from the probability estimate, F(s) = −γ 1/β log pWT(s), or via importance sampling reweighting, e.g., using a weighted kernel density estimation

In standard OPES, these two estimates are equivalent, while in OPES-explore (similarly to MetaD), they can differ significantly in the first part of the simulation until they eventually converge to the same estimate.

In Figure , we contrast an OPES and OPES-explore simulation of alanine dipeptide in a vacuum, which has become a standard test for enhanced sampling methods. Both simulations have the same input parameters and use the adaptive bandwidth scheme described in the Supporting Information. The bias is initially quite coarse, but the width of the kernels reduces as the simulation proceeds and the details of the FES are increasingly better described. It can clearly be seen that the OPES-explore variant employs fewer kernels compared to the original OPES. This is due to the fact that in OPES-explore, the kernel density estimation is used for pWT(s) ∝ [P(s)]1/γ that is a smoothed version of P(s), and thus requires fewer details. This more compact representation can be useful, especially in higher dimensions, where the number of kernels can greatly increase despite the compression algorithm. However, as a drawback, it can result in a less accurate bias estimate, especially for large values of γ.

Figure 1

Time evolution of a typical simulation of alanine dipeptide in a vacuum using the two OPES variants with the dihedral angles ϕ and ψ as CVs. For each method, the compressed kernels are shown on the left with the point size indicating the adaptive bandwidth, and the corresponding free energy estimate F(ϕ, ψ) on the right. (a) In the original OPES, kernels make up the unbiased distribution estimate P(ϕ, ψ) and F(ϕ, ψ) = −1/β log P(ϕ, ψ), while (b) in OPES-explore, kernels make up the sampled distribution estimate pWT(ϕ, ψ) and F(ϕ, ψ) = −γ 1/β log pWT(ϕ, ψ). All F(ϕ, ψ) are shifted to have zero minimum. Notice how OPES-explore requires fewer kernels and visits higher FES regions.

Fewer Transitions Can Lead to Better Convergence

The difference in performance between OPES and OPES-explore cannot be judged from the alanine dipeptide example because in this case, the CVs chosen are extremely efficient. To highlight the difference between the two methods, we study a simple two-dimensional model potential that is known as the Müller potential,[20] see Figure a, using the x coordinate as a collective variable. This is a clear example of suboptimal CV since it can discriminate the metastable states but not the transition state.

Figure 2

(a) Müller potential energy surface, U(x, y). (b) Free energy surface along the x coordinate, F(x), with and without the addition of the bias potential V(x) = −(1 – 1/γ)F(x), where γ = 20. (c) Potential energy modified by the bias potential, U(x, y) + V(x). All the energies reported are shifted to have zero minimum. It can be seen that despite the almost flat profile along x, the transition region between the states remains at high energy.

For two-dimensional systems, the free energy along the CV, F(x), can be computed precisely with numerical integration, Figure b. From F(x), the free energy difference between the two metastable states can be calculated as

While it is possible to distinguish better the two states using also the y coordinate, this does not result in a significant difference in the ΔF value (see Supporting Information, Sec. S3). On the other hand, x does a poor job of identifying the transition state, which is around x ≈ −0.7 and y ≈ 0.6, and not at x ≈ 0 as it would seem from F(x). As a consequence, it is not possible to significantly increase the transition rate between the states using a static bias that is a function of x only.

To show this, we consider the effect of adding to the system the converged well-tempered bias V(x) = −(1 – 1/γ)F(x), with γ = 20. In Figure b, we can see the effect of the bias on the FES along x, which becomes almost completely flat. However, when we consider the full 2D landscape, as shown in Figure c, we can see that such a bias does not really remove the barrier between the two states. From the height of the barrier at the transition state, one can roughly estimate that adding V(x) improves the transition rate by about one order of magnitude. Nevertheless, transitions remain quite rare, around one of every 106 uncorrelated samples (see Supporting Information, Sec. S3).

We want to compare the two OPES variants and well-tempered metadynamics in this challenging setting, where CVs are suboptimal and the total simulation time is not enough to reach full convergence. This type of situation is not uncommon in practical applications, and it is thus of great interest. Given enough time, all of the methods considered will converge to the same bias potential and sample the same target distribution, but we shall see that before reaching this limit, they behave very differently.

Figure a shows a typical run of the Müller potential obtained by biasing the x coordinate with OPES, OPES-explore, or MetaD. As a simple way to visualize the evolution of the bias, we also report in Figure b, the ΔF estimate obtained directly from the applied bias, using F(x) = −(1 – 1/γ)−1V(x) in eq . We can see a qualitative difference between OPES and the other two methods.

Figure 3

Typical simulations of the Müller potential using different methods for biasing the x coordinate. Given more time, the three methods will converge to the same bias potential and will sample the same target distribution. (a) Shows the trajectory along the CV and (b) shows the corresponding ΔF, eq , calculated using the FES estimate obtained directly from the applied bias, F(x) = −(1 – 1/γ)−1V(x). The correct ΔF value is highlighted by a blue stripe 1 kBT thick.

OPES reaches a quasi-static bias that is very close to the converged one but samples a distribution that is far from the well-tempered one, where the two basins would be about equally populated. On the other hand, the x distribution sampled by OPES-explore is closer to the target well-tempered one, but its bias is far from converged and makes ample oscillations around the correct value. Metadynamics behaves similarly to OPES-explore. This is the exploration–convergence tradeoff described in Section . Since the CV is suboptimal, even when using the converged bias V(x) to see a transition occur, one has to wait for an average number of steps τ ≈ 106, which is more than the total length of the simulation. However, it is possible to greatly accelerate transitions using a time-dependent bias that forces the system into higher energy pathways, which are not accessible at equilibrium.

In OPES-explore, the bias is based on the estimate of the sampled probability pWT(s) and pushes to make it similar to the almost flat well-tempered target. This means that to have a quasi-static bias, pWT(s) should both be almost flat and not change significantly as the simulation proceeds. Clearly, this cannot happen unless the simulation is longer than τ, otherwise, most of the time would be spent in the same basin and pWT(s) would be far from flat. On the contrary, in OPES, the bias is based on the reweighted estimate P(s), and thus it can reach a quasi-static regime even before properly sampling the target distribution.

In Figure a, we show the ΔF estimate averaged over 25 independent runs, all starting from the main basin x < 0. We can see that on average, OPES provides the best ΔF estimate at any n in spite of the fact that it induces far fewer transitions. In fact, most of the time, only one full back-and-forth transition is observed (see Supporting Information, Figure S5a). One should notice that after a single transition, the ΔF estimate is far from being accurate (see Figure b) but, since the bias quickly becomes quasi-static, it is possible to collect equilibrium samples and reliably reweight them, and the average estimate becomes more accurate the more simulations are run. Instead, in OPES-explore and MetaD, despite starting from independent initial conditions, the runs are highly correlated due to the transitions being mostly driven by the strong changes in the bias rather than the natural fluctuations of the system. As a consequence of this, a systematic error is present in the average estimate, even if ΔF is further averaged over time, to remove the oscillatory behavior of OPES-explore and MetaD. Such a systematic error depends on the characteristic of the system and the chosen CVs and is hard to predict whether it will be relevant or small. Nevertheless, one can be sure that it reduces over time as the bias converges.[24]

Figure 4

Estimate of the free energy difference ΔF for the Müller potential obtained by averaging 25 independent runs for each biasing method. The standard deviation is also shown for each estimate. Given more time, all of these estimates will converge to the correct ΔF. All simulations start from the main basin, x < 0, but with different initial conditions. (a) Shows the estimate obtained directly from the applied bias, as shown in Figure b, while (b) shows the corresponding estimate obtained via reweighting. For metadynamics, two different reweighting schemes are considered, bias-offset[21,22] and last-bias reweighting.[19,23] The correct ΔF value is highlighted by a blue stripe 1 kBT thick.

Estimates of ΔF using different reweighting schemes are shown in Figure b. For OPES and OPES-explore, eq has been used, while for MetaD, we consider two of the most popular reweighting schemes, namely last-bias reweighting[19,23] and bias-offset reweighting.[21,22] As expected, the reweighting estimate of OPES is virtually identical to the direct estimate obtained from the bias, while for the other two methods, the two estimates differ. The reweighting of OPES-explore has very small statistical uncertainty, which further highlights the presence of a systematic error in the free energy difference estimate. Like others before us,[22,25,26] we observe empirically that the last-bias reweighting for MetaD tends to be in agreement with the direct estimate, even when the simulation is far from converged, while the bias-offset reweighting provides a very unreliable estimate if the MetaD bias has not reached a quasi-static regime and the initial part of the simulation is not discarded. Once again, it must be noted that the simulations considered here are not fully converged, otherwise all of the different estimates of the various methods would have yielded the correct result, without systematic errors. However, for most practical purposes, they behave very differently; thus, it is important to choose between an exploration-focused or a convergence-focused enhanced sampling method, depending on the specific aim of the simulation.

Sometimes Exploration Is What Matters

In the example of the previous section, it was shown that OPES converges to a quasi-static bias faster than OPES-explore and provides more accurate FES estimates. However, FES estimation is not the only goal of an enhanced sampling simulation. In complex systems where good CVs are not available, convergence can remain out of reach; still, one might be interested in exploring the phase space and finding all of the relevant metastable basins. In such a situation, OPES-explore can be a useful tool.

We consider here alanine tetrapeptide in a vacuum as a test system, as in ref (4). It has three ϕ dihedral angles, each of them can change from positive to negative values and vice versa with a relatively low probability. This leads to 23 = 8 distinct metastable basins, each corresponding to a different combination of ϕ angles signs, as shown in Figure . Here, we are not interested in estimating the FES, but rather we want to compare the ability of different methods to explore this space and discover all metastable states.

Figure 5

Eight metastable basins of alanine tetrapeptide in a vacuum sampled via OPES-explore by biasing the three ψ angles, a suboptimal set of CVs. Each basin is identified by the sign of the three ϕ angles for a total of 23 possible combinations. The most stable basin has ϕ1, ϕ2, ϕ3 < 0, while the least stable basin has ϕ1, ϕ2, ϕ3 > 0.

Figure shows the number of explored basins averaged over 10 independent simulations for each enhanced sampling method. The simulations in the top panel (Figure a) bias the ϕ angles, V = V(ϕ1, ϕ2, ϕ3), which are good CVs, while in the bottom (Figure b), the suboptimal ψ angles are used, V = V(ψ1, ψ2, ψ3). In all methods, the exploration time increases approximately by two orders of magnitude when suboptimal CVs are used (please note the horizontal logarithmic scale). As expected, OPES and OPES-explore have similar exploration speeds when using good CVs, while with suboptimal CVs, OPES struggles to find all of the metastable basins. This is because the same region of CV space might correspond to two different metastable basins or to a basin and a transition state, as for the Müller potential.[18,19] In this situation, the previously estimated bias must change considerably for the simulation to escape the current metastable state quickly.

Figure 6

Exploration time of the eight metastable basins of alanine tetrapeptide over 100 ns. The lines are an average of over 10 independent runs for each method, showing the total number of visited basins. (a) Shows bias as a function of the three ϕ angles, V = V(ϕ1, ϕ2, ϕ3), while (b) shows the three ψ angles are used, V = V(ψ1, ψ2, ψ3). See the Supporting Information for results with different input parameters and other MetaD variants, such as parallel-bias MetaD.[11]

The exploration speed of MetaD depends critically on the input parameters and requires a trial-and-error tuning. We report here only the outcome of MetaD simulations in which a standard choice of the input parameters has been made. As can be seen in Figure , in these simulations, the exploration speed is roughly one order of magnitude slower than that of OPES-explore. However, the performance of MetaD simulations can be improved using different settings, as shown in the Supporting Information. In the Supporting Information, we also report and briefly discuss results obtained with nontempered metadynamics,[2] adaptive-Gaussians metadynamics,[23] and parallel-bias metadynamics.[11] None of these MetaD variants significantly improve the exploration speed, and some make it even worse.

Finally, in the Supporting Information, Sec. S5, we show how a preliminary OPES-explore run can be combined with a multithermal OPES simulation[5,6] to efficiently sample alanine tetrapeptide and reach a converged FES, even without explicitly biasing the ϕ angles.

Conclusions

With the help of model systems, we show that there is an exploration–convergence tradeoff in adaptive-bias methods when suboptimal CVs are used. This tradeoff should not be confused with the exit time problem, which is present also with optimal CVs, and is discussed in Section and refs (16, 17). Contrary to the exit time problem, the exploration–convergence tradeoff cannot be solved. It is an intrinsic limitation of CV-based adaptive-bias methods, which comes as a consequence of suboptimal CVs. We believe that the best way to handle this tradeoff is to have separate methods that clearly focus on one or the other aspect so that they can be used depending on the application. In a convergence-focused method, the bias soon becomes quasi-static to allow for accurate reweighting and free energy estimation. However, with suboptimal CVs, this leads to a slow transition rate and a long time is required to sample the target distribution. As discussed, even if one knows the true F(s) and directly applies the converged bias, one would not obtain a faster exploration. In an exploration-focused method, it is possible to improve the exploration speed by letting the bias change substantially even in a CV region that has already been visited. While this may increase the number of transitions, it comes at the cost of a slower convergence.

The original OPES method focuses on fast convergence to provide an accurate estimate of the free energy surface and reweighted observables. As a consequence, it is very sensitive to the quality of the CVs (see e.g., Figure a), and any improvement in the CVs results in a clear acceleration of the transition rate. This is a particularly useful property when developing machine learning-based CVs, and in fact, OPES has already been used several times in this context.[27−32]

In other situations, improving the CVs may first require a better exploration of the phase space.[28,33−35] Furthermore, one may be interested simply in exploring the metastable states of a system rather than estimating an accurate FES.[36−39] For this reason, we have introduced a variant of the OPES method, OPES-explore, which focuses on quickly sampling the target distribution and exploring the phase space.

We have shown that also well-tempered metadynamics is an exploration-focused method. One of the main advantages of OPES-explore over MetaD is that it is easier to use since it requires fewer input parameters and it has a more straightforward reweighting scheme (but more advanced ones can also be used[25,40]). Another important difference between the two methods is that OPES-explore, similarly to OPES, by default provides a maximum threshold to the applied bias potential, thus it avoids unreasonably high free energy regions. To obtain the same effect with MetaD, one typically has to define some ad hoc static bias walls by trial and error. This last feature of OPES-explore has been recently leveraged by Raucci et al. to systematically discover reaction pathways in chemical processes.[41]

Finally, we should clarify that OPES-explore, just as metadynamics, might not be able to exit any metastable state if the CVs are too poor,[19,42] and its improved exploration capability can only be harnessed if the CVs are close enough to the correct ones to make such transitions possible. The speed and the small number of input parameters of OPES-explore are extremely helpful for quickly testing several candidate CVs to find out which can drive transitions and discard the bad ones.

We believe that OPES-explore is an important addition to the OPES family of methods and will become a useful tool for researchers as it pushes forward the trend for more robust and reliable enhanced sampling methods.