Metadata-Version: 2.4
Name: matmalib
Version: 0.3.7
Summary: Add your description here
Requires-Python: >=3.12
Description-Content-Type: text/markdown
Requires-Dist: numpy
Requires-Dist: scipy
Requires-Dist: matplotlib
Requires-Dist: colormaps
Requires-Dist: networkx
Requires-Dist: py3dmol
Dynamic: requires-dist

# 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](https://pratiman-91.github.io/colormaps/docs/Sequential)

## Dependencies
You need the following python packages to use this package:  

Python >= 3.12  

Numpy  
Matplotlib  
Colormap <br>
Networkx <br>
py3Dmol <br>

In your conda environment, run: <br>
```
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](https://github.com/EChung0111/CarPpy/tree/4589b228cf8772b0f2bd6ccae8bb2f41ad1e7d70).
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.

```python
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)`:

```python
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:

```python
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 conformer counts 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, counts in 5° × 5° bins produce
up to eight pastel green contour bands and a colorbar labelled "Number of
conformers" on the right:

```python
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)
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].
```

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. These are raw counts per
angular cell, without density smoothing or normalization by spherical area.
The contours interpolate between bin centers; regions below half a count
are left white, and isolated one-conformer bins are still drawn. `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.

`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:

```python
# 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](https://doi.org/10.1021/ja00839a011), 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](tests/data/README.md).

## SNFG figures and sugar stereochemistry

```python
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](https://www.ncbi.nlm.nih.gov/glycans/snfg.html),
with reducing ends on the right and alpha/beta linkage labels. They are
Matplotlib vector drawings and 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

```python
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

```python
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](examples/analyze_glycan.py).
