Rare events, free energies and enhanced sampling

The timescale gap

A molecular dynamics simulation advances in steps of about one femtosecond, $10^{-15}$ s, because it must resolve the fastest motions of atoms, bond vibrations with periods of 10 to 100 fs. With a good GPU code and a few hundred thousand atoms, a day of computing yields perhaps a few hundred nanoseconds. Many of the processes that matter in chemistry and materials science are vastly slower. A vacancy in a metal at room temperature hops once every microsecond to millisecond; a small protein folds in microseconds to seconds; a crystal nucleates from a supercooled liquid after seconds to years; a catalytic turnover takes milliseconds; ions in a battery cathode move between sites on timescales from nanoseconds to hours. Between the femtosecond time step and these events lie nine to fifteen orders of magnitude.

The reason is not that these processes are intrinsically complicated. It is that they are rare. The system spends almost all its time vibrating inside a basin of the energy landscape, a stable or metastable state, and only very occasionally gathers enough energy in the right coordinates to cross the barrier into another basin. The crossing itself is fast, typically picoseconds; the waiting time before it is long, because the probability of reaching the top of a barrier of height $\Delta E$ is proportional to $e^{-\Delta E/k_BT}$. At room temperature, $k_BT \approx 0.026$ eV, so a barrier of 0.5 eV suppresses the crossing probability by a factor of about $10^{-8}$, and 1 eV by about $10^{-17}$. A brute-force simulation spends nearly all its effort re-sampling the same basin and almost never sees the event of interest.

This subject is about the methods that overcome this problem. Some locate barriers and paths on the energy surface directly, without dynamics; some compute rates from the properties of the barrier; some compute free energies by sampling cleverly; some bias the dynamics so that barriers are crossed often, and then remove the effect of the bias; and some, more recently, use machine learning to find the right coordinates or to generate samples directly. All of them start from the geometry of the energy landscape, which is the subject of this topic.

Basins and stationary points

The potential energy surface $U(\mathbf r_1, \dots, \mathbf r_N)$, introduced in the chapter on potential energy surfaces, is a function on a $3N$-dimensional configuration space. Its minima are the stable and metastable structures: a crystal, a vacancy at a particular site, a conformer of a molecule, a reactant or a product. Around each minimum lies its basin of attraction, the set of configurations from which steepest descent leads to that minimum. A dynamical trajectory at low temperature spends long periods inside one basin, rattling near its minimum, and occasionally passes into a neighbouring basin.

As the chapter on gradients showed, stationary points, where the gradient vanishes, are classified by the eigenvalues of the Hessian matrix of second derivatives $H_{ij} = \partial^2 U / \partial x_i \partial x_j$. At a minimum, all eigenvalues are positive (apart from zero eigenvalues associated with overall translation and rotation). At a first-order saddle point, exactly one eigenvalue is negative: the point is a maximum along one direction and a minimum along all others. Saddles of higher order, with several negative eigenvalues, are maxima in several directions and are not important for transitions, because there is always a lower route around them.

The first-order saddle connecting two basins is the transition state of the transition between them, the mountain pass in the hiking-map picture. Its energy above the starting minimum is the activation energy that enters the rate, and its vibrational frequencies (the positive Hessian eigenvalues) enter the prefactor, as the chapter on reaction rates explains. Locating transition states is therefore the first quantitative task in studying a rare event.

Minimum energy paths

Between two minima connected by a transition state runs a minimum energy path (MEP): the path that, at every point, is a minimum of the energy in all directions perpendicular to it. Equivalently, the force on the path, $-\nabla U$, has no component perpendicular to the path:

In mass-weighted coordinates, the steepest-descent path from the transition state down to the two minima is called the intrinsic reaction coordinate (IRC), a concept introduced by Kenichi Fukui in 1970 and central to computational chemistry. The MEP is the natural one-dimensional description of a transition at low temperature: the energy profile along it shows the barrier, any intermediate minima, and the shape of the reaction.

The Müller–Brown surface, a two-dimensional function introduced by Klaus Müller and Leo Brown in 1979 as a test case for path methods, makes these ideas concrete. It has three minima: a deep one at upper left (energy −146.7), a second at lower right (−108.2), and a shallow intermediate in between (−80.8), connected by two saddle points at −40.7 and −72.2. The MEP from the deepest minimum to the lower-right one is curved and passes through the intermediate; a straight line between the endpoints would climb high over a ridge.

Finding the path: elastic bands and strings

Double-ended methods find the MEP between two known minima by optimising a discretised path, a chain of $N$ configurations called images. The simplest idea, the elastic band of Ron Elber and Martin Karplus (1987), minimises the sum of the image energies plus springs between neighbouring images. It fails in an instructive way: the spring forces have components perpendicular to the path that pull it off the MEP and make it cut corners, and the true forces have components along the path that drag images down into the minima, leaving the barrier region poorly resolved.

The nudged elastic band (NEB) method of Hannes Jónsson, Greg Mills and Karsten Jacobsen (1998) fixes both problems by projection. Each image feels only the perpendicular component of the true force, which moves it towards the MEP, and only the parallel component of the spring force, which keeps the images evenly spaced:

The climbing-image NEB of Graeme Henkelman, Blas Uberuaga and Hannes Jónsson (2000) adds one modification: the highest image is freed from the springs and has the parallel component of its force inverted, so that it moves uphill along the path and downhill in all other directions. It converges to the saddle point itself, giving the barrier exactly rather than interpolated between images.

The string method of Weinan E, Weiqing Ren and Eric Vanden-Eijnden (2002) reaches the same goal differently. The images evolve by the perpendicular component of the force alone, with no springs, and after every step they are redistributed at equal arc length along the curve through them. The string is a curve in configuration space parametrised by arc length, and the reparametrisation keeps the images evenly spaced; its finite-temperature version extends the idea to free-energy landscapes.

Try the plain elastic band with a stiff spring: the path is pulled straight across the ridge and the predicted barrier is far too high. With a weak spring, images slide down into the minima. NEB, with the same number of images, settles onto the curved valley floor through the intermediate minimum. Climbing-image NEB then lifts its highest image onto the first saddle, reproducing the exact barrier of $−40.67 − (−146.70) = 106.0$ energy units.

In practice NEB calculations are run with codes such as ASE, VASP (with the VTST tools of the Henkelman group) or LAMMPS. The following ASE script computes the barrier for a vacancy hop in copper with the EMT potential, with five intermediate images and a climbing image:

import numpy as np
from ase.build import bulk
from ase.calculators.emt import EMT
from ase.optimize import BFGS, FIRE
from ase.mep import NEB

base = bulk("Cu", "fcc", a=3.6, cubic=True).repeat((3, 3, 3))      # 108 atoms
d = base.get_distances(0, range(len(base)), mic=True)
j = int(np.argsort(d)[1])                     # a nearest neighbour of site 0
initial = base.copy(); del initial[0]         # vacancy at site 0
final = base.copy(); final.positions[j] = base.positions[0]; del final[0]   # neighbour has hopped in
for a in (initial, final):
    a.calc = EMT(); BFGS(a, logfile=None).run(fmax=0.01)

images = [initial] + [initial.copy() for _ in range(5)] + [final]
for im in images[1:-1]:
    im.calc = EMT()
neb = NEB(images, climb=True)
neb.interpolate(mic=True)                     # straight-line initial path
FIRE(neb, logfile=None).run(fmax=0.02)
E = [im.get_potential_energy() for im in images]
print(f"barrier: {max(E) - E[0]:.3f} eV")      # 0.776 eV with EMT

The result, 0.78 eV, is close to the measured migration energy of vacancies in copper, about 0.7 eV, although agreement at this level from a simple potential is partly fortunate. With DFT or an MLIP, the same script gives the corresponding barrier for that model; only the calculator changes.

Finding saddles from one end

When only the starting minimum is known, as in exploring which transitions are possible from a given defect configuration, single-ended methods search for saddles without a final state. Eigenvector following and its relatives, such as the partitioned rational-function optimisation of Jon Baker (1986), use the Hessian: they climb uphill along the eigenvector of the lowest eigenvalue while minimising along all others. Computing the full Hessian is expensive for large systems, so the dimer method of Graeme Henkelman and Hannes Jónsson (1999) estimates the lowest-curvature direction using only forces, from a pair of closely spaced images (the dimer) that is rotated to minimise its energy; the dimer then translates uphill along that direction. The activation–relaxation technique (ART nouveau) of Gerard Barkema and Normand Mousseau follows a similar strategy.

Single-ended searches launched in many random directions from a minimum map out the possible transitions and their barriers. This is the basis of adaptive kinetic Monte Carlo (Henkelman and Jónsson, 2001), which finds saddles on the fly, computes rates from them, and advances a simulation event by event over timescales of seconds or longer. Global methods such as basin hopping (David Wales and Jonathan Doye, 1997), which alternate random perturbations with local minimisation, explore landscapes to find the lowest minima, and Wales’s discrete path sampling builds networks of minima and saddles from which kinetics can be computed.

From energy to free energy

The landscapes and paths of this chapter are properties of the potential energy, and they describe transitions at low temperature well. At finite temperature, a system does not follow the MEP exactly; it fluctuates around it, and the relevant barrier is a free-energy barrier, which includes the entropy of all the motions transverse to the path. A narrow pass and a wide valley of the same height are crossed at different rates, because the wide valley offers more ways through. In liquids, biomolecules and soft matter, where many configurations of similar energy coexist, entropy often dominates, and the energy landscape alone can be misleading: there may be no single saddle point at all, but a broad region of configuration space through which transitions pass.

The free-energy landscape is defined along chosen collective variables, as the chapter on ensembles introduced: $F(\xi) = -k_BT\ln P(\xi)$. Choosing collective variables that capture the slow, rare motion is one of the central problems of this subject, discussed in the chapter on collective variables; computing free energies along them is the subject of the free-energy and enhanced-sampling chapters that follow.

References and further reading

Rates and waiting times

When a system sits in a basin A and occasionally escapes to basin B, and the escapes are rare compared with the time it takes to forget how it got into A, every moment spent in A carries the same small chance of escape. The process is then memoryless: the probability of still being in A after time $t$ decays exponentially, and the waiting times are exponentially distributed.

The rate constant is what experiments measure and what simulations must predict. Because events are rare, it is usually far too small to measure by watching a simulation, and the theories in this chapter estimate it from properties of the barrier instead.

Arrhenius and Eyring

In 1889 Svante Arrhenius, analysing the temperature dependence of the inversion of sucrose, proposed that rate constants follow

so that a plot of $\ln k$ against $1/T$ is a straight line whose slope gives the activation energy $E_a$. He interpreted $E_a$ as the energy a molecule must acquire to react, and the exponential as the fraction of molecules that have it, which is the Boltzmann factor of statistical mechanics. The prefactor $A$ remained empirical.

In 1935 Henry Eyring, and independently Meredith Evans and Michael Polanyi, derived the prefactor from statistical mechanics, founding transition-state theory (TST). Eyring’s result, in its thermodynamic form, is

With $\kappa = 1$, every 5.7 kJ/mol of activation free energy at room temperature reduces the rate tenfold. The equation underlies how chemists interpret reaction rates and how computational chemists turn computed barriers into rate predictions.

Transition-state theory from first principles

TST can be derived cleanly in one dimension, and the derivation shows exactly what it assumes. Let the reaction coordinate be $x$, with the reactant basin at $x < x^\ddagger$ and the dividing surface at $x^\ddagger$, the top of the free-energy barrier. The rate is the flux of trajectories crossing the dividing surface from left to right, divided by the population of the reactant:

