Analyze trajectories and conformer populations¶

cremer_pople_trajectory applies a fixed ring assignment to many geometries and returns a Data object with named numeric columns. The same interface works for molecular dynamics frames and unordered conformer ensembles. Only a genuine sampled trajectory has a physical time axis or supports a statistical interpretation of populations.

To generate a trajectory, start with running MLMD simulations. That workflow covers potentials, minimization, equilibration, production, and checkpoint restart; return here to analyze saved carbohydrate geometries.

1. Read an XYZ ensemble and classify every ring¶

The fixture below contains 38 canonical glucose conformers, one per reference pucker. It is an ensemble, not an MD trajectory. The numerical time column is only a placeholder for this exercise.

from collections import Counter
from matmalib import Data, Glycan

reference = Glycan.from_file("data/glucose.xyz")
result = Data.from_file(
    "data/glucose_puckers.xyz",
    software="cp2k",
    observable="cremer_pople",
    rings=reference.rings,
    timestep=1.0,
    classify=True,
)
print("Numeric shape:", result.data.shape)
print("Columns:", result.columns)
print("Classification shape:", result.pucker_labels.shape)
print("Pucker families:", dict(Counter(result.pucker_families[:, 0])))
for label, distance in zip(result.pucker_labels[:6, 0], result.pucker_distances[:6, 0]):
    print(label, "distance:", round(float(distance), 3), "degrees")

Expect numeric shape (38, 6) and classification shape (38, 1). The numeric columns are time, monomer_0.Q, monomer_0.theta, monomer_0.phi, monomer_0.q2, and monomer_0.q3. Classification preserves those numeric columns and attaches three additional arrays: pucker_labels, pucker_families, and pucker_distances.

This reader streams coordinates from XYZ/CP2K format. It retains the results table, not the full coordinate trajectory. software="cp2k" selects the analysis reader even for a plain .xyz file. The reader requires an explicit timestep; it does not extract frame times from XYZ comments. Its time unit is fs. Supply rings when you have a reliable reference, as above; otherwise it detects sugars on the first frame and attaches that reference as result.glycan.

With explicit rings, the standalone reader's glycan metadata is None: it has geometric assignments but did not perform chemical detection. Use the Glycan method below when you want to retain a chosen reference object.

2. Choose how to supply frames¶

Input	How it is handled
`Space`	Iterate over `space.states`. `Space.from_xyz(path)` stores every frame as a `State`.
Iterable or generator of `State` objects	Consume once; compare atom counts and element sequences across frames.
Array shaped `(n_frames, n_atoms, 3)`	Each slice is a complete molecular geometry. Supply explicit rings or elements for detection.
Iterable/generator of coordinate arrays	Consume one finite `(n_atoms, 3)` array at a time.

Use a reference glycan when all frames share its atom order:

import numpy as np
from matmalib import Glycan, Space, cremer_pople_trajectory

reference = Glycan.from_file("data/glucose.xyz")
space = Space.from_xyz("data/glucose_puckers.xyz")
result = reference.cremer_pople_trajectory(space, timestep=0.5, classify=True)
print("Retains reference:", result.glycan is reference)
print("First three times:", result.data[:3, 0])
print("First three labels:", result.pucker_labels[:3, 0].tolist())

# The standalone function accepts full coordinate arrays too.
frames = np.stack([state.xyz for state in space.states[:3]])
array_result = cremer_pople_trajectory(
    frames,
    rings=reference.rings,
    times=[2.0, 2.5, 4.0],
    time_unit="ps",
    degrees=False,
)
print("Explicit times:", array_result.data[:, 0])
print("Units:", array_result.time_unit, array_result.angle_unit)

The glycan method checks each State against the reference element sequence and requires arrays to have the reference atom count. The standalone function checks consistency among its frames. These checks cannot detect a permutation between atoms of the same element; preserving atom identity is your responsibility.

For array input without explicit rings, supply elements, one element symbol per atom. Detection then uses the first geometry:

import numpy as np
from matmalib import cremer_pople_trajectory

frames = np.load("coordinates.npy")  # (n_frames, n_atoms, 3), unwrapped
elements = np.load("elements.npy").tolist()
result = cremer_pople_trajectory(frames, elements=elements, timestep=0.5)

Parameters and result metadata¶

