Metadata-Version: 2.4
Name: matmalib
Version: 0.4.0
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
Requires-Dist: tornado
Provides-Extra: md
Requires-Dist: ase<4,>=3.27; extra == "md"
Provides-Extra: aimnet
Requires-Dist: ase<4,>=3.27; extra == "aimnet"
Requires-Dist: aimnet>=0.2.0; extra == "aimnet"
Requires-Dist: torch>=2.8; extra == "aimnet"
Provides-Extra: solvation
Requires-Dist: ase<4,>=3.27; extra == "solvation"
Requires-Dist: packmol>=21.2.1; extra == "solvation"
Provides-Extra: metadynamics
Requires-Dist: ase<4,>=3.27; extra == "metadynamics"
Requires-Dist: plumed>=2.9; extra == "metadynamics"
Provides-Extra: mlmd
Requires-Dist: ase<4,>=3.27; extra == "mlmd"
Requires-Dist: aimnet>=0.2.0; extra == "mlmd"
Requires-Dist: torch>=2.8; extra == "mlmd"
Requires-Dist: packmol>=21.2.1; extra == "mlmd"
Requires-Dist: plumed>=2.9; extra == "mlmd"
Provides-Extra: test
Requires-Dist: pytest>=8; extra == "test"

# matmalib

## Documentation and interactive tutorials

