Skip to content

Plot trajectories, comparisons, and energy levels

Plot wraps Matplotlib for the arrays and molecular objects used by matmalib. This tutorial starts with generated data so every interactive example can run without uploading a simulation. Use the same calls with a Data object returned by the file readers.

By the end, you will be able to inspect a trajectory and its distribution, compare observables, customize a figure, and export it. For two-dimensional populations, continue to free-energy maps; for frequencies, see IR spectra.

Understand the plotting data contract

Data is a lightweight container: it stores the keyword attributes you give it without validating them. Most plotting methods read data.data, a NumPy array, and some read data.path, the directory used by savefig().

Method Input Figure access
trajectory() Data or list of Data; time in column 0, observable in col; path and time_unit required plot.fig, plot.ax[0] and histogram plot.ax[1]
scatter() Data or list of Data; x in column 0, y series in remaining columns plot.fig, plot.ax
fes() Data; two selected columns of samples and path plot.fig, plot.ax
contour() Data or list of Data; two-dimensional matrices and path plot.fig, plot.ax (array for multiple panels)
foo_plot() Data; whitespace-separated text rows and path plot.fig, plot.ax
cremer_pople_sphere() Named six-membered-ring angle data Returns (fig, ax) and retains both
snfg() Glycan Returns (fig, ax) and retains both
cluster_dendrogram() ClusterResult Returns (fig, ax) and retains both
space() List of Space objects Current Matplotlib figure; not retained by Plot
ir() State with Freq and Ints arrays Current Matplotlib figure; not retained by Plot

The older plotting methods return None. Call the method first, then edit plot.ax or save plot.fig. An interactive block renders its open Matplotlib figures automatically.

Read a trajectory and its histogram

Here the first column contains time in femtoseconds, and the second contains a synthetic distance in ångströms. The signal moves between two regions, so the histogram should have two broad peaks.

import numpy as np
from matmalib import Data, Plot

rng = np.random.default_rng(17)
time_fs = np.arange(1200) * 2.0
distance = 2.6 + 0.45 * np.sin(np.linspace(0, 12, len(time_fs)))
distance += rng.normal(0, 0.07, len(time_fs))
trajectory = Data(
    data=np.column_stack((time_fs, distance)),
    time_unit="fs",
    path=".",
)

plot = Plot()
plot.trajectory(trajectory, col=1, var_name="Distance (Å)",
                title="A synthetic collective variable", headers=["run A"])
print("Mean distance:", round(float(distance.mean()), 3), "Å")

The left panel preserves the order of observations. The right panel counts values using Matplotlib's rice bin rule; it is a count histogram, not a probability density. A populated region does not tell you how often the system transitions into it—use the time trace for that.

Try increasing the noise standard deviation from 0.07 to 0.2. The two peaks become less distinct even though the slowly varying part of the trajectory stays the same.

trajectory() applies these time conversions; it does not convert the observable's units:

data.time_unit Displayed x axis
"fs" Time divided by 1000, in ps
"ps" Time divided by 1000, in ns
"ns" Time unchanged, in ns
"Residue" Index unchanged, labeled residue

Use at least two rows, a valid numeric col, and one of the supported time-unit strings. The stepsize label is calculated from the first two rows. hist=False removes the histogram panel; for residue-wise RMSF, set time_unit="Residue" and use residue indices in column 0.

Smooth and overlay independent runs

average is the number of samples in a simple moving window. A window of 25 samples produces N - 24 smoothed values. The implementation aligns each average to the start of its window rather than its midpoint.

import numpy as np
from matmalib import Data, Plot

rng = np.random.default_rng(8)
time_ps = np.arange(1000) * 10.0
datasets = []
for shift in (0.0, 0.35):
    signal = 1.5 + shift + 0.2 * np.sin(time_ps / 1300)
    signal += rng.normal(0, 0.12, len(time_ps))
    datasets.append(Data(data=np.column_stack((time_ps, signal)),
                         time_unit="ps", path="."))

plot = Plot()
plot.set_colors(["#2563eb", "#d97706"])
plot.trajectory(datasets, var_name="RMSD (Å)", average=[25, 50],
                overlap=True, alpha=[0.9, 0.8], headers=["run A", "run B"])

The thin faded lines are the raw data, the thicker lines are moving averages, and the histograms use the raw values. In this implementation the raw overlay is truncated to the smoothed line's length.

Set overlap=False to keep just the smoothed curves. When smoothing is enabled without overlap, the histogram is not populated; use hist=False in that case. Use integer window sizes from 2 through the number of rows when overlap=True: overlap expects a smoothed series and fails if no averaging was performed.

