Run molecular dynamics with `mlmd`¶

mlmd takes a geometry, attaches a potential, and performs minimization, equilibration, and production dynamics through ASE. Its default potential averages four AIMNet2-2025 B97-3c/D3 models. You can also supply an ASE calculator for another potential. Trajectories stream to disk; the interface retains diagnostics and the current geometry rather than every frame.

This tutorial first runs a small CPU example with ASE's EMT calculator to teach the workflow without downloading models. It then runs a molecular example with AIMNet. The Run buttons inspect real configuration objects in the browser; simulations themselves run in local Python with the optional dependencies installed.

1. Install the simulation support¶

Use the current source checkout for these newly added APIs; the previously published PyPI release may not include them. Cloning the source repository requires GitHub access to it. From the repository directory in your Python environment, choose one installation:

# A small environment for the EMT examples and other supplied ASE calculators:
python -m pip install -e '.[md]'

# Or add AIMNet and PyTorch for machine-learning potentials:
python -m pip install -e '.[aimnet]'

# Or install all simulation extras, including solvation and PLUMED bindings:
python -m pip install -e '.[mlmd]'

The package requires Python 3.12 or newer; the repository's complete environment selects Python 3.12. With Conda, the complete environment also supplies the native PLUMED kernel and command-line tools:

conda env update -n matmalib -f environment.yml
conda activate matmalib
python -m pip check
python -c "from matmalib import mlmd; import ase; print(ase.__version__)"

md installs ASE; aimnet adds AIMNet and PyTorch. solvation adds ASE and Packmol, and metadynamics adds ASE and the PLUMED Python bindings. Those bindings alone do not supply the native PLUMED installation. Choose a PyTorch build appropriate for your CPU/GPU environment. Importing mlmd and creating its configurations do not load models or import these optional simulation packages.

2. Inspect a run before executing it¶

Configurations are dataclasses nested under mlmd. They hold parameters; validation that depends on the system or operation happens when you call minimize, equilibrate, or run_md. Constructing a config is not a simulation or a complete validity check.

from dataclasses import asdict, replace
from matmalib import mlmd

config = mlmd.MDConfig(
    steps=2000,
    timestep_fs=0.5,
    ensemble="NVT",
    temperature_K=300.0,
    seed=17,
    trajectory_interval=20,
)
print("Production duration (ps):", config.steps * config.timestep_fs / 1000)
print("Regular frame spacing (fs):", config.trajectory_interval * config.timestep_fs)
print("Thermostat:", asdict(config)["thermostat"])
short = replace(config, steps=100)
print("Original steps:", config.steps, "Short run steps:", short.steps)

This plans 1 ps of dynamics with regular saved frames 10 fs apart. The initial and final frames are also saved, so the final interval can be shorter than the requested spacing. Increasing trajectory_interval reduces disk output, not the number of force evaluations.

3. Run a complete CPU example¶

Save the following as run_emt.py in a working directory, then run python run_emt.py. It creates a small periodic copper system and performs a short minimization, NVT equilibration, and NVE production run. EMT demonstrates the API and integrators; it is not an ML potential and this copper example should not be reused with the default AIMNet models.

from pathlib import Path
from ase.build import bulk
from ase.calculators.emt import EMT
from matmalib import mlmd

run = Path("emt_run")
run.mkdir(exist_ok=True)
simulation = mlmd(calculator=EMT())
initial = bulk("Cu", cubic=True)  # four atoms, cell and PBC included
initial.positions[0, 0] += 0.03

relaxed = simulation.minimize(
    initial,
    mlmd.MinimizationConfig(
        optimizer="FIRE", fmax=0.02, steps=200,
        trajectory=run / "minimization.xyz",
        log_file=run / "minimization.log",
    ),
    output=run / "relaxed.xyz", output_format="state",
)
optimization = simulation.data

ready = simulation.equilibrate(
    relaxed,
    mlmd.EquilibrationConfig(
        steps=20, timestep_fs=0.2, temperature_K=100,
        seed=17, log_interval=5,
        trajectory=run / "equilibration.xyz",
        log_file=run / "equilibration.log",
    ),
    output=run / "equilibrated.xyz", output_format="state",
)
equilibration = simulation.data

