Lesson 7.3 · 16 min

Enhanced sampling: metadynamics

Many important events (a reaction, a conformational change, an ion leaving its solvation shell) are separated by free-energy barriers that ordinary MD almost never crosses. Metadynamics adds a history-dependent bias along a few collective variables that fills the free-energy wells, so the system escapes; the accumulated bias then gives the free-energy surface.

Why rare events need enhanced sampling, and what a collective variable (CV) is
How metadynamics builds a bias from Gaussian hills and why the bias estimates the free energy
Well-tempered metadynamics, convergence, and the importance of choosing a good CV

The idea

MD samples states with probability ∝ e^−F/kT. A barrier of ΔF slows transitions by roughly e^ΔF/kT. For a 10 kT barrier, that is about 20 000 times slower than a barrier-free motion, and 20 kT barriers (about 50 kJ/mol at room temperature) are effectively never crossed in a few-nanosecond run (code tab, Part 1). Enhanced sampling methods change the simulation so that rare events happen, and then correct for the change.

Collective variables. We first choose one or two collective variables s(r): low-dimensional coordinates that distinguish the states and describe the slow motion. Examples: a distance (Li⁺–anion), a coordination number (how many oxygens surround Li⁺), a dihedral, a path coordinate between two structures.

Metadynamics (Laio and Parrinello 2002) then adds a bias V(s) that grows in time: every few hundred steps a small Gaussian "hill" is placed at the current value of s. Visited regions get filled, so the system is pushed to explore new values, crosses barriers and eventually diffuses freely along s. When the wells are filled the bias mirrors the free energy, F(s) ≈ −V(s) + constant (Part 2).

Well-tempered metadynamics (Barducci, Bussi and Parrinello 2008) makes the hills smaller where the bias is already large: w = w₀ e^−V(s)/kΔT. The bias converges smoothly instead of fluctuating, and F(s) = −(T + ΔT)/ΔT · V(s). The bias factor γ = (T + ΔT)/T sets how strongly the barriers are flattened (Part 3).

The CV is everything. If the chosen CV misses the real slow motion, the bias pushes the system in an irrelevant direction and the barrier is never crossed, or only with hysteresis (Part 4 and the widget with CV = y). Signs of a poor CV: few back-and-forth transitions, an estimate that keeps drifting, and different answers from different starting points.

Related methods: umbrella sampling (fixed biases in windows + WHAM/MBAR), adaptive biasing force, replica exchange (temperature or Hamiltonian), and OPES (on-the-fly probability enhanced sampling), the newer method from the Parrinello group.

Key idea. Add Gaussian hills along a good collective variable until the wells are filled; the negative bias is the free energy. Use well-tempered metadynamics, check convergence over time, and question your CV.

The algorithm

choose CV s(r) and its range; hill height w0 (~0.5-1 kT), width σ (~1/3 of the CV fluctuation), stride (~500 steps)
V(s) = 0
loop over MD steps:
    F_total = F_physical - ∂V/∂s · ∇s(r)          # bias force through the chain rule
    integrate
    every stride steps:
        s_now = s(r)
        w = w0 · exp(-V(s_now) / (kB ΔT))           # well-tempered (ΔT = (γ - 1) T); plain: w = w0
        V(s) += w · exp(-(s - s_now)² / (2σ²))
F(s) = -(γ / (γ - 1)) · V(s)   (well-tempered)   or   -V(s) (plain, once the wells are filled)
check: F(s) vs time (converged?), many transitions, reweighting for other observables

In practice use PLUMED, which plugs into LAMMPS, GROMACS, CP2K, OpenMM, AMBER, NAMD and others. A two-line input defines a CV (d: DISTANCE ATOMS=1,2) and the bias (METAD ARG=d SIGMA=0.05 HEIGHT=1.0 PACE=500 BIASFACTOR=10).
Use one or two CVs; the cost of filling grows quickly with dimension.
Compute other averages by reweighting each frame with e^V/kT (with the time-dependent correction for well-tempered runs), or run unbiased simulations in each basin.
Multiple walkers sharing one bias speed up convergence on parallel machines.