matmalib

Documentation and interactive tutorials

Read the documentation and interactive tutorials, 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 for validation and publishing. The Documentation workflow validates changes on main and publishes the generated site to 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

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:

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:

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:

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

Machine-learning simulations (mlmd)

Start with the MLMD simulation tutorial, then follow the guides for solvation, metadynamics, checkpoint restart, and free-energy analysis. The configuration reference 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:

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 does not itself install the native kernel; it is needed for metadynamics and HILLS reconstruction. Packmol's Python wheels include the executable. No library function installs packages, re-executes Python, or changes global compilation settings.

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

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

Ensembles and integrators follow ASE's MD implementations:

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

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

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:

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:

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. The script feature map records the original settings, their library equivalents, and deliberate differences.