matmalib

About

This is a Python package designed by the Marianski Lab to easily create and design figures for publication, do data analysis from reading the output files of several commonly used QChem software such as Gaussian 16, ORCA, CP2k FHI-AIMS, and GROMMACS, create input files for these pieces of ssoftware and help automate several different computational chemistry workflows.

Documentation for the colors can be found here

Dependencies

You need the following python packages to use this package:

Python >= 3.12

Numpy
Matplotlib
Colormap
Networkx
py3Dmol

In your conda environment, run:

conda install numpy matplotlib colormap networkx py3Dmol

Install Instructions

To install the package you can run

pip install matmalib

To update the package:

pip install matmalib --upgrade

Version bumps

In GitHub Actions, select Bump version, choose Run workflow, select the target branch, and choose patch, minor, or major. The workflow uses GitHub's built-in token to update the package and bump-tool versions in pyproject.toml, create a commit and version tag, and push both atomically. For example, a patch bump from 0.2.6 creates tag 0.2.7.

Publish a GitHub release for that tag to trigger the existing PyPI publishing workflow. Branch protection rules still apply to the version commit.

Carbohydrate puckering and trajectories

For an executable walkthrough of the CarPpy-derived features and matmalib extensions, open examples/carppy_features.ipynb and run the cells in order from the project environment. It uses bundled benchmark structures to demonstrate sugar detection and numbering, stereochemistry and SNFG diagrams, all 38 canonical puckers, trajectories, torsions, RMSD clustering, energy filtering, duplicate removal, and exports. Generated figures and tables go into examples/carppy_features_output/. The sphere section uses the local Br-GalA / DCM / Ac donor-k1 / AI trajectory configured in the notebook and caches its derived puckering coordinates.

Glycan extends State, the library's molecular system class. Its monomers (also available as fragments) are Fragment objects that share the system's atoms. Sugar detection and backbone atom numbering are adapted from CarPpy. Detection covers substituted furanose and pyranose rings, with aldose and ketose numbering. It uses an existing connectivity matrix or atom graph when available; otherwise it infers bonds from coordinates.

from matmalib import Glycan, State

state = State.from_file("calculation/input.log", software="gaussian")
glycan = Glycan.from_state(state)
# Or load directly (XYZ/CP2K files use the last saved frame):
# glycan = Glycan.from_file("geometry.xyz")

puckering = glycan.cremer_pople()  # One dictionary per monomer
print(puckering[0]["Q"], puckering[0]["theta"], puckering[0]["phi"])
print(glycan.monomers[0].atom_assignment)  # {'O': ..., 'C1': ..., ..., 'C6': ...}
print(glycan.monomers[0].named_atoms["C1"].coord)

atom_assignment maps chemical names to zero-based positions in glycan.atoms, regardless of each atom's id. ring_indices records the ordered ring positions; atom_indices includes the monomer's substituents and hydrogens. named_atoms contains the corresponding Atom objects. Automatic ring order starts with the ring oxygen (O), then the anomeric carbon (C1 for aldoses, C2 for ketoses), followed by the remaining ring carbons. Monomers are ordered by their smallest ring atom index.

Glycosidic bonds separate the monomer fragments. The bridging oxygen belongs to the acceptor; for an anomeric-to-anomeric linkage, the lower monomer index is treated as donor. glycan.glycosidic_bonds records donor, acceptor, and linking atom indices. Covalently attached substituents remain with their monomer. Disconnected solvent and ions are in glycan.unassigned_atoms. For glycoconjugates, explicit monomer memberships can delimit the carbohydrate from a covalently attached protein or other aglycone.

Use the same assignments for a Space, an iterable of State objects, or a coordinate array of shape (n_frames, n_atoms, 3):

trajectory_data = glycan.cremer_pople_trajectory(frames, timestep=0.5)
print(trajectory_data.columns)
# ('time', 'monomer_0.Q', 'monomer_0.theta', 'monomer_0.phi',
#  'monomer_0.q2', 'monomer_0.q3')

Alternatively, stream an XYZ/CP2K trajectory directly through Data.from_file. This detects sugars on the first frame and reuses those assignments for every subsequent frame, without retaining the full coordinate trajectory:

from matmalib import Data, Plot