For a rectilinear trajectory, alpha must be a list with one value per dataset, such as [0.8] for one run; its default is 0.8 per run. average accepts either a single integer or one integer per run. headers needs one label per dataset, and the palette needs enough colors for every run.

The historical calc_qa option

calc_qa=True is a specialized legacy calculation, not a general bound/unbound-state estimator. It discards the first 500 values, constructs 50 histogram bins, selects the least-populated bin among positions 15–22, and divides the counts on either side. It also derives an annotation position using bin_index / 10, which is not generally the observable's bin coordinate. A zero upper-region count gives the sentinel value 1000.0.

Only use that option when this exact convention matches your analysis. For an explicit physical threshold, calculate populations from the array directly and annotate the plot yourself:

import numpy as np
from matmalib import Data, Plot

rng = np.random.default_rng(12)
values = np.r_[rng.normal(2.0, 0.15, 700), rng.normal(3.0, 0.15, 300)]
data = Data(data=np.column_stack((np.arange(len(values)), values)),
            time_unit="ns", path=".")
threshold = 2.5
below = np.count_nonzero(values < threshold)
above = np.count_nonzero(values >= threshold)
plot = Plot()
plot.trajectory(data, var_name="Distance (Å)")
plot.ax[0].axhline(threshold, color="black", linestyle="--")
print(f"Below / at-or-above {threshold} Å: {below}/{above} = {below / above:.3f}")

This reports a population ratio for the chosen sample set. An equilibrium constant requires an appropriate definition of states and sampling assumptions.

Compare multiple observables with scatter

For a single Data array, column 0 is x and every later column is a separate y series. Supply one header and one color per y column.

import numpy as np
from matmalib import Data, Plot

rng = np.random.default_rng(4)
x = np.linspace(1, 5, 30)
data = Data(data=np.column_stack((
    x,
    0.8 * x + rng.normal(0, 0.25, len(x)),
    1.1 * x + rng.normal(0, 0.25, len(x)),
)))
plot = Plot()
plot.set_config({"xlabel": "Reference value", "ylabel": "Prediction",
                 "title": "Compare two models"})
plot.set_colors(["#2563eb", "#dc2626"])
plot.scatter(data, headers=["model A", "model B"], format="o", figsize=(6, 4))
plot.ax.plot([1, 5], [1, 5], "k--", linewidth=1, label="agreement")
plot.ax.legend()

Points close to the dashed identity line agree with the reference. Change the 0.8 slope to 1.0 to remove that model's systematic slope error.

Passing a list of Data objects instead assigns labels and colors by dataset; use two-column arrays for an unambiguous legend. format is a Matplotlib marker string (for example "o" or "x"), not a line format such as "r--". The implementation computes a linear fit internally but does not display or return it; calculate a fit explicitly with NumPy if you need one.

Style axes and choose colors

Plot.Config holds configuration, and set_config() constructs it from a dictionary. These settings currently work through set_axes():

Key Meaning
xrange, yrange Two-element axis limits; not available on the Mollweide sphere
xlabel, ylabel, title Axis labels and figure-panel title
font Updates Matplotlib's global sans-serif font configuration
headers Default legend labels for trajectories and sphere titles; see scatter note below

set_config() fills missing keys from the original defaults, so each call replaces the configuration rather than merging with your previous call. It also adds those defaults to the dictionary passed in; use a fresh dictionary or a copy if you need to preserve the original.

The accepted xtick, ytick, title_fontsize, axis_fontsize, and tick_fontsize fields currently have no effect because the applying code checks different names. xextend and yextend are also unused. Set ticks and text sizes directly on the Matplotlib axes. For scatter(), pass headers explicitly; configured headers are only substituted when that argument is already non-None.

import numpy as np
import matplotlib.pyplot as plt
from matmalib import Plot

plot = Plot()
plot.set_config({"xlabel": "Angle (degrees)", "ylabel": "Cosine",
                 "title": "Apply matmalib styling to a Matplotlib axes",
                 "xrange": [0, 360], "yrange": [-1.1, 1.1]})
plot.cmap(color_num=3, map="ice", reverse=True, offset=0.05)
fig, ax = plt.subplots(figsize=(7, 3))
angles = np.linspace(0, 360, 200)
ax.plot(angles, np.cos(np.radians(angles)), color=plot.colors[0])
plot.set_axes(ax)
ax.set_xticks([0, 90, 180, 270, 360])
ax.tick_params(labelsize=10)
ax.set_title(ax.get_title(), fontsize=12)
print("Colors in the configured palette:", len(plot.colors))

