Free-energy and population maps¶

For one-dimensional MLMD/metadynamics profiles from FES, COLVAR, or HILLS files, use the MLMD free-energy tutorial and Plot.fes_1d. This page covers the existing two-dimensional population, XPM, and carbohydrate-sphere plotting interfaces.

A trajectory can visit several regions of a coordinate space. A population map shows where samples accumulate; a relative free-energy map expresses those counts on an energy scale. This tutorial covers Plot.fes(), Plot.contour(), and Plot.cremer_pople_sphere() and explains what their colors represent.

All examples use synthetic data to make the plotting mechanics reproducible. For a physical interpretation, use equilibrated samples with the correct weights. Neither fes() nor the sphere method removes a metadynamics bias or performs statistical reweighting.

Turn two sampled observables into a landscape¶

For a bin with population N, fes() computes relative values equivalent to:

ΔF = -R × T × ln(N / Nmax)

The most populated bin is the reference with zero free energy. Empty bins have no finite estimate and are masked. This is a count-based estimate on a fixed 72 × 72 grid; bins and weights are not configurable through this method.

import numpy as np
from matmalib import Data, Plot

rng = np.random.default_rng(25)
basin_a = rng.multivariate_normal([-1.0, -0.6], [[0.30, 0.12], [0.12, 0.35]], 24000)
basin_b = rng.multivariate_normal([1.0, 0.8], [[0.25, -0.08], [-0.08, 0.25]], 12000)
samples = np.vstack((basin_a, basin_b))
data = Data(path=".", data=np.column_stack((np.arange(len(samples)), samples)))

plot = Plot()
plot.cmap(color_num=8, map="ice", reverse=True)
plot.set_config({"xlabel": "Coordinate 1", "ylabel": "Coordinate 2",
                 "title": "Two sampled basins at 300 K"})
with np.errstate(divide="ignore"):
    plot.fes(data, cols=[1, 2], temp=300, num_levels=8)
plot.ax.set_facecolor("#f1f5f9")
plot.ax.set_aspect("equal", adjustable="box")
plot.fig.axes[-1].set_ylabel("Relative free energy (kJ/mol)")

The darker/lighter regions follow the selected palette: read the colorbar to determine which are low energy. The low-energy regions should lie near the two input Gaussian centers. The gray background marks areas without finite plotted estimates. Changing the temperature changes the energy scale, while changing the sample populations changes the relative basin depths.

np.errstate suppresses only the expected log(0) warning from empty bins. It does not fill them with a guessed energy. The implementation uses R = 8.314 J/(mol K) and divides by 1000, so values are kJ/mol, although its original colorbar label says only kJ. The example makes the molar unit explicit.

cols=[1, 2], temp=300, and num_levels=8 are the defaults. Column 0 has no special significance to fes() itself; we use it for a sample index by convention. Supply finite numeric samples with variation in both selected columns. num_levels is passed to Matplotlib's contour locator; it is not guaranteed to produce exactly that many bands.

The contour coordinates use the lower bin edges rather than their centers. If you need custom bins, weights, uncertainty estimates, or exact bin-center placement, calculate the histogram and free energy with NumPy, then plot the result directly with Matplotlib.

Plot matrices already containing energies¶

contour() displays a matrix or a list of matrices on a shared color scale. It does not convert populations to energy. This is useful for values loaded from supported GROMACS XPM files, or from your own calculation.

import numpy as np
from matmalib import Data, Plot

y, x = np.mgrid[-2:2:45j, -2:2:45j]
energy_a = 2.0 * ((x + 0.6) ** 2 + 0.7 * (y - 0.2) ** 2)
energy_b = 2.0 * ((x - 0.5) ** 2 + 0.7 * (y + 0.3) ** 2)
matrices = [Data(data=energy_a, path="."), Data(data=energy_b, path=".")]

plot = Plot()
plot.cmap(color_num=8, map="ice", reverse=True)
plot.set_config({"xlabel": "Column index", "ylabel": "Row index"})
plot.contour(matrices, limit=12)
for ax, title in zip(plot.ax, ["Surface A", "Surface B"]):
    ax.set_title(title)
    ax.set_aspect("equal")
for ax in plot.fig.axes[2:]:
    ax.set_ylabel("Input energy units")

Both panels share nine contour boundaries from zero through limit (eight bands). Matrix positions are row and column indices, not coordinates from the original file. The same scale allows direct color comparisons between panels when the inputs use the same energy units and reference.

The current masking rule hides values greater than limit - 1, even though the colorbar extends to limit. With the default limit=16, values above 15 are masked. Use Matplotlib directly if you need another clipping rule, negative contour levels, or explicit physical x/y coordinates.

For an existing supported XPM file, the equivalent local workflow is:

from matmalib import Data, Plot

matrix = Data.from_file("free-energy.xpm", software="gromacs")
plot = Plot()
plot.contour(matrix, limit=16)
plot.fig.savefig("free-energy.svg", bbox_inches="tight")

The XPM reader handles a narrow single-character, square-matrix format. Check the resulting matrix and units against the input; see file formats and limitations.

Map a six-membered ring onto the Cremer–Pople sphere¶

For a six-membered ring, the Cremer–Pople angles theta and phi locate a puckering direction. cremer_pople_sphere() displays their distribution with a Mollweide projection. It accepts the output of cremer_pople_trajectory(), or any correctly named Data object with the same contract.

This standalone example creates angle samples around two regions. These are illustrative populations, not an inferred molecular trajectory.