trajectory = simulation.run_md(
    ready,
    mlmd.MDConfig(
        steps=40, timestep_fs=0.2, ensemble="NVE",
        initialize_velocities=False,
        trajectory=run / "production.xyz", trajectory_interval=10,
        log_file=run / "production.log", log_interval=10,
    ),
)
production = simulation.data
print("Minimization converged:", optimization.converged)
print("Equilibration steps:", equilibration.steps)
print("Production:", production.status, production.steps, "steps")
print("Production duration (fs):", production.simulation_time_fs)
print("Trajectory:", trajectory)
print("Checkpoint:", production.checkpoint_file)
print("Mean temperature (K):", production.averages["temperature_K"])

Expect a converged minimization, 20 equilibration steps, and a complete 40-step production run lasting 8 fs. production.xyz has five frames at steps 0, 10, 20, 30, and 40. Each phase has its own output files; production starts its own step/time count at zero. This very short run checks the workflow, not thermodynamic equilibration or sampling convergence.

The equilibration result carries masses and momenta into production. initialize_velocities=False explicitly requires them. The production temperature_K setting does not thermostat NVE or rescale existing velocities. To repeat the example, use a new output directory, remove its earlier outputs yourself, or deliberately set overwrite=True in each configuration whose files you want to replace. Existing outputs are protected by default.

Check the energy trace¶

Run this after the example, from the same directory:

from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt

records = np.loadtxt(Path("emt_run") / "production.log", ndmin=2)
time_fs = records[:, 1]
total_energy_eV = records[:, 8]
print("Total-energy range (eV):", np.ptp(total_energy_eV))
fig, ax = plt.subplots()
ax.plot(time_fs, total_energy_eV - total_energy_eV[0], marker="o")
ax.set(xlabel="Time (fs)", ylabel="Total-energy change (eV)")
fig.savefig("emt_run/energy_change.png", dpi=150, bbox_inches="tight")
plt.close(fig)

The log's header identifies its columns. Assess an NVE timestep using energy behavior over a suitably long test and compare smaller timesteps. Temperature fluctuations in a four-atom example are not representative of a bulk sample. ASE explains how timestep choice affects numerical stability in its MD guide.

4. Use the default AIMNet ensemble¶

With the aimnet extra installed, save this separate example as run_aimnet.py. It creates an isolated neutral water molecule and uses the actual default ML potential. The first energy/force evaluation can download the model weights and take longer than subsequent steps.

from pathlib import Path
from ase.build import molecule
from matmalib import mlmd

run = Path("aimnet_water_run")
run.mkdir(exist_ok=True)
simulation = mlmd(device="cpu")  # use "cuda" with a suitable GPU installation
trajectory = simulation.run_md(
    molecule("H2O"),
    mlmd.MDConfig(
        charge=0, multiplicity=1,
        steps=200, timestep_fs=0.25, temperature_K=300,
        ensemble="NVT", thermostat="langevin", seed=23,
        minimization=mlmd.MinimizationConfig(fmax=0.03, steps=1000),
        equilibration=mlmd.EquilibrationConfig(
            steps=100, timestep_fs=0.25, temperature_K=300, seed=23,
            log_file=run / "equilibration.log",
        ),
        trajectory=run / "production.xyz",
        log_file=run / "production.log",
        trajectory_interval=10, log_interval=10,
    ),
)
report = simulation.data
print("Trajectory:", trajectory)
print("Models:", report.models)
print("Preparation phases:", list(report.phases))
print("Final temperature (K):", report.final["temperature_K"])
print("Mean ensemble energy spread (eV):", report.averages["energy_std_eV"])

This requests 25 fs of equilibration and 50 fs of production. The numbers make a small initial run, not a recommended research sampling duration. Use your own XYZ/PDB path or State in place of molecule("H2O") after choosing a suitable model, charge, multiplicity, cell, and simulation protocol.