The derivation makes two assumptions. First, the reactant is in equilibrium, so its population and the density at the barrier top follow the Boltzmann distribution; this holds when relaxation within the basin is fast compared with escape. Second, every trajectory that crosses the dividing surface forwards is counted as a reaction, even if it turns round and recrosses a moment later. The second assumption means that TST is an upper bound on the true rate.

For a solid or a molecule, where the barrier is a saddle point on the potential energy surface, expanding the energy to second order around the minimum and around the saddle gives harmonic transition-state theory, derived by George Vineyard in 1957 for diffusion in crystals:

The vacancy hop in copper computed with NEB in the previous chapter has a barrier of 0.78 eV with the EMT potential. With a prefactor of $5 \times 10^{12}$ s⁻¹, harmonic TST predicts about 0.4 hops per second at 300 K, $1.4 \times 10^6$ at 600 K and $2 \times 10^8$ at 900 K: from essentially frozen to rapid diffusion over a few hundred kelvin, the hallmark of an activated process. No direct MD simulation could see the room-temperature hops, but a single NEB calculation predicts their rate.

Recrossing and the transmission coefficient

Real trajectories that reach the dividing surface do not all continue to the product. Collisions with surrounding atoms can turn them round, so that they recross the surface one or more times before settling on either side. The true rate is the TST rate multiplied by a transmission coefficient $\kappa \le 1$, the fraction of forward crossings that lead to a committed reaction. The Bennett–Chandler method computes it from simulation: start many short trajectories at the dividing surface with Maxwell–Boltzmann velocities, follow each for a short time, and measure the fraction that ends in the product, weighted by the initial velocity. The rate is then $k = \kappa\,k_{\text{TST}}$, with $k_{\text{TST}}$ obtained from a free-energy calculation of $F(x)$. This two-step “reactive flux” approach was introduced by Charles Bennett (1977) and David Chandler (1978) and is still a standard route to rates in solution.

Because TST is an upper bound, the dividing surface can be moved to minimise the predicted rate; this variational transition-state theory, developed from ideas of Eugene Wigner and James Keck and systematised by Donald Truhlar and others, gives the best TST estimate and reduces the recrossing correction.

Kramers’ theory: the role of friction

In 1940 Hendrik Kramers asked how a particle escapes from a well when it is coupled to a heat bath through friction $\gamma$ and random forces, the Langevin dynamics described in the chapter on thermostats. His answer shows that the rate depends on friction in a non-monotonic way. At very low friction, the particle exchanges energy with the bath too slowly to reach the barrier energy, and the rate is limited by energy diffusion, increasing with friction. At high friction, the particle reaches the barrier easily but its motion across the barrier top becomes diffusive, and it recrosses many times; the rate falls as $1/\gamma$. In between, the rate peaks near the TST value: the Kramers turnover. In the moderate-to-high friction regime the rate is

The figure below measures rates directly. Three hundred independent particles move by Langevin dynamics in a double well with a barrier of 5 (in units where the temperature can be varied); each time a particle passes from the left well to the right one, a transition is counted, and the rate is the number of transitions divided by the total time spent in the left well. Record points at several temperatures to build an Arrhenius plot and compare with TST (dashed) and Kramers (solid). Then change the friction: at low friction the measured rate approaches TST; at high friction it falls well below, as Kramers predicted.

The straight lines on the Arrhenius plot have a slope equal to the barrier height, and their intercepts differ by the transmission coefficient, which the readout gives as $\kappa = k_{\text{Kr}}/k_{\text{TST}}$. At very low friction, below the turnover, Kramers’ high-friction formula itself breaks down and the measured rate falls again, a regime that matters for gas-phase reactions and for weakly damped degrees of freedom.

Rates from simulation in practice

When barriers are low enough that transitions occur within simulation time, rates can be counted directly, as in the figure, or computed from mean first-passage times. When they are not, the strategies of the following chapters apply. A free-energy profile along a good collective variable gives the TST rate, and short trajectories from the barrier top give $\kappa$. Path-sampling methods (the chapter on paths) sample the rare reactive trajectories themselves and compute rates without a dividing surface. Markov state models (their own chapter) build rates from many short simulations. And for crystals, where transitions are hops between well-defined sites, harmonic TST with NEB barriers feeds kinetic Monte Carlo, which advances a lattice of sites event by event, choosing each event with probability proportional to its rate and advancing the clock by an exponentially distributed waiting time. This is how diffusion, ordering and growth in materials are simulated over seconds to hours.

The quantity that connects all these approaches is the committor $p_B(\mathbf x)$: the probability that a trajectory started at configuration $\mathbf x$, with random velocities, reaches B before A. It is 0 in A, 1 in B, and $\tfrac12$ on the true transition state ensemble. The committor is the ideal reaction coordinate, and much of modern rare-event theory, from transition path theory to machine-learned collective variables, can be seen as attempts to approximate it.

References and further reading

Describing a transition with a few numbers

A transition in a system of thousands of atoms involves, in principle, all $3N$ coordinates. Yet chemists describe a reaction by a bond that forms and another that breaks, a biophysicist describes folding by how many native contacts have formed, and a materials scientist describes crystallisation by the size of the largest crystalline cluster. Each of these is a collective variable (CV): a function $\xi(\mathbf x)$ of the atomic coordinates, chosen to capture the progress of the process while ignoring the irrelevant details.

Collective variables are the language of almost every enhanced-sampling method. Free energies are computed as functions of them; biases are applied along them; and the quality of the result depends on how well they describe the transition. Choosing them is the most important, and often the most difficult, decision in a rare-event study.

A catalogue of collective variables

The simplest CVs are geometric. A distance between two atoms or centres of mass describes binding, dissociation and bond breaking; an angle or dihedral describes conformational change, as the backbone dihedrals $\phi$ and $\psi$ do for the classic test system alanine dipeptide. A coordination number counts how many atoms of one group are within a cutoff of another, made smooth and differentiable with a switching function,

and describes solvation, adsorption and the number of bonds an atom has. The root-mean-square deviation (RMSD) from a reference structure measures how far a molecule is from, say, its folded state.

For more complex transitions, CVs are built from many atoms. Path collective variables (Davide Branduardi, Francesco Gervasio and Michele Parrinello, 2007) describe progress along a reference path of frames, $s$, and distance from that path, $z$. Steinhardt bond-order parameters (Paul Steinhardt, David Nelson and Marco Ronchetti, 1983), built from spherical harmonics of the bond directions around each atom, distinguish liquid-like from crystal-like environments and different crystal structures; averaged or counted over atoms, they track nucleation and growth. Coordination-based order parameters, the number of hydrogen bonds, the electric dipole of a region, the volume of a cavity, and many others are in use. The PLUMED library, a plugin that adds CVs and enhanced-sampling methods to most major MD codes (GROMACS, LAMMPS, AMBER, CP2K, OpenMM and others), implements hundreds of them.

The free-energy surface

Once a CV is chosen, the free-energy surface along it is defined from the probability distribution of the CV in equilibrium, as introduced in the chapter on ensembles:

Because $F(s)$ is a marginal of the full distribution, it depends on the choice of CV in a way the energy does not. Two consequences are worth knowing. First, a coordinate transformation contributes a Jacobian term: the free energy as a function of the distance $r$ between two particles includes $-2k_BT\ln r$, because the volume of the shell at distance $r$ grows as $r^2$; even two non-interacting particles have a free energy that decreases with distance. Comparing free energies along different CVs, or converting a distance profile into a binding free energy, requires accounting for such terms. Second, free-energy differences between states are independent of the CV only when the states are defined by regions that do not depend on it; barrier heights depend on the CV and are physically meaningful only for good CVs.

The derivative of the free energy is the mean force, the average force along the CV with all other coordinates in equilibrium. For a CV that is a combination of Cartesian coordinates it is simply the average of $-\partial U/\partial \xi$ at fixed $\xi$; in general, as Emily Carter, Giovanni Ciccotti, James Hynes and Raymond Kapral showed in 1989 for simulations with the CV held fixed by a constraint (the blue-moon ensemble), it contains geometric correction terms. Integrating the mean force along the CV gives $F(s)$: this is thermodynamic integration along a CV, and it is also the principle behind the adaptive biasing force method of the chapter on adaptive biasing.

What makes a collective variable good

A good CV must at least distinguish the states of interest: its values in the reactant, the product and the transition region must be different. But that is not enough. A good CV should also capture the slow degrees of freedom of the transition, so that once the CV is fixed, everything else relaxes quickly. If an important slow coordinate is left out, the free-energy surface along the CV will look smooth while the real barrier hides in the coordinate that was ignored, and biasing along the CV will fail to accelerate the transition, a problem known as hidden barriers or orthogonal degrees of freedom.

The figure shows this in two dimensions. Two states, A and B, are separated by a barrier along a diagonal direction. The collective variable $s$ is a projection onto a direction you choose. When $s$ points along the true reaction direction, the free-energy profile shows the barrier clearly and the two states project onto separate regions of $s$. Rotate $s$ away and the barrier shrinks and disappears, while the projections of A and B overlap until they become indistinguishable: a free-energy profile along that CV would suggest there is no barrier at all, although nothing about the system has changed.

The rigorous test of a reaction coordinate uses the committor $p_B$ introduced in the chapter on rates. If $s$ is a good reaction coordinate, then all configurations with the value $s^\ddagger$ at the top of the barrier should have $p_B \approx \tfrac12$. The committor test (Peter Bolhuis, Christoph Dellago, Phillip Geissler and David Chandler, around 1999–2000) samples configurations at fixed $s^\ddagger$, shoots many trajectories from each, and histograms the resulting committor values: a histogram peaked at $\tfrac12$ means a good coordinate, a flat or bimodal one means that $s$ misses something important. The test is expensive but decisive, and it motivates the machine-learning methods of the chapter on learned collective variables that approximate the committor or the slowest dynamical modes directly.

Practical choices

In practice, CVs are chosen from physical insight, tested, and refined. Start from the simplest coordinates that distinguish the states; check whether biased simulations along them actually produce transitions and whether the resulting free-energy surface converges; use two-dimensional surfaces when one CV is not enough; look at the transitions that occur for degrees of freedom that change but are not in the CV; and, where the transition is complex, use path CVs or data-driven CVs learned from simulations. The examples throughout the following chapters use one-dimensional CVs for clarity, but real applications often need two or more, and the cost of most methods grows quickly with the number of biased CVs.

References and further reading

The problem: a ratio of partition functions

Free energies are not averages of anything a simulation records. The free energy of a state is $-k_BT\ln Z$, and the partition function $Z$ is an integral over all of configuration space, which a simulation never samples completely. What can be computed are free-energy differences between two states, A and B, defined by two potential energy functions $U_A$ and $U_B$ on the same configuration space, because their ratio of partition functions can be written as an average:

The states can be two conformations, a molecule in water and in vacuum, a ligand bound and unbound, or two different molecules altogether in an alchemical calculation, where atoms are gradually switched from one element to another or made to disappear. These calculations are the basis of solvation free energies, binding affinities in drug design, the relative stability of crystal polymorphs and the chemical potentials of solids and liquids.

Free-energy perturbation

In 1954 Robert Zwanzig noticed that the ratio of partition functions is an equilibrium average over state A:

The formula is exact but treacherous. The exponential average is dominated by configurations with low $\Delta U$, which are the configurations typical of state B. If those configurations are rare in state A, the average depends on a few samples in the tail of the distribution and converges extremely slowly, with a large bias for any finite sample. The method works well only when the two states overlap in configuration space: when the configurations important for B are also visited, at least occasionally, by A. The same formula applied in the reverse direction, sampling B, gives a second estimate, $\Delta F = k_BT\ln\langle e^{\beta\Delta U}\rangle_B$, and the two usually disagree when overlap is poor.

