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

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)

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.