The default AIMNet2-2025 family covers H, B, C, N, O, F, Si, P, S, Cl, As, Se, Br, and I. Its training domain is molecular systems; a periodic cell does not make it a potential for metallic copper or arbitrary bulk materials. Open-shell systems need a model designed for their spin treatment. Review the AIMNet model selection guide for your chemistry before running it.

Select or supply a potential¶

Construction	Model selection
`mlmd()`	All four `aimnet2-b973c-2025-d3_0` through `_3` members; energy, forces, and requested stress are averaged.
`mlmd(model_path="aimnet2-2025")`	The same four-member family.
`mlmd(model_path="aimnet2-b973c-2025-d3_0")`	One registry member. Ensemble disagreement is unavailable for a single member.
`mlmd(model_path="/models/inference.pt")`	One existing AIMNet inference file.
`mlmd(model_path=["/models/a.pt", "/models/b.pt"])`	The mean of the supplied inference models.
`mlmd(calculator=calculator)`	Use the supplied ASE calculator instance; do not also set `model_path`.

For the default models, the resolver first checks the working directory and its checkpoints/ subdirectory for matching files, including aimnet2_2025_b973c_d3_0.pt through _3.pt. Otherwise AIMNet handles registry downloads. A missing explicit file path raises an error; it does not silently choose a different potential. Supply inference files rather than raw training checkpoints. Other supported family aliases include aimnet2, aimnet2-b973c, aimnet2-nse, aimnet2-pd, and aimnet2-rxn; matmalib family aliases expand to four members.

device and compile_model are explicit constructor arguments. Additional AIMNet calculator options go in model_kwargs. compile_model=True raises when TORCH_COMPILE_DISABLE=1 or TORCH_DISABLE_NATIVE_JIT=1 conflicts with that request; the library does not change those environment variables. An injected calculator must provide energy and forces, and stress for NPT. The interface clears input calculator attachments and uses the calculator you supplied to mlmd.

5. Choose the ensemble and coupling¶

Configuration	Integrator	Requirements and interpretation
`ensemble="NVE"`	ASE velocity Verlet	No thermostat or barostat. Temperature is only an initialization parameter when drawing velocities.
`ensemble="NVT", thermostat="langevin"`	Langevin	Default NVT choice. Friction is the reciprocal of `thermostat_time_fs`; accepts the supported constraints.
`ensemble="NVT", thermostat="berendsen"`	NVT Berendsen	Temperature coupling uses `thermostat_time_fs`.
`ensemble="NVT", thermostat="nose-hoover"`	Nosé–Hoover chains	Uses `thermostat_time_fs` and `thermostat_chain`; constraints are rejected.
`ensemble="NPT", barostat="mtk"`	Isotropic MTK with chains	Default NPT choice. Full 3D PBC, a finite-volume cell, stress, and no constraints required.
`ensemble="NPT", barostat="berendsen"`	NPT Berendsen	Same geometry/constraint requirements; additionally uses `compressibility_per_bar`.

NPT integrators include their thermostat; the NVT thermostat selector does not choose a separate NPT thermostat. pressure_bar defaults to 1 bar, thermostat_time_fs to 100 fs, and barostat_time_fs to 1000 fs in MDConfig. Chain lengths default to three. EquilibrationConfig instead defaults to a 10 fs thermostat time. Berendsen coupling suppresses equilibrium fluctuations, so its output should not be treated as exact canonical or isothermal-isobaric sampling. See ASE's ensemble documentation.

Here is an independent short NPT mechanics example using EMT and a periodic copper cell:

from ase.build import bulk
from ase.calculators.emt import EMT
from matmalib import mlmd

simulation = mlmd(calculator=EMT())
trajectory = simulation.run_md(
    bulk("Cu", cubic=True),
    mlmd.MDConfig(
        ensemble="NPT", barostat="mtk", pressure_bar=1.0,
        steps=20, timestep_fs=0.2, temperature_K=100, seed=31,
        thermostat_time_fs=100, barostat_time_fs=1000,
        trajectory="npt_run/trajectory.xyz", log_file="npt_run/md.log",
    ),
)
print("Final volume (angstrom^3):", simulation.data.final["volume_A3"])

