How to plot a free-energy surface from an MSM#

Goal#

Render a 2D free-energy surface on any pair of projected dimensions (typically the slowest TICA components) from a fitted Model, optionally overlaying the macrostate assignments so you can see which basin each macrostate occupies.

Minimal example#

model.plotFES(0, 1, temperature=360)

Free-energy surface on TIC1 / TIC2 for the villin folding system

plotFES() bins the projected data along dimX and dimY and accumulates the stationary distribution there: each microstate marks the bins it visits and contributes its equilibrium population pi_i to those bins (so heavily-sampled but kinetically rare microstates are correctly down-weighted). The accumulated mass p is then converted to a free-energy surface via -kT log(p) and drawn as a filled contour map. Above is the FES from the villin folding tutorial on TIC1 / TIC2 - the deep central basin around (0.5, -0.5) is the folded state; the surrounding higher-energy contours are unfolded / partially-folded conformations.

Parameters that matter#

Parameter	What it does
`dimX`, `dimY`	Indices of the projected dimensions to use as axes. For a TICA-reduced model, `(0, 1)` is the two slowest components.
`temperature`	Simulation temperature in K. Used in the `kT` factor for converting populations to kcal/mol. Match the production temperature.
`states`	`False` (default) for a clean FES; `True` to overlay microstate scatter points coloured by macrostate.
`cmap` / `fescmap` / `statescmap`	Matplotlib colormaps. `cmap` sets both at once; the split lets you override the FES and state colour scales independently.
`levels`	Number of contour levels in the FES. Default 7.
`s`	Scatter-marker size when `states=True`. Default 10. Larger when you have a few macrostates and want their assignments to read clearly.
`plot`, `save`	`plot=True` (default) opens the figure; `save="fig.png"` writes to disk. With both `plot=False` and `save=None` the function returns `(fig, ax)` so you can post-process.
`data`	Pass a different `MetricData` than the one used to build the model - e.g. project the FES onto RMSD even though the model was built on contacts.

Common variations#

Overlay macrostate assignments#

model.plotFES(0, 1, temperature=360, states=True)

Free-energy surface with macrostate scatter overlay

Each microstate is drawn as a small marker coloured by its macrostate. This is the right diagnostic for “did PCCA produce sensible basins?” - macrostates should cluster within the deep wells, not spray across the saddle regions.

Save the figure to disk without opening a window#

model.plotFES(0, 1, temperature=298, save="./fes.png", plot=False)

Useful in scripted pipelines where you want the PNG but don’t want to pop a matplotlib window.

Plot the FES on a different feature space than the model was built on#

from moleculekit.projections.metricrmsd import MetricRmsd
from moleculekit.projections.metricgyration import MetricGyration

# Build the model on contacts, but plot the FES on RMSD-to-native vs radius of gyration
ref = Molecule("crystal.pdb")
m = Metric(model.data.simlist)
m.set([
    MetricRmsd(ref, "protein and name CA"),
    MetricGyration("protein"),
])
alt_data = m.project()

model.plotFES(0, 1, temperature=298, data=alt_data)        # dim 0 = RMSD, dim 1 = Rg

Two distinct 1-D projections stacked through a single Metric give you a (n_frames, 2) feature array - exactly what plotFES(0, 1, ...) needs. data= must have the same simlist as the model (same number of frames, same ordering); only the per-frame feature values differ. Lets you keep a single Markov model while exploring its energy landscape in different coordinate systems.

Custom colour map and contour levels#

import matplotlib.pyplot as plt

model.plotFES(0, 1, temperature=298,
              fescmap=plt.cm.gray, levels=15)

More contour lines reveal finer basin substructure; greyscale is sometimes clearer in printed figures.

Post-process the figure#

fig, ax = model.plotFES(0, 1, temperature=298, plot=False, save=None)
ax.set_title("My system - 298 K free-energy surface")
ax.axhline(0, color="white", linewidth=0.5)
fig.savefig("./fes_annotated.png", dpi=300, bbox_inches="tight")

With both plot=False and save=None, plotFES returns the (fig, ax) pair instead of displaying. Hand back to matplotlib for any further customisation.

Gotchas#

The FES is computed from microstate populations, not raw frame populations. Sparse microstates (very few visits) contribute noisy bins - their free energies can look artificially high. Cluster with enough microstates that each one is hit by many frames before reading absolute energy values off the FES.
temperature= matters because the FES is in kcal/mol via -kT log(p). Use the production temperature - mixing 300 K and 360 K runs into one model and asking for one temperature makes the y-axis numerically meaningless.
Axis labels come from data.description if present, otherwise they’re generic Dimension N strings. After TICA the projection description does carry useful labels (TIC1, TIC2, …) but custom projections often have empty descriptions.
dimX, dimY must be valid indices into the projected data. For a TICA-reduced model with tica.project(3), valid indices are 0, 1, 2.
Don’t read kinetics off the FES. Basin depths tell you equilibrium populations; the transition timescales between them come from plotTimescales() and Kinetics. A shallow well can have a slow kinetic exit if the saddle is high.