The figure below makes this concrete with two harmonic wells, for which the exact answer is known. Each frame draws samples from both states and plots the distributions of $\Delta U$ seen in A and in B. With the states close together, the two histograms overlap and all estimates converge quickly to the exact value. Move them apart and the histograms separate: the forward and reverse estimates converge slowly from opposite sides, with large systematic errors at any feasible number of samples.

The practical remedy is to break a large change into many small ones. Introduce a coupling parameter $\lambda$ that mixes the two potentials, for example $U(\lambda) = (1 - \lambda)U_A + \lambda U_B$, run simulations at a series of intermediate values $\lambda_1, \lambda_2, \ldots$, and add up the small free-energy differences between neighbours, each of which has good overlap. This staging is used in essentially every alchemical calculation.

Thermodynamic integration

An alternative, older than perturbation, was given by John Kirkwood in 1935. Differentiating $F(\lambda) = -k_BT\ln Z(\lambda)$ with respect to $\lambda$ gives an equilibrium average, and integrating it gives the free-energy difference:

Thermodynamic integration avoids exponential averages, so it does not suffer from the same bias, but its accuracy depends on the integrand being smooth enough for quadrature. For alchemical transformations in which atoms appear or disappear, the integrand can diverge near the endpoints, because a nascent atom with a tiny interaction can overlap other atoms; soft-core potentials, which cap the interaction at short range when the atom is partly decoupled, remove this problem and are standard in GROMACS, AMBER and OpenMM.

Free-energy calculations are almost always organised in thermodynamic cycles. Free energy is a state function, so the sum of free-energy changes around a closed cycle is zero, and an unphysical but computable path can replace a physical but intractable one. The relative binding free energy of two ligands to a protein, for example, is computed by alchemically transforming ligand 1 into ligand 2 once in the binding site and once in water, rather than by simulating the slow physical binding and unbinding of each: $\Delta\Delta G_{\text{bind}} = \Delta G_{1\to2}^{\text{site}} - \Delta G_{1\to2}^{\text{water}}$. Calculations of this kind, run routinely in pharmaceutical research, typically reach errors of about 1 kcal/mol (4 kJ/mol) relative to experiment.

Optimal estimators: BAR and MBAR

Using samples from both states together gives a much better estimator than either one-sided average. Charles Bennett derived in 1976 the estimator that minimises the variance for given samples from A and B, the Bennett acceptance ratio (BAR). With $n_A$ forward values $w_F = \Delta U$ sampled in A and $n_B$ reverse values $w_R = -\Delta U$ sampled in B, the free-energy difference solves

In the figure, the BAR estimate (rose) converges fastest and most reliably, even when the one-sided estimates disagree. BAR is still limited by overlap, but it extracts the most information from whatever overlap there is. The multistate Bennett acceptance ratio (MBAR) of Michael Shirts and John Chodera (2008) generalises BAR to any number of states sampled at once, such as all the $\lambda$ windows of an alchemical calculation or the windows of umbrella sampling, and it reduces to the WHAM equations of the next chapter when the data are histogrammed. MBAR, available in the pymbar library and its successors, is the standard analysis tool for multi-state free-energy calculations.

Non-equilibrium work: Jarzynski and Crooks

All the methods above require equilibrium sampling at each state. In 1997 Christopher Jarzynski proved a remarkable result: free-energy differences can be extracted from irreversible, non-equilibrium processes. Start the system in equilibrium in state A, then change the parameter $\lambda$ from A to B at a finite rate, for example by pulling a molecule with a moving harmonic trap, and record the work $W = \int (\partial H/\partial\lambda)\,\mathrm d\lambda$ done on the system. Repeat many times. The second law says that the average work is at least the free-energy difference, $\langle W\rangle \ge \Delta F$, with equality only for infinitely slow (reversible) changes. Jarzynski showed that an exponential average of the work gives $\Delta F$ exactly, at any speed:

Gavin Crooks extended the result in 1999 to a fluctuation theorem relating the work distributions of the forward process and of its time reverse: $P_F(W)/P_R(-W) = e^{\beta(W - \Delta F)}$, so the two distributions cross exactly at $W = \Delta F$. Both results were confirmed experimentally by stretching single RNA molecules with optical tweezers (Liphardt and colleagues, 2002; Collin and colleagues, 2005).

The figure pulls a particle through a double-welled landscape with a moving harmonic trap, many times over, and histograms the work. At low speed the work distribution is narrow and centred near $\Delta F$. At high speed it shifts to larger values and broadens: much of the work is dissipated as heat. The Jarzynski estimate remains close to the exact $\Delta F$, but only because a few rare trajectories with unusually low work dominate the exponential average; at very high speeds those trajectories become so rare that a feasible number of pulls cannot find them, and the estimate is biased upwards. This is the same overlap problem as in free-energy perturbation, which is indeed the instantaneous-switching limit of Jarzynski’s equality.

In molecular simulation, Jarzynski’s equality underlies steered molecular dynamics calculations of unbinding and unfolding free energies, and non-equilibrium switching is used in alchemical calculations as an alternative to equilibrium staging. Its practical limit is always the same: the rarer the low-work trajectories, the more realisations are needed.

References and further reading

The idea: make the improbable probable

A free-energy profile $F(s) = -k_BT\ln P(s)$ along a collective variable can be computed from a histogram of $s$, but only where the simulation goes. Across a barrier of 10 or 20 $k_BT$, the probability density at the top is $e^{-10}$ to $e^{-20}$ times that in the well, so an unbiased simulation never samples it. In 1977 Glenn Torrie and John Valleau proposed adding a bias potential that makes the improbable regions probable, and then removing the effect of the bias analytically. They called the bias an umbrella because it spreads sampling over a region that the natural distribution would leave uncovered.

The commonest form today uses a series of windows: separate simulations, each with a harmonic restraint that holds the collective variable near a different centre $s_k$,

Unbiasing one window

In window $k$, the biased distribution of $s$ is related to the unbiased one by the Boltzmann factor of the known bias:

where $f_k = -k_BT\ln\langle e^{-\beta w_k}\rangle$ is an unknown constant, the free energy of switching on the bias. Each window therefore gives the free-energy profile over its own narrow range, but shifted by an unknown constant. The whole problem of analysis is to determine the constants $f_k$ so that the pieces join up, which is possible only where neighbouring windows overlap: where both sample the same values of $s$.

The weighted histogram analysis method

The weighted histogram analysis method (WHAM), developed by Shankar Kumar, John Rosenberg, Peter Kollman and colleagues in 1992 from earlier work of Alan Ferrenberg and Robert Swendsen on combining Monte Carlo histograms, finds the best combination of all windows. It writes the unbiased distribution as a weighted combination of the windows’ histograms, with weights chosen to minimise the statistical error, which leads to two equations that are solved by iteration:

The figure runs eleven windows across an asymmetric double-well profile in parallel, each sampling its biased distribution by Langevin dynamics. The bottom panel shows each window’s histogram. As soon as every window has collected some samples, WHAM is solved and the reconstruction (orange) is drawn over the true profile (blue).

With enough windows and moderate springs, the histograms overlap, WHAM converges in a few hundred iterations, and the reconstruction matches the true profile within a fraction of $k_BT$ across the barrier. Remove windows or stiffen the springs and gaps open between histograms. WHAM then has no information to connect the pieces on either side of a gap: the reconstructed profile on one side can be shifted by any amount relative to the other, and the error grows. In the opposite direction, springs that are too soft cannot hold the windows on the steep parts of the profile, where the walkers slide down into the wells and the barrier region is undersampled. The spring constant must be large enough to overcome the local slope of $F(s)$, roughly $\kappa \gtrsim |\mathrm dF/\mathrm ds|/\Delta s$ for window spacing $\Delta s$, and small enough that neighbouring windows still overlap; a useful rule of thumb is that each window’s distribution should have a width comparable to the window spacing.

MBAR, introduced in the previous chapter, solves the same problem without histogramming: each sample is used individually, with no binning error, and the equations reduce to WHAM in the limit of fine bins. Modern analyses usually use MBAR, or WHAM implementations that provide bootstrap error estimates.

Practical umbrella sampling

A typical umbrella-sampling study proceeds as follows. Generate starting configurations along the collective variable, usually from a short steered-MD pull or from an interpolated path. Run each window long enough for its distribution to converge, discarding an initial equilibration period. Check the overlap of neighbouring histograms and add windows where it is poor. Combine with WHAM or MBAR, and estimate errors by bootstrapping or by block averaging over time (trajectory analysis), remembering that consecutive samples are correlated. Finally, check convergence by comparing profiles from the first and second halves of the data.

The most dangerous pitfall is not visible in the histograms: slow degrees of freedom orthogonal to the collective variable. If, at fixed $s$, the system has two configurations separated by a barrier the window never crosses, each window samples only one of them, every histogram looks fine, and the profile is wrong. This is the hidden-barrier problem of the chapter on collective variables. Running windows from different starting configurations, or combining umbrella sampling with replica exchange between windows, helps detect and reduce it.

Umbrella sampling remains among the most widely used free-energy methods, because it is simple, robust and easy to parallelise: every window is an independent simulation. Its limitations are that the collective variable and the window placement must be chosen in advance, and that the cost grows quickly with the number of collective variables. The adaptive methods of the next chapters, metadynamics and its relatives, construct the bias on the fly instead.

References and further reading

Filling the wells

A bias built from history

Umbrella sampling needs the bias to be designed in advance. A more adaptive idea is to let the simulation build its own bias from its history: wherever the system has already been, add a little repulsive potential, so that it is gradually pushed to explore somewhere new. Once the bias has filled the well the system started in, it spills over the barrier into the next well, fills that too, and eventually the whole free-energy landscape is levelled. The accumulated bias is then a mirror image of the free energy.

Origins

Versions of this idea appeared several times. Thomas Huber, Andrew Torda and Wilfred van Gunsteren’s local elevation (1994) penalised conformations already visited in conformational searches; Helmut Grubmüller’s conformational flooding (1995) filled a potential-energy well with a Gaussian to accelerate escape; and the Wang–Landau algorithm (Fugao Wang and David Landau, 2001) built up a flat histogram in energy by adaptively modifying the density of states in Monte Carlo. In 2002 Alessandro Laio and Michele Parrinello formulated metadynamics for collective variables in molecular dynamics, in a form that has since become one of the most widely used enhanced-sampling methods in chemistry, biophysics and materials science. Laio and Parrinello described it as filling the free-energy wells with “computational sand”.

The algorithm

Depositing Gaussian hills

At regular intervals $\tau_G$, metadynamics adds a small Gaussian hill to a bias potential at the current value of the collective variable:

Metadynamics on a one-dimensional profile

The figure is a laboratory for all the adaptive methods of this part. A single collective variable diffuses by Langevin dynamics on a free-energy profile (blue). With plain MD it stays in its starting well. With metadynamics, hills accumulate (violet fill, the profile plus the bias), the system escapes, and the negative of the bias (orange dashed) approaches the true profile.

Watch the first well fill: the system rattles in its starting basin while the violet region rises; when the bias reaches the barrier height, the system crosses, and the time series shows the first transition. Thereafter it moves back and forth freely while both wells continue to fill. The estimate converges to the true profile, but not exactly. Because hills keep being added at full height forever, the bias never stops changing: it fluctuates around $-F(s)$ with an error that depends on the hill height and the deposition rate, and it overfills regions where the dynamics is slow to respond. Increase the hill height and the errors grow; make it small and convergence takes longer. Laio, Rodriguez-Fortea, Gervasio, Ceccarelli and Parrinello analysed this error in 2005, showing that it scales as the square root of the deposition rate.

Well-tempered metadynamics