cmap(map="ice") reads the external colormaps package. Its map names are separate from Matplotlib colormap names accepted by cremer_pople_sphere(cmap=...). reverse=True reverses the palette. offset subtracts that amount from each RGB channel and clips to [0, 1]; positive values darken the colors. It does not rotate or shift positions along a colormap.

Currently, cmap(color_num=n) selects n + 1 colors. Omitting color_num keeps the full palette. set_colors(["#2563eb", "#d97706"]) gives precise control of the line colors. Both methods configure the palette for a future Cremer–Pople sphere; calling cmap() after drawing a sphere also recolors its heatmap. Use actual color lists for general plotting; set_colors(None) is not a safe universal reset for all methods.

Draw categorical comparison marks with foo_plot

Despite its historical name, foo_plot() is useful for displaying several values per category as short horizontal segments. It expects a list of text lines, not a numeric CSV matrix: the first line contains category names, and each later line begins with the series identifier followed by one number per category.

from matmalib import Data, Plot

data = Data(path=".", data=[
    "Glucose Mannose Galactose",
    "SCR1 1.2 2.1 1.8",
    "SCR2 2.0 1.6 2.7",
    "SCR3 0.9 1.4 2.2",
])
plot = Plot()
plot.set_colors(["#2563eb", "#d97706", "#059669"])
plot.set_config({"ylabel": "Response (arbitrary units)",
                 "title": "Responses across sugars"})
plot.foo_plot(data, SCR=["SCR1", "SCR3"], w=0.55)

SCR=None draws every row; an explicit list selects rows and controls their color order. w=0.5 is the default segment width in category coordinates. The plot does not add a series legend, and its default y axis starts at zero. Set an explicit yrange when negative values matter. To load this text layout from a file, use Data(data=path.read_text().splitlines(), path=str(path.parent)); the numeric CSV reader is intended for a different layout.

Compare conformer energy ladders

space() calls Space.diff(quant=...) for each conformer collection. Each collection therefore has its own minimum set to zero. These heights compare the spread within collections, not absolute energies between different collections.

import matplotlib.pyplot as plt
from matmalib import Plot, Space, State

space_a = Space(states=[State(E=e) for e in (-100.000, -99.998, -99.994)])
space_b = Space(states=[State(E=e) for e in (-120.000, -119.999, -119.997)])
plot = Plot()
plot.space([space_a, space_b],
           names=[["ensemble A"] * 3, ["ensemble B"] * 3], quant="E")
fig, ax = plt.gcf(), plt.gca()
ax.set_ylabel("Relative electronic energy (kcal/mol)")
ax.set_xlabel("Ensemble")
ax.set_title("Each ensemble is referenced to its own minimum")
print("Ensemble A relative energies:", space_a.diff())

Use quant="E", "F", or "H" for electronic energies, free energies, or enthalpies; every state must have the selected attribute. names contains one list per space, with one category label per state; repeat a label to put the whole collection in one column, or omit names for numeric collection labels. Values on states are assumed to be in hartree, and Space.diff() multiplies differences by 627.5095 to report kcal/mol. space() currently labels its axis Ha, so the example corrects that label. It does not retain fig and ax on the Plot; obtain them from Matplotlib as shown, or save the fig directly.

For structural similarity rather than energy alone, follow conformer clustering and dendrograms. To draw carbohydrate symbols, use Plot.snfg() and the SNFG tutorial.

Save a finished figure

The most explicit export is Matplotlib's fig.savefig(). Plot.savefig(filename="fig", format="png") is a convenience wrapper that saves plot.fig at 300 dpi with a tight bounding box under plot.path.

from pathlib import Path
import numpy as np
from matmalib import Data, Plot

plot = Plot()
plot.scatter(Data(data=np.array([[0, 0], [1, 1], [2, 4], [3, 9]])),
             headers=["squared"], format="o")
plot.path = "."
plot.savefig(filename="squared", format="png")
print("Saved PNG bytes:", Path("squared.png").stat().st_size)

In the browser this file is written to the Python runtime's temporary filesystem and appears as a Download squared.png link below the output. Download it before leaving the page. On your computer it is written to the specified directory.

Some readers set data.path to the input file, while the convenience writer expects a directory. Set plot.path = "figures" after plotting, create that directory first, or call plot.fig.savefig("figures/plot.svg", bbox_inches="tight") directly. scatter(), snfg(), and cluster_dendrogram() do not initialize an output path, so set it yourself before using the convenience writer.

Choose the next plot