Parameter	Default	Meaning
`rings`	`None`	Ordered, disjoint five-/six-membered rings expressed as zero-based atom positions. Without rings, detect once on the first frame.
`timestep`	`1.0`	Interval between saved frames, finite and positive when `times` is absent.
`times`	`None`	One finite value per frame; replaces the generated time column. Available on the standalone function and `Glycan` method.
`time_unit`	`"fs"`	Label for supplied/generated times; no numeric conversion.
`degrees`	`True`	CP angular units; false produces radians.
`classify`	`False`	Attach per-frame, per-ring label/family/distance arrays. Uses classification's default tolerances.
`elements`	`None`	Element symbols for automatic detection from coordinate-array input.
`conn_mat`	`None`	Trusted reference connectivity for automatic detection.
`distXX`, `distXH`	`1.65`, `1.15`	Bond inference cutoffs in coordinate units; defaults assume angstroms.
`path`, `software`	`""`, `None`	Metadata attached to the returned `Data`.

Glycan.cremer_pople_trajectory supplies its own rings. The Data.from_file(..., software="cp2k", observable="cremer_pople") adapter exposes rings, timestep, degrees, classify, conn_mat, distXX, and distXH. For explicit nonuniform times or another time unit, use the standalone or glycan method.

Results include data, columns, rings, time_unit, angle_unit, path, software, and glycan. Columns are grouped by monomer. A five-membered ring contributes Q, phi, and q2; a six-membered ring contributes Q, theta, phi, q2, and q3. Mixed ring sizes are supported in one trajectory. Select by columns.index("monomer_0.theta") instead of assuming an offset.

times need not be evenly spaced; the implementation validates length and finiteness rather than monotonicity. Use increasing times for meaningful time-series plots. classify=True does not expose a trajectory-wide max_distance; to use custom rejection, classify individual frames with classify_ring_pucker(..., max_distance=...).

3. Follow a controlled shape change in time¶

This synthetic trajectory isolates an oscillating theta and phi, with constant amplitude. It demonstrates data flow and plotting; its coordinates are not a molecular simulation.

import numpy as np
from matmalib import cremer_pople_trajectory, Plot

time_fs = np.arange(201) * 10.0
theta = np.radians(45 + 30 * np.sin(2 * np.pi * time_fs / 1000))
phi = np.radians((time_fs / 5) % 360)
Q = 0.55
j = np.arange(6)
angle = 2 * np.pi * j / 6
frames = []
for polar, phase in zip(theta, phi):
    q2, q3 = Q * np.sin(polar), Q * np.cos(polar)
    z = np.sqrt(2 / 6) * q2 * np.cos(2 * angle + phase)
    z += q3 / np.sqrt(6) * (-1.0) ** j
    frames.append(np.column_stack((1.4 * np.cos(angle), -1.4 * np.sin(angle), z)))
result = cremer_pople_trajectory(frames, rings=[range(6)], times=time_fs)

plot = Plot()
plot.trajectory(result, col=result.columns.index("monomer_0.theta"),
                var_name="theta (degrees)", title="Synthetic ring oscillation")
print("Amplitude range:", np.round(result.data[:, 1].min(), 4),
      np.round(result.data[:, 1].max(), 4))

The plot shows theta oscillating from 15° to 75° and a companion histogram; amplitude stays 0.55. The ordinary Plot.trajectory view converts fs to ps, ps to ns, and leaves ns unchanged for display. Its col selects the measured value, while column zero supplies time. It stores the figure as plot.fig; the time-series method does not return (fig, ax).

Experiment: halve the period in the sine expression, rerun, and observe the extra oscillations. Changing saved-frame spacing changes the plotted time axis, not the geometry measured at each frame. Avoid arithmetic smoothing of phi across its 0°/360° seam; ordinary averages of 359° and 1° produce a misleading 180°.

4. Map conformer counts on the CP sphere¶

A six-membered ring has polar theta and azimuthal phi, so its conformer distribution can be drawn on a Mollweide map. Here the source is a reference library, so counts are the appropriate descriptive display.

from matmalib import Glycan, Space, Plot

reference = Glycan.from_file("data/glucose.xyz")
space = Space.from_xyz("data/glucose_puckers.xyz")
result = reference.cremer_pople_trajectory(space)
plot = Plot()
fig, ax = plot.cremer_pople_sphere(
    result,
    monomer=0,
    quantity="counts",
    bins=(24, 12),
    pucker_labels=True,
    cmap="viridis",
    title="38 reference conformers: counts, not equilibrium populations",
)
print("Conformers counted:", int(plot.conformer_counts.sum()))
print("Stored histogram shape:", plot.conformer_counts.shape)

All 38 valid conformers contribute one count. The map shows theta 0° at the top (4C1 in the standard ring order), theta 180° at the bottom (1C4), and phi decreasing from 360° on the left to 0° on the right. The 0°/360° seam wraps. Pure chairs have arbitrary phi and are counted at the poles using phi 180°; planar frames and other undefined-angle frames are omitted.