Shrinking hills

In 2008 Alessandro Barducci, Giovanni Bussi and Michele Parrinello solved the convergence problem with well-tempered metadynamics (WTMetaD). The height of each new hill is reduced according to how much bias has already accumulated at that point:

The bias factor

The well-tempered bias does not flatten the landscape completely. It flattens it to a fraction $1/\gamma$ of its original height, as if the collective variable were at an effective temperature $T + \Delta T$ while all other degrees of freedom stay at $T$. The free energy is recovered as $F(s) = -\frac{\gamma}{\gamma - 1}V(s)$. The bias factor sets the compromise: large enough that barriers of a few tens of $k_BT$ become crossable (a bias factor of 10 turns a 20 $k_BT$ barrier into 2 $k_BT$), small enough not to waste time sampling extremely high free energies. Select well-tempered metadynamics in the figure and compare with the plain version: the hills shrink as the profile fills, the violet surface converges to a tilted copy of the profile rather than a flat line, and the estimate stops fluctuating.

Metadynamics of a real cluster

Two collective variables for seven atoms

The one-dimensional profile hides the hardest part of a real calculation: choosing collective variables and watching the bias act on actual atoms. The two-dimensional Lennard-Jones heptamer has four stable shapes. The hexagon, a ring of six atoms around a central one, is the ground state at $E = -12.53\,\varepsilon$; three other isomers, here called B, C and D, lie at $-11.50$, $-11.48$ and $-11.40\,\varepsilon$. A good pair of collective variables, introduced for this system in the PLUMED tutorials, describes how the coordination numbers of the atoms are distributed. Each atom’s coordination number is a smooth count of its neighbours,

The bias acts on the atoms through the chain rule. Writing $\partial\mu_2/\partial c_i = 2(c_i - \bar c)/N$ and $\partial\mu_3/\partial c_i = 3[(c_i - \bar c)^2 - \mu_2]/N$, and differentiating each switching function with respect to the distance it depends on, gives the force of the bias on every atom. The figure computes exactly this at every step, and the implementation was checked against finite differences to six significant figures.

Watching the bias work

The simulation starts in the hexagon at $k_BT = 0.12\,\varepsilon$. Turn the bias off first: the cluster vibrates, the point on the map rattles around the hexagon minimum, and in long unbiased runs the cluster spends about 98 % of its time there. Turn the bias on and the map begins to colour from the hexagon outwards. Within a few hundred $\tau$ the bias has filled the hexagon basin enough to push the cluster into one of the other isomers, and the atoms rearrange in front of you: the central atom is squeezed out, its coordination number drops from six to four or five, and $\mu_3$ falls. The point then wanders between B, C and D, the basins fill in turn, and the free-energy estimates in the lower plot settle.

How accurate are they? Four independent simulations were run to check. Two long well-tempered runs ($4\times10^4\,\tau$ each) gave free energies relative to the hexagon of $0.57$, $0.54$ and $0.71\,\varepsilon$ for B, C and D, agreeing with each other to $0.02\,\varepsilon$. Three unbiased runs of $2\times10^5\,\tau$ each, five times longer still, gave $0.54$–$0.56\,\varepsilon$ for B and $0.52$–$0.56\,\varepsilon$ for C, in agreement, but for D anything from $0.62$ to $0.99\,\varepsilon$: even $2\times10^5\,\tau$ of plain molecular dynamics visits isomer D too rarely to measure its free energy, while metadynamics converged all three in a fifth of that time. That is the case for enhanced sampling in one sentence.

The same calculation in PLUMED

The PLUMED input for this calculation, with the MD engine running the Lennard-Jones cluster in reduced units, is

# plumed.dat (reduced Lennard-Jones units)
UNITS NATURAL
c1: COORDINATIONNUMBER SPECIES=1-7 MOMENTS=2-3 SWITCH={RATIONAL R_0=1.5 NN=8 MM=16}
metad: METAD ARG=c1.moment-2,c1.moment-3 SIGMA=0.04,0.04 HEIGHT=0.012 PACE=250 BIASFACTOR=10 TEMP=0.12 FILE=HILLS
PRINT ARG=c1.moment-2,c1.moment-3,metad.bias STRIDE=100 FILE=COLVAR

COORDINATIONNUMBER with MOMENTS=2-3 computes the central moments of the coordination numbers, and the rational switching function with NN=8 and MM=16 is the one in the equation above.

Reweighting and rates

Unbiased averages

The bias distorts the sampled distribution, so ordinary averages from a metadynamics run are biased averages. Because the bias is known, they can be reweighted. For well-tempered metadynamics in its quasi-stationary regime, each frame receives a weight $e^{\beta(V(s(t), t) - c(t))}$, where $c(t)$ is a time-dependent offset that accounts for the growth of the bias; Pratyush Tiwary and Michele Parrinello gave a practical estimator for it in 2015. Reweighting allows free energies along variables that were not biased to be computed from the same simulation, for example the free energy as a function of a second order parameter in a nucleation study.

Rates from infrequent metadynamics

Metadynamics can also give rates. In infrequent metadynamics (Tiwary and Parrinello, 2013), hills are deposited slowly enough that the transition state region is never biased. The system then escapes from a well faster than in an unbiased simulation, but by passing through the same unbiased transition states. Each escape time can be rescaled by the acceleration factor $\langle e^{\beta V}\rangle$, accumulated while waiting in the well, to give an unbiased escape time, and the distribution of rescaled times should be exponential if the assumptions hold, which provides a statistical check. This approach has been used to compute residence times of drugs bound to proteins on timescales of seconds.

Variants and software

Variants

Several variants address the limitations of the basic method. Multiple-walker metadynamics (Paolo Raiteri, Laio, Gervasio, Cristian Micheletti and Parrinello, 2006) runs many replicas that deposit hills into a shared bias, parallelising the filling. Bias-exchange metadynamics (Stefano Piana and Laio, 2007) runs several replicas, each biasing a different collective variable, and exchanges configurations between them, which allows many CVs to be used without the cost of a high-dimensional bias. Parallel-bias metadynamics (Jim Pfaendtner and Massimiliano Bonomi, 2015) applies several one-dimensional biases simultaneously in a single simulation. Funnel metadynamics (Vittorio Limongelli, Bonomi and Parrinello, 2013) adds a funnel-shaped restraint for ligand-binding calculations, so that the unbound state is sampled in a finite volume and the standard binding free energy can be computed.

Running metadynamics with PLUMED

Metadynamics is implemented in PLUMED, which works with most major MD and ab initio codes. A well-tempered run biasing the distance between two atoms needs only a few lines:

# plumed.dat
d1: DISTANCE ATOMS=1,10
metad: METAD ARG=d1 SIGMA=0.05 HEIGHT=1.2 PACE=500 BIASFACTOR=10 TEMP=300 FILE=HILLS
PRINT ARG=d1,metad.bias STRIDE=100 FILE=COLVAR

HEIGHT is in the energy units of the MD engine (kJ/mol in GROMACS), PACE in MD steps and SIGMA in the units of the CV. After the run, plumed sum_hills --hills HILLS reconstructs the bias on a grid, and from it the free energy.

When to use metadynamics

Metadynamics is at its best when one or two good collective variables are known and the free-energy landscape along them is to be explored without prior knowledge of where the minima and barriers are. It discovers unexpected intermediate states, it is easy to set up, and its well-tempered form converges to a well-defined answer. Its weaknesses are those of every CV-based method: a poor collective variable leaves hidden barriers, and the cost grows steeply with the number of biased variables, since hills must fill a space of growing dimension. The next chapter describes OPES, which reformulates the same goal in a way that converges faster, together with variationally enhanced sampling and the adaptive biasing force.

References and further reading

Target distributions instead of hills

Well-tempered metadynamics, viewed from the end, drives the collective variable towards a definite target distribution, $P(s)^{1/\gamma}$, by depositing hills until it gets there. Two later methods start from that target and ask directly what bias produces it. If the unbiased distribution of the collective variable is $P(s)$ and we want the system to sample a target $p_{\text{tg}}(s)$, the required bias is

because adding $V$ multiplies the distribution by $e^{-\beta V}$. The difficulty is that $P(s)$ is unknown; it is what we want to compute. OPES estimates it on the fly from the simulation itself; VES finds the bias by minimising a functional that does not require knowing $P(s)$ in advance.

OPES: on-the-fly probability enhanced sampling

Michele Invernizzi and Michele Parrinello introduced OPES in 2020. As the simulation runs, it builds a kernel density estimate $\tilde P_n(s)$ of the unbiased distribution, by depositing a Gaussian kernel at each sampled point weighted by $e^{\beta V}$ to undo the current bias, and it sets the bias at every step from that estimate:

The essential difference from metadynamics is what is accumulated. Metadynamics accumulates the bias itself, so every early mistake stays in the bias until later hills compensate. OPES accumulates an estimate of a probability distribution, which converges like any reweighted histogram, and computes the bias from it afresh at every step. As a result the bias becomes quasi-static quickly: after the first few transitions it hardly changes, which makes reweighting straightforward and convergence much faster. In practice OPES needs one main parameter, the expected barrier height BARRIER, from which sensible defaults for the bias factor and kernel width follow.

In the figure, OPES typically makes its first crossing a little later than an aggressive metadynamics run, because the bias is capped and grows from a probability estimate rather than from a pile of hills, but once both states have been visited its bias settles and its estimate converges smoothly, without the persistent fluctuations of plain metadynamics.

Several variants extend the method. OPES-explore (2022) targets the well-tempered distribution of the biased rather than the unbiased simulation, which explores faster at the cost of slower convergence, a useful first step for mapping a landscape. OPES expanded ensembles (Invernizzi, Pablo Piaggi and Parrinello, 2020) apply the same idea to temperatures or Hamiltonians instead of collective variables, sampling several thermodynamic states at once in a single simulation, a CV-free alternative to replica exchange. OPES flooding (Dhiman Ray and Parrinello, 2023) builds the bias only in the reactant basin and keeps the transition state unbiased, giving rates in the same spirit as infrequent metadynamics. In PLUMED, a basic OPES run is

d1: DISTANCE ATOMS=1,10
opes: OPES_METAD ARG=d1 PACE=500 BARRIER=40
PRINT ARG=d1,opes.bias STRIDE=100 FILE=COLVAR

with the barrier in the engine’s energy units (kJ/mol in GROMACS).

Variationally enhanced sampling

Omar Valsson and Michele Parrinello introduced variationally enhanced sampling (VES) in 2014. They showed that the bias producing a chosen target distribution minimises a convex functional:

The bias is expanded in a basis set, for example Legendre polynomials, Chebyshev polynomials or wavelets, or represented by a neural network, and its coefficients are optimised by stochastic gradient descent during the simulation, using averages collected over short intervals. VES turned the construction of a bias into an optimisation problem, which opened the way to neural-network biases and to the machine-learning methods of the last part of this subject. Its practical difficulty is the noise of the stochastic optimisation, which is managed with averaged stochastic gradient schemes.

The adaptive biasing force

A different route, based on forces rather than energies, is the adaptive biasing force (ABF) method of Eric Darve and Andrew Pohorille (2001). It estimates the mean force along the collective variable, $-\mathrm dF/\mathrm ds$, by averaging the instantaneous force acting on the CV in bins of $s$, and applies a biasing force equal and opposite to the current estimate:

Because the mean force is a local quantity, ABF converges well once each bin has enough samples; a ramp that applies only a fraction of the estimate until a bin has a minimum number of samples avoids acting on noise. Select ABF in the figure: the CV diffuses across the whole range almost at once, and the free energy is built from the integrated mean force. The original form requires CVs whose instantaneous force can be computed, which restricts it; the extended-system ABF (eABF) of Adrien Lesage, Tony Lelièvre, Gabriel Stoltz and Jérôme Hénin (2017) couples the CV to a fictitious particle by a stiff spring and applies ABF to that particle, which works with any differentiable CV. ABF is implemented in the Colvars module of NAMD, LAMMPS and GROMACS, and in PLUMED.

Choosing an adaptive method

All of these methods, metadynamics, well-tempered metadynamics, OPES, VES and ABF, solve the same problem along the same kind of collective variables, and in their converged limits they give the same free-energy surface. They differ in how fast they explore, how fast and how smoothly they converge, and what they need from the user. Plain metadynamics explores aggressively but does not converge cleanly; well-tempered metadynamics converges and is the best documented; OPES converges fastest to a quasi-static bias with few parameters and is increasingly the default choice in PLUMED; VES is the most flexible in the form of the bias and the target; ABF is efficient for one or two CVs whose mean force can be estimated and makes no assumption about the shape of the free energy. Every one of them shares the fundamental weakness of collective-variable methods: if the CV misses a slow degree of freedom, the result is wrong in ways that may not be obvious. Methods that avoid collective variables altogether, replica exchange and its relatives, are the subject of the next chapter.

References and further reading

Tempering

Heating without losing the low-temperature ensemble

The collective-variable methods of the previous chapters accelerate sampling along chosen coordinates. A blunter but more general way to cross barriers is to raise the temperature: at high $T$, barriers that are many $k_BT$ high at room temperature become small, and the system explores freely. The difficulty is that the high-temperature ensemble is not the one we want. The idea of tempering is to run at high temperature to explore, and at low temperature to sample correctly, letting configurations move between the two in a way that preserves the correct distribution at every temperature.

Parallel tempering

The method

Parallel tempering, or replica exchange, runs $M$ copies (replicas) of the system simultaneously at temperatures $T_1 < T_2 < \cdots < T_M$, each sampling its own canonical ensemble. Periodically, a swap of configurations between two replicas is attempted and accepted with a probability that keeps every replica’s ensemble exact. The method was introduced for spin glasses by Robert Swendsen and Jian-Sheng Wang (1986), generalised by Charles Geyer (1991) and by Koji Hukushima and Koji Nemoto (1996) as exchange Monte Carlo, and brought to molecular dynamics by Yuji Sugita and Yuko Okamoto in 1999 as replica-exchange molecular dynamics (REMD).

The acceptance rule

The acceptance rule follows from detailed balance in the extended ensemble of all replicas, whose probability is a product of Boltzmann factors. Swapping configurations $\mathbf x_i$ and $\mathbf x_j$ between temperatures $T_i$ and $T_j$ changes that product by a factor that depends only on the two energies:

A random walk in temperature

Each configuration therefore performs a random walk in temperature: it heats up, crosses barriers easily, cools down, and is sampled at the target temperature in a new basin. The figure runs a ladder of walkers by Metropolis Monte Carlo on a rugged one-dimensional landscape with several wells separated by barriers of several $k_BT$ at the lowest temperature.

Without exchanges, the lowest-temperature walker stays in its starting well, and its histogram covers one well only, a large error in the sampled distribution. With exchanges, configurations visited at high temperature flow down the ladder, and the lowest-temperature histogram converges to the full Boltzmann distribution across all wells. Reduce the number of replicas while keeping $T_{\max}$ fixed and the acceptance falls: neighbouring temperatures are now so different that their energy distributions do not overlap, and swaps are almost always rejected.

Replica exchange in a simulation box

A one-dimensional landscape shows the principle; a real system shows what it looks like in a simulation. The figure below runs replica-exchange molecular dynamics on the two-dimensional Lennard-Jones heptamer, LJ7: seven atoms interacting through the Lennard-Jones potential, in reduced units ($\varepsilon = \sigma = m = 1$, time in $\tau = \sigma\sqrt{m/\varepsilon}$), with Langevin dynamics. Christoph Dellago, Peter Bolhuis and David Chandler used this cluster in 1998 to develop transition path sampling, and it has been a standard test of rare-event methods since. Its ground state is a hexagon around a central atom with six neighbours, at $E = -12.53\,\varepsilon$. Three other isomers lie about $1\,\varepsilon$ higher, separated from the hexagon by barriers that a cluster at $k_BT = 0.06\,\varepsilon$ crosses only rarely.

Every replica starts in the same higher isomer, a trapezoid at $E = -11.50\,\varepsilon$. A separate box on the left runs ordinary molecular dynamics at the lowest temperature as a control. The ladder boxes keep their temperatures; the atoms carry the colour of their replica, so when a swap is accepted you see the two configurations change places, and over time each colour wanders up and down the ladder. A green flash between two boxes is an accepted swap, a red one a rejected attempt.

Watch three things. First, the hot boxes change shape constantly, the cold ones rarely; the coldest replica finds the hexagon only after a configuration that reached it at high temperature has walked down the ladder. In the runs used to prepare this chapter, five replicas from $k_BT = 0.06$ to $0.30\,\varepsilon$ with a swap attempt every 100 steps brought the hexagon into the coldest box after $47\,\tau$ (the figure, with a different random seed, typically does so within a few tens of $\tau$), while five independent plain MD runs at $k_BT = 0.06\,\varepsilon$ stayed in the trapezoid for between $1400\,\tau$ and more than $2000\,\tau$. Second, the energy distributions: neighbouring curves overlap, and the overlap is largest at the cold end, where the measured acceptances were 0.50 and 0.46, falling to 0.26 and 0.25 between the hotter pairs. Third, turn the exchanges off: every replica is then an independent simulation, and the coldest one is as trapped as the control. Reduce the number of replicas to three with the same highest temperature, and the gaps between the distributions open up; swaps are rejected, the colours stop moving and the method fails, which is the subject of the next section.

Designing the ladder

The acceptance between neighbouring replicas depends on the overlap of their potential-energy distributions. For a system with $N$ degrees of freedom, the energy fluctuates by $\sim\sqrt N$ while its mean changes with temperature by $\sim N$, so the temperature spacing must shrink as $1/\sqrt N$ and the number of replicas needed to span a given range grows as $\sqrt N$, as David Kofke analysed in 2002. A solvated protein of tens of thousands of atoms needs dozens to hundreds of replicas to span 300–450 K, which is the main cost of the method. Temperatures are usually spaced geometrically, giving roughly uniform acceptance, with target acceptance rates around 20–30 %; the highest temperature must be high enough to cross the relevant barriers quickly, but not so high that the system spends its time in irrelevant high-energy states, such as an unfolded protein that takes long to refold.

Variants of tempering

Hamiltonian replica exchange and REST

For large solvated systems, most of the energy fluctuation comes from the solvent, which is why so many replicas are needed. Hamiltonian replica exchange (Hiroaki Fukunishi, Osamu Watanabe and Shoji Takada, 2002) instead varies the Hamiltonian across replicas, for example by scaling particular interactions, while keeping the temperature fixed; the acceptance rule involves the energies of each configuration under both Hamiltonians. Replica exchange with solute tempering (REST, Pu Liu, Byungchan Kim, Richard Friesner and Bruce Berne, 2005) and its improved form REST2 (Lingle Wang, Friesner and Berne, 2011) scale only the interactions involving a solute, which is equivalent to heating the solute while the solvent stays cold. Because only a small part of the system is effectively heated, far fewer replicas are needed. REST2 is widely used for protein–ligand binding and peptide conformations, and as a component of alchemical free-energy protocols.

Simulated tempering and expanded ensembles

Simulated tempering (Enzo Marinari and Giorgio Parisi, 1992) achieves the same random walk in temperature with a single simulation, by treating the temperature as a dynamical variable that changes by Monte Carlo moves. To make the walk uniform it needs weights for each temperature, essentially the free energy at each temperature, which must be learned; adaptive schemes such as Wang–Landau-style updating or the OPES expanded-ensemble method of the previous chapter provide them. The single-simulation approach avoids running many replicas but is harder to tune.

Combinations and boost potentials

Parallel-tempering metadynamics

Tempering and collective-variable biasing are complementary: tempering helps all degrees of freedom a little, while CV biasing helps a few degrees of freedom a lot. Parallel-tempering metadynamics (Giovanni Bussi, Francesco Gervasio, Alessandro Laio and Michele Parrinello, 2006) combines them, and many protocols use replica exchange between umbrella windows to reduce hidden barriers in orthogonal coordinates.

Accelerated molecular dynamics

A different CV-free strategy adds a boost potential that raises the energy of low-lying regions without changing the high ones, flattening the landscape everywhere at once. Accelerated MD (Donald Hamelberg, John Mongan and Andrew McCammon, 2004) adds, whenever the potential energy $U$ is below a threshold $E$, the boost $\Delta U = (E - U)^2/(\alpha + E - U)$. Gaussian-accelerated MD (GaMD, Yinglong Miao, Victoria Feher and McCammon, 2015) uses a harmonic boost $\Delta U = \tfrac12 k(E - U)^2$, chosen so that the distribution of boost energies is nearly Gaussian, which makes reweighting by a cumulant expansion accurate. Both are implemented in AMBER and other codes and are popular for exploring protein conformational changes, but reweighting becomes noisy for large boosts and large systems, which limits their use for accurate free energies.

References and further reading

Sampling paths instead of states

The methods of the previous part bias the dynamics along a chosen collective variable and then correct for the bias. They work well when a good collective variable is known. When it is not, when the mechanism itself is the question, a different philosophy is attractive: leave the dynamics unbiased and instead sample trajectories, concentrating effort on the rare trajectories that actually make a transition. Path-sampling methods do this. Because the trajectories they produce are true dynamical trajectories, they give rates and mechanisms directly, without assuming a reaction coordinate.

The central object is the transition path ensemble: the set of all trajectories of a fixed length $L$ that start in state A and end in state B, each weighted by its probability under the true dynamics. For a trajectory $\mathbf X = (\mathbf x_0, \mathbf x_1, \ldots, \mathbf x_L)$,

Transition path sampling

Transition path sampling (TPS), developed by Christoph Dellago, Peter Bolhuis, Félix Csajka and David Chandler in 1998, samples this ensemble by a Monte Carlo random walk in the space of trajectories. The most important move is the shooting move: choose a time slice of the current path at random, perturb its momenta slightly (or draw new ones from the Maxwell–Boltzmann distribution), and integrate the equations of motion forwards to the end and backwards (with reversed momenta) to the beginning, creating a new trajectory that shares one point with the old one. With symmetric generation of the perturbation and equilibrium momenta, the new path is accepted simply if it still connects A and B; otherwise the old path is kept. Shifting moves, which translate a path in time, improve sampling of the endpoints.

Because each new path shares a point with the previous one, TPS explores the transition path ensemble gradually, like a random walk; it never needs to find a transition from scratch. Starting from any connecting path, even an unphysical one such as a straight line between the states, the first accepted shooting moves replace it with a genuine dynamical trajectory, and the sampled paths soon forget the initial guess. In the figure, the initial path is a straight line; watch it relax into the dynamical channel through which transitions really occur, and watch accepted paths explore the width of that channel.

The accepted paths constitute a sample of how the transition really happens: which degrees of freedom move first, through which bottleneck the system passes, how long the transition takes. Analysing the ensemble, for example by computing committors along the paths or by fitting a model that predicts the outcome of shooting moves from collective variables (the likelihood maximisation of Baron Peters and Bernhardt Trout, 2006), identifies reaction coordinates from data. TPS has been used to study the mechanisms of ion dissociation in water, autoionisation of water, nucleation, protein conformational changes and chemical reactions. The OpenPathSampling library implements TPS and the methods below.

Rates from path sampling

