Why — and how — we compute a free energy difference that is, on its face, impossible to compute.
A free energy difference is the log of a ratio of partition functions, and you can estimate that ratio — without ever computing either partition function — by sampling one state and asking what its energy would be in the other.
A system held at temperature \(T\) does not sit still in its lowest-energy configuration. Thermal noise constantly bumps it between configurations, and the fraction of time it spends in each one is set by a single weight: the Boltzmann factor \(e^{-\beta E}\), where \(\beta = 1/k_B T\). The lower the energy, the more time spent there. The lower the temperature, the sharper that preference becomes.
Turning these weights into probabilities just needs a normalizer. That normalizer has a name — the partition function:
Don't let the name make it sound deeper than it is. \(Z\) is the sum of all the weights, the thing we divide by so that probabilities add to one. But this innocuous-looking sum will turn out to be the doorway through which every thermodynamic question gets asked — and the thing that will give us trouble.
The widget runs the actual mechanism: 200 particles, each at one of five energy levels you can drag, hopping between levels with Metropolis dynamics. Every frame, each particle proposes moving to another level and accepts with probability \(\min(1, e^{-\beta \Delta E})\). After a few seconds, the populations stop changing — they have found the Boltzmann distribution. The histogram you watch settle into place is one realization of \(P_i = e^{-\beta E_i}/Z\), now arrived at by living the dynamics rather than by writing down the formula.
For a system in thermal contact with a bath at temperature \(T\), the stationary distribution over microstates \(\{i\}\) is
\[ P_i \;=\; \frac{e^{-\beta E_i}}{Z}, \qquad Z \equiv \sum_i e^{-\beta E_i}, \qquad \beta \equiv \frac{1}{k_B T}. \]This is the canonical ensemble. It is derived by treating the bath plus system as a microcanonical (fixed-energy) ensemble and using \(S = k_B \ln \Omega\) and a first-order Taylor expansion of \(\ln \Omega_{\text{bath}}\) in the small system's energy. The Boltzmann factor is what's left after the bath has been marginalized out.
Metropolis dynamics is one way to sample \(P_i\): a Markov chain whose transition rates satisfy detailed balance \(P_i \,W_{i\to j} = P_j \,W_{j\to i}\) with the choice \(W_{i\to j} = \min(1, e^{-\beta(E_j - E_i)})\). The widget's particles are independent walkers on this chain; their pooled occupations form an empirical estimate of \(P_i\) whose variance scales as \(1/\sqrt{N_{\text{particles}}}\).
The Boltzmann distribution from the last step is a complete description of how the system at temperature \(T\) apportions its time across configurations. But "complete" is the problem — we want a single number that summarizes the whole ensemble, something we can compare across systems and watch change with temperature. Energy alone won't do: at finite \(T\) the system isn't sitting at its ground state, and "lower \(\langle E \rangle\) wins" makes the wrong prediction the moment entropy is involved. The right summary quantity must know about both energy and the spread of the distribution.
There is a clean criterion for "which way does the system want to go." The combined system-plus-bath is closed, and a closed system spontaneously increases its total entropy. We need to express that total in terms of quantities the system's distribution actually controls.
In the canonical setup (weak coupling between system and bath), the joint distribution over (system microstate \(i\), bath microstate \(b\)) factorizes as \(P(i, b) = P_i \cdot P_{b\mid i}\), with \(P_{b\mid i}\) uniform over the bath microstates compatible with the system being in state \(i\). The joint Shannon entropy splits cleanly:
where \(S_{\text{sys}} = -k_B \sum_i P_i \ln P_i\) is the system's Shannon entropy and \(\langle S_{\text{bath}}\rangle = \sum_i P_i \, S_{\text{bath}}(E_{\text{tot}} - E_i)\) is the bath's entropy averaged over the system's state. So far this is exact.
The bath's term we get from Boltzmann's counting formula \(S_{\text{bath}}(E) = k_B \ln \Omega_{\text{bath}}(E)\), where \(\Omega_{\text{bath}}(E)\) is the number of bath microstates at energy \(E\). When the system is in microstate \(i\) with energy \(E_i\), the bath holds the rest, \(E_{\text{tot}} - E_i\). Because the bath is huge, \(E_i \ll E_{\text{tot}}\), and a first-order Taylor expansion of the log-density of states around \(E_{\text{tot}}\) is accurate:
The derivative on the right is the bath's inverse temperature by definition: \(\beta \equiv 1/(k_B T) = \partial \ln \Omega_{\text{bath}} / \partial E\). Multiplying through by \(k_B\),
Averaging this over the system's distribution \(\{P_i\}\),
where \(\langle E_{\text{sys}}\rangle \equiv \sum_i P_i E_i\) is the system's average energy. Substituting into the split for \(S_{\text{tot}}\):
The bath baseline \(S_{\text{bath}}(E_{\text{tot}})\) depends only on the (fixed) total energy and is therefore the same for every choice of system distribution. When we optimize over \(\{P_i\}\), it drops out. What remains to be maximized is
Multiply through by \(-T\). Since \(T > 0\), multiplication by a negative number reverses the optimization direction: maximizing the line above is identical to minimizing
That quantity is the free energy: the function that decreases as the system relaxes and sits at its minimum at equilibrium. The minus sign on the entropy term is not a convention — it's the bath's entropy contribution, \(-\langle E_{\text{sys}}\rangle/T\), after multiplication by \(-T\) flipped the entropy maximization into a free-energy minimization.
To see what \(F\) actually evaluates to, plug the Boltzmann distribution \(P_i = e^{-\beta E_i}/Z\) into the definition. The two energy terms cancel themselves:
The equilibrium value of "what the system minimizes" is determined entirely by the partition function. Two properties make \(F = -k_B T \ln Z\) the form worth carrying around. Additivity: independent systems multiply their partition functions, \(Z_{12} = Z_1 Z_2\), and the log turns that product into a sum, \(F_{12} = F_1 + F_2\) — properly extensive, the way an energy must be. Soft minimum: \(-k_B T \ln \sum_i e^{-\beta E_i}\) is a temperature-smoothed minimum of the energies. At \(T \to 0\) it returns the ground-state energy \(\min_i E_i\) exactly; at finite \(T\) it drops below the ground state by an amount that grows with how many states are accessible. Free energy is the system's effective lowest energy, with entropy folded in — both "how deep can I go" and "how many ways can I get there."
Now watch what that quantity actually does. Give a system two places it could live — say, two compartments with different sets of energy levels — and let particles hop freely between them. They do not distribute equally, and they do not distribute by ground-state energy. They distribute by free energy:
Lower \(F\) wins. And the winner can be the compartment with the higher ground state if it has enough accessible states above that ground state. Entropy literally pulls particles in: when a compartment offers many configurations, that compartment's \(F\) is dragged down, and population follows.
This is the working definition: \(F\) is the quantity whose minimum determines equilibrium at constant temperature. It is the only number that matters for "which side wins" when energy and entropy fight.
Both compartments start with the same number of accessible levels, but with different spectra. Drag the lines up and down. Try the energy wins preset: A has a single deep level, B has a single shallow one, so A captures the particles even at moderate \(T\). Then try entropy wins: A still has a slightly lower ground state, but B has a stack of nearby states above its ground. At low \(T\) the ground state still wins; turn \(T\) up and watch the population cross over to B as the entropy of B's stack starts to count.
Treating each particle as independent and distinguishable, with stationary distribution \(P(X,i) \propto e^{-\beta E_i^X}\) over (compartment \(X\), level \(i\)) pairs, the marginal over compartments is
\[ P(X) \;=\; \frac{Z_X}{Z_A + Z_B}, \qquad Z_X \;=\; \sum_i e^{-\beta E_i^X}. \]So the equilibrium ratio of expected populations is \(\langle N_A \rangle / \langle N_B \rangle = Z_A/Z_B = e^{-\beta(F_A - F_B)}\). The two-compartment system is just a Boltzmann distribution over a discrete combined state space, and the marginalization recovers the "competing free energies" picture.
The dynamics implemented in the widget propose, per particle, a uniform random new state across all (compartment, level) pairs and accept Metropolis-style on the energy difference. This satisfies detailed balance with respect to the target distribution.
This is also exactly the structure of two competing thermodynamic states — folded vs. unfolded protein, ligand bound vs. unbound, two ligand conformations — whenever they share the same temperature and exchange a conserved population. Lower \(F\) wins.
We have a target: \(\Delta F = F_B - F_A = -k_B T \ln (Z_B/Z_A)\). Both \(Z\)'s are sums over all microstates of a protein-sized system, and there are roughly \(10^{6000}\) of them. Direct computation is finished as a strategy. But notice that we don't actually want \(Z_A\) or \(Z_B\). We want their ratio. And a ratio admits a trick that a single \(Z\) does not.
Write the ratio and force a factor of the \(A\) distribution to appear inside the integral:
where \(\Delta U(x) = U_B(x) - U_A(x)\) and \(\langle \cdot \rangle_A\) is the expectation under the Boltzmann distribution of state \(A\). Take a log and you have it:
This is the Zwanzig identity (1954). Read what it instructs you to do. Sample the world of state \(A\). At every sample, evaluate the energy of that same configuration as if it lived in state \(B\). The exponentiated negative difference is the weight. Average those weights. Take a log. That is the free energy difference.
You never enumerated microstates. You never computed \(Z_A\) or \(Z_B\). You computed their ratio from a single running ensemble in state \(A\). The "perturbation" in free energy perturbation is exactly this: you treat \(B\) as a perturbation of \(A\), and the perturbation shows up only as the energy difference \(\Delta U\) on the configurations you already have.
Both potentials are harmonic wells: \(U_A(x) = \tfrac{1}{2}x^2\) and \(U_B(x) = \tfrac{1}{2}(x-d)^2 + \Delta\). The walker explores only \(U_A\) — that's the only state we sample. Each visited configuration is evaluated in both wells, and the exponentiated energy gap becomes its weight. The running mean of those weights, log‑converted, drifts onto the true \(\Delta F = \Delta\). Notice something worth pausing on: the true answer is \(\Delta\) no matter what \(d\) is — the wells have the same curvature, so the shift never changes \(Z_B/Z_A\). This isn't a coincidence of the widget; it falls out of the integral itself. Expanding \(\Delta U = U_B - U_A = -dx + \tfrac{1}{2}d^2 + \Delta\) and pulling the \(x\)-independent pieces out of the average, \(Z_B/Z_A = e^{-\beta(\Delta + d^2/2)}\,\langle e^{\beta d x}\rangle_A\). The remaining piece, \(\langle e^{\beta d x}\rangle_A\), is the moment‑generating function of the Gaussian \(p_A\) (mean 0, variance \(1/\beta\)), which evaluates in closed form to exactly \(e^{\beta d^2/2}\) — canceling the \(e^{-\beta d^2/2}\) prefactor precisely, for any \(d\), and leaving \(Z_B/Z_A = e^{-\beta\Delta}\). The shift \(d\) never had anywhere to hide in the final answer. Yet slide \(d\) out toward 2 and the estimate gets visibly noisier even though the target it chases hasn't moved: the cancellation is exact only in the infinite-sample average; with finitely many samples, a large \(d\) means the handful of samples that land far out in \(A\)'s tail (where \(x\) is large) carry disproportionately huge weights, and it's those rare samples the mean is struggling to average over. The shift touches only the variance, not the answer. Hold that thought — it is the whole subject of the next step.
The Zwanzig identity is exact for any pair of states \(A, B\) with the same configuration space. The proof is the line of algebra above — one multiplication by \(1 = e^{-\beta U_A}/e^{-\beta U_A}\) and a relabel of the resulting integrand as an expectation.
Statistically, it is an instance of importance sampling: the proposal distribution is \(p_A\), the target distribution is \(p_B\), and the importance weight at configuration \(x\) is
\[ w(x) \;=\; \frac{p_B(x)\,Z_B}{p_A(x)\,Z_A} \cdot \frac{Z_A}{Z_B} \;=\; e^{-\beta \Delta U(x)} \cdot \frac{Z_A}{Z_B}. \]So \(\langle w \rangle_A = 1\) (correctly normalized weights) and \(\langle e^{-\beta\Delta U} \rangle_A = Z_B/Z_A\), which is the Zwanzig identity. The exponentially-weighted form is what you actually compute; the partition function ratio is what you get out.
The same mathematical structure underlies the ELBO of variational inference, score-based density-ratio estimation, and noise-contrastive estimation. In every case you have two distributions, you can sample one, you want a property of the other, and you accept the bias-variance trade in exchange for tractability.
The Zwanzig identity is exact. But "exact in the infinite-sample limit" is not what you get from a 10‑ns simulation. With finite samples the estimator has variance, and the variance depends on a property of the two distributions that has nothing to do with how clever the math is: their overlap.
If \(p_A\) and \(p_B\) live on roughly the same configurations, samples from \(A\) get weights of order one, and the running mean converges fast. But if \(p_B\) puts substantial probability on configurations where \(p_A\) is vanishingly small — the world of \(B\) is somewhere the walker in \(A\) just never goes — the estimator becomes a lottery. Almost every weight is microscopic. A handful of weights are astronomical, and they carry the entire estimate. The mean of a heavy-tailed distribution converges painfully slowly, if at all, in finite samples.
The standard diagnostic is the effective sample size:
If your run has \(N = 10^6\) samples but ESS = 3, you are computing your free energy from three effective samples regardless of how long you wait. The error bar does not shrink with more sampling. You are not sample-limited — you are overlap-limited, and the wall is geometric.
Start with \(d = 1\) and the estimator looks reasonable. Slide \(d\) past 2 and the overlap region shrinks; past 3 it nearly vanishes. The weight histogram becomes a spike at zero with a far tail; the running estimate stops tracking the truth and jumps every time a rare large-weight sample lands. This is the failure mode of FEP applied directly between two real ligands. Two ligands that differ by even one atom can have non-overlapping conformational ensembles — and direct Zwanzig between them just doesn't work.
The variance of the Zwanzig estimator is bounded below by a quantity that measures the dissimilarity of \(p_A\) and \(p_B\):
\[ \text{Var}\big[\widehat{Z_B/Z_A}\big] \;\ge\; \frac{1}{N}\bigg(\!\int \frac{p_B^2}{p_A}dx \;-\; \Big(\!\int p_B\,dx\Big)^{\!2}\bigg). \]The integral \(\int p_B^2/p_A\,dx\) is finite only when the support of \(p_B\) is contained in the support of \(p_A\) with bounded ratio. In high dimensions, this is almost never the case for two distinct ligands. The integral diverges in the configurations where \(p_A\) is essentially zero but \(p_B\) is not — exactly the configurations that should have huge weights and that the \(A\)-walker, by construction, almost never visits.
The \(\chi^2\)-divergence between \(p_A\) and \(p_B\) bounds this variance and motivates the use of overlap-aware estimators (BAR, MBAR) as well as overlap-aware path design (the next step). The same divergence governs the variance of importance-weighted ELBOs in deep variational inference.
The fix is in the geometry of the path, not in the estimator. If the two distributions don't overlap, build a sequence of distributions between them that do, and walk the staircase.
Introduce a coupling parameter \(\lambda \in [0,1]\) and an interpolated potential:
At \(\lambda = 0\) you have state \(A\); at \(\lambda = 1\), state \(B\); in between, a continuum of synthetic intermediate worlds. Pick a ladder \(0 = \lambda_0 < \lambda_1 < \cdots < \lambda_K = 1\) such that each adjacent pair \((\lambda_i, \lambda_{i+1})\) does overlap. The total free energy difference is a sum:
Notice what that sum is actually doing. Each window's contribution is itself a free energy difference, \(\Delta F_{\lambda_i \to \lambda_{i+1}} = F(\lambda_{i+1}) - F(\lambda_i)\), so the per-window pieces fit together like a chain. Write the \(K = 4\) case out longhand, with \(F_i \equiv F(\lambda_i)\):
The intermediate \(F_1, F_2, F_3\) each appear once with a \(+\) and once with a \(-\) and vanish in pairs:
This is a telescoping sum. The synthetic intermediate worlds drop out of the final answer; only the two real endpoints survive. So the windowing trick is not introducing any new approximation — algebraically it is exactly the same \(\Delta F_{A \to B}\). What it changes is purely the variance of how we estimate it: instead of one impossible long jump, we compute a chain of short, well-overlapped jumps whose estimator errors are small and which sum back up to the same exact quantity.
Each term is a small Zwanzig step between two distributions that do overlap, so each one converges. The impossible long jump has been factored into a sequence of possible short ones. This is the FEP that actually runs in software: a series of MD simulations, one per \(\lambda\) window, with reweighting between adjacent rungs.
The shift is held at \(d = 6\) — the regime where single-shot FEP collapsed in the last widget. At \(K = 1\), you see the same failure: one impossible jump, runaway variance. Raise \(K\): the chain of Gaussians fills in along the path, each pair overlapping its neighbor, and the per-window contributions sum back up to the true \(\Delta F\). The estimator recovers because the path was repaired, not because the estimator was.
The choice of estimator on the ladder is itself a research subject. Forward FEP per window (used here) is the simplest. BAR (Bennett, 1976) combines forward and backward Zwanzig within each window:
\[ \Delta F_i \;=\; k_B T \ln \frac{\big\langle f(\beta(U_{i+1}-U_i) + C)\big\rangle_{i+1}}{\big\langle f(\beta(U_i - U_{i+1}) - C)\big\rangle_{i}} + C, \quad f(x) = (1+e^x)^{-1}, \]and chooses \(C\) self-consistently to minimize variance. MBAR (Shirts & Chodera, 2008) generalizes to all windows at once via reweighting across the full ladder.
Thermodynamic integration takes the same ladder but a different identity: \(\Delta F = \int_0^1 \langle \partial_\lambda U \rangle_\lambda\,d\lambda\), computed by quadrature over the windows. TI is a derivative-based method; FEP/BAR/MBAR are ratio-based. The path is the same.
Modern variants: replica exchange across the ladder (HREX, Hamiltonian replica exchange), nonequilibrium switching (Crooks fluctuation theorem, Jarzynski), and learned \(\lambda\)-paths from generative models. All inherit this fundamental architecture of "bridge non-overlapping distributions with a chain."
Now the payoff. Nothing in the construction required \(\lambda\) to interpolate between two physical states. The whole machinery still works if \(U_A\) and \(U_B\) differ by an arbitrary, even unphysical, transformation — one ligand's atoms smoothly morphing into another's, charges turning on and off, methyl groups appearing from vacuum. These nonphysical paths are alchemical, and they are what makes FEP useful for drug design.
In practice, the medicinal chemist's question is not "what is the absolute binding free energy of compound X?" but "does this proposed modification \(A \to B\) bind tighter?" That answer requires a difference of differences: \(\Delta\Delta G_{\text{bind}}^{A \to B} = \Delta G_{\text{bind}}^B - \Delta G_{\text{bind}}^A\). The physical legs of this calculation — ligand wandering through bulk water until it docks — require timescales no MD simulation will ever cover.
So we close a cycle and let the impossible legs cancel.
Stage through the cycle. Four states sit at the corners: ligand \(A\) free in water, ligand \(B\) free in water, ligand \(A\) bound to the protein pocket, ligand \(B\) bound. The horizontal arrows are physical binding events — the thing we actually care about, and the thing we cannot simulate. The vertical arrows are alchemical transformations: \(A \to B\) in water, and \(A \to B\) in the pocket. At stages 4 and 5 the widget actually runs those two transformations — each is a live windowed-FEP calculation exactly like Widget 5, a ladder of six \(\lambda\)-windows with a Metropolis walker in each, reweighting between neighbors. The two ΔG numbers you watch converge are not scripted; they are running estimates accumulating from the sampling. (The pocket is modeled as a stiffer well than bulk water, which changes the per-window variance but — as you saw in Widget 3 — not the leg's true ΔG.) Because free energy is a state function, going around the cycle returns zero, which means
The impossible legs cancel. Two alchemical calculations — the same ligand transformation done in two different environments — give the experimentally meaningful difference in binding affinity, often within 1 kcal/mol of measurement. This is, mechanically, how every relative binding free energy calculation in production drug discovery works.
Free energy is a thermodynamic state function: the integral of \(dG\) around any closed cycle vanishes. Concretely, going around the box A‑water → A‑complex → B‑complex → B‑water → A‑water gives
\[ \Delta G_{\text{bind}}^A + \Delta G_{\text{complex}}^{A\to B} - \Delta G_{\text{bind}}^B - \Delta G_{\text{solvent}}^{A\to B} \;=\; 0, \]which rearranges to the boxed identity. The "alchemical" legs are simulated using the FEP/BAR/MBAR pipeline on a dual-topology or single-topology Hamiltonian where the molecular mechanics force field smoothly interpolates between \(A\) and \(B\). To avoid singularities as atoms appear or vanish, soft-core potentials replace the bare Lennard-Jones interaction near the endpoints (Beutler 1994; Steinbrecher 2007).
In practice the windows are run with replica exchange (HREX), trajectories are checked for convergence by overlap matrices, and the cycle is closed by comparing forward and reverse runs. Modern pipelines (Schrödinger FEP+, OpenFE, Amber TI) automate the dual-topology setup and ladder design; on standard benchmark sets they reach RMSE ≈ 1 kcal/mol against experimental \(\Delta\Delta G\), enough to guide medicinal-chemistry decisions on a series of analogs.
Six steps that arrived where they began.
Free energy is the quantity that decides equilibrium at fixed temperature: \(F = -k_B T \ln Z\). Differences in it are log‑ratios of partition functions, \(\Delta F = -k_B T \ln(Z_B / Z_A)\). A ratio of partition functions admits a trick no single \(Z\) does — it equals a re‑weighted average under either endpoint, the Zwanzig identity. That estimator works when the two distributions overlap and dies when they don't. The cure is geometric: bridge them with a chain of intermediate distributions, each adjacent pair overlapping, so the impossible jump factors into possible steps. Once that machinery exists, \(\lambda\) does not have to interpolate between physical states — it can interpolate between any two Hamiltonians. Alchemical transformations let us close a thermodynamic cycle in which the impossible physical legs of binding cancel out, leaving the experimentally meaningful \(\Delta\Delta G_{\text{bind}}\) as a difference of two computable alchemical legs.
The whole subject in one sentence again: a free energy difference is the log of a ratio of partition functions, and you can estimate that ratio — without ever computing either — by sampling one state and asking what its energy would be in the other.