Lab: free energies by hand

The setup

Why a special method? The free energy F = −k_BT ln Z counts how many states are accessible, not just how low they are. An ordinary simulation gives averages like ⟨U⟩; the free energy difference between two states can be written as an average too, but of an exponential, which brings its own problems.

Step 1

Potential of mean force from g(r)

Why: g(r) says how much more often two atoms are found at distance r than in an ideal gas. Turning this probability into an energy, w(r) = −k_BT ln g(r), gives the free energy of the pair at that distance, averaged over everything else: the potential of mean force (PMF). For ion pairs it tells whether contact pairs or solvent-separated pairs are preferred.

What the numbers say: In a gas the PMF is just the pair potential. In the liquid it is not: the well is shallower (−0.71 instead of −1.00), and at 1.675, where the bare potential is still attractive, the PMF is repulsive (+0.35), because the pair at that distance has to push away the neighbours that would like to sit between them. The PMF contains the solvent. Merit: it costs nothing beyond g(r). Limitation: it only works where g(r) is well sampled; at high PMF values (rare distances) g is tiny and noisy, which is why barriers need the biased methods of steps 5–6.

Step 2

Free-energy perturbation by hand

Why: Zwanzig's formula gives the free-energy change from state 0 to state 1 using samples of state 0 only: ΔF = −k_BT ln ⟨e^(−ΔU/k_BT)⟩₀, where ΔU = U₁ − U₀ is evaluated on each sample. Test it on a spring whose stiffness doubles, where the exact answer is known.

What the numbers say: Six samples already give 0.365 against the exact 0.347. The free energy rises although the average energy does not change at all (½k_BT in both states, by equipartition): the stiffer spring confines the particle to a narrower range, and the whole ΔF is entropy. This is why free energies cannot be computed from energy averages alone. Note the sample at x = 1.90 contributes only 0.16 to the average: FEP is dominated by the samples where U₁ is low.

Step 3

The overlap problem

Why: FEP only works if the samples of state 0 include the configurations that matter for state 1. Make the spring 100 times stiffer and run FEP in both directions: forward (sample the soft spring, evaluate the stiff one) and backward (sample the stiff spring, evaluate the soft one).

What the numbers say: Forward works: the soft spring visits the narrow region the stiff one needs. Backward fails badly, and more samples barely help (1.35 after 100 000 samples, still 40 % low), because the stiff spring never visits |x| > 1, where the soft spring has a third of its weight: the missing part is never seen, so the estimate is biased, not just noisy. The cures: go in the direction from the broad to the narrow distribution, split a large change into stages so that neighbours overlap (2.18 with 4 stages), and use both directions together (the Bennett acceptance ratio, BAR/MBAR, the standard in GROMACS and alchemical tools).

Step 4

Charging an ion: FEP fails, TI works

Why: Solvation free energies are computed by slowly switching on an ion's charge. A simple model of the solvent: the electrostatic potential φ at the ion site fluctuates like a Gaussian with spread s; the energy of charge q is qφ. Linear response gives the exact answer, so we can test one-step FEP and thermodynamic integration (TI), which integrates the average slope ⟨∂U/∂λ⟩ along a switching parameter λ.

What the numbers say: One-step FEP is off by 300 kJ/mol even with a million samples: the charged ion needs solvent configurations 20 standard deviations away from what the neutral solute ever sees. TI with just three λ points gets −500.2 against −501.2, because it only needs ordinary averages at each λ, where the samples are typical. Staged FEP works for the same reason. In this model ⟨φ⟩ is a straight line in λ, so three points suffice; for real ions (and especially for switching off LJ, where soft-core potentials are needed) the curve bends and needs 10–20 λ windows. A value of −500 kJ/mol is the right order of magnitude for the hydration of a small ion such as Na⁺ (about −400 kJ/mol), and a real calculation also needs the charged-box correction of Lab 5.5.

Step 5

Umbrella sampling by hand

Why: A barrier top is visited too rarely to measure its probability. Umbrella sampling adds a spring w(s) = ½k(s − s₀)² that holds the system near s₀, measures the histogram there, and removes the known bias afterwards: F(s) = −k_BT ln P(s) − w(s) + constant.

What the numbers say: Without the spring, the barrier top would be visited a fraction e^(−5) of the time; with it, 10⁵ samples per bin. After removing the bias, the window reproduces the shape of the barrier (the F + C₀ column falls from 1.93 to 0.80 as the exact U falls from 4.98 to 3.85), but only up to an unknown constant C₀. Neighbouring windows are connected through bins they both sample. By hand we matched them at one bin and got −3.94 against −4.02. Limitation: the matching relies on the overlap; a bin seen rarely by one of the windows makes the shift noisy.

Step 6

WHAM: all windows at once

Why: Matching windows one bin at a time wastes data and accumulates errors. The weighted histogram analysis method (WHAM) finds the constants of all windows together, so that every bin uses all the windows that visited it, weighted by how much they saw. How many windows, and how long, does a 5 k_BT barrier need?

What the numbers say: Umbrella sampling with WHAM recovers the 5 k_BT barrier, which a plain run of the same length crosses only three times (far too few to measure anything). With 5 windows spaced 0.6 apart the overlap is poor and the result is high by 0.3; with 17 windows the error drops to ±0.14 for the same total cost as 9 long windows: overlap matters more than length. Merits: every window is an ordinary, independent simulation (easy to run in parallel), and the method is the same as metadynamics' result in Lab 7.4 (5.00), from a different route. Limitations: the collective variable must be chosen in advance, as in metadynamics, and slow motions orthogonal to it make windows converge slowly. In practice: PLUMED for the biases, gmx wham or WHAM by Grossfield for the analysis.

Summary