TPS samples transitions but not, by itself, how often they happen. Transition interface sampling (TIS, Titus van Erp, Daniele Moroni and Peter Bolhuis, 2003) computes the rate by placing a series of interfaces $\lambda_0 < \lambda_1 < \cdots < \lambda_n$ between A and B, defined by an order parameter, and writing the rate as a product of the flux out of A through the first interface and a series of conditional probabilities:

Each conditional probability is moderate, say 0.1 to 0.3, even when their product is $10^{-20}$, so each can be estimated with a reasonable number of trajectories. Unlike a free-energy method, the order parameter only defines the interfaces; it does not have to be a good reaction coordinate, although a better one makes the calculation more efficient. Replica-exchange TIS swaps paths between interface ensembles to improve sampling.

Forward-flux sampling (FFS, Rosalind Allen, Patrick Warren and Pieter Rein ten Wolde, 2005) uses the same factorisation but generates trajectories only forwards in time: launch many trajectories from the configurations that first crossed interface $i$, record those that reach interface $i+1$, and use their endpoints as starting points for the next stage. Because it never integrates backwards, FFS works for non-equilibrium systems without detailed balance, such as driven systems and biochemical networks, and it is widely used for crystal nucleation, where the order parameter is the size of the largest crystalline cluster.

Milestoning and weighted ensemble

Milestoning (Anton Faradjian and Ron Elber, 2004) also places interfaces, called milestones, along the transition, but it runs only short trajectories between neighbouring milestones and records the times and probabilities of reaching each neighbour. Combining these local statistics, under the assumption that the system loses memory between milestones, gives the mean first-passage time over the whole path. Its strength is that the short trajectories are cheap and independent, which makes it efficient for very slow processes such as ligand unbinding.

The weighted-ensemble method (Gary Huber and Sangtae Kim, 1996) runs many trajectories in parallel, each carrying a statistical weight. Configuration space is divided into bins along a progress coordinate, and at regular intervals trajectories in under-populated bins are split into several copies, sharing the parent’s weight, while trajectories in crowded bins are merged. Splitting and merging keep the total weight exact, so the method is unbiased, yet it concentrates computing effort on the rare trajectories that advance towards B. Weighted ensemble gives rates and pathways and has been applied to large biomolecular processes such as protein–protein binding and viral capsid dynamics, notably with the WESTPA software.

Strings at finite temperature and transition path theory

The string method of the chapter on energy landscapes finds minimum energy paths at zero temperature. Its finite-temperature version (Weinan E, Weiqing Ren and Eric Vanden-Eijnden, 2005) evolves the string according to averages computed in restrained simulations around each image, so that it converges to the principal curve of the transition tube in collective-variable space, accounting for entropy. The string method with swarms of trajectories (Albert Pan, Deniz Sezer and Benoît Roux, 2008) moves each image by the average drift of many short unrestrained trajectories launched from it, which avoids the need for restraints and is popular for large conformational changes in proteins.

All these methods can be understood within transition path theory (E and Vanden-Eijnden, 2006), which describes the statistics of reactive trajectories in terms of the committor $q(\mathbf x)$, the probability of reaching B before A introduced in the chapter on rates. The density of reactive trajectories, the probability current of transitions and the rate can all be written in terms of $q$; the committor’s isosurfaces are the natural interfaces, and the committor itself is the ideal reaction coordinate. Much of the recent work on machine-learned collective variables aims to learn $q$ from data (the chapter on learned collective variables).

References and further reading

Long kinetics from short trajectories

Enhanced-sampling methods change the dynamics to obtain thermodynamics; path sampling targets particular transitions. Markov state models (MSMs) take a third approach that has become central in biomolecular simulation: run many ordinary, unbiased simulations, each much shorter than the processes of interest, and stitch their statistics together into a model of the long-time dynamics. Because the simulations can be run independently and in parallel, on many GPUs or on distributed computing projects such as Folding@home, MSMs reach timescales of milliseconds and beyond from simulations of microseconds.

The idea is to divide configuration space into many small regions, microstates, and to model the dynamics as jumps between them, observed at a fixed time interval $\tau$, the lag time. If $\tau$ is long enough that the system forgets where within a microstate it came from, the jumps are memoryless, and the whole dynamics is described by a matrix of transition probabilities.

The transition matrix

Let $c_{ij}(\tau)$ count how often a trajectory in microstate $i$ at time $t$ was in microstate $j$ at time $t + \tau$, summed over all trajectories and all times. The maximum-likelihood transition matrix is the row-normalised count matrix, often symmetrised (or estimated with a reversibility constraint) to respect detailed balance in equilibrium:

The power of the model is that it can be propagated far beyond the length of any single trajectory: $\mathbf T(\tau)^n$ predicts the evolution after $n$ lag times. Everything about the long-time dynamics is contained in the eigenvalues and eigenvectors of $\mathbf T$.

Eigenvalues and implied timescales

The largest eigenvalue of a transition matrix is exactly 1, and its left eigenvector is the stationary distribution. The other eigenvalues $\lambda_2 \ge \lambda_3 \ge \ldots$ are smaller than 1, and each describes a relaxation process: a pattern of population flowing between groups of microstates that decays by a factor $\lambda_k$ every lag time. Writing $\lambda_k = e^{-\tau/t_k}$ defines the implied timescale of that process:

The eigenvectors give a natural way to find metastable states: groups of microstates that exchange population quickly among themselves and slowly with other groups. Clustering microstates by their values in the slowest eigenvectors, as the PCCA+ algorithm of Peter Deuflhard and Marcus Weber does, lumps hundreds of microstates into a few macrostates with a clear physical meaning.

Validating the model

The Markov assumption is never exact; it is a good approximation only when the lag time is long enough. The standard test, introduced by William Swope, Jed Pitera and Fred Suits in 2004, plots the implied timescales as a function of $\tau$. If the dynamics between microstates is Markovian, the timescales are independent of $\tau$; if not, they increase with $\tau$, because memory makes the short-lag model look faster than the real dynamics. One chooses the smallest lag time at which the slowest timescales have levelled off. A stricter test is the Chapman–Kolmogorov test: compare the populations predicted by $\mathbf T(\tau)^n$ with those estimated directly at lag $n\tau$.

The figure builds an MSM live. A particle diffuses on a landscape with three wells; the plane is divided into a grid of microstates; as the trajectory grows, the count matrix and its eigenvalues are recomputed. The left panel colours each microstate by the slowest eigenvector, whose sign separates the two groups of states that exchange most slowly. The right panel shows the implied timescales against lag time.

At short lag times the implied timescales are too small: the coarse grid of microstates hides fast motions within each cell, the model has memory, and the transitions look faster than they are. As the lag increases, the timescales rise and level off, the sign that the model has become Markovian. Once the lag approaches the timescale itself, statistics run out. A good MSM is built at the smallest lag on the plateau, and the plateau values are its predictions for the slow relaxation times.

Better microstates: TICA and the variational principle

The quality of an MSM depends on how well the microstates resolve the slow processes. Microstates are usually defined by clustering conformations, but clustering in raw coordinates such as atomic positions emphasises large, fast fluctuations and may split slow processes across many clusters. Since 2013, the standard practice has been to first project the data onto slow collective coordinates found by time-lagged independent component analysis (TICA, introduced to molecular simulation by Guillermo Pérez-Hernández and colleagues and by Christian Schwantes and Vijay Pande), which finds linear combinations of input features with maximal autocorrelation at a lag time, the coordinates that change most slowly. The chapter on learned collective variables shows TICA finding a slow direction that principal component analysis misses.

TICA is a special case of a general variational principle. Frank Noé and Feliks Nüske showed in 2013 that, for reversible dynamics, the eigenvalues of the transfer operator are bounded from above by those estimated with any finite set of basis functions, just as the Rayleigh–Ritz principle bounds energies in quantum mechanics. Better features give higher eigenvalues and longer timescales, so models can be compared and optimised by a score: the variational approach to conformational dynamics (VAC), and its generalisation to non-reversible processes, the variational approach for Markov processes (VAMP, Hao Wu and Noé). These scores turned MSM construction into an optimisation problem and led directly to VAMPnets (Andreas Mardt, Luca Pasquali, Hao Wu and Noé, 2018), neural networks that learn the optimal soft assignment to metastable states end to end.

Fluxes, committors and pathways

An MSM also gives the mechanisms of transitions through transition path theory (Weinan E and Eric Vanden-Eijnden, applied to MSMs by Philipp Metzner, Christof Schütte and Vanden-Eijnden, and by Noé and colleagues for protein folding in 2009). The forward committor $q_i$ of each microstate, the probability of reaching B before A, solves a linear system built from $\mathbf T$; from it follow the reactive flux between microstates, the dominant transition pathways, and the rate $k_{AB}$. MSM studies of protein folding, ligand binding and conformational change in kinases and receptors have used these tools to show how multiple parallel pathways contribute to a single observed process.

In practice

Building an MSM involves choosing features (distances, dihedrals, contacts), reducing them with TICA or a VAMP-trained network, clustering into a few hundred microstates with k-means, choosing the lag time from implied timescales, estimating a reversible transition matrix, validating with the Chapman–Kolmogorov test, and coarse-graining into macrostates with PCCA+. Errors are estimated by Bayesian sampling of transition matrices or by bootstrapping over trajectories. The PyEMMA and deeptime libraries, from the Noé group, and MSMBuilder, from the Pande group, implement this workflow. The main requirement is enough sampling: every transition the model describes must have been observed, in both directions, enough times in the short trajectories, which is often achieved by adaptive sampling, starting new trajectories from poorly sampled microstates.

References and further reading

Learning the coordinates

Every collective-variable method depends on choosing the right variables, and the chapter on collective variables showed how a poor choice hides barriers. For simple processes, physical intuition suffices: a distance, a dihedral, a coordination number. For complex ones, such as the folding of a protein, the nucleation of a crystal or a reaction in a disordered solvent, the relevant combination of hundreds of possible coordinates is not obvious. Since about 2010, a family of methods has learned collective variables from simulation data instead. They differ in what they ask the variable to capture: the largest fluctuations, the difference between known states, the slowest dynamics, or the ability to predict the future.

All of them share a practical pattern. Choose a large set of physically meaningful input features, such as distances, contacts, dihedral angles or symmetry functions; collect data, from unbiased simulations, from short simulations in each known state, or from a first round of enhanced sampling; train a model that maps features to a few collective variables; use those variables to bias a new simulation, which samples transitions and produces better data; and repeat until the results converge.

Variance is not slowness

The oldest and simplest approach is principal component analysis (PCA): find the linear combinations of features with the largest variance. Its weakness is fundamental. Large fluctuations are often fast and irrelevant, like the floppy motion of a loop or the vibration of a side chain, while the rare transition that matters may involve small, slow changes. Nonlinear unsupervised methods, such as diffusion maps (applied to molecules by Mary Rohrdanz, Cecilia Clementi and colleagues, 2011), sketch-map (Michele Ceriotti, Gareth Tribello and Michele Parrinello, 2011) and autoencoders (used for enhanced sampling by Wei Chen and Andrew Ferguson, 2018), capture more complex shapes of the data but share the same blind spot when trained on static configurations: they see where the system is, not how slowly it moves.

Time-lagged independent component analysis (TICA) uses time information instead. It seeks the linear combination of features $\mathbf v^{\mathsf T}\boldsymbol\chi(\mathbf x)$ whose autocorrelation after a lag time $\tau$ is largest, that is, the coordinate that changes most slowly:

The figure runs a two-dimensional system with two basins separated along a narrow direction, with large but fast fluctuations across it. PCA picks the direction of the large fluctuations; TICA picks the direction along which the system jumps between the basins. The time series on the right make the difference plain: the projection on the PCA direction is noise, while the projection on the TICA direction shows clean switching between two states.