trajectory_data = Data.from_file(
    "trajectory.xyz",
    software="cp2k",
    observable="cremer_pople",
    timestep=0.5,  # fs between saved frames
)
Plot().trajectory(trajectory_data, col=2, var_name="theta (degrees)")
print(trajectory_data.glycan.monomers[0].atom_assignment)

For six-membered rings, plot relative free energy on a Mollweide projection of the Cremer–Pople sphere using Plot.cremer_pople_sphere. Like Plot.fes and Plot.contour, this uses filled color bands on a white background with a light dashed coordinate grid. By default, populations in 5° × 5° bins are converted to free energy at 300 K, with eight pastel green contour bands and a colorbar labelled "Free energy (kcal/mol)" on the right:

sphere_plot = Plot()
sphere_plot.cmap(color_num=16, map='tealgrn')  # Optional custom palette.
fig, ax = sphere_plot.cremer_pople_sphere(trajectory_data, monomer=0, temp=300)
fig.savefig("puckering-mollweide.svg", bbox_inches="tight")
# Set pucker_labels=True to add the 38 canonical pucker labels.
# Use bins=(144, 72) for a finer 2.5° × 2.5° grid (phi bins, theta bins).
# Use num_levels=6 to change the number of filled contour bands.
# Pool conformers from multiple runs with [data_a, data_b].
# Use quantity="counts" to display populations instead of free energy.
# Use energy_unit="kJ/mol" for kilojoules instead of the default kcal/mol.

Every valid frame contributes one count, regardless of time spacing or frame order. Counts are pooled when multiple datasets are supplied; a time column and matching time units are not required. Following Plot.fes, the conversion is ΔG = −RT ln(N/Nmax), using R = 8.314 J/(mol K) and temp in kelvin, then converting to the selected energy unit (1 kcal = 4184 J). Set temp to the simulation temperature. The most populated bin has zero free energy; unsampled bins are masked in the stored grid. This is free energy per angular cell, without spherical-area correction or trajectory reweighting. A finer drawing grid interpolates energies using occupied neighbors and rounds the boundary at half occupancy. Isolated sampled regions remain visible and empty bin centers stay white. This interpolation changes only the rendering; the histogram and calculated energies retain their original resolution. quantity="counts" retains the original population contours. bins accepts a positive integer for both dimensions or a (phi, theta) pair of bin counts. Plot.cmap() configures the palette using its map, reverse, offset, and color_num options; Plot.set_colors() supplies custom colors, as in the FES and contour plots. Passing cmap="plasma" (or a Matplotlib colormap object) to cremer_pople_sphere overrides the palette for that call. Calling sphere_plot.cmap(...) after plotting also recolors the existing sphere and its colorbar. Redisplay fig to see the change in a notebook's static output. The default figure size is (8, 4) inches; use figsize= to adjust it. Use a shared alpha value to change the heatmap opacity (default 1). Optional headers are joined into a default title. The filled contours and colorbar are available as plot.heatmap and plot.colorbar. Exact integer counts are available as plot.conformer_counts, in displayed (latitude, longitude) order; plot.latitude_edges and plot.longitude_edges contain the corresponding bin boundaries in radians. In free-energy mode, plot.free_energy holds the masked grid in the selected energy_unit ("kcal/mol" by default, or "kJ/mol") in the same order.

monomer is a zero-based ring index. Both degree and radian trajectory data are supported; displayed angles are degrees. The map has φ decreasing from 360° on the left to 0° on the right, θ=0° (⁴C₁) at the top, and θ=180° (¹C₄) at the bottom. Labels use O for the ring oxygen (position 6 in the classification results) and assume ring order O, C1, C2, C3, C4, C5. This is the automatic assignment order; as in classification, ketose labels describe ring positions rather than chemical carbon numbers. Labels mark ideal CP coordinates, independently of whether classify=True was requested.

Angles wrap at the 0°/360° seam. Pure chairs are counted at the poles using φ=180°, since their phase is arbitrary; planar frames and other undefined angles are omitted. Five-membered rings have no θ and cannot use this projection. Configured xrange/yrange are unsupported for the sphere. Plot.trajectory(data, projection="mollweide", ...) remains a shortcut to cremer_pople_sphere; its col, var_name, and hist options do not apply, and smoothing, overlap, and calc_qa are unsupported. The ordinary Plot.trajectory(data, ...) call keeps the time-series plot and histogram.

