Computing Free Energy Surfaces

Estimating $p(\xi)$ from Samples: Ergodicity¶

Computing an FES from Eq. eq:FES requires estimating the marginal equilibrium distribution $p(\xi)$, typically estimated from molecular dynamics or Monte Carlo trajectories. This estimate assumes ergodicity, hence implying that sufficiently long trajectories sample phase space according to equilibrium weights, making time averages equivalent to ensemble averages. Formally, for an observable $O(\mathbf{r})$:

where $f(\mathbf{r}) = \frac{e^{-\beta U(\mathbf{r})}}{Q_{NVT}}$ is the canonical distribution, $U(\mathbf{r})$ is the potential energy, and $Q_{NVT}$ the configurational partition function. In the specific case of a collective variable $\xi(\mathbf{r})$, the ergodic hypothesis ensures that the empirical distribution of $\xi$ obtained from sampling, converges to equilibrium in the long-time limit:

When systems feature metastable states separated by high barriers, unbiased simulations on practical timescales often fail to achieve ergodic sampling, leading to incomplete estimates of $p(\xi)$; enhanced sampling methods address this by accelerating rare transitions, ensuring thorough exploration of relevant states, and yielding reliable FESs that capture both thermodynamic stability and kinetic accessibility.

Recovering Ergodicity via Biased Sampling: Foundations¶

Biased enhanced sampling methods aim to efficiently achieve ergodicity by introducing an artificial bias potential $V(\xi)$ that modifies the system Hamiltonian,

to alter the equilibrium distribution $p(\xi)$ into a biased one $p_b(\xi)$. Reweighting techniques, often rooted in Zwanzig’s free energy perturbation Zwanzig, 1954, enable the recovery of unbiased distributions. Building on these principles, methods such as Umbrella Sampling, Blue Moon Ensemble, and Metadynamics enable the computation of FESs from biased ensembles; for a comprehensive overview, see Ref. Hénin et al., 2022.

Free Energy Perturbation and Thermodynamic Integration¶

The term Free Energy perturbation (FEP) describes a series of methods originating in the thermodynamic perturbation (TP) one, which evaluates the free-energy change when switching from a reference Hamiltonian $H_0$ to a perturbed one $H_1$ via Zwanzig’s formula Zwanzig, 1954,

where the subscript 0 is a reminder that the average has to be taken using the unperturbed Hamiltonian $H_0$ to integrate the equations of motions. The formula is exact for any value of $H_1$, so the term “perturbation” is somewhat misleading. However, the sampling is efficient only when configurations sampled under $H_0$ overlap well with those favored by $H_1$. Stratifying the jump into neighboring states $\xi_i\to \xi_{i+1}$ gives multistep TP,

The principle of TP is at the core of numerous methods described in the following sections.

Thermodynamic Integration¶

Thermodynamic integration (TI) is in effect an application of multistep TP in the limiting case of infinitesimally small changes:

Eq. eq:TI can be rigorously derived from Eq. eq:multiTP using the first term of the cumulant expansion Kubo, 1962

$\ln\left\langle e^{\epsilon \delta s}\right\rangle \simeq \left\langle \epsilon \right\rangle \delta s + O(\delta s^2)$, or by integrating $\partial F/\partial s = -(kT/Z_{NVT}) \partial Z_{NVT}/\partial s$.

The TI formula has a precise interpretation as the integral of the generalized mean force; however, one needs to take care to ensure the quasi-adiabatic limit. In fact, the early single-configuration TI (SCTI) implementation replaced the ensemble average by a single instantaneous value taken while changing $s$ continuously but showed hysteresis Mitchell & McCammon, 1991 at finite switching rates $ds/dt$, a clear sign of non-equilibrium. Multiconfiguration TI (MCTI) Straatsma & McCammon, 1991 addresses this issue by calculating proper ensemble averages at each value of $s$, yielding per-window statistical errors, and is also an embarrassingly parallel algorithm. Despite the development of more efficient methods, the flaw of implementing a non-equilibrium sampling has sparked new interest in SCTI in the context of Jarzynski’s identity and steered molecular dynamics Hummer, 2001Park & Schulten, 2004, see Section sec:single-molecule-spectroscopy on single molecule spectroscopy.