The figure labels mark ideal CP positions independently of whether classify=True was requested. Labels use O for the ring oxygen (position 6 in the classification strings), with the oxygen-first positional convention explained in Puckering. Five-membered rings have no theta and cannot use this map.

Map options¶

Option	Default	Effect
`monomer`	`0`	Zero-based six-membered ring index.
`quantity`	`"free_energy"`	Choose relative population-derived energy or `"counts"`.
`bins`	`(72, 36)`	`(phi_bins, theta_bins)`; default cells are 5° × 5°. A positive integer uses that count for both axes.
`temp`	`300`	Positive temperature in kelvin for population-to-energy conversion.
`energy_unit`	`"kcal/mol"`	`"kcal/mol"` or `"kJ/mol"`.
`num_levels`	`8`	Positive number of filled color bands.
`cmap`	Configured palette or pastel green	Matplotlib colormap name or object for this call.
`pucker_labels`	`False`	Annotate 38 ideal pucker positions.
`figsize`	`(8, 4)`	Figure size in inches.
`alpha`	`1`	Shared opacity, between zero and one.
`title`	`None`	Explicit title.
`headers`	`None`	One header per dataset, joined to form a default title.

Plot.cmap(color_num=..., map=..., reverse=..., offset=...) or Plot.set_colors(...) configures the palette. An explicit cmap on the sphere call overrides that palette for the call. Calling plot.cmap(...) after plotting updates the existing sphere and colorbar. Redisplay a notebook figure to see its changed static output.

The method returns (figure, axes) and stores them as plot.fig, plot.ax. Further results are plot.heatmap, plot.colorbar, exact integer plot.conformer_counts, plot.latitude_edges, and plot.longitude_edges. Histograms are stored in displayed (latitude, longitude) order; bin edges are radians. In free-energy mode, plot.free_energy is the masked energy array in that same order.

Pass a list such as [run_a, run_b] to pool counts. Matching time units and even a time column are unnecessary for this map; the required columns are the selected monomer's named theta/phi and a valid angle_unit. Every valid row carries equal weight, regardless of spacing or order. Opacity must be common across pooled datasets. Configured xrange/yrange are unsupported because the map has fixed limits.

Plot.trajectory(result, projection="mollweide", ...) is a shortcut. In that view col, var_name, and hist have no effect, and smoothing, overlap, and calc_qa are unsupported. Prefer cremer_pople_sphere when you want its complete set of map options.

5. Interpret a population-derived energy map¶

The default quantity="free_energy" converts each occupied bin count with

ΔG[i] = -R * T * ln(N[i] / max(N))

Here R=8.314 J/(mol K) and temp is in kelvin. Energies are converted using 1 kcal = 4184 J. The most populated bin is zero; empty bins are masked. The calculation uses counts per angular cell, without spherical-area correction or trajectory reweighting. Interpolated contour rendering smooths the appearance of occupied neighborhoods; it does not increase the resolution of stored counts or energies.

Use a thermodynamic interpretation only for an appropriately sampled ensemble. A canonical-conformer library, energy-filtered search, biased simulation, or irregularly saved trajectory is not automatically an equilibrium distribution. More angular bins can reveal detail but also produce sparse, noisy populations.

For an actual trajectory in a local workflow:

from matmalib import Data, Plot

result = Data.from_file(
    "production.xyz", software="cp2k", observable="cremer_pople",
    timestep=5.0,  # fs between saved frames, not necessarily the MD integrator step
)
plot = Plot()
fig, ax = plot.cremer_pople_sphere(result, temp=310, energy_unit="kJ/mol")
fig.savefig("puckering-populations.svg", bbox_inches="tight")

6. Export a labeled table¶

import csv
from matmalib import Data

result = Data.from_file(
    "production.xyz", software="cp2k", observable="cremer_pople",
    timestep=5.0, classify=True,
)
with open("puckering.csv", "w", newline="") as stream:
    writer = csv.writer(stream)
    label_columns = [f"monomer_{i}.pucker" for i in range(len(result.rings))]
    writer.writerow([*result.columns, *label_columns])
    for values, labels in zip(result.data, result.pucker_labels):
        writer.writerow([*values, *labels])

All frames must have the same atom count and atom order and contain whole, unwrapped rings. No periodic-boundary correction or topology update is performed. The reader reports frame and atom context for malformed XYZ input; the analysis reports frame and ring context for invalid geometry. Empty trajectories, invalid indices, changing atom counts, nonfinite coordinates, and mismatched time arrays are rejected.

Next, use conformer clustering to choose representative geometries from an ensemble while preserving original frame indices.