For Berendsen NPT, supply a compressibility appropriate to your system. Its default 4.57e-5 per bar is a water-like value, not a universal material constant. Geometry minimization optimizes atomic positions at the supplied cell; it is not a cell relaxation. A common liquid workflow is minimization, NPT equilibration to adjust volume, then production in the desired ensemble. A combined run preserves the cell produced by NPT equilibration.

6. Control geometry, velocities, and constraints¶

Inputs may be a .xyz/.extxyz/.pdb/.ent path, a matmalib State (including Glycan), or ASE Atoms. A multi-frame path is read sequentially and its last frame is used. Input objects are copied. Config values for cell, pbc, charge, and multiplicity override input metadata; otherwise input values are preserved, with charge 0 and multiplicity 1 as fallbacks. Plain XYZ files do not supply a periodic cell.

Coordinates are in angstrom, energies in eV, forces in eV/angstrom, time in fs, temperature in K, and pressure in bar. State.velocities uses angstrom/fs, while State.momenta uses ASE's native momentum units. If a state has both arrays, momenta take precedence when reading it. Returned states preserve cell, PBC, charge, multiplicity, masses, and supported constraints.

from dataclasses import asdict
from matmalib import mlmd

config = mlmd.MDConfig(
    ensemble="NVT", thermostat="langevin",
    fixed_atoms=(0,),          # zero-based atom positions
    fix_com=True,
    initialize_velocities=None,
    max_velocity_angstrom_fs=0.1,
)
for key in ("ensemble", "thermostat", "fixed_atoms", "fix_com",
            "initialize_velocities", "max_velocity_angstrom_fs"):
    print(key, "=", asdict(config)[key])

The snippet inspects settings; whether atom 0 exists and leaves enough free degrees of freedom is checked against the real system when running.

Setting	Behavior
`initialize_velocities=None`	Preserve existing momenta; otherwise draw velocities at `temperature_K` using the seeded NumPy generator.
`initialize_velocities=True`	Draw new velocities, even if the input has momenta.
`initialize_velocities=False`	Require input momenta/velocities; missing values raise.
`remove_com_motion=True`	Remove center-of-mass motion only when initializing an unconstrained system with more than one atom. It does not constrain later dynamics.
`fix_com=True`	Add a persistent ASE center-of-mass constraint.
`fixed_atoms=(0, 2)`	Fix selected zero-based atom positions; existing fixed-atom constraints are retained.
`shake_hydrogens=True`	Add fixed-length constraints to nearby pairs involving H, using `shake_cutoff=1.3` angstrom and `shake_tolerance=1e-6`; also add persistent center-of-mass fixation.
`hydrogen_mass=...`	Replace all hydrogen masses with the given amu value.
`heavy_mass_multiplier=...`	Multiply non-hydrogen masses by this positive factor.

SHAKE pair selection uses the distance cutoff, not a chemical bond-order model. Inspect the geometry and choose the cutoff deliberately. NPT and NVT Nosé–Hoover reject all constraints, including inherited constraints and fix_com. Langevin and NVE support the documented fixed atoms/SHAKE workflow. A system with no remaining free degrees of freedom is rejected.

Mass modification is scaling, not mass-conserving hydrogen mass repartitioning. Existing velocities are preserved when masses change, so kinetic energies and momenta change accordingly. This alters the dynamics; validate your timestep and interpretation. A combined equilibration/production call applies mass settings once and carries velocities into production. If you run phases as separate method calls, avoid requesting the same heavy-atom multiplication again on an already scaled state.

7. Minimize and combine preparation phases¶

MinimizationConfig defaults to FIRE, fmax=0.05 eV/angstrom, at most 10,000 steps, and maxstep=0.1 angstrom. Optimizer names are matched case-insensitively to ASE's exported optimizer classes, for example BFGS, LBFGS, and BFGSLineSearch. Additional settings go in optimizer_kwargs; use the dedicated fields for maxstep, trajectory, and log paths. maxstep is passed when the chosen optimizer accepts it explicitly. FIRE enables its downhill check by default.