Deep-TICA (Luigi Bonati, GiovanniMaria Piccini and Parrinello, 2021) replaces the linear features with a neural network trained to maximise the TICA eigenvalues of its outputs, giving nonlinear slow variables. The same objective, in its general variational form, drives VAMPnets (Andreas Mardt, Luca Pasquali, Hao Wu and Frank Noé, 2018) and state-free reversible VAMPnets (Wei Chen, Hythem Sidky and Andrew Ferguson, 2019), which learn slow coordinates or soft metastable-state assignments directly from trajectories. Deep-TICA was designed to be trained on data from a first metadynamics or OPES run, with the bias reweighted so that the learned variables describe the unbiased dynamics, and then used as collective variables in a second, better-targeted run.

Discriminating known states

When the metastable states are known but the path between them is not, a supervised approach trains a variable that separates them. Linear discriminant analysis finds the linear combination of features that best separates two or more labelled sets of configurations relative to their internal spread; its harmonic variant, HLDA (Dan Mendels, GiovanniMaria Piccini and Parrinello, 2018), was proposed specifically for enhanced sampling. Deep-LDA (Bonati, Valerio Rizzi and Parrinello, 2020) applies the same criterion to the output of a neural network. Training needs only short unbiased simulations in each state, which are cheap, and the learned variable can then drive metadynamics or OPES between them. Its limitation is that separating the states is not the same as describing the transition between them: the variable can be excellent at telling A from B and still poor at the transition state.

Predicting the future: RAVE and SPIB

Pratyush Tiwary and colleagues approached the problem through information theory. A good reaction coordinate should be a compressed description of the present that retains the information needed to predict the future, an information bottleneck. The reweighted autoencoded variational Bayes method (RAVE, João Marcelo Lamim Ribeiro, Pablo Bravo, Yihang Wang and Tiwary, 2018) trained a variational autoencoder to predict configurations a time $\Delta t$ ahead from a low-dimensional latent variable. Its successor, the state predictive information bottleneck (SPIB, Dedi Wang and Tiwary, 2021), predicts not the future configuration but the future state label, which is both easier and more relevant:

The state labels are not fixed in advance. SPIB starts from an initial coarse labelling, for example by clustering, and iteratively relabels each configuration by the state its latent variable predicts, so that the number and definition of metastable states are learned along with the coordinate. The time delay $\Delta t$ sets the timescale of interest: short delays resolve fast processes, long delays keep only the slowest. The learned latent variable is a natural candidate for biasing with metadynamics, and SPIB has been used for protein conformational changes, ligand unbinding and crystal nucleation, often in iterative loops of biased sampling and retraining.

Learning the committor

The ideal reaction coordinate is the committor $q(\mathbf x)$, the probability of reaching B before A. Several methods learn it directly. Ao Ma and Aaron Dinner (2005) fitted committors estimated from shooting trajectories with neural networks; more recent approaches train networks to satisfy the variational principle that the committor minimises a weighted average of $|\nabla q|^2$ subject to boundary conditions, using data from biased simulations, which gives the committor and a characterisation of the transition state ensemble without running large numbers of shooting trajectories. Committor-based variables are attractive because they measure progress through the transition itself, but learning them requires data from the transition region, which is exactly what is rare; most practical schemes therefore alternate between learning and biased sampling.

Using learned variables in practice

Learned collective variables are only useful if they can be differentiated and evaluated inside a simulation. The mlcolvar library (Bonati, Enrico Trizio, Andrea Rizzi and Parrinello, 2023) provides PyTorch implementations of autoencoders, Deep-LDA, Deep-TICA, committor models and others, and exports trained models in a form that PLUMED loads (through its PYTORCH_MODEL collective variable) to drive metadynamics or OPES in any supported MD engine. SPIB and VAMPnet implementations are available from their authors and in the deeptime library.

Three cautions apply. A learned variable can only describe processes present in its training data: if a transition has never been observed, no amount of learning will find it, which is why iterative schemes are standard. Neural-network variables can behave unpredictably outside the region where they were trained, producing large biasing forces in unexplored regions, so it is wise to check them on new data and to cap biases. And interpretability matters: a variable that works but cannot be understood teaches less than one that can, so it is worth analysing which input features a learned variable depends on most.

References and further reading

Sampling without dynamics

Why trajectories are slow

Every method in the previous chapters generates samples of the Boltzmann distribution $p(\mathbf x) \propto e^{-U(\mathbf x)/k_BT}$ by moving through configuration space one small step at a time. A molecular dynamics step changes positions by a fraction of an ångström; a Monte Carlo move is accepted only if it does not raise the energy too much. Consecutive samples are therefore strongly correlated, and the time needed to pass from one metastable state to another grows exponentially with the barrier between them, as the chapter on rates showed. Enhanced sampling methods attack this by biasing the dynamics, heating the system or pulling it along a collective variable, but they all keep the basic picture of a walker moving through the landscape.

A different question

In 2019 Frank Noé, Simon Olsson, Jonas Köhler and Hao Wu asked a different question. Suppose a neural network could be trained to produce configurations directly: feed it random numbers, and it returns a molecular configuration, independent of the previous one, drawn from a distribution $q(\mathbf x)$ close to the Boltzmann distribution. Configurations from distant metastable states would appear in the same batch, with no transitions to wait for. They called such a network a Boltzmann generator. Two problems must be solved before the idea is useful. A learned $q$ is never exactly $p$, so its errors must be corrected; and the network must be trainable when the very thing we lack is a good sample of $p$. This chapter develops both solutions, then trains a Boltzmann generator in front of you.

Importance sampling

Reweighting to the Boltzmann distribution

The errors of an approximate sampler can be removed exactly if its density can be evaluated. The idea is older than computers: Herman Kahn and others used it in the Monte Carlo neutron calculations of the early 1950s. Draw samples $\mathbf x_i$ from $q$, and give each a weight equal to the ratio of the target density to the proposal density:

The estimate converges to the exact Boltzmann average for any $q$ that is non-zero wherever $p$ is. A sample in a region the generator over-populates receives a small weight; a sample in a region it under-populates receives a large one.

The effective sample size

How well $q$ matches $p$ decides how many samples are needed. If the weights are nearly equal, every sample counts; if one sample carries most of the total weight, the average is effectively computed from that single sample. Leslie Kish’s effective sample size, introduced for survey statistics in 1965, measures this:

The ratio $n_{\text{eff}}/n$ is the standard diagnostic of a Boltzmann generator: close to 1, the generator is nearly exact; close to 0, a few samples dominate and the estimates are noisy.

What reweighting cannot fix

Reweighting has a blind spot. If $q$ is practically zero in a region where $p$ is large, samples from that region essentially never appear, and no weight can be attached to samples that do not exist. The estimate is then wrong, and, worse, the effective sample size can still look excellent, because it is computed only from the samples that were drawn. The figure shows this with a proposal built by hand from Gaussian modes: shrink the modes and a few samples carry enormous weights; switch to mode collapse and the population of the missing state is estimated as zero with no warning. This limitation will return when we train a real generator.

Normalizing flows

Changing variables

To reweight, the generator must report the density $q(\mathbf x)$ of every configuration it produces. Most neural generators cannot do this, but one family can by construction. Start from random numbers $\mathbf z$ with a simple density, a standard Gaussian $p_z$, and transform them with an invertible function $\mathbf x = f(\mathbf z)$. Probability is conserved under the transformation: where $f$ stretches space, the same probability is spread over a larger volume and the density falls; where it compresses space, the density rises. The local stretching factor is the Jacobian determinant, and the result is the change-of-variables formula of calculus:

A transformation built from a sequence of simple invertible steps that turns a simple density into a complicated one is called a normalizing flow: the density “flows” through the steps and stays normalised throughout. Esteban Tabak and Eric Vanden-Eijnden formulated density estimation this way in 2010; Danilo Rezende and Shakir Mohamed brought the name and the idea into deep learning in 2015.

Coupling layers

A general neural network is neither invertible nor equipped with a cheap determinant: for $d$ coordinates the determinant costs of order $d^3$ operations. Laurent Dinh, David Krueger and Yoshua Bengio found a way around both problems in 2014 (NICE), extended by Dinh, Jascha Sohl-Dickstein and Samy Bengio in 2017 (RealNVP). Split the coordinates into two groups, $\mathbf x_a$ and $\mathbf x_b$. Leave the first group unchanged and transform the second by a scale and a shift computed from the first:

The Jacobian of this map is triangular, because $\mathbf y_a$ does not depend on $\mathbf x_b$, so its determinant is the product of the diagonal entries $e^{s_j}$. The networks $\mathbf s$ and $\mathbf t$ can be as expressive as we like; invertibility is guaranteed by the structure, not by the networks.

Stacking layers

One coupling layer cannot do much, because half of the coordinates pass through it unchanged. Alternating which group is transformed, and stacking many layers, produces a flexible map whose log-determinant is simply the sum over layers. The figure below uses eight such layers on two coordinates. Its slider shows the samples after each layer: watch the Gaussian cloud being sheared one coordinate at a time, folded into a curved valley, and finally split into two separate states.

Training a Boltzmann generator

Training by energy

How should the flow be trained when no good samples of $p$ exist? The answer of Noé and colleagues was to minimise a distance between $q$ and $p$ that requires only the energy function. The Kullback–Leibler divergence of $q$ from $p$ is

The average is over Gaussian noise, so its gradient with respect to the network parameters is computed by drawing $\mathbf z$, pushing it through the flow and backpropagating through the energy and the layers, the reparameterisation trick of Diederik Kingma and Max Welling (2014). The force $-\nabla U$ appears directly in this gradient: the generator is trained by the same forces that drive molecular dynamics.

Training by example

The second way to train is by example, from configurations $\mathbf x_n$ produced by short simulations. Here the flow is used in reverse: each configuration is mapped back to noise, and its likelihood is maximised,

which minimises the opposite divergence, $\mathrm{KL}(p\,\|\,q)$, as far as the data represent $p$. The short simulations do not need the correct relative populations of the states, only reasonable coverage of each; the populations are corrected later by reweighting.

Mode-seeking and mode-covering

The two divergences fail in opposite ways, and the difference is the most important practical fact about Boltzmann generators. The reverse divergence $\mathrm{KL}(q\|p) = \int q\ln(q/p)$ is an average over the generator’s own samples, so it heavily penalises samples in high-energy regions but pays nothing for regions the generator never visits. Starting from a Gaussian, the flow finds one well, fits it well, and has no incentive to cross the barrier: it is mode-seeking and liable to mode collapse. The forward divergence $\mathrm{KL}(p\|q) = \int p\ln(p/q)$ is an average over the data, so it heavily penalises any data point the generator assigns low probability: it is mode-covering, but only for the modes the data contain. Noé and colleagues combined the two: likelihood on short simulations in each known state keeps all states in the generator, and the energy term sharpens the generator where data are sparse.

A Boltzmann generator in your browser

The figure trains a real normalizing flow, eight affine coupling layers with 24 hidden units each, on a two-dimensional double well with a curved valley,

whose barrier of about $8\,k_BT$ keeps ordinary Langevin dynamics in the well where it starts. The right-hand state lies $0.97\,k_BT$ higher in free energy and holds 27.4 % of the population. All of it, the coupling networks, the backpropagation through every layer and the Adam optimiser, runs in JavaScript in this page, about a hundred training steps per second.

Work through the four training modes in order.