Read the [documentation and interactive tutorials](https://EChung0111.github.io/matmalib-docs/), with editable Python
examples, real browser calculations, figure and file downloads, and a complete
API reference. Tutorials cover molecular geometry, file I/O, ensembles and
energies, glycans, puckering, trajectories, conformer clustering, and plotting.

To preview it locally, install `requirements-docs.txt` alongside the package,
then run `python -m mkdocs serve`. See
[the documentation guide](docs/contributing.md) for validation and publishing.
The `Documentation` workflow validates changes on `main` and publishes the
generated site to [matmalib-docs](https://github.com/EChung0111/matmalib-docs),
the public GitHub Pages hosting repository. Documentation source stays here
alongside the library.

## 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
Python **3.12 or newer** is required. The base installation automatically installs
NumPy, SciPy, Matplotlib, `colormaps`, NetworkX, py3Dmol, and Tornado.

Simulation packages are optional dependencies declared in `pyproject.toml`:

| Extra | Additional packages | Enables |
| --- | --- | --- |
| `md` | ASE >=3.27,<4 | Minimization, equilibration, MD and restarts with a supplied ASE calculator |
| `aimnet` | ASE, AIMNet >=0.2.0, PyTorch >=2.8 | Default 2025 AIMNet ensemble or local AIMNet model files |
| `solvation` | ASE, Packmol >=21.2.1 | Solvation boxes |
| `metadynamics` | ASE, PLUMED Python bindings >=2.9 | Biased MD with a supplied calculator; native PLUMED also required |
| `mlmd` | All of the above | Complete simulation workflow |
| `test` | pytest >=8 | Test suite |

These packages are imported only when their features are used. Importing
`matmalib` or `mlmd`, creating configs, and analyzing/plotting an existing FES
do not require ASE, AIMNet, PyTorch, Packmol, or PLUMED. Reconstructing an FES
from HILLS does require the native `plumed sum_hills` command.

## Install Instructions
To install the package you can run  
```
pip install matmalib
```

To update the package:
```
pip install matmalib --upgrade
```

For a source checkout, activate the desired environment and choose extras:

```bash
conda activate matmalib
python -m pip install -e '.[mlmd]'          # full simulation support
# python -m pip install -e '.[solvation]'  # solvation only
# python -m pip install -e '.[md]'         # use your own ASE calculator
```

PLUMED's Python package does not include its native kernel or command-line
tools. For metadynamics or HILLS reconstruction, also install:

```bash
conda install -n matmalib -c conda-forge 'plumed>=2.9' 'py-plumed>=2.9' fftw
```

To install the complete workflow and test dependencies in one environment,
run `conda env update -n matmalib -f environment.yml` from this repository,
then `conda run -n matmalib python -m pip check`. The environment file selects
Python 3.12 and the `mlmd,test` extras. Choose a PyTorch build appropriate for
your machine if you need a specific GPU runtime.

## 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](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](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 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:

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

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

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

## Machine-learning simulations (`mlmd`)

Start with the [MLMD simulation tutorial](https://echung0111.github.io/matmalib-docs/learn/mlmd/), then follow the guides for [solvation](https://echung0111.github.io/matmalib-docs/learn/mlmd-solvation/), [metadynamics](https://echung0111.github.io/matmalib-docs/learn/mlmd-metadynamics/), [checkpoint restart](https://echung0111.github.io/matmalib-docs/learn/mlmd-restarts/), and [free-energy analysis](https://echung0111.github.io/matmalib-docs/learn/mlmd-analysis/). The [configuration reference](https://echung0111.github.io/matmalib-docs/reference/mlmd-config/) lists every configuration field and current default.

`mlmd` integrates the AIMNet project's MD, solvation, and metadynamics analysis
scripts. Configurations are dataclasses nested under the class. Molecular input
can be an XYZ path or a matmalib `State` (including a `Glycan`); PDB paths and ASE
`Atoms` are also accepted. Inputs are copied. Simulations return one final state or a disk
path, and never collect trajectories into a `Space` or a list in RAM.

Install in the `matmalib` conda environment from the repository directory:

```bash
conda env update -n matmalib -f environment.yml
conda run -n matmalib python -m pip check
```

The optional `mlmd` extra includes ASE, AIMNet, PyTorch, Packmol, and PLUMED
bindings. Smaller extras are listed under Dependencies above. `environment.yml` additionally
installs the native PLUMED kernel, its command-line tools, and FFTW. The
[PLUMED Python package](https://pypi.org/project/plumed/) does not itself install
the native kernel; it is needed for metadynamics and HILLS reconstruction.
[Packmol's Python wheels](https://github.com/m3g/packmol) include the executable.
No library function installs packages, re-executes Python, or changes global
compilation settings.

```python
from matmalib import mlmd, Plot

simulation = mlmd()  # ensemble of all four AIMNet2-2025 B97-3c/D3 members
# simulation = mlmd(model_path="/models/custom_inference.pt", device="cuda")
# simulation = mlmd(model_path=["/models/member0.pt", "/models/member1.pt"])

box = simulation.solvate(
    "solute.xyz", "water.xyz",
    mlmd.SolvationConfig(cell=30, density_g_ml=1.0),
)
packing = simulation.data                  # Data object with packing diagnostics
relaxed = simulation.minimize(box)
optimization = simulation.data            # convergence, steps, energy, forces
ready = simulation.equilibrate(
    relaxed,
    mlmd.EquilibrationConfig(steps=2000, ensemble="NPT", pressure_bar=1.0),
    output="equilibrated.xyz", output_format="state",
)
trajectory = simulation.run_md(
    ready,
    mlmd.MDConfig(steps=100000, ensemble="NVT", trajectory="trajectory.xyz"),
)
statistics = simulation.data              # statistics, timings, output paths
```

`model_path` accepts AIMNet inference weights, not raw training checkpoints.
Without paths, the [four 2025 registry members](https://isayevlab.github.io/aimnetcentral/models/guide/)
are used, first checking for the project's `aimnet2_2025_b973c_d3_0.pt` through
`_3.pt` in the working directory or its `checkpoints` subdirectory. Otherwise AIMNet resolves/downloads its registry
weights. Custom paths must exist and never silently fall back to another model.
Known suffixless model aliases are also accepted: for example,
`model_path="aimnet2-2025"` selects its four-member family, while
`model_path="aimnet2-b973c-2025-d3_0"` selects one member. An explicit missing
`.pt` filename raises an error. Compilation requests fail clearly if conflicting
PyTorch environment settings disable compilation.
Other potential architectures can be supplied as `mlmd(calculator=my_ase_calculator)`;
NPT requires that calculator to implement stress. `model_kwargs` passes additional
settings to `AIMNet2Calculator`; `device` and `compile_model` are explicit arguments.

### Inputs, outputs, and diagnostics

| Method | Default return | File output |
| --- | --- | --- |
| `solvate(solute, solvent, config)` | final `State` | `output="box.xyz"` |
| `minimize(state_or_xyz, config)` | final `State` | `output="relaxed.xyz"` |
| `equilibrate(state_or_xyz, config)` | final `State` | `output="equilibrated.xyz"` |
| `run_md(state_or_xyz, config)` | trajectory `Path` | `config.trajectory` or `trajectory=...` |
| `analyze_metadynamics(config)` | FES/summary `Data` | `config.output_dir`; return its `Path` with `output_format="files"` |

For the first three methods, supplying `output` returns the XYZ `Path` by default;
`output_format="state"` writes the file and returns a `State` instead. Outputs
use extended XYZ, with cell, PBC, masses, charge, multiplicity and momenta. A
plain XYZ has no periodic cell: set `cell` and `pbc` explicitly when needed.
Multi-frame XYZ input is read forward one frame at a time, selecting the last
frame. Output files are protected unless `overwrite=True` is set. An input
geometry cannot be used as its own output path.

`simulation.data` contains the latest operation's `Data` object. Save that
reference before another operation to retain earlier diagnostics. Its `summary`
is a dictionary. Minimization also exposes `converged`, `steps`,
`initial_energy_eV`, `final_energy_eV`, and `final_fmax`. Solvation exposes
`atom_count`, `molecules`, and `achieved_density_g_ml`. Dynamics exposes `status`,
`steps`, `start_step`, `last_step`, `time_fs`, `elapsed_seconds`, `averages`,
`trajectory`, `log_file`, `checkpoint_file`, and `checkpoint_step`;
combined runs keep preparation summaries in `phases`. Running averages include
step zero and samples at `log_interval`, plus the final step. Nothing from the
original scripts' status reports is printed by this API; diagnostics are returned
as data. Third-party libraries can still emit their own startup messages.

Dynamics statistics include potential, kinetic and total energies, temperature,
maximum force, ensemble energy/force disagreement, volume, and step wall time.
`finite_samples` records the count used for each average; unavailable disagreement
statistics have `None` averages. `initial`, `final`, `performance`, `system`,
`models`, and `settings` retain run details. Combined runs also expose
`workflow_averages`, `workflow_finite_samples`, `total_elapsed_seconds`, and
`total_simulation_time_fs`. Failed runs retain their completed-step diagnostics.

Coordinates are angstrom, forces eV/angstrom, energies eV, time fs, temperature K,
and pressure bar. Returned states have `cell`, `pbc`, `charge`, `multiplicity`,
atom masses, and (after dynamics) `momenta` in ASE units and `velocities` in
angstrom/fs. Existing momenta are preserved by default when starting a new run.
XYZ frames record `step` and `time_fs`. Use the checkpoint API below to retain
random state, thermostat/barostat variables, and metadynamics history when restarting.

### Simulation configuration

| Dataclass | Principal settings |
| --- | --- |
| `SystemConfig` | `cell`, `pbc`, `charge`, `multiplicity`, `min_distance`, `fixed_atoms`, `fix_com`, `shake_hydrogens`, `shake_cutoff`, `shake_tolerance` |
| `MinimizationConfig` | system settings; ASE `optimizer` name (for example `FIRE`, `BFGS`, `LBFGS`, `BFGSLineSearch`), `fmax`, `steps`, `maxstep`, `optimizer_kwargs`, `require_convergence`, optional `trajectory`/`log_file` |
| `MDConfig` | system settings; `steps`, `timestep_fs`, `temperature_K`, `ensemble`, thermostat/barostat settings, `seed`, `initialize_velocities`, mass settings, trajectory/log intervals, `checkpoint_interval`, `checkpoint_file`, `mode`, `metadynamics`, optional preparation configs, `continue_equilibration_bias` |
| `EquilibrationConfig` | MD settings; defaults to 2000 NVT steps with a 10 fs thermostat time, no trajectory, and `equilibration.log` |
| `SolvationConfig` | `cell` (cube or three orthorhombic lengths), `density_g_ml`, optional `counterion`, `counterion_gap`, `packmol`, `tolerance`, `nloop`, `timeout`, `seed`, `volume_samples`, `vdw_radii`, `work_dir`, `dry_run`, total `charge`/`multiplicity` |

Ensembles and integrators follow [ASE's MD implementations](https://docs.ase-lib.org/ase/md.html):

- `ensemble="NVE"`: velocity Verlet.
- `ensemble="NVT"`: `thermostat="langevin"` (default), `"berendsen"`, or `"nose-hoover"`.
- `ensemble="NPT"`: `barostat="mtk"` (default) or `"berendsen"`; each includes its thermostat.

`thermostat_time_fs` and `barostat_time_fs` control coupling; NPT uses
`pressure_bar`. Berendsen pressure coupling additionally uses
`compressibility_per_bar` (default is approximately liquid water; choose the
value for your system). Berendsen coupling is useful for equilibration and does
not reproduce the exact equilibrium ensemble fluctuations. `thermostat_chain`
and `barostat_chain` control chain lengths. NPT requires full 3D periodicity and
stress support. NPT and Nose–Hoover chains reject constraints; use Langevin or
Verlet for fixed atoms/SHAKE. `fixed_atoms` uses zero-based indices, while
PLUMED atom indices are one-based.

`hydrogen_mass` sets hydrogen masses; `heavy_mass_multiplier` scales heavy atom
masses. This is mass scaling, not mass-conserving repartitioning. It is applied
once across a combined equilibration/production workflow. `initialize_velocities`
can force initialization (`True`), require existing momenta (`False`), or
preserve them when available (`None`, default). `remove_com_motion` applies only
to newly initialized velocities; `fix_com=True` constrains the center of mass
throughout a run. SHAKE also enables this persistent constraint, as in the
original script. `max_force_eV_A` aborts production above 150 eV/angstrom by
default; set it to `None` to disable this guard. Equilibration defaults to no
finite-force threshold. `max_velocity_angstrom_fs` optionally aborts unstable
dynamics; it does not clip velocities. Nonfinite energies/forces also stop a
run. Partial trajectories remain readable on disk.

Supply `minimization=mlmd.MinimizationConfig(...)` and/or
`equilibration=mlmd.EquilibrationConfig(...)` in `MDConfig` to run preparation
before production. Equilibration is unbiased unless explicitly configured
otherwise. A failed minimization raises by default; setting
`require_convergence=False` returns the final state and marks it unconverged.
Trajectories include the initial frame, configured intervals, and final frame;
output wrapping affects only a copy, preserving live force/bias caches.

### Checkpoints and restarts

Every MD and equilibration run writes a checkpoint **every 1,000 steps and at
completion** by default. This applies to NVE, NVT and NPT, with or without PLUMED,
including custom `plumed_file` and `input_lines`. Saving positions and velocities
does not depend on HILLS/COLVAR files or on whether bias history can be restored
automatically. Ordinary MD checkpointing and restart do not require PLUMED to
be installed.
The latest snapshot is replaced atomically; a failed write leaves the previous
snapshot intact. Default filenames follow the trajectory, for example
`trajectory.checkpoint.npz`; equilibration without a trajectory uses the log
filename, or `equilibration.checkpoint.npz` if logging is disabled.
Set `checkpoint_file` to choose a path, `checkpoint_interval` to change the
frequency, or `checkpoint_file=None` to disable checkpointing.

```python
simulation.run_md("initial.xyz", mlmd.MDConfig(
    steps=10000, trajectory="production.xyz",
    checkpoint_file="production.checkpoint.npz", checkpoint_interval=1000,
))
checkpoint = simulation.data.checkpoint_file

# Also works in a new Python process with the same MLIP/calculator.
simulation.restart_md(
    checkpoint, steps=20000,              # additional integration steps
    trajectory="production_part2.xyz", log_file="production_part2.log",
)
print(simulation.data.last_step)          # 30000 after successful completion
```

Restart returns a new streamed XYZ segment with continuous step/time labels.
If `trajectory` is omitted, its name includes the checkpoint step. Existing
output files are protected. Checkpoints contain positions, velocities
(angstrom/fs), momenta, masses, cell, periodicity, constraints, charge/spin,
simulation settings, the NumPy RNG state, and ASE thermostat/barostat variables.
Restart retains the checkpoint interval saved in the original simulation settings.
They require the same ASE version and the same model files/calculator settings;
model weights are not embedded. The format uses JSON and NumPy arrays without
pickle. Coordinates for only the current frame are stored in memory.

For **generated** `MetadynamicsConfig` input, the checkpoint also contains
flushed HILLS/COLVAR snapshots. Restart restores these into a new `.plumed`
directory next to the new trajectory and continues the bias with the original
step count. Bias files are copied in chunks, so growing history is not loaded
into RAM. Restarted bias trajectories can differ slightly because PLUMED writes
HILLS values with finite precision. Automatic bias restart is not supported for
custom `plumed_file`/`input_lines`, whose actions may require additional state;
`restart_md` raises instead of silently discarding that state. Those runs still
produce molecular-state checkpoints on the same schedule. Ordinary MD and
generated PLUMED configurations both support automatic restart.

Solvation retains the script's total-density counting, Monte Carlo vdW-volume
estimates, fixed initial solute/counterion placement, periodic collision checks,
and packing retry. It includes one solute and optionally one counterion, then
adds whole solvent molecules. `charge=None` sums component charges;
`multiplicity` must be supplied for the combined system when it is not a singlet.
Packmol's intermediate PDB format limits a box to 99,999 atoms. Input isotope
and custom masses are retained and used for density counting. Intermediate files
live in `work_dir` when supplied, or a temporary directory otherwise. Reusing
generated files requires `overwrite=True`; unrelated files are preserved.
Reports are retained in `.data` and `State.solvation_report`.

For the original script's dry-run workflow, use
`SolvationConfig(dry_run=True, work_dir="packing_plan")`. This writes persistent
Packmol inputs and the density/count report without invoking Packmol. It returns
a `Data` object (also available as `simulation.data`), since no packed geometry
exists yet; omit `output` and `output_format` for a dry run.

### Metadynamics, analysis, and 1D FES plotting

```python
meta = mlmd.MetadynamicsConfig(
    atoms=(3, 11),                 # one-based atom indices; DISTANCE by default
    sigma=0.02,                    # nm for DISTANCE
    height_kj_mol=2.0, pace=300, biasfactor=10,
    calc_rct=True, grid_auto=True, # periodic, fixed-cell DISTANCE only
    hills_file="HILLS", colvar_file="COLVAR",
    save_input="plumed.dat", log_file="plumed.log",
)
simulation.run_md(ready, mlmd.MDConfig(
    steps=100000, mode="metadynamics", metadynamics=meta,
    ensemble="NVT", trajectory="metad.xyz", log_file="metad.log",
))

fes = simulation.analyze_metadynamics(mlmd.AnalysisConfig(
    plumed_input="plumed.dat", start_time=10.0, bins=100,
    output_dir="metad_analysis", smooth_window=11,
))
# fes is also simulation.data; fes.summary holds the script's reported results.
plot = Plot()
fig, ax = plot.fes_1d(fes, show_raw=True, show_basins=True)
fig.savefig("fes.png", dpi=180)
```

Generated CVs support `DISTANCE`, `ANGLE`, and `TORSION`; angle/torsion `sigma` is
in radians. `cv_label`, `metad_label`, and `wall_label` set action labels;
`metad_arg`, `wall_arg`, and `colvar_arg` configure their arguments and PRINT
columns. Unspecified arguments follow the selected CV and available bias fields.
PRINT lists omit disabled wall and normalized-bias fields.
`wall_type="UPPER_WALLS"` or `"LOWER_WALLS"` enables a wall with
`wall_at`, `wall_kappa`, `wall_exp`, `wall_eps`, and `wall_offset`. Alternatively,
set `plumed_file` or `input_lines` for arbitrary PLUMED actions/multiple CVs;
those inputs control their own output paths and units. `save_input` optionally
records the resolved input, which is also available in `simulation.data.plumed_input`.
Generated PLUMED uses nm, kJ/mol and ps. Bias updates occur once per integration
step, and the PLUMED virial is included for NPT. Metadynamics adds a time-dependent
bias, so an NVE integrator does not conserve unbiased energy in a biased run.

To reproduce the original script's biased equilibration followed by production
with different thermostat damping, set `continue_equilibration_bias=True` and
use the same metadynamics configuration in both phases:

```python
equilibration = mlmd.EquilibrationConfig(
    steps=2000, timestep_fs=0.5, thermostat_time_fs=0.5,
    mode="metadynamics", metadynamics=meta, log_file="biased_equilibration.log",
)
simulation.run_md(ready, mlmd.MDConfig(
    steps=100000, timestep_fs=0.5, thermostat_time_fs=100,
    mode="metadynamics", metadynamics=meta, equilibration=equilibration,
    continue_equilibration_bias=True, trajectory="production.xyz",
))
```

This retains the live PLUMED bias, HILLS/COLVAR history, velocities and continuous
step/time labels across the phase boundary. It also works with custom PLUMED
input in a single process. Both phases must use the same ensemble, timestep,
`temperature_K`, potential, cell/PBC/charge/spin overrides and metadynamics
configuration; thermostat damping can change. The RNG continues across phases:
use matching seeds or leave one seed unset. NPT's equilibrated cell is retained.
The default remains independent phases with unbiased equilibration unless that
phase explicitly requests metadynamics.

`calc_rct=True` records normalized bias for population reweighting and needs a
grid. The default `calc_rct=None, grid_auto=None` enables this automatically for
a fixed-cell periodic DISTANCE, matching the original script, or when a complete
explicit grid is supplied. Other cases require HILLS analysis or an explicit
grid for normalized COLVAR analysis.
Automatic bounds work only for a periodic DISTANCE in a fixed cell. For NPT,
use explicit `grid_min`, `grid_max`, and `grid_bins` with suitable bounds, or
set `calc_rct=False` and reconstruct from HILLS. Explicit `True` requires a valid
grid. Read `.data.summary` and assess sampling before interpreting a FES.

`AnalysisConfig` preserves the analysis script's three FES sources: normalized
COLVAR populations (`colvar`, `rbias_column` or `bias_column` plus `rct_column`),
PLUMED reconstruction (`hills`), or an existing grid (`fes_file`, `energy_column`).
`plumed_input` can infer the CV, units, temperature, output paths and walls.
For example, to reconstruct without normalized COLVAR data:

```python
fes = simulation.analyze_metadynamics(mlmd.AnalysisConfig(
    hills="HILLS", cv="cv", temperature=300, cv_unit="nm",
))
fig, ax = Plot().fes_1d(fes, smooth_window=0, energy_unit="kcal/mol")
```

Other analysis options include `grid_range`, `start_time`/`end_time`,
`additional_bias`, `upper_wall`/`keep_wall`, explicit `reactant`/`product` ranges,
`kappa`, and `infrequent` with `acceleration_column`. `colvar` can contain multiple
independent infrequent-metadynamics runs for first-passage kinetics. Raw
metadynamics bias alone is never treated as a normalized population weight.
Matching HILLS plus a configured grid can recover missing normalization as in
the original script. Multidimensional and periodic-boundary-crossing FES
analysis remains unsupported; supply a suitable 1D marginal grid instead.

FES output is a `Data` object with `x`, `free_energy` (kJ/mol),
`free_energy_smoothed`, two-column `data`, `summary`, and `warnings`.
`Plot.fes_1d` retains unsampled gaps, optionally shades inferred/explicit basins,
and marks the projected barrier. Its local cubic smoothing matches the script
and affects only the displayed curve; `smooth_window=0` disables it. Free-energy
and rate calculations always use the raw profile. `Plot.fes` remains the existing
2D population plotting method. `plot=True` with an analysis `output_dir` writes
`free_energy.png` using `Plot.fes_1d` alongside CSV/JSON files. Projected-barrier
Eyring estimates are labeled as proxies, not measured dynamical rates.

A complete example is in [examples/mlmd_workflow.py](examples/mlmd_workflow.py).
The [script feature map](docs/reference/mlmd-script-parity.md) records the
original settings, their library equivalents, and deliberate differences.