When the perturbation affects an internal CV rather than an interaction parameter, the system has to be kept at or in the vicinity of the chosen value of $\xi=s$. If the configuration change is driven by a $\xi$-dependent perturbation in the form of a restraint $V(q,\xi)$, typically a harmonic potential, the free energy of the unrestrained system is recovered by unbiasing the TI result with TP of the restraint Straatsma & McCammon, 1991:

Another way of looking at the TI when the RC is a function of atomic coordinates instead of a parameter of the Hamiltonian, for example $\xi(\mathbf{r}) = \left| \mathbf{r}_2 - \mathbf{r}_1 \right|$. In this case the PMF along the RC in terms of the conditional probability $p(s)$ is written as

The free energy difference takes the form

where the conditional average is defined by

This can easily be recovered as the stiff limit of a harmonic restraint. In the case of holonomic constraints being imposed rather than restraints, correct unbiasing requires more complex approaches described in Sec.sec:BM.

Umbrella Sampling¶

Umbrella sampling (US), introduced by Torrie and Valleau in the 1970s Torrie & Valleau, 1977Kästner, 2011, is one of the earliest enhanced-sampling methods to compute an FES along a CV $\xi$. US is a static bias method, i.e., the biasing potential does not change in time. The central idea of US is to overcome the poor sampling of high-energy regions by introducing a ---usually harmonic--- restraint, which localizes the sampling around a reference value $\xi_i^{\mathrm{ref}}$:

In the context of TP, $V_i(\xi)$ can be seen as the perturbation introduced in the physical Hamiltonian of the system in each of the $i$ simulations.

In US, a series of such biased simulations (“windows”, Fig. fig:biasedC) is performed to span the full range of $\xi$ (see Fig. fig:biasedD). Each window produces a biased distribution $p_i^b(\xi)$. The unbiased distribution in window $i$ is formally recovered as

from which the free energy (or potential of mean force, PMF) is obtained using Eq. eq:FES. Because each window only samples a narrow range of $\xi$, the problem of stitching the windows together arises: the free-energy offsets between windows $F_i$ are not known directly.