import numpy as np
from matmalib import Data, Plot

rng = np.random.default_rng(31)
theta = np.r_[rng.normal(35, 12, 2200), rng.normal(110, 15, 1000)]
theta = np.clip(theta, 0, 180)
phi = np.r_[rng.normal(70, 20, 2200), rng.normal(240, 25, 1000)] % 360
data = Data(
    path=".",
    data=np.column_stack((np.arange(len(theta)), theta, phi)),
    columns=("time", "monomer_0.theta", "monomer_0.phi"),
    angle_unit="degrees",
)

plot = Plot()
fig, ax = plot.cremer_pople_sphere(
    data, monomer=0, bins=(36, 18), temp=300,
    energy_unit="kcal/mol", cmap="viridis", num_levels=8,
    title="Synthetic ring-pucker populations",
)
print("Samples counted:", plot.conformer_counts.sum())
print("Stored grid shape:", plot.free_energy.shape)
print("Minimum relative free energy:", float(plot.free_energy.min()), "kcal/mol")

The output should count all 3200 rows, store an 18 × 36 grid, and have minimum energy 0. The two islands correspond to the generated angle regions. Change bins=(36, 18) to (72, 36) to use the default five-degree cells: smaller bins resolve more detail but have fewer samples each.

The method expects columns named monomer_0.theta and monomer_0.phi for monomer=0. Select another zero-based monomer index only when matching column names exist. angle_unit must be "degrees" or "radians"; theta must lie in [0, 180] degrees or [0, pi] radians. Phi wraps around its full circle. Undefined angles are omitted; a pure chair at a pole is counted with an arbitrary common phase of 180 degrees.

Read the map correctly¶

The horizontal and vertical labels show Cremer–Pople angles in degrees, although the geographic projection uses longitude 180 - phi and latitude 90 - theta. Free energy is relative to the most populated angular cell, using the same -RT log(N/Nmax) convention as fes().

The default energy_unit is "kcal/mol", unlike fes(); "kJ/mol" is also supported. These are counts per equal angular cell. There is no spherical-area correction: cells near the poles occupy different areas on a sphere than cells near its equator. Do not interpret this automatically as a probability density per unit solid angle.

Contours interpolate between occupied bins for display. White regions indicate unsampled locations; interpolated colors do not create additional observations. Numerical analysis should use the exact stored arrays:

Attribute	Meaning
`conformer_counts`	Integer histogram in displayed `(latitude, longitude)` order
`free_energy`	Masked energy array in the same order; `None` in count mode
`latitude_edges`, `longitude_edges`	Bin edges in projection coordinates, in radians
`heatmap`, `colorbar`	Matplotlib contour object and colorbar
`fig`, `ax`	Figure and Mollweide axes

Pool runs and show counts¶

Multiple input datasets are pooled into one histogram, not overlaid as separate colored surfaces. Every valid row contributes one observation regardless of its timestamp or the source trajectory's sampling interval. Equalize sampling or account for weights outside this method if trajectories should contribute differently.

import numpy as np
from matmalib import Data, Plot

rng = np.random.default_rng(6)
runs = []
for center in (120, 200):
    theta = rng.uniform(45, 120, 500)
    phi = rng.normal(center, 20, 500) % 360
    runs.append(Data(
        data=np.column_stack((theta, phi)), path=".",
        columns=("monomer_0.theta", "monomer_0.phi"), angle_unit="degrees",
    ))

plot = Plot()
plot.set_colors(["#e0f2fe", "#7dd3fc", "#0284c7", "#0c4a6e"])
plot.cremer_pople_sphere(runs, quantity="counts", bins=(24, 12),
                         headers=["run A", "run B"], alpha=0.9,
                         pucker_labels=True)
print("Pooled samples:", int(plot.conformer_counts.sum()))

The colorbar now reports numbers of conformers rather than energy. pucker_labels=True annotates 38 ideal puckers. Those labels assume the ring atom order O, C1, C2, C3, C4, C5; “O” corresponds to ring position 6 in classification results. For another ordering, do not interpret the label names as if they used that convention. See puckering coordinates and classification.

headers must have one entry per dataset and becomes a combined title unless you pass title. alpha accepts one shared number or a list of identical values; pooled datasets cannot use different opacities. Temperature must be finite and positive even in count mode. bins accepts a positive integer for both dimensions or a (phi_bins, theta_bins) pair. num_levels must be a positive integer.

The map has fixed geographic limits: remove xrange and yrange from the plot configuration. figsize=(8, 4) is the default. A Matplotlib cmap name or object supplied directly overrides the palette previously selected with Plot.cmap() or Plot.set_colors() for that draw.

Use the trajectory shortcut¶

The following call delegates to the sphere method and returns (fig, ax):

fig, ax = plot.trajectory(
    pucker_data,
    projection="mollweide",
    monomer=0,
    bins=(72, 36),
    quantity="free_energy",
    energy_unit="kcal/mol",
    pucker_labels=True,
)

This shortcut forwards monomer, pucker_labels, cmap, bins, num_levels, quantity, temp, energy_unit, title, alpha, and headers. The usual col, var_name, and hist arguments do not apply. Moving-average smoothing, overlap=True, and calc_qa=True are unsupported in this view. projection=None or "rectilinear" selects the usual time-series plot instead.

Continue with coordinate-based trajectory analysis to replace the synthetic angles with calculated ring puckers, or plot styling and saving to finish a figure.