require_convergence=True raises if minimization ends above the requested force threshold. The failed optimization report remains in simulation.data, including converged, initial_fmax, final_fmax, and target_fmax. With require_convergence=False, the final unconverged state is returned; check its report before using it. Minimization preserves explicitly supplied momenta, so request initialize_velocities=True afterward if you intend a fresh thermal initialization.

MDConfig(minimization=..., equilibration=...), as used in the AIMNet example, runs preparation in order. Each phase must use distinct output paths. The production report keeps preparation summaries in phases, and adds total_simulation_time_fs, total_elapsed_seconds, and pooled workflow_averages/workflow_finite_samples. Individual production statistics remain under summary, averages, initial, and final. Ordinary equilibration is unbiased unless its own mode and metadynamics configuration request bias.

An EquilibrationConfig supports the MD settings and an optional minimization, but cannot contain another equilibration. It defaults to 2,000 NVT steps, no trajectory, equilibration.log, and no finite-force threshold. Both standalone equilibration and production write checkpoints even when no PLUMED calculation is involved.

8. Read output and diagnose a run¶

minimize and equilibrate return a State by default. Adding output="final.xyz" writes extended XYZ and returns its Path; adding output_format="state" writes the file and still returns the state. run_md returns a trajectory Path, and requires a trajectory filename. Its trajectory= keyword overrides config.trajectory for that call.

The extended XYZ stream contains geometry, cell, PBC, masses, charge/spin metadata, momenta, constraints, and step/time_fs. Force/energy calculator results are deliberately omitted from XYZ; use the log and report for energetics. wrap_output=True wraps a copied periodic output frame; it leaves the live simulation coordinates unchanged. Use unwrapped molecular geometries for analyses that assume molecules are intact across boundaries.

simulation.data is replaced by each operation. Save it as report = simulation.data before the next call to keep that report. MD reports expose:

Field	Meaning
`status`, `failure`	Completion status and, for errors caught during integration, exception type/message.
`steps`, `start_step`, `last_step`	Completed integration steps in this segment and its step range.
`time_fs`, `simulation_time_fs`	Absolute final step time and this segment's duration.
`trajectory`, `log_file`, `checkpoint_file`, `checkpoint_step`	Output paths and the last checkpoint written during this segment.
`averages`, `finite_samples`, `samples`	Online sample means, finite counts for each quantity, and total sampled records.
`initial`, `final`	Energies, temperature, maximum force, volume, ensemble spread, and step wall time at the endpoints.
`performance`	Wall-clock milliseconds per step and estimated ns/day for this run.
`system`, `models`, `settings`	Geometry/constraints/masses, calculator metadata, and principal run settings.

Averages sample the initial step, every log_interval, and the final step, whether or not a log file is enabled. They are sample means, not time-weighted integrals; a shorter final interval and duplicated phase boundary samples affect interpretation. Ensemble energy_std_eV is the population standard deviation across model energies, and force_std_eV_A is the maximum componentwise force standard deviation across atoms/components. An injected calculator or single model without these values yields nan in logs and None averages with zero finite counts. Disagreement is a diagnostic, not a calibrated error guarantee.

min_distance=0.5 angstrom rejects severe initial overlaps. Nonfinite energies, forces, coordinates, or velocities stop integration. Production also aborts above max_force_eV_A=150 eV/angstrom by default; None disables that threshold. max_velocity_angstrom_fs optionally aborts on excessive atomic speed; it never clips velocities. Partial trajectory/log output remains on disk, and the last successfully saved checkpoint remains usable. Validation or calculator errors before integration may not have the full completed-step report.

Continue with solvation to build a solvent box, metadynamics to add bias, FES analysis to interpret biased samples, or checkpoints and restarts to resume a run. The configuration reference lists the fields and defaults together.

Use trajectory analysis for saved geometries. Supply the interval between saved frames, rather than the MD integration timestep; use explicit frame times when intervals differ.

Run molecular dynamics with mlmd¶