The standalone cremer_pople_trajectory(frames, ...) function supports the same inputs. For coordinate arrays, provide elements=["C", "O", ...] for automatic detection, or supply explicit rings. Its times argument can supply one time per frame; otherwise time starts at zero with the given timestep (default 1.0). time_unit defaults to "fs".

Distance-based detection uses CarPpy's distXX=1.65 and distXH=1.15 angstrom cutoffs. Supply a conn_mat indexed by atom-list position, or adjust the cutoffs, when the input uses different length units. A candidate sugar has one ring oxygen, otherwise carbon ring atoms, and exocyclic O/N/S substituents on at least two ring carbons. Ambiguous anomeric positions and overlapping rings raise ValueError rather than guessing an atom order.

For unusual chemistry, supply ordered rings and complete monomer memberships:

# Example positions only: replace with your system's assignments.
rings = [[0, 1, 2, 3, 4, 5]]
monomers = [list(range(24))]
glycan = Glycan.from_state(state, rings=rings, monomers=monomers)

Each ring must follow cyclic order; use the same starting atom and traversal direction in every frame. Monomers must include their ring atoms and cannot share atoms. Supplying only rings still infers monomer membership. Explicit rings can also be passed to the trajectory function or CP2K reader.

cremer_pople(coordinates) also accepts a single ordered ring or a batch with shape (..., 5, 3) or (..., 6, 3). The calculation follows the Cremer–Pople definition, with normal = R_sine × R_cosine and the negative-sine phase convention. It returns Q, q2, and phi; six-membered rings also return signed q3 and theta. Angles are degrees by default (degrees=False selects radians). Amplitudes retain the coordinate length unit. phi is in [0, 360) and theta in [0, 180]. Undefined angles are NaN: phi for planar rings and pure chairs, and theta for planar rings. A degenerate mean plane raises ValueError.

Trajectory frames must keep the same atom order and contain whole, unwrapped rings. These functions do not apply periodic-boundary corrections. The XYZ reader uses the supplied timestep rather than parsing times from comments.

The CarPpy adaptation and benchmark provenance are documented in tests/data/README.md.

SNFG figures and sugar stereochemistry

from matmalib import Glycan, Plot

glycan = Glycan.from_file("geometry.xyz")
print(glycan.identify_sugars())
# [{'name': 'Glc', 'configuration': 'D', 'anomer': 'beta',
#   'ring_type': 'pyranose'}, ...]

fig, ax = glycan.snfg(show_names=True)
fig.savefig("glycan.svg", bbox_inches="tight")
# Also available as Plot().snfg(glycan) or draw_snfg(glycan).

Figures use the official SNFG shapes and colors, with reducing ends on the right and alpha/beta linkage labels. Alpha linkages are dashed; beta linkages are solid. Automatic bond directions follow the acceptor carbon position: 2 down, 3 lower left, 4 left, 6 upper left, and 8 up. These directions are fixed even for single links and nested branches. bond_length=1.2 sets the horizontal and vertical center-to-center length in drawing units. Diagonal bonds are 20% longer to separate nearby symbols; diagonal_bond_scale=1.0 restores equal lengths. Reduce bond_length for more compact diagrams. Residue abbreviations are hidden unless show_names=True. Symbol borders and linkage lines default to 2-point strokes; customize them with symbol_linewidth= and linkage_linewidth=. Unknown or other acceptor positions and overlapping layouts require explicit positions=, which retain straight connections between the coordinates you supply. These Matplotlib vector drawings support SVG, PDF, and PNG export. Automatic naming recognizes the eight aldohexopyranose stereochemistries, C2 amino and N-acetyl variants, C6 deoxy variants, and uronic acids. It works without explicit hydrogens. Names describe base residues, not every protecting group or modification.

Unsupported or ambiguous chemistry stays Unknown (a white hexagon). For other sugars, including furanoses and ketoses, supply standard SNFG abbreviations: glycan.snfg(labels={0: "Fru"}), or one name per monomer in a list. Explicit names override automatic names; detected stereochemistry is retained. The renderer includes pentoses, sialic acids, deoxy sugars, amino sugars, and other standard symbol families. Nondefault D/L forms and furanose rings are annotated. Disconnected glycans are supported; cyclic or multiply linked topologies require a positions={monomer_index: (x, y), ...} mapping covering every monomer. ax= draws into an existing axes.

