Paracetamol: free energy along each torsion¶
Three torsions, three PMFs, from the 1 us cMD run and the cold REST2 rung.
Which torsions, and why these three. Paracetamol has one amide, one ring-to-amide bond and one
hydroxyl. The REST2 scaler protects the amide omega and all six aromatic ring bonds; what it heats
are the bonds OUTSIDE the ring. So aryl and phenol are the two coordinates REST2 is actually
working on, and omega is the control -- it should look the same in both runs because neither
run heats it.
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
PARA_TORSIONS = {
# name: (A, B, C, D) -- indices into the solute-only topology AND the TYL.sdf, which this
# notebook verifies are the same atom ordering before using either.
"omega": (0, 1, 3, 4), # C1-C2-N1-C3 the amide, protected by the REST2 convention
"aryl": (1, 3, 4, 5), # C2-N1-C3-C4 the acetamido group vs the ring. ORDINARY torsion:
# its central bond 3-4 is NOT part of the ring
"phenol": (6, 7, 8, 17), # C5-C6-O2-H7 the hydroxyl rotation
}
cmd = td.solute("PARA-1us", stride=1)
quads = list(PARA_TORSIONS.values())
names = list(PARA_TORSIONS)
tors = md.compute_dihedrals(cmd, quads)
print(f"{cmd.n_frames} frames, {tors.shape[1]} torsions")
fig, axes = plt.subplots(1, 3, figsize=(14, 4.0), layout="constrained")
for ax, name, column in zip(axes, names, range(tors.shape[1])):
centres, F = td.pmf(tors[:, column])
ax.plot(np.rad2deg(centres), F, lw=2)
ax.set_title(f"{name}\n{np.isfinite(F).sum()}/{F.size} bins, span {np.nanmax(F):.1f} kJ/mol")
ax.set_xlabel(f"{name} (degrees)"); ax.grid(alpha=0.3)
axes[0].set_ylabel("F (kJ/mol)")
100000 frames, 3 torsions
Text(0, 0.5, 'F (kJ/mol)')
omega is the narrow one. It occupies a fraction of the circle because an amide is not a
rotatable bond on this timescale -- which is exactly why the REST2 convention leaves it unscaled.
A torsion with no second state to reach gains nothing from being heated and would only distort.
aryl and phenol populate the whole circle, with two-fold structure: the ring's two faces are
chemically identical, so the molecule has a symmetry that makes the lobes equivalent. That is the
subject of the clustering and symmetry notebooks.
rest2 = td.solute("PARA-REST2", trajectory="REST2-run1/solute_state0_prod1.nc")
tors_r = md.compute_dihedrals(rest2, quads)
fig, axes = plt.subplots(1, 3, figsize=(14, 4.0), layout="constrained")
for ax, name, column in zip(axes, names, range(tors.shape[1])):
c1, F1 = td.pmf(tors[:, column]); c2, F2 = td.pmf(tors_r[:, column])
ax.plot(np.rad2deg(c1), F1, lw=2, label=f"cMD ({cmd.n_frames})")
ax.plot(np.rad2deg(c2), F2, lw=2, ls="--", label=f"REST2 state 0 ({rest2.n_frames})")
ax.set_title(name); ax.set_xlabel(f"{name} (degrees)"); ax.grid(alpha=0.3)
axes[0].set_ylabel("F (kJ/mol)"); axes[0].legend(fontsize=8)
<matplotlib.legend.Legend at 0x7e88bbeded20>
The cold rung of a REST2 ladder samples the UNSCALED Hamiltonian, so its PMF is the physical one and should agree with cMD. It has far fewer frames here, so expect its tails to be noisier -- that is sample size, not a disagreement about the free energy.