ALA: the free energy along phi¶
The potential of mean force along the backbone torsion phi, from the 1 us cMD run and from the cold rung of the REST2 ladder.
What a PMF here is. F(phi) = -kT ln p(phi), with p the histogram of the sampled angle.
It is a free energy along ONE coordinate, with everything else averaged over -- not an energy
barrier, and not a rate. Two runs that sample the same ensemble must give the same PMF; that is
the comparison this page is for.
Empty bins stay empty. A bin nothing visited has no estimate, and tutorial_data.pmf leaves
it NaN rather than substituting a large number. Filling it would draw a barrier whose height is
set by the sample size.
import sys, pathlib
sys.path.insert(0, str(pathlib.Path.cwd().parents[1] / "_analysis")) # docs/tutorial/_analysis
import numpy as np
import matplotlib.pyplot as plt
import mdtraj as md
import tutorial_data as td
# Where the runs are. NOT defaulted: see docs/tutorial/_analysis/tutorial_data.py.
# The resolved path is deliberately NOT printed. These notebooks are committed WITH
# their outputs, so an echoed absolute path would put one machine's layout into a
# published file -- which this project forbids, and which executing a notebook turns
# from a parameter into a stored string.
print("tutorial runs:", "found" if td.runs_root().is_dir() else "MISSING")
tutorial runs: found
The cMD run¶
One microsecond, 100 000 frames. compute_phi returns the torsion for every residue that has
one; alanine dipeptide has exactly one.
cmd = td.solute("ALA-1us", stride=1)
phi_atoms, phi = md.compute_phi(cmd)
psi_atoms, psi = md.compute_psi(cmd)
print(f"{cmd.n_frames} frames, phi over atoms {phi_atoms[0]}, psi over atoms {psi_atoms[0]}")
centres, F = td.pmf(phi[:, 0])
print(f"span {np.nanmax(F):.1f} kJ/mol over {np.isfinite(F).sum()} populated bins of {F.size}")
100000 frames, phi over atoms [ 4 6 8 14], psi over atoms [ 6 8 14 16] span 23.5 kJ/mol over 52 populated bins of 72
The REST2 cold rung¶
State 0 of the ladder is the PHYSICAL ensemble -- tau = 0, unscaled -- so its PMF is directly comparable with the cMD one. The hotter states are not: they sample a scaled Hamiltonian and their marginals are not the physical free energy.
rest2 = td.solute("ALA-REST2", stride=1,
trajectory="REST2-run1/solute_state0_prod1.nc")
_, phi_r = md.compute_phi(rest2)
centres_r, F_r = td.pmf(phi_r[:, 0])
print(f"REST2 state 0: {rest2.n_frames} frames, span {np.nanmax(F_r):.1f} kJ/mol")
REST2 state 0: 5000 frames, span 16.2 kJ/mol
fig, ax = plt.subplots(figsize=(7.5, 4.4), layout="constrained")
ax.plot(np.rad2deg(centres), F, lw=2, label=f"cMD, 1 us ({cmd.n_frames} frames)")
ax.plot(np.rad2deg(centres_r), F_r, lw=2, ls="--",
label=f"REST2 state 0 ({rest2.n_frames} frames)")
ax.set_xlabel("phi (degrees)"); ax.set_ylabel("F (kJ/mol)")
ax.set_title("Alanine dipeptide: free energy along phi")
ax.legend(); ax.grid(alpha=0.3)
Reading it¶
The two curves should agree wherever BOTH runs sampled. Where they diverge, look at which one has the emptier histogram there before concluding anything about the physics -- a PMF is only as good as the occupancy of its bins, and the shorter run's tail is noise long before it is a discrepancy.
This notebook makes no claim about convergence. Showing that a PMF has converged needs repeats or a block analysis, which is a separate question from how to compute one.