Several post-processing strategies exist, all rooted in FEP and TI concepts Roux, 1995. The Weighted Histogram Analysis Method (WHAM) Kumar et al., 1992 is the de facto standard post-processing tool in umbrella sampling and determines the free-energy offsets self-consistently by minimizing statistical error, merging all biased histograms into a single distribution (see also Souaille & Roux, 2001Hub et al., 2010. The Multistate Bennett Acceptance Ratio (MBAR) method provides an equivalent, statistically optimal alternative to WHAM Shirts & Chodera, 2008Tan et al., 2012.

An alternative that avoids explicit offset estimation is umbrella integration (UI) Kästner & Thiel, 2005, which reconstructs the mean force directly from the biased distribution,

yielding the PMF by integration over $\xi$. This way, UI makes explicit the proximity between US and thermodynamic integration methods, effectively implementing an FES estimator similar to Eq. eq:TI_gen We note that UI inspired the development of Mean Force Integration, which can be applied to the post-processing of sampling gathered with both static and adaptive biasing methods Marinova & Salvalaglio, 2019Bjola & Salvalaglio, 2024F. et al., 2024.

Overall, umbrella sampling can be viewed as a restrained-sampling framework analyzed either by histogram reweighting (WHAM/MBAR) or by mean-force integration (UI); both approaches converge to the same PMF given sufficient sampling.

Over the years US has inspired many extensions. Adaptive US iteratively refines the bias toward uniform sampling Mezei, 1987, while local elevation US introduces a history-dependent bias reminiscent of early metadynamics Huber et al., 1994Hansen & Hünenberger, 2010Laio & Parrinello, 2002. Multidimensional formulations Bartels & Karplus, 1997Kästner, 2009 treat coupled coordinates, and hybrid schemes, such as self-learning adaptive variants Wojtas-Niziurski et al., 2013 or replica-exchange US Jiang et al., 2012, enable efficient, massively parallel implementations. These developments confirm US remains a versatile and widely used framework for free-energy calculations.

Constrained reaction coordinates: the Blue Moon Ensemble¶

Often, holonomic constraints, where the constraining force acts along $\partial\xi/\partial \mathbf{r}$ Goldstein et al., 2014Van Kampen & Lodder, 1984, are used as an alternative to harmonic restraints to enforce a specific value of the CV $\xi(\mathbf {r})=s$ using algorithms like SHAKE Ryckaert et al., 1977. In this case, one needs to apply an unbiasing factor $Z_{\xi}= \sum_i \left|\partial {\xi}/\partial \mathbf r_i\right|^2 / m_i$ to correct for momentum loss along the constraint Carter et al., 1989. Combining this factor with the contribution from the integration over (all, including the unbiased) kinetic degrees of freedom yields the Blue Moon ensemble configurational formula (so named because it helps sample events that happen “once in a blue moon”)

where $M$ is the mass-metric tensor appearing in Eq. eq:metric-factor-origin. In this formulation, the RC must be one of the generalized coordinates used to describe the configurations. Sprick and Ciccotti derived an equivalent expression that requires only the constraint force magnitude along with a curvature correction Sprik & Ciccotti, 1998. The case of multidimensional reaction constants is discussed in Ref.Ciccotti, 1991.

Adaptive Bias Methods: Metadynamics and derived methods¶

A key ingredient of US is the use of a fixed set of windows chosen a priori, with convergence depending on their overlap. A complementary class of methods instead builds the bias adaptively as the simulation progresses. Metadynamics (MetaD), introduced by Laio and Parrinello Laio & Parrinello, 2002, is one of the most influential adaptive-bias algorithms for reconstructing FESs Barducci et al., 2011Laio & Parrinello, 2002Bussi & Laio, 2020Valsson et al., 2016. Its central idea is to discourage a molecular simulation from revisiting previously explored regions of the space of CVs $\xi$; in doing so, metadynamics enhances the fluctuations along $\xi$, thus speeding up the sampling Valsson et al., 2016. A MetaD simulation evolves under a time-dependent bias potential $V(\xi,t)$ that is incrementally constructed as the trajectory progresses. At regular time intervals, a small repulsive Gaussian hill of height $w$ and width $\sigma$ is deposited at the instantaneous CV value $\xi(t)$:

This “computational sand-filling” Barducci et al., 2011 progressively raises the free-energy of ensembles of configurations visited during sampling, allowing the system to escape local minima and visit new regions of phase space. In the long-time limit, the accumulated bias offsets the underlying free-energy surface $F(\xi)$ up to an additive constant, so that $V(\xi,t \to \infty)\approx -F(\xi)$. Once this condition is reached, the biased dynamics samples a uniform probability distribution in CV space, effectively restoring ergodicity (see Sec. sec:Theory).

MetaD is conceptually related to other history-dependent approaches, such as the local elevation method Huber et al., 1994. In particular, MetaD shares with local elevation the general principle of discouraging revisits to previously explored regions of collective variable space, while differing in its formulation and bias-update protocol. Moreover, compared with earlier mean-force-based schemes such as the Adaptive Biasing Force (ABF) method Comer et al., 2015Darve & Pohorille, 2001Darve et al., 2002, metadynamics constructs the bias from local visitation history rather than explicit force estimates. Despite its simplicity and general applicability, MetaD suffers from a few drawbacks: the continual deposition of Gaussians can lead to systematic overshooting of the free-energy surface; convergence depends on the choice of Gaussian height and width; and the time dependence of $V(\xi,t)$ complicates rigorous reweighting. These issues motivated the development of statistically controlled variants that address these shortcomings Barducci et al., 2008Valsson et al., 2016Bussi & Laio, 2020.

Well-tempered metadynamics Barducci et al., 2008, is the most popular variant of metadynamics that addresses several of these issues. For instance, WTMetaD introduces a smooth tempering of the bias deposition rate to achieve self-limiting convergence. This is achieved by imposing a dynamic evolution of the Gaussian heights that decreases exponentially with the local value of the bias already accumulated:

where $\gamma > 1$ is the bias factor. At the beginning of a simulation, the Gaussian bias height is $w_0$; as $V(\xi,t)$ grows, the added bias diminishes, so that $V$ asymptotically approaches a smoothed version of the underlying FES as $V(\xi,t\to\infty) = -({1-1/\gamma}) F(\xi)$. The bias factor $\gamma$ tunes the trade-off between exploration (large $\gamma$) and accuracy (small $\gamma$). The stationary distribution sampled by WTMetaD in the long-time limit is no longer flat but well-tempered, i.e. $p(\xi)\propto {e^{-\beta F(\xi)/\gamma}}$. The sampling along $\xi$ therefore occurs at an effectively elevated temperature $T_\text{eff}=\gamma T$, which thus enhances barrier crossing while maintaining a known, analytically recoverable bias. WTMetaD has, de facto, become the standard formulation implemented in modern software (e.g., PLUMED Tribello et al., 2014, GROMACS Lindahl et al., 2022, LAMMPS Thompson et al., 2022) because it provides controlled convergence, improved statistical efficiency, and the possibility to monitor the flattening of the free-energy landscape on-the-fly.

One can combine WTMetaD with static biases (e.g., harmonic restraints, walls, or custom static biases) in a straightforward, additive manner @Awasthi2016, Limongelli et al., 2013, @MFI_Bjola. If the WTMetaD bias is deposited sufficiently slowly - so that transition states remain effectively bias-free - the infrequent metadynamics framework allows barrier-crossing times to be rescaled to recover physical rate constants Tiwary & Parrinello, 2013Salvalaglio et al., 2014Palacio-Rodriguez et al., 2022. The bias potential accumulated during a MetaD or WTMetaD simulation modifies the underlying probability distribution, so that direct estimates of the free energy $F(\xi)$ from the bias $V(\xi,t)$ are inherently time-dependent. Several formulations have been proposed to obtain time-independent free-energy estimators, which recover $F(\xi)$ or the distribution of other observables from a trajectory sampled under the effect of a time-dependent bias Tiwary & Parrinello, 2015Bonomi et al., 2009Marinova & Salvalaglio, 2019Giberti et al., 2019Ono & Nakai, 2020.

Recent developments have further generalized ideas originating from metadynamics by recasting bias construction in probabilistic or variational terms, giving rise to methods such as Variationally Enhanced Sampling (VES) Valsson & Parrinello, 2014, On-the-fly Probabilty Enhanced Sampling (OPES) Invernizzi et al., 2020, and GAMBES Debnath & Parrinello, 2020 which differ in how the target distribution is enforced. VES does so by minimizing a KL-based variational functional. In contrast, OPES achieves it via an explicit on-the-fly estimate of the marginal biased probability density and a reweighting-based update of the bias; GAMBES, instead, reconstructs the bias from a Gaussian Mixture approximation of the sampled probability distribution. We refer the reader interested in additional details to the comprehensive review of Henin et al. Hénin et al., 2022.

Free Energy Profiles and Equilibrium Constants¶

Free energy profiles and PMFs link probabilities to partition functions, allowing computation of equilibrium constants for binding reactions. Care must be taken when connecting to standard free energies, as equilibrium constants in statistical mechanics depend on the reference concentration and are dimensionless, while experimental $\Delta G^\circ$ values refer to $c^\circ=1$ M. To avoid confusion, here we refer to the (concentration-based) mass action equilibrium constant using the symbol $K_c$. For association reactions $R+L\rightleftharpoons RL$, the mass action law gives

but this holds only in the thermodynamic limit. In finite simulation boxes containing one receptor and one ligand, the equilibrium constant is instead expressed through the ratio of Boltzmann probabilities for bound and unbound states De Jong et al., 2011:

where $c$ is the box concentration, and volume change effects are neglected. Writing $p(RL)/p(R+L)=Q_{RL}/Q_{R+L}$, the corresponding mass action equilibrium constant becomes $K_c=VQ_{RL}/Q_{R+L}$. Roux proposed an equivalent form using a delta function to fix the ligand in bulk Woo & Roux, 2005Roux et al., 2008. Gilson et al. derived the same constant from activities for a rigid ligand,

where the $8\pi^2$ term arises from rotational integration, $\sigma_i$ are symmetry numbers, and $Q_S$ the solvent partition function. The contribution coming from volume changes $\Delta \bar V_{LR} = V_{LR} -V_L -V_R$ is spelled out explicitly, and the infinite dilution identity $Q_RQ_L/ Q_S = Q_{R+L}$ has been used. This formulation underpins the double--decoupling method (DDM), in which the ratio of partition functions is obtained by decoupling the ligand at the binding site and in the bulk. Boresch et al. introduced a minimal set of six restraints (one distance $r$, two angles $\theta_i$, and three dihedral angles $\phi_i$) to restore the free energy cost of decoupling. The correction is expressed in terms of the respective equilibrium distance $r_0$, angles $\theta_{i,0}$, and force constants $K_r$, $K_{\theta_i}$ and $K_{\phi_i}$, as

The DDM thus provides a general alchemical framework - encompassing confine-and-release variants - whose rigor relies on the proper application and unbiasing of restraints Gilson et al., 1997Bian et al., 2025.

Reaction Paths on Free Energy Landscapes: Transition Path Sampling¶

As mentioned in Sec.sec:kinetics, FEPs can be used to extract useful information about reaction kinetics, but care must be taken in presence of complex landscapes. A complementary class of approaches focuses on determining representative transition pathways by optimizing geometric or variational principles. Across these families of methods, algorithms inspired by minimum‐action geometry - such as Nudged Elastic Band and its extensions [Jónsson et al. (1998); Henkelman & Jónsson (2000);wang2016automated; Ásgeirsson et al. (2021)] as well as finite-temperature string approaches [weinan2002string; weinan2005finite; E et al. (2007)] and Onsager-Machlup formulations Onsager & Machlup, 1953Adib, 2008Faccioli et al., 2006 - contrast with flux-based kinetic formalisms such as the Bennett-Chandler theory [Ruiz-Montero et al. (1997); Maragliano et al. (2006); vanden2010transition] in that the former determine representative transition pathways by extremizing geometric or variational principles, whereas the latter quantify reaction rates through time-correlation functions and dynamical recrossing statistics.

These geometric and variational approaches provide valuable insight into likely pathways, but a complete statistical description of rare events requires sampling full reactive trajectories, as achieved by Transition Path Sampling (TPS). For a Markovian dynamics, the probability density of a trajectory $\{\mathbf{r}\}=(\mathbf{r}_0,\mathbf{r}_1,\dots,\mathbf{r}_t)$ is

where $\rho(\mathbf{r}_0)$ is the equilibrium distribution and $p(\mathbf{r}_{i+1}|\mathbf{r}_i)$ the conditional propagator Bolhuis et al., 2002. The probability that a trajectory initiated in $A$ reaches $B$ at time $t$ is

with

where $\mathcal{Z}_A$ is the probability of a trajectory to start from $\mathbf{r}_0$, regardless of where it ends at time $t$, and $\mathcal{Z}_{AB}$ that to start in $r_0$ and end in $r_t$. These probabilities are expressed in terms of path integrals, effectively functional integrals over the set of all possible paths Feynman, 1948Kleinert, 2009Seifert, 2012. The rate constant, for example, follows from the long-time derivative Bolhuis et al., 2002,

and the problem of evaluating the rate $k_{AB}$ reduces to sampling the path ensemble in the same conceptual way as equilibrium properties are obtained from the Boltzmann distribution.

One can generate trajectories that sample the biased probability of starting in region $A$ and ending in region $B$ without recrossing in a straightforward manner. Starting from one reactive path, new trial trajectories $\{\mathbf{r}'\}$ are generated by small stochastic modifications and integrated forward and backward in time, accepting or rejecting the new path $\{\mathbf{r}'\}$ using a Metropolis scheme Metropolis & Ulam, 1949Kalos & Whitlock, 2008Frenkel & Smit, 2023. In the microcanonical ensemble or other generalized ensembles with extended Hamiltonians, the acceptance probability becomes remarkably simple Dellago et al., 2002 as all reactive proposals are accepted. Quantities like the rate $k_{AB}$ are not directly accessible as they are computed with a normalization factor over the path starting in $A$ and ending at any point. However, TPS provides direct access to the transmission function $\kappa(t) = k(t)/k_{tst}$ via $\kappa(t) = \langle \dot h_B(t) \rangle_{AB} / \langle \dot h_B(0) \rangle_{AB}$. The rate constant can be accessed by supplementing TPS with a US in the TPE Dellago et al., 1999. Sampling the TPE enables a straightforward computation of the committor function $p_B(\mathbf r)$ (see Sec.sec:committor). With TPS one can easily identify transition-state configurations by launching short trajectories from frames along the generated paths and selecting those for which $p_B\simeq{0.5}$.

Related path-sampling approaches include Transition Interface Sampling (TIS) Van Erp et al., 2003Van Erp & Bolhuis, 2005 Forward Flux Sampling (FFS) Allen et al., 2009, and Markov State Models (MSMs) Prinz et al., 2011Chodera & Noé, 2014.

TIS introduces a hierarchy of interfaces between $A$ and $B$, estimating the rate as the product of conditional crossing probabilities and the initial flux. FFS instead propagates trajectories forward across interfaces, making it suitable for non-equilibrium or irreversible systems. MSMs reconstruct long-time kinetics from short unbiased trajectories by estimating transition probabilities between metastable states, with rates and committors derived from the transition matrix.

Further refinements include Precision, S-, and Aimless Shooting Grünwald et al., 2008Menzl et al., 2016Peters & Trout, 2006, and hybrid schemes combining TPS with metadynamics or replica exchange Borrero & Dellago, 2016Van Erp, 2007. A detailed overview is given in Bolhuis & Swenson, 2021.

FESs from experiments: single molecule force spectroscopy¶

Single-molecule force spectroscopy (SMFS) methods---such as AFM and optical or magnetic tweezers---enable experimental probing of molecular free-energy landscapes. A molecule is stretched between a surface and a probe applying controlled force or displacement, with the extension $s$ serving as the reaction coordinate. The resulting mechanical response encodes the potential of mean force $F(s)$, recoverable from nonequilibrium work measurements via statistical--mechanical relations. A central theoretical foundation of this connection is Jarzynski’s equality Jarzynski, 1997, which relates nonequilibrium work to equilibrium free-energy differences:

where $W(t)$ is the external work performed along the pulling trajectory, and $\langle \dots \rangle$ denotes an average over many repetitions of the process. Remarkably, Eq. eq:Jarzynski holds even for transformations driven arbitrarily far from equilibrium, providing a formal bridge between dynamical experiments and equilibrium thermodynamics. Hummer and Szabo Hummer & Szabo, 2005, extended this result to the reconstruction of a full free-energy profile along the molecular extension coordinate $s$, obtaining

where $V(s,t)$ is the time-dependent external potential applied during pulling. For instance, in optical-tweezer or AFM experiments, a harmonic trap of stiffness $k_s$ is displaced at constant velocity $v$, $V(s,t)=\tfrac12 k_s(s-vt)^2$.

Equation eq:HummerSzabo enables reconstruction of the equilibrium free-energy profile $F_0(s)$ from nonequilibrium pulling trajectories. In the adiabatic limit, the work $W$ equals the reversible $\Delta F$, but real experiments require averaging over many trajectories to account for stochastic fluctuations. Liphardt et al. Liphardt et al., 2002 first validated this approach, showing RNA hairpin unfolding reproduced equilibrium free energies with sub-$kT$ accuracy. These combined experimental--computational approaches reveal intermediates and barriers, extending free-energy surface reconstruction into the experimental domain through nonequilibrium work theorems.