Canonical puckers and glycosidic torsions

puckers = glycan.identify_ring_puckers()
print(puckers[0]["label"], puckers[0]["family"], puckers[0]["distance"])
# e.g. 4C1, chair, distance from the ideal reference

links = glycan.glycosidic_torsions()  # phi, psi, and omega where applicable
backbone = glycan.backbone_torsions()  # hydroxymethyl omega per monomer
print(links[0]["atom_indices"])  # Exact quadruplets defining the angles

trajectory_data = glycan.cremer_pople_trajectory(frames, timestep=0.5, classify=True)
print(trajectory_data.pucker_labels)  # shape (n_frames, n_monomers)

Six-membered rings are compared with CarPpy's 38 canonical chair, boat, skew, half-chair, and envelope references. This uses three transformed dihedrals, stored as carppy_theta, independently of the Cremer–Pople coordinates also returned. distance is Euclidean distance in that three-angle space, in degrees; it is not a CP polar angle. Labels use CarPpy's ring positions 1–5 after oxygen, with 6 denoting oxygen. For ketoses, these positions differ from chemical carbon numbering. Five-membered rings use the nearest of 20 ideal envelope/twist CP phases, with positions 0–4 in input ring order and wrapped phase distance in degrees.

Planar rings are labeled planar. Passing max_distance= to identify_ring_puckers() rejects distant matches as unassigned, retaining the closest candidate in nearest_label. These labels describe proximity to reference shapes. The low-level classify_ring_pucker(coordinates) and ring_dihedrals(coordinates) functions are also public.

classify=True works with the standalone trajectory function and the CP2K Data.from_file interface. It adds pucker_labels, pucker_families, and pucker_distances arrays without changing the numeric data columns.

Torsions use heavy-atom definitions: phi is donor ring O–anomeric C–linking O–acceptor C; psi continues through the preceding acceptor backbone carbon (or ring O for an anomeric acceptor). Exocyclic omega is linking O–terminal C–ring C–ring O. Hydroxymethyl omega reverses that order. Results include atom indices and units; degrees=False returns radians. Missing or ambiguous terminal oxygens produce omega=None in backbone results. These functions do not unwrap periodic coordinates.

Conformer clustering and duplicate removal

from matmalib import Space, Plot

space = Space.from_xyz("conformers.xyz")  # All frames; glycans=True builds Glycans
clusters = space.cluster(threshold=0.5, method="average")
print(clusters.labels)                   # Zero-based cluster per input conformer
print(clusters.members)                  # Original conformer indices per cluster
print(clusters.representative_indices)
representatives = clusters.representatives  # New Space sharing selected States

fig, ax = Plot().cluster_dendrogram(clusters)
fig.savefig("clusters.svg", bbox_inches="tight")
unique = space.remove_duplicates(threshold=0.1)

Clustering uses aligned heavy-atom RMSD by default, with Kabsch rotation and translation, and single/complete/average/weighted hierarchical linkage. The threshold is in coordinate length units, normally angstroms. Supply atom_indices= to select part of a system, or heavy_only=False to include hydrogens. Atom correspondence is positional: all conformers must have the same elements and atom order. Symmetry permutations and reflections are not used. rmsd(first, second), rmsd_matrix(states), and space.rmsd_matrix() are available separately; pairwise matrices require quadratic memory in the number of conformers.

Representatives are lowest-energy members when every input State has a finite E; otherwise they are medoids (smallest total within-cluster RMSD). energy="F" selects another energy attribute; energy=None forces medoids. energy_window=5.0 excludes conformers more than 5 kcal/mol above the minimum; energy_unit describes stored energies ("hartree" by default, "kcal/mol", or "kJ/mol"). Filtering requires complete energies. Excluded inputs have label -1, while included_indices maps linkage-tree leaves back to the original input. XYZ comments are not parsed for energies; assign them on States or load electronic-structure calculations first.

Duplicate removal uses complete linkage, so every pair within a duplicate group is within the specified threshold. Neither clustering nor duplicate removal reorders or modifies the input Space or its coordinates.

For a runnable example that exports figures, pucker tables, and cluster representatives, see examples/analyze_glycan.py.