1. Likelihood, one state. The training data come from one simulation started in the left well, which never crosses the barrier. The flow learns the left well very well, the effective sample size is high, and the estimated population of the right state is zero. Nothing in the output signals that anything is missing.
2. Energy only. No data at all. The flow finds a well and collapses into it; which one depends on the random initial network. Again the effective sample size looks healthy and the free-energy difference is wrong.
3. Likelihood, both states. One short simulation in each well, equally long, so the data have the wrong ratio of 50 : 50. The raw samples inherit that ratio, but reweighting corrects it, and the free-energy difference converges to the exact $0.97\,k_BT$ within statistical error.
4. Likelihood + energy. The combination used by Noé and colleagues: the data hold both states in place and the energy term refines the shape of each well. The raw samples are already closer to the Boltzmann distribution, and the reweighted estimate is exact.

Press Play the transformation after training to see what the network has learned: the first layers shear and stretch the Gaussian, the middle layers bend it along the curved valley, and the final layers tear it into two pieces with a gap at the barrier, where the flow places almost no probability.

The same model in PyTorch

The browser version is written by hand to show every step. In practice a flow is written in a framework with automatic differentiation. The program below trains the same model in PyTorch with both losses and reweights 200 000 samples; it runs in about two minutes on a laptop CPU.

import math
import torch
from torch import nn

torch.manual_seed(0)
kT = 1.0

def energy(x):                                   # x: (n, 2)
    a, b = x[:, 0], x[:, 1]
    return 8 * (a**2 - 1)**2 + 0.5 * a + 2 * (b - 0.6 * torch.sin(1.8 * a))**2

class Coupling(nn.Module):
    """Affine coupling: keep x[:, keep], transform x[:, move] by a scale and shift
    computed from x[:, keep]. The Jacobian is triangular, so log|det J| = sum of s."""
    def __init__(self, keep, hidden=32):
        super().__init__()
        self.keep, self.move = keep, 1 - keep
        self.net = nn.Sequential(nn.Linear(1, hidden), nn.Tanh(),
                                 nn.Linear(hidden, hidden), nn.Tanh(),
                                 nn.Linear(hidden, 2))
        nn.init.zeros_(self.net[-1].weight); nn.init.zeros_(self.net[-1].bias)   # start as the identity

    def st(self, a):
        s, t = self.net(a).chunk(2, dim=1)
        return 2 * torch.tanh(s / 2), t                      # bounded scale for stability

    def forward(self, x):                                    # z -> x direction
        a, b = x[:, [self.keep]], x[:, [self.move]]
        s, t = self.st(a)
        out = torch.empty_like(x)
        out[:, [self.keep]], out[:, [self.move]] = a, b * torch.exp(s) + t
        return out, s.squeeze(1)

    def inverse(self, y):                                    # x -> z direction
        a, yb = y[:, [self.keep]], y[:, [self.move]]
        s, t = self.st(a)
        out = torch.empty_like(y)
        out[:, [self.keep]], out[:, [self.move]] = a, (yb - t) * torch.exp(-s)
        return out, -s.squeeze(1)

class Flow(nn.Module):
    def __init__(self, n_layers=8):
        super().__init__()
        self.layers = nn.ModuleList([Coupling(keep=k % 2) for k in range(n_layers)])

    def sample(self, n):                                     # returns x and log q(x)
        z = torch.randn(n, 2)
        logq = -0.5 * (z**2).sum(1) - math.log(2 * math.pi)
        x = z
        for layer in self.layers:
            x, ld = layer(x)
            logq = logq - ld
        return x, logq

    def log_prob(self, x):
        logdet = torch.zeros(len(x))
        for layer in reversed(self.layers):
            x, ld = layer.inverse(x)
            logdet = logdet + ld
        return -0.5 * (x**2).sum(1) - math.log(2 * math.pi) + logdet

def langevin(x0, steps=20000, dt=0.004, every=20):
    """Overdamped Langevin dynamics: the 'short simulations' that provide training data."""
    x = torch.tensor([[x0, 0.6 * math.sin(1.8 * x0)]])
    frames = []
    for i in range(steps):
        x.requires_grad_(True)
        f = -torch.autograd.grad(energy(x).sum(), x)[0]
        x = (x + f * dt + math.sqrt(2 * kT * dt) * torch.randn_like(x)).detach()
        if i % every == 0:
            frames.append(x)
    return torch.cat(frames)

data = torch.cat([langevin(-1.0), langevin(1.0)])           # one short run in each state
print(f"training frames: {len(data)}, fraction in right state: {(data[:, 0] > 0).float().mean():.2f}")

flow = Flow()
opt = torch.optim.Adam(flow.parameters(), lr=2e-3)
for step in range(3001):
    batch = data[torch.randint(len(data), (256,))]
    loss_ml = -flow.log_prob(batch).mean()                  # training by example
    x, logq = flow.sample(256)
    loss_kl = (energy(x) / kT + logq).mean()                # training by energy
    loss = loss_ml + loss_kl
    opt.zero_grad(); loss.backward(); opt.step()
    if step % 1000 == 0:
        print(f"step {step:5d}   NLL {loss_ml.item():7.3f}   KL {loss_kl.item():7.3f}")

with torch.no_grad():                                        # reweight to the Boltzmann distribution
    x, logq = flow.sample(200_000)
    logw = -energy(x) / kT - logq
    w = torch.exp(logw - logw.max())
    ess = w.sum()**2 / (w**2).sum() / len(w)
    right = x[:, 0] > 0
    p_raw = right.float().mean()
    p_rw = w[right].sum() / w.sum()
    dF = -kT * torch.log(p_rw / (1 - p_rw))

grid = torch.linspace(-3, 3, 200_001)
bz = torch.exp(-(8 * (grid**2 - 1)**2 + 0.5 * grid))
p_exact = bz[grid > 0].sum() / bz.sum()
print(f"effective sample size: {100 * ess:.1f} %")
print(f"P(right): raw {p_raw:.3f}, reweighted {p_rw:.3f}, exact {p_exact:.3f}")
print(f"dF(right - left): {dF:.3f} kT, exact {-math.log(p_exact / (1 - p_exact)):.3f} kT")

Its output:

training frames: 2000, fraction in right state: 0.50
step     0   NLL   2.616   KL  13.500
step  1000   NLL   1.038   KL   0.179
step  2000   NLL   0.996   KL   0.135
step  3000   NLL   1.030   KL   0.140
effective sample size: 89.8 %
P(right): raw 0.357, reweighted 0.273, exact 0.274
dF(right - left): 0.980 kT, exact 0.973 kT

The training data hold the two states in equal amounts, the raw generator puts 35.7 % of its samples in the right state, and reweighting brings the population to 27.3 %, against the exact 27.4 %. The KL column is the reverse divergence minus the constant $\ln Z$. For this potential $\ln Z = -0.098$, so the final value of $0.140$ corresponds to a divergence of about $0.04$: the generator is close to, but not exactly, the Boltzmann distribution, and reweighting removes the remaining difference.

Free energies and exact sampling

Free-energy differences by reweighting

With weights in hand, the free-energy difference between two states $A$ and $B$ follows from their reweighted populations,

This is the calculation the figure performs. It needs no path between the states and no collective variable: the generator connects them through its latent space, where the two states are neighbouring regions of the same Gaussian. The same idea, mapping one ensemble onto another with a learned invertible transformation, underlies targeted free-energy perturbation, proposed by Christopher Jarzynski in 2002 and realised with normalizing flows by Peter Wirnsberger and colleagues in 2020.

A generator as a Monte Carlo proposal

Reweighting is exact only in the limit of many samples. A generator can also be made exact sample by sample, by using it as the proposal of a Monte Carlo chain. In the independence Metropolis sampler, analysed by Luke Tierney in 1994, each new configuration $\mathbf x'$ is drawn from $q$ independently of the current one $\mathbf x$, and accepted with probability

the ratio of importance weights. The chain samples the Boltzmann distribution exactly, and its acceptance rate is a direct measure of the quality of the generator. Generators can also propose local moves combined with molecular dynamics, the strategy of Timewarp below.

Beyond the first Boltzmann generators

Symmetry and coordinates

A molecule’s energy does not change if it is rotated or translated or if identical atoms are exchanged; a flow on raw Cartesian coordinates must learn these symmetries from data, wastefully. Two remedies are used. The first transforms to internal coordinates, bond lengths, angles and torsions, in which rigid-body motion has been removed; the original Boltzmann generators did this for proteins. The second builds the symmetry into the network: the equivariant flows of Jonas Köhler, Leon Klein and Frank Noé (2020) give densities that are invariant under rotations and permutations by construction.

Continuous flows and flow matching

Instead of a stack of discrete layers, the transformation can be a continuous flow defined by a learned velocity field, $\mathrm d\mathbf x/\mathrm dt = \mathbf v_\theta(\mathbf x, t)$, integrated from $t = 0$ to $1$. These continuous normalizing flows (Ricky Chen and colleagues, 2018) need no special layer structure, and the density follows from the divergence of the velocity field,

They were expensive to train until flow matching (Yaron Lipman and colleagues, 2023) replaced likelihood training by a simple regression of $\mathbf v_\theta$ onto the velocities of straight-line paths between noise and data. Equivariant flow matching (Klein, Andreas Krämer and Noé, 2023) applied this to molecules and is the basis of most current Boltzmann generators.

Transferable generators

A Boltzmann generator trained for one molecule is useless for another. Timewarp (Leon Klein and colleagues, 2023) learns instead the distribution of configurations a long time step ahead in molecular dynamics, conditioned on the current configuration, across many small peptides, and uses Metropolis acceptance to keep sampling exact on peptides it was not trained on. Transferability across chemistry, rather than speed on one system, is now the main goal of the field.

Diffusion models for molecular ensembles

The score is the force

Diffusion models generate samples by learning to reverse a gradual noising process; what they learn is the score, the gradient of the log-density, at every noise level. For the Boltzmann distribution the score has a familiar meaning:

A diffusion model trained on equilibrium configurations therefore learns, at its lowest noise level, an approximation to the force field. Marloes Arts and colleagues (2023) used this to obtain coarse-grained force fields from diffusion models trained on simulation data.

Conformers and proteins

Diffusion models are now the dominant generative architecture for molecular structures. GeoDiff (Minkai Xu and colleagues, 2022) and torsional diffusion (Bowen Jing and colleagues, 2022) generate conformer ensembles of small molecules; equivariant diffusion models generate new molecules and crystals; and models trained on protein structures and large amounts of simulation data, such as BioEmu (Sarah Lewis and colleagues, 2025), generate ensembles of protein conformations orders of magnitude faster than molecular dynamics.

The density problem

Diffusion models do not provide the density of a generated sample cheaply. It can be computed by integrating the associated probability-flow ordinary differential equation, as for a continuous flow, but at considerable cost. Without it, a generated ensemble cannot be reweighted and must be validated against simulation or experiment. Whether a generated ensemble is the Boltzmann distribution, or only resembles it, is the central question of generative sampling.

What generative sampling can and cannot do

Generative samplers deliver what enhanced sampling has always sought: independent configurations from widely separated states, without waiting for transitions. They are strongest when tied to physics, trained on simulation data, refined with the energy and corrected by reweighting or Metropolis acceptance. Their weaknesses are equally clear. A generator cannot be trusted in regions it has never been shown, and mode collapse can hide a state while every diagnostic looks healthy, as the first two modes of the figure demonstrate. Scaling is the open problem: the overlap between $q$ and $p$, and with it the effective sample size, falls roughly exponentially with the number of degrees of freedom, so the success on two dimensions here, or on small proteins in vacuum, does not carry over directly to a solvated system of ten thousand atoms. The practical frontier combines generative proposals with exact acceptance steps and with molecular dynamics, using the generator for speed and statistical mechanics for correctness.

References and further reading