Introduction

HeddleMD is a molecular dynamics engine for simulating condensed-phase systems — liquids, solutions, and solvated biomolecules — on a single NVIDIA GPU. You describe a system with a force field, a starting configuration, and a thermodynamic ensemble, and the engine propagates it in time and writes a trajectory you can analyse.

Out of the box you can:

• Run NVE, NVT, and NPT dynamics, with a choice of thermostats (Nosé-Hoover chains, CSVR, Andersen, Berendsen) and barostats (stochastic cell-rescale, Monte-Carlo, Berendsen).
• Build a force field from Lennard-Jones van der Waals terms, smooth particle-mesh Ewald (SPME) electrostatics, and bonded interactions (harmonic and Morse bonds, harmonic angles, periodic torsions) with automatic 1-2/1-3 exclusions and 1-4 scaling.
• Keep water and other rigid groups rigid with SETTLE or SHAKE/RATTLE constraints, so you can take the longer timesteps those force fields assume.
• Minimize a starting structure before sampling, then chain multiple phases (equilibration followed by production) in a single run.
• Post-process trajectories in place with heddlemd analyze — the radial distribution function ships built-in.

Every quantity that affects the physics is set in a human-readable TOML file, in SI units by default, and inputs can be checked on a login node with heddlemd lint before a job reaches the queue.

Reproducibility

A distinguishing feature of HeddleMD is bit-wise reproducibility: rerun the same input on the same GPU and you get a byte-identical trajectory and log, every time. This makes debugging, regression testing, and sharing an exact result with a collaborator straightforward. The scope and limits of the guarantee are spelled out in Reproducibility.

What’s in this book

• Getting Started walks through installation and running the bundled 10,000-atom Lennard-Jones argon example.
• User Guide covers building your own simulation: the TOML config, the extended-XYZ starting-structure file, output formats, trajectory analysis, and what reproducibility does and does not guarantee.
• Reference documents the command-line interface.

This book is the user-facing guide; its Developer Guide is for contributors working on the engine itself. It frames the engine’s internal design — the GPU compute pipeline and the deterministic-reduction strategy — and, in its Extending HeddleMD sub-section, walks through adding new physics such as thermostats, integrators, barostats, and pair or bonded potentials. The exhaustive design reference lives in docs/architecture.md in the repository.

Features

This chapter is a capability map: the ensembles, force-field terms, and methods HeddleMD supports today, and — at the end — what it does not yet do. Every method below is selected and parameterised in the TOML config; see the Configuration Reference for the fields.

Ensembles and integrators

Choose an ensemble by composing an integrator with an optional thermostat and barostat in the phase’s config.

• velocity-verlet — symplectic NVE time-stepping. Add a [phase.thermostat] for NVT and/or a [phase.barostat] for NPT.
  • lossless = true — an opt-in variant whose position and velocity updates are exactly invertible, so a run can be stepped backward to its exact starting state. Incompatible with constraint groups. See Reproducibility.
• langevin-baoab — stochastic NVT via the Leimkuhler-Matthews BAOAB splitting. Thermostats through the friction term; needs no separate [phase.thermostat].
• mtk-npt — deterministic extended-system NPT (Martyna-Tobias-Klein, isotropic). Integrates its own Nosé-Hoover chains and cell degree of freedom, so it needs neither a separate thermostat nor barostat.

Initial velocities, when the starting-structure file omits them, are drawn from a Maxwell-Boltzmann distribution at the configured temperature, with centre-of-mass momentum removed and the ensemble rescaled to hit the target temperature exactly.

Thermostats

Compose any of these with velocity-verlet to sample the canonical ensemble (or to equilibrate).

• nose-hoover-chain — deterministic canonical sampling (Martyna-Klein-Tuckerman 1992). Configurable chain length, Suzuki-Yoshida sub-stepping, and RESPA sub-cycling; reports a conserved-quantity column (nhc_conserved) you can watch for drift.
• csvr — stochastic canonical-sampling velocity rescaling (Bussi-Donadio-Parrinello 2007).
• andersen — stochastic per-particle Maxwell-Boltzmann resampling (Andersen 1980), with a configurable collision rate.
• berendsen — weak-coupling temperature control. Equilibration only — it relaxes to the target temperature quickly but does not sample the canonical ensemble.

Barostats

Compose with velocity-verlet (and a thermostat) for NPT.

• c-rescale — stochastic isotropic cell-rescaling (Bernetti-Bussi 2020). Samples the NPT distribution exactly.
• monte-carlo — Metropolis Monte-Carlo barostat: periodic trial volume moves that rigidly scale molecular centres of mass and accept/reject against the NPT weight. Needs no instantaneous virial, and because it scales whole molecules it composes cleanly with rigid-water constraints.
• berendsen — weak-coupling pressure control. Equilibration only.

Energy minimization

• steepest-descent — adaptive-step steepest descent: each iteration takes a trial step down the force, accepting and growing the step on an energy decrease or rejecting and shrinking it on an increase. Convergence criteria (max force per atom, relative energy change, iteration cap) are all configurable. Declared as a [[minimization]] phase that interleaves freely with MD [[phase]] blocks; velocities and the box pass through untouched, so a subsequent MD phase starts from a relaxed configuration. Writes a per-iteration log (energy, max force, step size, accept/reject).
• Constraint-aware. When the topology declares rigid groups, every trial configuration is projected back onto the constraint manifold before its energy is evaluated.

Non-bonded interactions

Van der Waals

• Lennard-Jones — per type-pair (σ, ε), with a switching function that ramps the interaction smoothly to zero over [r_switch, cutoff] (continuous in value and first derivative, so energy and forces have no discontinuity at the cutoff), or a hard cutoff when r_switch = cutoff.

Electrostatics

• Smooth particle-mesh Ewald (SPME) — the full long-range treatment: a real-space screened-Coulomb term plus a reciprocal-space mesh contribution. Per-particle charges come from [[particle_types]].charge. You control the Ewald splitting parameter α, the real-space cutoff, the FFT grid, and the B-spline interpolation order (4–8). Requires the cell-list neighbor mode.

  A short determinism self-test runs at startup when SPME is configured, so a GPU whose FFT library would break reproducibility is caught before any dynamics run rather than after.

Bonded interactions

Bonded terms are declared in an optional topology file and parameterised by named type entries in the config.

• Bonds — harmonic (U = ½ k (r − r₀)², the usual stiff-spring bond) and Morse (U = D_e (1 − e^{−a(r−r_e)})²). A system may mix both; each bond is routed to its potential by type.
• Angles — harmonic (U = ½ k_θ (θ − θ₀)²).
• Dihedrals — periodic torsions (U = k_φ (1 + cos(n·φ − φ₀)), multiplicity n ∈ [1, 6]). Express a multi-term torsion by naming the same quadruple once per Fourier term.
• Exclusions — 1-2 and 1-3 neighbours are excluded automatically from the non-bonded interaction, with explicit per-pair overrides available. Periodic dihedrals additionally introduce scaled 1-4 LJ and Coulomb interactions with per-type scale factors.

Rigid constraints

Constraint groups let you hold selected internal coordinates rigid — the standard way to run rigid-water models and to take the longer timesteps that come with removing the fastest bond vibrations.

• settle — analytic (non-iterative) constraint for a symmetric three-atom rigid water molecule (Miyamoto & Kollman 1992). The fast path for the common rigid three-site water models (such as SPC/E and TIP3P): it holds the three intramolecular distances rigid with a closed-form solve, no iteration, no tolerance.
• shake — iterative SHAKE (position) paired with RATTLE (velocity) for arbitrary rigid clusters up to 8 atoms / 12 constraint pairs per group, with per-pair target distances on the constraint type.

Both work with the standard velocity-verlet integrator and the steepest-descent minimizer. The Langevin, MTK-NPT, and lossless velocity-Verlet integrators do not currently accept constraint groups.

Boundary conditions

• Orthorhombic and triclinic periodic boxes (the triclinic box takes the six lower-triangular lattice parameters lx, ly, lz, xy, xz, yz).
• Minimum-image convention for pair distances.
• Wrapped positions in the trajectory, with an optional per-particle integer image triple so you can reconstruct unwrapped coordinates exactly for diffusion and other unwrapped analyses.

Finding neighbours

Two schemes enumerate the short-range non-bonded pairs:

• cell-list (default) — spatial cell list with a skin distance, rebuilt only when an atom has drifted far enough to need it. The r_skin margin trades rebuild frequency against pairs-per-step and is configurable. Buffers size themselves to the actual pair count, so there is no per-particle neighbour cap to tune.
• all-pairs — the O(N²) direct sum, with no neighbor list. Convenient for small systems and as a reference.

Before any dynamics run, the engine checks that the box is large enough for the chosen cutoff (minimum perpendicular width ≥ 3 · (cutoff + r_skin)), so a too-small box is reported up front.

Input and output

• Inputs: a TOML config (*.in.toml), an extended-XYZ starting-structure file (*.in.xyz), and an optional topology file. Output filenames default from the config’s name, so a directory listing cleanly separates inputs (*.in.*) from generated outputs (*.out.*).
• Trajectory (*.out.<phase>.xyz) — self-describing extended-XYZ frames, each a valid starting-structure file, so any frame can be lifted out to restart a run.
• Log (*.out.<phase>.log) — a CSV of step, time, kinetic energy, and temperature, plus any conserved-quantity columns the chosen thermostat or barostat contributes. Loads directly into pandas.
• Timings (*.out.<phase>.timings) — a per-stage wall-clock performance summary.
• Overwrite protection — the runner refuses to start if any output file already exists, so a run never silently clobbers previous results.
• Units — SI by default; set units = "atomic" to work in Hartree atomic units. The choice applies uniformly to the config, the starting structure, and every output file.

Trajectory analysis

• heddlemd analyze post-processes a trajectory in place, driven by a small .in.analysis file. It is CPU-only, cheap enough to run on a login node, and its output is byte-identical across runs.
• The radial distribution function (rdf) between any two particle types ships built-in; see Analysis. Custom builds can register additional analyses alongside it.

Reproducibility

• Bit-wise reproducibility on the same GPU. Two runs of the same config produce byte-identical trajectory and log files. This covers the trajectory, the log, initial-velocity generation, and any stochastic thermostat or barostat (each seeded explicitly in the config). The full scope and its limits — cross-hardware differences, the intentionally non-deterministic timings file — are in Reproducibility.

Requirements

• NVIDIA GPU, with the CUDA Toolkit 11.8 or newer available at build time. The engine is GPU-only: there is no CPU or non-NVIDIA execution path.
• Single precision (f32) for positions, velocities, and forces throughout. This is well suited to short-range condensed-phase workloads and maximises throughput; a double-precision build path is reserved for the future.
• Built from source. See Installation.

Not yet supported

The following are outside the engine’s current scope:

• Double-precision storage and force kernels (reserved for a future build option).
• Pair and bond potentials beyond those above (e.g. Buckingham, FENE, Coulomb-Wolf).
• Ryckaert-Bellemans and improper dihedrals (only the periodic torsion form exists today).
• Anisotropic or fully flexible-cell NPT.
• Bit-wise reproducibility across different hardware — the guarantee is same-GPU only (see Reproducibility).
• Binary trajectory formats (NetCDF, HDF5) and compressed output.
• CPU or non-NVIDIA execution.

Installation

HeddleMD is built from source. There are no pre-built binaries.

Prerequisites

• NVIDIA GPU with a recent driver. The engine is CUDA-only; there is no CPU or non-NVIDIA fallback.
• CUDA Toolkit 11.8 or newer. nvcc must be on PATH so the build can compile the engine’s CUDA kernels. Verify with:

  nvcc --version
  

• Rust (Cargo edition 2024). Install via rustup.

The kernels are compiled once at cargo build time and embedded in the binary; nothing extra needs to be installed at runtime.

Build

From the repository root:

cargo build --release

This produces the binary at target/release/heddlemd. A debug build (cargo build, no flag) lives at target/debug/heddlemd and is suitable for development but several times slower per timestep.

Verify the install

Run the bundled example (described in detail in Your First Simulation):

./target/release/heddlemd run examples/lj-10000-argon/argon.in.toml

A successful run finishes in roughly a second on a recent GPU and prints a single [heddlemd] complete: ... line on stdout.

Container build

The repository ships a Containerfile for Podman/Docker that pins a known-good toolchain. Use it when you do not want to install CUDA on the host or when running the project under an AI-assistant sandbox.

Your First Simulation

The repository ships a complete 10,000-atom Lennard-Jones argon system at examples/lj-10000-argon/. It is the fastest way to confirm a working install and the easiest reference when you start building your own input decks.

What’s in the bundle

examples/lj-10000-argon/
├── argon.in.toml       # simulation config
├── argon.in.xyz        # initial particle state (10,000 atoms)
└── generate_init.py    # regenerates argon.in.xyz deterministically

The .in.{toml,xyz} suffixes follow the file-naming convention the loader uses to distinguish inputs from outputs: every input file ends in .in.<ext>, and the runner derives the output filenames from the config’s root by replacing .in. with .out..

argon.in.toml runs 100 integration steps at dt = 1 fs and a target temperature of 100 K. Atoms sit on a 20 × 20 × 25 simple-cubic lattice with 4.0 Å spacing, centred at the origin in an 8 × 8 × 10 nm box. The pair potential uses standard LJ-argon parameters (σ = 3.4 Å, ε ≈ 120 k_B) and the O(N²) all-pairs kernel.

Run it

From the repository root, with a release build available (cargo build --release):

./target/release/heddlemd run examples/lj-10000-argon/argon.in.toml

cargo run --release -- run examples/lj-10000-argon/argon.in.toml works the same way. On a recent NVIDIA GPU the run finishes in roughly a second and prints one line on stdout:

[heddlemd] complete: 100 steps in <N> ms (frames: 11, log rows: 21)

An exit code of 0 means every requested step ran and every requested output flushed. Non-zero exit codes are documented in the CLI reference.

What you should see on disk

Three files appear in examples/lj-10000-argon/ alongside the config:

• argon.out.run.xyz — 11 extended-XYZ frames: the initial state plus one frame every 10 steps. Each frame is self-describing (lattice vectors, column layout, step index, simulation time) and is itself a valid init file. Format details in Output Files and Init Files.
• argon.out.run.log — 22 lines of CSV: one header plus 21 rows logged every 5 steps. Columns are step,time,kinetic_energy,temperature:

  step,time,kinetic_energy,temperature
  0,0.000000000e0,2.070973475e-17,9.999999878e1
  5,5.000000000e-15,2.070582761e-17,9.998113260e1
  ...
  

  The realised step-0 temperature is exactly 100 K (within f32 storage round-off) because the initial velocities were sampled from a Maxwell-Boltzmann distribution and rescaled to the configured target.
• argon.out.run.timings — fixed-width text table with one row per instrumented kernel and host stage. Wall-clock measurements vary every run by design; the trajectory and log do not.

Did it work?

Two cheap sanity checks:

1. argon.out.run.log’s kinetic_energy column should hold steady to about four significant figures across all 21 rows (energy drift is small over 100 steps of NVE).
2. Re-run the command after deleting the three output files. The regenerated argon.out.run.xyz and argon.out.run.log must be byte-identical to the previous run:

   diff argon.out.run.xyz <prior copy>
   diff argon.out.run.log <prior copy>
   

   That is the load-bearing reproducibility guarantee. argon.out.run.timings will differ — that file is intentionally non-deterministic.

Re-running

The runner refuses to overwrite existing outputs. Delete (or move) argon.out.run.xyz, argon.out.run.log, and argon.out.run.timings between runs, or set explicit output.*_path fields in argon.in.toml.

Validating without running (heddlemd lint)

Before queueing a long job on shared GPU hardware, run

heddlemd lint examples/lj-10000-argon/argon.in.toml

to check the config, init file, output-path collisions, and box-vs-cutoff geometry on a login node without touching the GPU. Add --with-gpu to extend the lint through GPU initialisation and force-field allocation when you do have a GPU available. The full reference lives in the CLI Reference.

Next steps

• Writing a Simulation — build a config and init file for your own system.
• Configuration Reference — every TOML field.

Writing a Simulation

A simulation is fully specified by two files:

• a TOML config that pins everything affecting the trajectory — timestep, step count, integrator, particle types, pair potentials, output cadences;
• an extended-XYZ init file that carries the simulation box, the particles, and their starting positions (and optionally velocities).

Run them with heddlemd run <config>. The runner resolves the init path relative to the config’s directory.

This chapter walks through writing both files from scratch for a small custom system. For exhaustive field reference, see Configuration Reference and Init Files.

Worked example: 8 argon atoms on a 2×2×2 lattice

Goal: a 2 nm³ box containing 8 argon atoms on a simple-cubic lattice, run for 500 fs at 100 K with the standard Lennard-Jones argon parameters. Everything stays in SI units.

We will keep the two files side by side in a directory. The file-naming convention requires the config to end in .in.toml; we’ll use the matching .in.xyz form for the init file too so a future ls of the directory keeps inputs and outputs visually grouped:

my-run/
├── argon.in.toml
└── argon.in.xyz

Step 1: design choices to make before writing anything

Step 2: write argon.in.xyz

A 2 nm³ box centred at the origin means each position lies in [-1.0e-9, 1.0e-9) per axis. The four lattice corners at fractional coords (±0.25, ±0.25, ±0.25) give a 1 nm spacing:

8
Lattice="2.0e-9 0 0 0 2.0e-9 0 0 0 2.0e-9" Properties=species:S:1:pos:R:3
Ar -5.0e-10 -5.0e-10 -5.0e-10
Ar  5.0e-10 -5.0e-10 -5.0e-10
Ar -5.0e-10  5.0e-10 -5.0e-10
Ar  5.0e-10  5.0e-10 -5.0e-10
Ar -5.0e-10 -5.0e-10  5.0e-10
Ar  5.0e-10 -5.0e-10  5.0e-10
Ar -5.0e-10  5.0e-10  5.0e-10
Ar  5.0e-10  5.0e-10  5.0e-10

Things to remember:

• Lattice is a 9-component row-major list. For an orthorhombic box, the three off-diagonal pairs are all 0.
• Properties is a fixed enum — see Init Files for the four accepted forms. Here we omit velo:R:3 so the runner will sample velocities for us.
• Positions are in metres and must satisfy pos ∈ [-L/2, L/2) per axis. A particle at exactly +L/2 is rejected (upper bound is exclusive). The eight lattice points above sit safely inside the box.
• Particle IDs are implicit (row order). No ID column is supported.

Step 3: write argon.in.toml

schema_version = 1
init = "argon.in.xyz"

[simulation]
seed = 1            # RNG seed for the initial Maxwell-Boltzmann sampling
temperature = 100.0 # K — used because argon.in.xyz has no `velo:R:3`

[[phase]]
name = "run"
n_steps = 500
dt = 1.0e-15        # 1 fs

[phase.integrator]
kind = "velocity-verlet"
lossless = false

[phase.output]
trajectory_every = 50
log_every = 10
include_velocities = true

[[particle_types]]
name = "Ar"
mass = 6.6335e-26   # kg

[[pair_interactions]]
between = ["Ar", "Ar"]
potential = "lennard-jones"
sigma = 3.40e-10    # m
epsilon = 1.65e-21  # J  (≈ 120 K · k_B)
cutoff = 8.5e-10    # m  (2.5 σ)

Why these particular values:

• seed = 1 is sufficient — any u64 is fine, but the same value gives the same byte-identical trajectory on the same GPU.
• temperature = 100.0 is only used because argon.in.xyz lacks velocities. If you supply them, the field is still required and validated but ignored at runtime.
• cutoff = 8.5e-10 is well below L/2 = 1.0e-9, so the minimum-image convention is safe. For very small boxes like this, always sanity-check cutoff + r_skin < L/2 along the shortest box vector. The default cell-list configuration adds a 0.3·cutoff skin, which keeps you comfortably under the limit here.
• We let phase.output.trajectory_path, phase.output.log_path, and phase.output.timings_path default to argon.out.run.xyz, argon.out.run.log, argon.out.run.timings — derived from the config’s root and the phase’s name by the filename convention.

Step 4: run it

cd my-run
heddlemd run argon.in.toml

On success the runner prints one line per phase plus a final aggregate:

[heddlemd] phase `run`: 500 steps in <T> ms (frames: 11, log rows: 51)
[heddlemd] complete: 500 steps in <T> ms

and writes argon.out.run.xyz, argon.out.run.log, and argon.out.run.timings next to the config. Format details for each file live in Output Files.

Step 5: re-run reproducibility check

Move the outputs aside, re-run, and diff — argon.out.run.xyz and argon.out.run.log should match byte-for-byte. argon.out.run.timings will differ; that file is intentionally non-deterministic. See Reproducibility.

Common variations

Supply explicit velocities

Replace Properties=species:S:1:pos:R:3 with Properties=species:S:1:pos:R:3:velo:R:3 and append three velocity columns (m/s) to every data row. simulation.temperature is still required and validated, but the runner will not generate or rescale velocities.

Switch to NVT (thermostat composition)

Add a [phase.thermostat] section to the phase. For a deterministic canonical ensemble:

[phase.thermostat]
kind = "nose-hoover-chain"
temperature = 100.0
tau = 1.0e-13       # 100 fs coupling time

For a stochastic alternative:

[phase.thermostat]
kind = "csvr"
temperature = 100.0
tau = 1.0e-13
seed = 7            # independent of simulation.seed

For an integrator that owns its own thermostat use kind = "langevin-baoab" under [phase.integrator] instead, and omit [phase.thermostat] — combining the two is rejected at load time.

Reach NPT (barostat composition)

Add a [phase.barostat] section, or switch the integrator to kind = "mtk-npt" (which owns both thermostat and barostat). See Configuration Reference for the field reference.

Minimize before sampling

Add a [[minimization]] phase before the [[phase]] block. The minimizer relaxes positions along the negative energy gradient until the maximum per-atom force or the relative energy change drops below a tolerance; velocities and the box pass through untouched, so any subsequent MD phase starts from a clean low-strain configuration.

[[minimization]]
name = "min"

[minimization.algorithm]
kind = "steepest-descent"
# Optional — all of these have sensible defaults:
#   initial_step = 1.0e-12     # m
#   max_step = 1.0e-10         # m
#   step_increase = 1.2        # multiplier on accept
#   step_decrease = 0.2        # multiplier on reject
#   force_tolerance = 1.0e-10  # N
#   energy_tolerance = 1.0e-7  # relative
#   max_iterations = 1000

[[phase]]
name = "run"
n_steps = 500
dt = 1.0e-15
...

Phases run in source-document order, so the minimization above executes before the MD [[phase]]. The minimization writes <root>.out.min.minlog (per-iteration energy, max-force, step size, accept/reject flag); see Output Files.

Non-convergence at max_iterations is a hard error (exit code 2); subsequent phases do not run. If you are tuning a new system, inspect the .minlog to see whether energy is still decreasing at the cap and either raise max_iterations or relax one of the tolerances.

Equilibrate, then sample (multi-phase)

Append a second [[phase]] to the config. Particle state (positions, velocities, box) carries across the boundary; the integrator, thermostat, and barostat slots reset. Common pattern: a short NVT equilibration that suppresses trajectory output, followed by a longer NPT production phase that writes a frame per picosecond. See examples/spc-water-8192/water.in.toml for a fully worked version.

Bonded forces, exclusions, rigid constraints

Add a topology field at the top level pointing to a .topology file, then declare the matching parameter arrays:

• [[bond_types]] — harmonic (potential = "harmonic", k/r0) or Morse (potential = "morse") two-atom bonds; a system may mix both.
• [[angle_types]] — harmonic three-atom angles (k_theta/theta_0).
• [[dihedral_types]] — periodic torsions (potential = "periodic", k_phi/n/phi_0); repeat a [dihedrals] row per Fourier term for a multi-term torsion.
• [[constraint_types]] — rigid groups. Use kind = "settle" (d_OH/d_HH) for rigid three-atom water (the fast analytic path), or kind = "shake" (atoms + a constraints list) for a general rigid cluster. Constraints require the plain velocity-verlet integrator (not lossless, Langevin BAOAB, or MTK NPT).

See the Configuration Reference for the full field tables. Worked bundles:

• examples/spc-water-256/ — flexible SPC water (harmonic + Morse bonds, harmonic angles).
• examples/ethane-216/ — liquid ethane exercising the periodic dihedral force.
• examples/chignolin/ — a solvated all-atom protein combining bonds, angles, dihedrals, and SETTLE-constrained rigid water.

Restart from a previous frame

Trajectory frames are themselves valid init files. Extract a single frame into its own file, point a fresh config’s init field at it, and run. To keep the velocity field exactly, write the original trajectory with include_velocities = true and ensure the new config’s simulation.temperature is set to some value (it will be ignored because the frame supplies velocities).

Pre-flight checklist

Before launching a long run:

• Every [[particle_types]] name appears in the init file’s species column, and vice versa.
• You supplied a [[pair_interactions]] entry for every unordered pair, including same-type self-pairs. For N types that is N · (N+1) / 2 entries.
• All positions lie in [-L/2, L/2) per axis (the upper bound is exclusive).
• cutoff + r_skin is at most 1/3 of the minimum perpendicular width of the box (the runner enforces this — easier to fix in your editor than to debug at startup).
• None of the resolved output files already exist; the runner refuses to overwrite.

Configuration Reference

A simulation is specified by a TOML config file. The config pins every parameter that affects the trajectory; everything else (positions, velocities, box) lives in the init file. Run a config with:

heddlemd run path/to/argon.in.toml

This chapter is the user-facing rendering of the field reference. The canonical, exhaustive schema lives at rqm/io/config-schema.md in the repository — consult it for the deserialiser/validator error catalogue and any field this chapter elides.

Config filename convention

The config-file path passed to heddlemd run must end in .in.toml. The loader rejects any other name (InvalidConfigFilename) before it opens the file. The convention has two purposes:

• An ls of a simulation directory groups inputs (*.in.*) and outputs (*.out.*) into two visible blocks, so it is obvious which files were generated and can safely be deleted before a re-run.
• The loader derives every default output path from the <root> of the config filename without an extra config field: strip the .toml extension, then strip one trailing .in. The result is the config’s root.

Examples:

The suffix match is case-sensitive on the whole .in.toml string. argon.IN.toml is rejected. Filenames whose <root> derivation would yield the empty string (e.g. .in.toml) are likewise rejected.

The init-file path, the topology-file path, and any explicit output.*_path field are not subject to the convention — they are arbitrary user-supplied paths.

Units

The top-level units field selects the unit system used by the config file, the referenced init (.in.xyz) file, and every output file the run produces. Two values are accepted:

• units = "si" — the default (equivalent to omitting the field). SI throughout:

  Quantity	SI unit
  length	metres
  mass	kilograms
  time	seconds
  energy	joules
  charge	coulombs
  temperature	kelvin
  pressure	pascals

• units = "atomic" — Hartree atomic units (length in Bohr radii, mass in electron masses, time in atomic time units, energy in Hartrees, temperature as k_B · T in Hartrees, and so on).

The comparison is case-sensitive; any other value is rejected (UnknownUnits). The engine stores and computes in atomic units internally regardless — the loader converts every unit-bearing value at the I/O boundary, and threads the chosen system back through to every output writer so your view stays consistent end to end. All field descriptions below name their SI form for clarity; in atomic mode you write (and read) the corresponding atomic-unit quantity instead. No unit suffixes are supported in either mode.

Path resolution

All file paths in the config (init, topology, phase.<name>.output.*) are interpreted relative to the directory containing the config file, not the current working directory. Absolute paths are honored as-is. After resolution, every supplied path must be pairwise distinct.

Phases

A simulation is composed of one or more phases. Phases come in two kinds:

• MD phases declared as [[phase]] array elements. Each MD phase carries its own n_steps, dt, and slot blocks ([phase.integrator], optional [phase.thermostat], optional [phase.barostat], optional [phase.output]).
• Minimization phases declared as [[minimization]] array elements. Each minimization phase carries an algorithm selector (steepest-descent in v1), algorithm parameters, convergence criteria, and an optional [minimization.output] block.

Particle state (positions, velocities, box) carries from one phase into the next regardless of phase kind; per-phase slot state resets at every boundary. Velocities and the simulation box are unchanged by a minimization phase — only positions move.

The two arrays may be freely interleaved: phases execute in source-document order across [[phase]] and [[minimization]] combined. Typical pattern is one minimization phase followed by one or more MD phases:

[[minimization]]
name = "min"

[minimization.algorithm]
kind = "steepest-descent"

[[phase]]
name = "equil"
n_steps = 5000
dt = 1.0e-15
...

[[phase]]
name = "prod"
n_steps = 10000
dt = 1.0e-15
...

Each phase produces its own output files derived from <root> (the config file stem with one trailing .in stripped) and the phase’s name:

• MD phases write <root>.out.<phase>.xyz, <root>.out.<phase>.log, <root>.out.<phase>.timings.
• Minimization phases write <root>.out.<phase>.minlog and <root>.out.<phase>.timings (always), plus <root>.out.<phase>.xyz when trajectory_every > 0.

Per-phase paths can be overridden via the phase’s output block ([phase.output] for MD, [minimization.output] for minimization).

Phases must have distinct names across the union of both arrays. Stochastic MD slots of the same kind across phases (e.g. two CSVR thermostats) must use distinct seeds.

Worked example

A complete single-phase config for 10,000 argon atoms in NVE:

Saved as argon.in.toml:

schema_version = 1
init = "argon.in.xyz"

[simulation]
seed = 1
temperature = 100.0 # K

[[phase]]
name = "run"
n_steps = 100
dt = 1.0e-15        # 1 fs

[phase.integrator]
kind = "velocity-verlet"
lossless = false

[phase.output]
trajectory_every = 10
include_velocities = true
log_every = 5

[[particle_types]]
name = "Ar"
mass = 6.6335e-26   # kg

[[pair_interactions]]
between = ["Ar", "Ar"]
potential = "lennard-jones"
sigma = 3.40e-10    # m
epsilon = 1.65e-21  # J  (~120 K · k_B)
cutoff = 8.5e-10    # m  (2.5σ)

This is the bundled examples/lj-10000-argon/argon.in.toml file. The Configuration sections below walk through every field.

A multi-phase NPT example (equilibration then production) is shown in the bundled examples/spc-water-8192/water.in.toml.

Configuration sections

Top level

When topology is omitted no bonded forces are evaluated and the LJ and SPME real-space kernels see no exclusion list.

[simulation]

Per-step settings (timestep, step count, integrator/thermostat/barostat composition, output cadences) live in the [[phase]] array, not here.

[[phase]]

At least one [[phase]] array element is required. Each phase carries:

A phase is then specified by its [phase.integrator] (required), [phase.thermostat] (optional), [phase.barostat] (optional), and [phase.output] (optional) sub-tables. The schema of each of those slot blocks is described in the sections that follow.

[phase.integrator]

A required kind field selects the integrator; all other fields are kind-specific. Configs that pair an integrator that owns its own thermostat (langevin-baoab, mtk-npt) with a [phase.thermostat] section are rejected at load time.

kind = "velocity-verlet"

Symplectic NVE time-stepping. Compose with [phase.thermostat] for NVT.

kind = "langevin-baoab"

Stochastic NVT via the BAOAB splitting (Leimkuhler–Matthews). Owns its own thermostat; incompatible with [phase.thermostat].

kind = "mtk-npt"

Deterministic extended-system NPT (Martyna-Tobias-Klein, isotropic). Owns both its thermostat and its barostat; incompatible with both [phase.thermostat] and [phase.barostat].

[phase.thermostat] (optional)

Omit for NVE composition with the integrator. A required kind field selects the thermostat; per-kind fields are listed below. Configs that combine [phase.thermostat] with an integrator that owns its own thermostat are rejected.

kind = "nose-hoover-chain"

Deterministic NVT. Adds nhc_conserved to the log.

kind = "csvr"

Stochastic canonical sampling velocity rescaling (Bussi-Donadio-Parrinello).

kind = "andersen"

Per-particle stochastic Maxwell-Boltzmann resampling.

kind = "berendsen"

Deterministic weak-coupling thermostat. Equilibration only — does not sample the canonical ensemble.

[phase.barostat] (optional)

Omit for constant-volume composition. A required kind selects the barostat. Configs that combine [phase.barostat] with an integrator that owns its own barostat (currently mtk-npt) are rejected.

kind = "berendsen"

Deterministic isotropic weak-coupling barostat. Equilibration only.

kind = "c-rescale"

Stochastic isotropic cell-rescaling barostat (Bernetti-Bussi). Samples the canonical NPT distribution exactly.

[[minimization]]

A minimization phase advances positions along the negative gradient of the potential energy. Each entry carries:

A minimization entry has a required [minimization.algorithm] block and an optional [minimization.output] block. Unlike MD phases, minimization rejects n_steps, dt, [minimization.integrator], [minimization.thermostat], and [minimization.barostat] — those concepts do not apply.

[minimization.algorithm]

A required kind field selects the algorithm; all other fields are optional with documented defaults.

kind = "steepest-descent"

Adaptive-step steepest descent. Each iteration proposes a trial position x_trial = x + step · F / F_max, where F_max = max_i ||F_i||; if the trial energy is lower, the step is accepted and step is multiplied by step_increase (capped at max_step); otherwise positions are restored and step is multiplied by step_decrease.

The phase terminates with success when either F_max <= force_tolerance or the energy-tolerance criterion fires on an accepted iteration.

[minimization.output] (optional; all fields have defaults)

Velocities never appear in minimization trajectory frames; setting include_velocities under [minimization.output] is rejected as an unknown field. See Output Files for the .minlog format.

[[particle_types]] (array, ≥ 1 entry)

One entry per species. Names must be unique.

[[pair_interactions]] (array)

One entry per unordered pair of declared types. For N types you must supply exactly N · (N + 1) / 2 entries (including same-type self-pairs). ["A", "B"] and ["B", "A"] refer to the same pair.

Common fields:

Fields for potential = "lennard-jones":

[spme] (optional)

Activates smooth particle-mesh Ewald — the engine’s electrostatics model. Per-particle charges come from [[particle_types]].charge. Requires the cell-list neighbor mode ([spme] combined with mode = "all-pairs" is rejected).

The top-level [coulomb] truncated-Coulomb table used by earlier versions of the config has been retired. A config that declares [coulomb] is rejected at load time (CoulombRetired); use [spme] for electrostatics.

[[bond_types]], [[angle_types]], [[dihedral_types]], [[constraint_types]]

These four optional arrays declare parameter sets referenced by name from the optional .topology file’s [bonds], [angles], [dihedrals], and [constraints] sections respectively. Each is empty by default; declared-but-unused entries are permitted. Every entry has a name (unique within its array) and — for the force potentials — a potential selector; constraint types use a kind selector instead.

[[bond_types]]

Two-atom bonded potentials. A system may freely mix bond potentials; each bond is routed to the matching slot by its type.

potential = "harmonic" — U = ½ k (r − r₀)² (harmonic stiff spring):

potential = "morse" — U = D_e (1 − e^{−a(r−r_e)})²:

[[angle_types]]

Three-atom angle potentials. potential = "harmonic" is the only supported value — U = ½ k_θ (θ − θ₀)²:

[[dihedral_types]]

Four-atom torsion potentials. potential = "periodic" is the only supported value — U = k_φ · (1 + cos(n·φ − φ₀)). A multi-term torsion on one (i, j, k, l) quadruple is expressed by declaring one dihedral type per Fourier term and naming each from its own [dihedrals] row on the same quadruple.

[[constraint_types]]

Rigid-constraint groups. Each entry’s kind selects the algorithm that processes any group declared with that type in the topology’s [constraints] section. The presence of any constraint group is incompatible with the Langevin BAOAB, MTK NPT, and lossless velocity-Verlet integrators (see Validation).

kind = "settle" — analytic three-atom rigid water (fast path). A [constraints] row for this type lists three atoms in canonical order: oxygen, then the two hydrogens.

kind = "shake" — general iterative SHAKE/RATTLE for arbitrary rigid clusters. A [constraints] row lists this type’s atoms global atom indices.

The exhaustive validator error catalogue for all four arrays lives in rqm/io/config-schema.md, with the physics in the matching pages under rqm/forces/ and rqm/integration/.

[neighbor_list] (optional)

Selects the algorithm short-range pair-force slots use to enumerate non-bonded pairs. Defaults to the cell-list pipeline.

There is no max_neighbors field: the packed-neighbour pair-force architecture sizes its buffers to the actual interaction count plus a growth margin, so supplying max_neighbors is a load-time error. For mode = "all-pairs", r_skin is likewise rejected as an unknown field, and combining all-pairs with [spme] is rejected. The simulation box’s minimum perpendicular width must satisfy >= 3 · (cutoff_max + r_skin); the runner validates this once the box is known.

[phase.output] (optional; all fields have defaults)

<root> is the config root derived from the filename per the Config filename convention, and <phase> is the phase’s name. See Output Files for the full format spec.

Validation

The loader rejects configs that:

• live at a path whose final filename component does not end in .in.toml, or whose <root> derivation is empty (InvalidConfigFilename; see the Config filename convention);
• carry an unrecognised top-level or section field;
• choose an unknown kind for the integrator, thermostat, barostat, or any constraint type;
• omit a required field, or supply a non-finite or out-of-domain value (negative dt, empty type name, r_switch > cutoff, etc.);
• declare a [[pair_interactions]] entry referencing an undeclared type, omit a required pair, or duplicate a pair;
• resolve two supplied file paths to the same location;
• declare the retired [coulomb] table (CoulombRetired), or combine mode = "all-pairs" with [spme];
• combine [phase.thermostat]/[phase.barostat] with an integrator that owns its own;
• declare zero phases (the union of [[phase]] and [[minimization]] must be non-empty), duplicate a phase name across either array, or share a seed between two stochastic slots of the same kind across phases;
• choose an unknown minimization kind, or supply out-of-domain values for the SD parameters (e.g. step_decrease >= 1.0, max_step < initial_step, max_iterations = 0);
• pair the constraint framework with an integrator that does not support constraints (Langevin BAOAB, MTK NPT, and lossless velocity-Verlet currently do not).

Errors are reported on stderr as error: <message> and exit with code 1 (no integration steps run). The runner only proceeds when the config validates cleanly.

Init Files (Extended XYZ)

The init file carries everything the runner needs about the state of the system at t = 0: the particle count, the simulation box, per-particle type names, positions, and (optionally) velocities and integer image flags. Per-type properties such as mass and charge live in the TOML config, not the init file.

The format is a restricted subset of the widely-used extended-XYZ convention — restricted because the runner is strict about what it will accept.

File structure

N
key1=value1 key2="value with spaces" ...
<row 1>
<row 2>
...
<row N>

• Line 1 — a non-negative integer N, the particle count.
• Line 2 — the comment line: space-separated key=value attributes. Values may be double-quoted to embed spaces. Unknown keys are ignored.
• Lines 3..N+2 — one data row per particle.

The file is UTF-8. Lines end in \n or \r\n. Blank lines after the last data row are tolerated; non-blank trailing content is an error.

Required attributes

Lattice

A nine-element space-separated f64 list in double quotes — the three Cartesian box vectors in row-major order:

Lattice="lx 0 0 xy ly 0 xz yz lz"

The parser accepts only lower-triangular matrices: the three upper-triangular slots (positions 2, 3, 6 — a_y, a_z, b_z) must be exactly 0.0. The orthorhombic case (xy = xz = yz = 0) is the common form:

Lattice="8.0e-9 0 0 0 8.0e-9 0 0 0 1.0e-8"

The three diagonal entries (lx, ly, lz) must be finite and strictly positive; the three tilt entries (xy, xz, yz) may be any finite value. All values are in metres.

Properties

A colon-separated column-layout spec. Exactly one of these four forms is accepted, and the column order is fixed:

Properties=species:S:1:pos:R:3
Properties=species:S:1:pos:R:3:image:I:3
Properties=species:S:1:pos:R:3:velo:R:3
Properties=species:S:1:pos:R:3:velo:R:3:image:I:3

The velo and image blocks are independent — either may be present or absent — but when present each one must supply all three components in every row.

Data rows

Whitespace separates columns. Non-finite (NaN, ±Inf) values in any real column are rejected. Integer columns that fail to parse as i32 are rejected.

Positions must lie in the primary image

Every position must satisfy s_a, s_b, s_c ∈ [-1/2, 1/2) in fractional coordinates. For an orthorhombic box this is pos_x ∈ [-lx/2, lx/2) (and similarly for y, z) — the lower bound is inclusive, the upper bound exclusive. A particle sitting exactly on +L/2 is rejected. The runner does not silently wrap out-of-cell positions; if your generator can produce them, wrap or wrap-and-record the image triple before writing the file.

Particle IDs

Particle IDs are implicit: row 1 is particle 0, row 2 is particle 1, through N-1. There is no explicit ID column.

Image flags

When image:I:3 is omitted, every particle’s image triple defaults to (0, 0, 0). When present, the unwrapped position is pos + n_a · a + n_b · b + n_c · c, which for an orthorhombic box reduces to pos + image · (lx, ly, lz).

Velocities: file vs. config

• Velocities present in the file — the runner uses them verbatim, cast from f64 to f32. simulation.temperature in the config is still required and validated, but it is ignored for velocity setup.
• Velocities absent from the file — the runner samples a Maxwell-Boltzmann distribution at simulation.temperature using a ChaCha8 RNG seeded by simulation.seed, subtracts the centre-of-mass momentum so the total momentum is zero, then rescales by a single scalar so the realised flat-3N temperature equals the configured target. The procedure is fully deterministic in (seed, temperature, masses), so two runs with the same config and init file produce byte-identical starting velocities.

Round-tripping with the trajectory

Trajectory frames written by the engine are valid init files. The frame carries an extra Step= and Time= on the comment line (which the init parser ignores as unknown attributes), but the lattice, properties, and data rows match the init format exactly. To restart from frame k of a previous run, extract that frame into its own file and point a fresh config’s init field at it.

A minimal example

Two argon atoms, no velocities, in a 1 nm cubic box:

2
Lattice="1.0e-9 0 0 0 1.0e-9 0 0 0 1.0e-9" Properties=species:S:1:pos:R:3
Ar  0.0      0.0  0.0
Ar  3.4e-10  0.0  0.0

A larger example lives at examples/lj-10000-argon/argon.in.xyz; the Python script next to it regenerates it deterministically.

What’s not supported

• The crystallographic a, b, c, α, β, γ lattice syntax. Convert to the 9-component matrix yourself.
• Column types other than species:S:1, pos:R:3, velo:R:3, image:I:3 — no masses, charges, forces, or arbitrary user columns in schema v1.
• Multi-frame init files. Only the first frame is read; trailing non-blank content after the last data row is an error.
• Gzip, bzip2, xz, or binary-format variants.

Output Files

Every successful run writes files per phase alongside the config. The exact file set depends on the phase kind:

• MD phases ([[phase]]) write <root>.out.<phase>.xyz, <root>.out.<phase>.log, and <root>.out.<phase>.timings.
• Minimization phases ([[minimization]]) write <root>.out.<phase>.minlog and <root>.out.<phase>.timings (always), plus <root>.out.<phase>.xyz when trajectory_every > 0. No .log file (the .minlog replaces it).

<root> is the config’s root derived per the Config filename convention, and <phase> is the phase’s name. So a single-phase config argon.in.toml with name = "run" writes argon.out.run.xyz, argon.out.run.log, and argon.out.run.timings. A two-phase config with phases equil and prod writes six files. Paths and cadences are controlled by the [phase.output] section (or [minimization.output]) of the TOML config and can be overridden per phase.

The runner refuses to start when any of these files already exists at the resolved path. Delete or move them before re-running. The check is done up front, before the init file is read, so the runner fails fast.

Trajectory file (*.out.<phase>.xyz)

Extended-XYZ frames concatenated into a single file. Each frame is fully self-describing — particle count, box, column layout, step, and time — and matches the format the init-file parser accepts, so any single frame can be lifted out and used to restart.

Frame layout

N
Lattice="lx 0 0 xy ly 0 xz yz lz" Properties=species:S:1:pos:R:3[:velo:R:3][:image:I:3] Step=<u64> Time=<f64>
<row 1>
...
<row N>

• Lattice repeats verbatim in every frame even though the box does not change (it makes single-frame extraction easy).
• Properties is fixed at writer construction by the phase’s output.include_velocities and output.include_images, and never varies within a file.
• Step is the integration-step index, phase-local (0 for each phase’s initial frame).
• Time is step * dt in seconds, where dt is the phase’s dt.

Positions are always written wrapped into the primary image. When include_images = true, the per-particle integer image triple (images_x, images_y, images_z) is appended to each row; the unwrapped position is

pos + images_x · a + images_y · b + images_z · c

which reduces to pos + image · (lx, ly, lz) for an orthorhombic box.

Cadence

A frame is written for the initial state (Step=0) plus one every trajectory_every steps up to and including n_steps. Total frame count when trajectory_every > 0 is floor(n_steps / trajectory_every) + 1. Setting trajectory_every = 0 disables trajectory output entirely (not even the step-0 frame is written, and no file is created).

Number formatting

Floats use Rust’s {:.9e} formatter — nine fractional digits and a lower-case e exponent, which round-trips every f32 value exactly.

Log file (*.out.<phase>.log)

A plain CSV with a fixed four-column header followed by one row per log interval. No comment characters, no quoting, no trailing summary — easy to load with pandas or grep.

step,time,kinetic_energy,temperature
0,0.000000000e0,2.070973475e-17,9.999999878e1
5,5.000000000e-15,2.070582761e-17,9.998113260e1
...

Extra columns

The integrator, thermostat, and barostat each append their own diagnostic columns, in that order (integrator first, then thermostat, then barostat). The header row is always self-describing, so the exact set depends on the phase’s slot composition:

For example, a velocity-verlet + nose-hoover-chain + c-rescale phase writes:

step,time,kinetic_energy,temperature,nhc_conserved,pressure,box_volume

The *_conserved columns report the slot’s extended-system conserved quantity — a useful drift check for a well-behaved run. When every slot declares no extras (e.g. plain velocity-verlet, or langevin-baoab which owns a thermostat but adds no column), the header is exactly the four base columns.

Temperature convention

The temperature column uses a flat-3N degrees-of-freedom convention. This is exact for Langevin-thermostatted runs and for the initial sampled velocities (which the runner rescales to this convention). For an NVE run with centre-of-mass momentum removed, the per-thermal-DOF equipartition temperature is N/(N-1) times this value — a difference that vanishes for non-trivial system sizes.

Cadence

A row is written for Step=0 and then every log_every steps. Total rows when log_every > 0 is floor(n_steps / log_every) + 1 plus the one header line. Setting log_every = 0 disables the log entirely (no header, no file).

Minlog file (*.out.<phase>.minlog)

Minimization phases write a per-iteration diagnostic CSV in place of the per-step .log file. Each row records the post-iteration accepted state plus whether the trial step was accepted.

iter,energy,max_force,step,accepted
0,4.123456789e-18,2.345678901e-10,0.000000000e0,1
1,4.012345678e-18,2.234567890e-10,1.000000000e-12,1
2,4.012345678e-18,2.234567890e-10,1.200000000e-12,0
3,3.998765432e-18,2.198765432e-10,2.400000000e-13,1
...

Number formatting is identical to the .log: floats use {:.9e}, the integer columns use unpadded base-10.

Cadence

A row is written for the pre-loop initial state (iter=0) plus one every minlog_every accepted-or-rejected iterations. The final convergence iteration always appears as the last row, even when its index is not a multiple of minlog_every. Setting minlog_every = 0 disables the file (no header, no file). See [minimization.output] for the configurable fields.

Trajectory frames during minimization

When [minimization.output].trajectory_every > 0, the runner writes a <root>.out.<phase>.xyz file alongside the .minlog. Frames are written for the initial state, every trajectory_every accepted iterations, and the final convergence iteration. Velocities never appear in minimization frames (they do not change during minimization); the file’s Properties string omits velo:R:3. The Time= attribute is fixed at 0.0 (minimization has no physical time).

Timings file (*.out.<phase>.timings)

A fixed-width text table with one row per instrumented stage that collected at least one sample. The runner times every kernel launch with a pair of CUDA events on the default stream and every host stage with std::time::Instant. There is no opt-out.

stage                             count       total_ms       mean_us      min_us      max_us
vv_kick_drift                       100          0.996          10.0         6.1       111.8
neighbor_displacement_squared        100          0.518           5.2         4.1        13.3
...
total_runtime                         1        460.234      460233.7    460233.7    460233.7

Each row is exactly 92 columns wide followed by \n. A value that does not fit its nominal width is written in full rather than truncated.

Stages that recorded zero samples are omitted from the file, so the exact set of rows depends on which integrator, force slots, and neighbor-list mode were active. The fixed row ordering is documented in rqm/performance-analysis.md.

The .timings file is written only on successful exit, once, just before the runner returns; the path is reserved at startup but no partial data is left behind on failure.

Why this file is not reproducible

Wall-clock measurements vary run-to-run for reasons that have nothing to do with the simulation: GPU clocks, OS scheduling, driver state. Mixing them into the deterministic outputs would silently break reproducibility checks. They live in their own file precisely so a diff of *.out.<phase>.xyz and *.out.<phase>.log against a reference run is a clean yes/no answer. See Reproducibility for the full guarantee.

stdout

On success the runner prints one line per phase plus a final aggregate. MD phases report a step count and frame/log totals; minimization phases report an iteration count and the convergence reason:

[heddlemd] phase `min`: 87 iters in 412 ms (converged: force_tolerance, frames: 0, log rows: 88)
[heddlemd] phase `prod`: 10000 steps in 5234 ms (frames: 101, log rows: 101)
[heddlemd] complete: 2 phases, 10087 steps in 5646 ms

Convergence reasons are force_tolerance, energy_tolerance, force_zero, or max_iterations. The last of these is a hard error and exits with code 2; see the CLI reference for exit codes.

On failure the runner prints error: <message> on stderr and exits non-zero.

Reproducibility

Bit-wise reproducibility is the engine’s primary design goal. This chapter spells out exactly what that means, what it does not mean, and how to configure the lossless integrator mode that takes the guarantee one step further (bit-exact time reversal).

The guarantee

Two runs of the same config and init file, on the same GPU, produce byte-identical *.out.<phase>.xyz and *.out.<phase>.log files.

“Byte-identical” is meant literally — diff and sha256sum both say zero. The guarantee covers:

• Every floating-point value written to the trajectory and log,
• The Maxwell-Boltzmann velocity-generation step (when the init file omits velocities), and
• Any thermostat or barostat that draws stochastic noise: each one is seeded by an explicit seed field in its TOML section, and the RNGs are counter-based or deterministic.

What is not covered

*.timings

The timings file is intentionally non-deterministic. Wall-clock measurements vary every run for reasons that have nothing to do with the simulation (GPU clocks, OS scheduling, driver state). They live in their own file so that comparing trajectories against a reference is a clean yes/no operation. See Output Files.

Cross-hardware bit-exactness

The guarantee is same GPU only. A run on an A100 will not produce the same bytes as a run on an H100, and a CPU reference implementation will not match either. The reason is fused multiply-add (FMA): IEEE-754 permits implementations to fuse a*b + c into a single rounded operation, and CUDA’s nvcc is free to make that decision differently on different hardware (and at different optimisation levels). The engine ships with nvcc’s default FMA contraction enabled because the performance benefit is real; cross-hardware bit-match is not promised.

For tests, the rule of thumb is:

• GPU vs. same GPU → exact equality.
• GPU vs. CPU reference, or GPU vs. different GPU → small relative tolerance (f32 round-off scale).

Driver and CUDA-toolkit changes

Reproducibility holds for a single binary against a single GPU. Changing the CUDA toolkit, the NVIDIA driver, or the Rust compiler can perturb the generated PTX or the kernel’s compiled SASS and break the bit-exact match. A reference trajectory is only meaningful when checked back against the binary that produced it.

Why it works, in one paragraph

Floating-point addition is not associative: (a + b) + c and a + (b + c) can differ in the last bits. On a GPU, many threads contribute to each atom’s force in an order that varies run to run, so simply summing them as they arrive would give a slightly different result each time. HeddleMD avoids this by fixing the summation order — every run adds the same contributions in the same order, with no order-dependent atomic operations — and no precision is sacrificed to do it (there is no fixed-point quantization). The same discipline applies to the sums that produce log quantities such as kinetic energy. The mechanics of how this is arranged on the GPU are an implementation matter that lives in docs/architecture.md; the Developer Guide frames that design, and the determinism rules a new kernel must preserve are summarized in its Extending HeddleMD sub-section.

Fast math and the guarantee

[simulation].fast_math (default true) builds the JIT kernels with nvcc’s --use_fast_math. It does not weaken the same-GPU guarantee: the deterministic force accumulation above does not depend on the per-pair division/transcendental precision, so a fast_math run is still byte-identical to itself run-to-run. What changes is which f32 trajectory you get — a fast_math = true run and a fast_math = false run follow different (each independently reproducible) trajectories, so compare a reference trajectory only against a binary built with the same setting.

Lossless integrator mode

The velocity-Verlet integrator has an opt-in lossless mode, enabled per run via:

[phase.integrator]
kind = "velocity-verlet"
lossless = true

What it does, and what it does not do:

What it does

lossless adds a compensated-summation (Kahan-style) low-part to every particle’s position and velocity. Each accumulator becomes an (f32, f64) pair: the f32 carries the running value, the f64 carries the rounding residual. The integrator’s x += v · dt and v += a · dt/2 updates become exactly invertible — running the simulation backward from any timestep recovers the bits of the starting state.

What it does not do

lossless is not a double-precision simulation. Forces, a = F/m, and dt itself are still single-precision regardless of the flag. A run with lossless = true and a run with lossless = false follow different f32 trajectories — both are reproducible against themselves; only the former is reversible.

If you want better physical accuracy, you want f64 storage and f64 force kernels. That is not yet implemented and will be a compile-time feature flag (the f32 build should pay no abstraction cost).

Cost

Per-particle memory roughly doubles for the integrator buffers. Per-step compute is somewhat higher because each update reads and writes the low-part as well. The exact margin depends on the workload; for most systems the integrator stage is not the bottleneck, so the end-to-end slowdown is modest.

Compatibility

lossless = true is currently incompatible with constraint groups in the topology file. Configs that combine the two are rejected at load time. The constraint-supporting integrator is the standard lossless = false variant.

Verifying reproducibility yourself

Run the bundled example, save the outputs, delete and re-run:

./target/release/heddlemd run examples/lj-10000-argon/argon.in.toml
mv examples/lj-10000-argon/argon.out.run.xyz     /tmp/traj.ref
mv examples/lj-10000-argon/argon.out.run.log     /tmp/log.ref
mv examples/lj-10000-argon/argon.out.run.timings /tmp/timings.ref
./target/release/heddlemd run examples/lj-10000-argon/argon.in.toml
diff examples/lj-10000-argon/argon.out.run.xyz /tmp/traj.ref
diff examples/lj-10000-argon/argon.out.run.log /tmp/log.ref

Both diff invocations should print nothing. A diff on the timings file will, by design, show differences.

Analysis (heddlemd analyze)

heddlemd analyze runs post-processing analyses on a trajectory written by an earlier heddlemd run. The work is declared in a <root>.in.analysis TOML file alongside the simulation inputs; each analysis writes a CSV next to the input. v1 ships one built-in analysis kind — the radial distribution function (rdf) — and an open registry that lets custom builds add more without touching the framework.

The subcommand is CPU-only, runs cheaply on a login node, and its outputs are byte-identical across runs on identical inputs.

Quick start

In a directory containing the simulation outputs from a previous heddlemd run (argon.in.toml, argon.in.xyz, argon.out.run.xyz), write a minimal argon.in.analysis:

schema_version = 1

[[analyses]]
name = "ar-ar"
kind = "rdf"
between = ["Ar", "Ar"]
r_max = 3.5e-9    # m; must satisfy r_max <= min_perp_width / 2
n_bins = 200

Then run:

heddlemd analyze argon.in.analysis

A successful invocation prints

[heddlemd] analyze complete: 1 analyses over 11 frames in <T> ms

and writes argon.out.ar-ar.csv next to the analysis file.

The <root>.in.analysis file

Filename convention

The path passed to heddlemd analyze must end in .in.analysis. The derived root (filename with .in.analysis stripped) is used to default the sibling simulation config and every analysis’s output path:

A filename whose <root> would be empty (e.g. .in.analysis alone) is rejected.

Implicit pairing

When the analysis file does not set simulation or trajectory explicitly, the runner pairs with the sibling <root>.in.toml (must exist) and selects a phase from that config. The phase field on the analysis file picks the phase by name; when omitted, the last phase in the config is used. The chosen phase’s resolved output.trajectory_path (which itself defaults to <root>.out.<phase>.xyz per the Configuration Reference) is then used as the trajectory.

You can override either default explicitly:

schema_version = 1

# Analyse the equilibration phase rather than the default last phase.
phase = "equil"

# Or, analyse a trajectory from a different simulation directory.
simulation = "../other-run/other.in.toml"
trajectory = "../other-run/other.out.prod.xyz"

[[analyses]]
name = "ar-ar"
kind = "rdf"
between = ["Ar", "Ar"]
r_max = 3.0e-9
n_bins = 100

Frame selection

Three optional top-level fields select which trajectory frames each analysis consumes:

Frame positions count frames in the file, not the Step= value on the comment line. last_frame >= file_frame_count is rejected at trajectory-open time with FrameOutOfRange; last_frame < first_frame is rejected at load time.

[[analyses]] array

Each entry is a TOML table. Required common fields:

Kind-specific fields follow on the same entry; see RDF parameters below.

Output naming

Default output is <root>.out.<name>.csv next to the analysis file. The framework refuses to overwrite a pre-existing file at the resolved output path — delete or move the old CSV between runs.

RDF parameters

The rdf kind computes the radial distribution function g_AB(r) between two particle types.

Algorithm

Per consumed frame, the analysis enumerates unordered type-pair distances in particle-index order (i < j for same-type; the Cartesian product A × B for cross-type), applies the minimum-image convention via the trajectory’s Lattice attribute, and increments the histogram bin containing each distance under r_max. After the trajectory pass it converts the integer histogram into g_AB(r) against an ideal-gas reference using exact shell volumes and the constant box volume from the first frame.

Output CSV

r,g_r,count
<r_0>,<g_0>,<count_0>
<r_1>,<g_1>,<count_1>
...

Exactly n_bins data rows after the one-line header.

Linting an .in.analysis file

heddlemd lint dispatches on file extension: pointing it at a .in.analysis runs the analyze lint pipeline, which performs the same four setup-phase stages heddlemd analyze would run but stops before the trajectory pass and writes no files:

$ heddlemd lint argon.in.analysis
[heddlemd lint] OK
  config       argon.in.analysis
  output paths none pre-exist
  trajectory   resolved, 10000 particles, box 8.0e-9 × 8.0e-9 × 1.0e-8 m
  analyses     1 analysis builders validated

A failure surfaces the standard FAIL — <reason> line on the offending stage plus an error: <message> line on stderr — same shape as the simulation lint. Useful for catching geometric mistakes (r_max too large) or output-path collisions on a login node before queueing the analysis job.

Reproducibility

Two heddlemd analyze runs on the same .in.analysis, .in.toml, and trajectory produce byte-identical output CSVs. The guarantee is unconditional on hardware in v1 (CPU-only). It rests on deterministic pair enumeration order, integer histogram accumulation, and fixed numeric formatting; see rqm/analysis/rdf.md for the full argument.

Out of scope

• GPU-accelerated analyses.
• Streaming output (each analysis writes its CSV once at end of pass).
• Variable-box trajectories (the first frame’s lattice is taken as canonical).
• Cross-trajectory analysis or restart/append mode.
• Selecting frames by Step= value rather than file position.

Command-Line Interface

The heddlemd binary has three subcommands: run (executes the simulation), lint (validates inputs without running), and analyze (post-processes a trajectory written by run).

Usage

heddlemd run     <config-path>
heddlemd lint    <config-path> [--with-gpu]
heddlemd analyze <analysis-path>

Invoking heddlemd with no arguments, an unknown subcommand, or a subcommand without its required path argument prints

usage: heddlemd run     <config-path>
       heddlemd lint    <config-path> [--with-gpu]
       heddlemd analyze <analysis-path>

to stderr and exits with code 1.

heddlemd run <config-path>

Loads the TOML config at <config-path> and runs the simulation it describes to completion.

• <config-path> is the path to a TOML config file. Relative paths are resolved against the current working directory.
• All output paths in the config (init, output.trajectory_path, output.log_path, output.timings_path, optional topology) are resolved relative to the config file’s directory, not the current working directory. Absolute paths are honored as-is.

What it does, in order

1. Parses the CLI arguments.
2. Loads and validates the TOML config.
3. Checks that none of the enabled output files already exist. Failing this check up front means the runner never starts a long run that would be unable to write its results.
4. Loads the init file and, when supplied, the topology file.
5. Initialises CUDA and uploads the particle state to the GPU.
6. Generates initial velocities (only when the init file omits them).
7. For each phase in source-document order across [[phase]] and [[minimization]]: opens the phase’s output files, runs either the timestep loop ([[phase]]) for phase.n_steps iterations writing frames and log rows at the configured cadences, or the steepest-descent outer loop ([[minimization]]) writing per-iteration .minlog rows until a convergence criterion fires (see Configuration → [[minimization]]). Flushes the files, then writes the phase’s .timings file before moving on.

On success

Prints one line per phase plus a final aggregate on stdout. MD phases report a step count; minimization phases report an iteration count and the convergence reason (force_tolerance, energy_tolerance, force_zero):

[heddlemd] phase `<name>`: <N> steps in <T> ms (frames: <F>, log rows: <R>)
[heddlemd] phase `<name>`: <N> iters in <T> ms (converged: <reason>, frames: <F>, log rows: <R>)
...
[heddlemd] complete: <total_phases> phases, <total_N> steps in <total_T> ms

and exits with code 0. For very short runs (< 10 ms) the elapsed time is shown in microseconds instead (<T> µs).

On failure

Prints a single line to stderr beginning with error: , followed by a human-readable description that names the offending file or field where applicable, and exits with the appropriate non-zero code.

Exit codes

The split between exit code 1 and 2 makes it cheap for a wrapper script to distinguish input mistakes (re-edit the config) from mid-flight failures (likely transient: re-run, check GPU health, or relax minimization tolerances).

heddlemd lint <config-path> [--with-gpu]

Validates the config and its referenced inputs against every error the runner can detect without executing the integration loop, then exits without writing any output files. Designed for HPC contexts where a long submission queue makes trial-and-error iteration expensive: run heddlemd lint on a login node before queueing a job, and fix any reported issues up front instead of after the queue eventually grants GPU time.

The subcommand dispatches on the file extension:

• <path>.in.toml runs the simulation lint pipeline (described in What it checks below).

• <path>.in.analysis runs the analyze lint pipeline (described under Analysis and rqm/analysis/framework.md).

• --with-gpu (optional) is accepted only for the simulation lint pipeline; it extends the lint to include device initialisation, the cuFFT determinism smoke test (when SPME is configured), the host-to-device upload, slot construction, and the force-field allocation. Passing --with-gpu together with a .in.analysis path is rejected (analysis is CPU-only in v1).

• Lint writes no files at any time. Pre-existing output files are detected with Path::exists(); the filesystem is otherwise unchanged.

• Stops at the first failed check (short-circuit). Subsequent stages appear as skipped (earlier check failed) in the per-stage report.

What it checks

Stages run in the order below. Each stage reports Ok, FAIL, Skipped, or not checked (re-run with --with-gpu) in the per-stage report.

1. config — TOML parse, the .in.toml filename convention, every per-field domain check, path-collision check, and per-kind registry dispatch.
2. output paths — checks each enabled output path with Path::exists(). A pre-existing file is reported as a FAIL with the same OutputExists payload run would surface.
3. init — loads the extended-XYZ init file, reports the particle count and box dimensions.
4. box/cutoff — for cell-list mode, verifies min_perpendicular_width ≥ 3 · (cutoff_max + r_skin). For all-pairs mode, reported as not applicable.
5. topology — loads the topology file (when supplied) and reports the bond, angle, dihedral, and constraint-group counts. When the config omits topology, reported as not supplied.
6. gpu — only attempted with --with-gpu. Runs init_device (including the cuFFT smoke test for SPME configs), allocates the particle buffers, constructs every slot, and constructs the force field. Without --with-gpu, reported as not checked.

On success

Prints the per-stage report on stdout (header [heddlemd lint] OK) and exits with code 0:

[heddlemd lint] OK
  config       /path/to/argon.in.toml
  output paths none pre-exist
  init         resolved, 10000 particles, box 8.0e-9 × 8.0e-9 × 1.0e-8 m
  box/cutoff   min perp width 8.00e-9 m ≥ required 3.32e-9 m
  topology     not supplied
  gpu          not checked (re-run with --with-gpu)

On failure

Prints the per-stage report on stdout (header [heddlemd lint] FAIL), followed by a single error: <message> line on stderr that matches what run would print for the same condition, and exits with code 1:

[heddlemd lint] FAIL
  config       /path/to/argon.in.toml
  output paths none pre-exist
  init         resolved, 10000 particles, box 2.0e-9 × 5.0e-9 × 5.0e-9 m
  box/cutoff   FAIL — min perp width 2.00e-9 m along `a` < required 3.30e-9 m
  topology     skipped (earlier check failed)
  gpu          not checked (re-run with --with-gpu)
error: simulation box perpendicular width along lattice direction `a` is 2e-9, below the required 3.3e-9

Exit codes

heddlemd analyze <analysis-path>

Runs every analysis declared in <analysis-path> (a <root>.in.analysis file) against the trajectory it points at, and writes one CSV per analysis. See Analysis for the input-file schema, the implicit pairing rule, the trajectory selection knobs, and the built-in rdf kind.

• <analysis-path> is the path to a <root>.in.analysis file. The filename must end in .in.analysis. Relative paths are resolved against the current working directory.
• Implicit pairing: when the analysis file does not set simulation or trajectory explicitly, the runner pairs with the sibling <root>.in.toml (the same <root> as the analysis filename) and reads the trajectory from that config’s resolved output.trajectory_path.
• v1 is CPU-only; no --with-gpu flag.
• Output files default to <root>.out.<name>.csv (where <name> comes from each [[analyses]] entry’s name field). Each per-analysis output_path field overrides.
• Pre-existing output files cause a hard error (OutputExists) at pre-flight, before any frame is read.

On success

Prints one line on stdout:

[heddlemd] analyze complete: <K> analyses over <F> frames in <T> ms

where <K> is the number of analyses and <F> is the number of frames consumed after first_frame, last_frame, and stride are applied.

Exit codes

What is not provided

Schema v1 deliberately keeps the CLI minimal. There are no environment variables or alternative config sources — every parameter affecting the trajectory lives in the TOML config so that two runs of the same file produce identical bits.

The following do not exist:

• A --seed flag (set simulation.seed and any thermostat/barostat seed field in the config).
• A --steps or --dt flag (set phase.n_steps and phase.dt in the relevant [[phase]]).
• A --output-dir flag (set the paths under [phase.output]).
• A --help flag. The only help is this book and the usage line above.

Developer Guide

This part of the book is for people who want to change the HeddleMD engine. It assumes that you know the basic ideas of molecular dynamics, including force fields, integrators, thermostats, and neighbor lists. It does not assume that you have programmed in Rust or CUDA.

This page introduces the vocabulary used throughout the developer guide. You do not need to learn Rust or GPU programming in full before you begin, but you should understand the host/device split and be able to recognize the Rust constructs described below. The extension guides explain the project-specific interfaces in more detail:

• Schedules and CUDA Graphs explains how HeddleMD represents and runs a timestep.
• Testing Extension Components explains the shared physical test harnesses and how they complement focused component tests.
• Extending HeddleMD explains how to add a thermostat, integrator, barostat, or potential.
• docs/architecture.md in the repository is the detailed design reference. It is useful when changing the framework itself, but it is not a prerequisite for following the extension guides.

The simulation loop

At a high level, HeddleMD performs the same loop as other molecular-dynamics programs:

1. Use the current positions to build a neighbor list of nearby atoms.
2. Evaluate the force field and sum the force on each atom.
3. Advance the positions and velocities by one timestep.
4. Apply any thermostat, barostat, or constraint operations required by the selected method.
5. Repeat.

The principal data flow is:

positions -> neighbor lists -> force evaluation -> net forces -> integration
     ^                                                               |
     +---------------------- next timestep ---------------------------+

HeddleMD runs its numerical work on a single NVIDIA GPU. On a fixed GPU, two runs of the same input are designed to produce byte-for-byte identical output. That reproducibility requirement affects how reductions, random numbers, and timestep operations are implemented.

Phases: the stages of a run

A single HeddleMD run is divided into one or more phases, each a self-contained stage that runs the simulation loop above for a fixed number of steps. A phase is the unit the user configures: every [[phase]] (and every [[minimization]]) block in the input file is one phase, and phases execute in the order they appear. Particle state — positions, velocities, and the box — carries over from one phase to the next, so the phases of a run form a continuous trajectory.

Each phase fixes its own method for its whole duration: one integrator, at most one thermostat, at most one barostat, its own timestep, and its own output cadences. A typical run uses several phases to do different jobs with different settings — for example, an energy minimization to relax the starting geometry, then an equilibration phase under a thermostat, then a production phase (perhaps adding a barostat for constant pressure) that writes the trajectory you analyze. The Configuration Reference describes the [[phase]] schema in full.

The developer guide refers to “a phase” constantly, and the reason it is the natural unit is exactly that its method is fixed for its whole span. Anything the engine only needs to prepare once — validating the integrator’s schedule (see Validating Timestep Dependencies), recording a CUDA graph (see Schedules and CUDA Graphs) — is set up once at the start of each phase and reused for every step of that phase, because the operations that setup depends on cannot change until the next phase begins.

Host code and device code

HeddleMD uses two languages with different responsibilities:

• Rust host code lives under src/ and runs on the CPU. It reads the configuration, allocates GPU memory, arranges operations in the correct order, launches GPU work, handles errors, and writes output.
• CUDA C device code lives under kernels/ and runs on the GPU. A CUDA function launched across many GPU threads is called a kernel.

The word host therefore means the CPU and its Rust code. The word device means the GPU and its memory. A device buffer is memory allocated on the GPU; the CPU cannot use its contents directly. Moving a value from the device to the host requires an explicit copy and usually forces the CPU to wait for earlier GPU work to finish.

At build time, CUDA sources are compiled to PTX, an intermediate GPU instruction format, and embedded in the HeddleMD executable. At run time, the Rust code loads the PTX and launches its kernels.

A useful mental model for a kernel

A kernel is similar to the body of a parallel loop. HeddleMD launches many instances of that body at once, usually assigning one particle or interaction to each GPU thread. Threads are grouped into blocks, and groups of 32 threads within a block execute together as a warp. An individual thread within a warp is sometimes called a lane.

GPU operations are normally asynchronous: launching a kernel asks the GPU to do work but does not necessarily wait for it to finish. A device-to-host copy of a result often introduces a synchronization point. This distinction is important for CUDA graphs and for performance.

Device buffers and scratch space

The extension guides frequently refer to a scratch buffer. This is temporary device memory used to hold an intermediate result, such as the total kinetic energy or a velocity scale factor. Allocate scratch buffers when the component is constructed and reuse them on every timestep. Allocating device memory inside a per-step method is both expensive and incompatible with some execution optimizations.

A “length-1 device buffer” is simply a GPU allocation containing one scalar. It is used when one GPU kernel produces a value that a later kernel consumes. Keeping the scalar on the device avoids copying it to the CPU between kernels.

Reading the Rust examples

The examples in this guide use a small set of recurring Rust constructs. This section is a notation guide, not a general Rust tutorial.

Structs and methods

A struct groups named fields, much like a C struct or a simple Python data class:

#![allow(unused)]
fn main() {
struct ThermostatState {
    target_temperature: f64,
    kinetic_energy: CudaSlice<Real>,
}
}

An impl block defines methods for a type. Self means the type named after impl, and self is the current value:

#![allow(unused)]
fn main() {
impl ThermostatState {
    fn target(&self) -> f64 {
        self.target_temperature
    }
}
}

References and mutable references

Rust passes borrowed access with references:

• &value borrows read-only access.
• &mut value borrows mutable access, allowing the called function to change the value.
• &self is a method that reads its object.
• &mut self is a method that may update its object.

For example, compute_kinetic_energy_on_device(buffers, &mut scratch) may write into the existing scratch buffer. It does not allocate a new buffer or transfer ownership of the old one.

Traits and trait implementations

A trait defines behavior that a type must provide. It is similar to a C++ interface or an abstract base class. The Thermostat trait, for example, defines the operations that every thermostat implementation may perform.

#![allow(unused)]
fn main() {
impl Thermostat for MyThermostat {
    fn apply_post(&mut self, /* arguments */) -> Result<(), ThermostatError> {
        // implementation
        Ok(())
    }
}
}

The first line means “implement the Thermostat interface for MyThermostat.” Traits can provide default methods, so an implementation only needs to override behavior that differs from the default.

A type such as Box<dyn Thermostat> is a heap-owned value used through the Thermostat trait. This lets the simulation hold whichever concrete thermostat the user selected without encoding every thermostat type in the runner.

Results, options, and the question-mark operator

Rust uses Result<T, E> for operations that can fail. Ok(value) represents success and Err(error) represents failure. A function returning Result can use ? after another fallible call:

#![allow(unused)]
fn main() {
launch_kernel()?;
}

If the launch succeeds, execution continues. If it fails, the error is returned to the caller. This replaces the unchecked error codes common in C.

Option<T> represents a value that may be absent. Some(value) contains a value and None means that no value is present. Potential builders use this distinction: Ok(Some(slot)) activates a potential, while Ok(None) means that the configuration did not request it.

Derive attributes and configuration types

An attribute such as #[derive(Debug, Clone, Deserialize)] asks the compiler to generate standard behavior for a struct. In this guide:

• Debug makes a value printable for diagnostics.
• Clone allows the registry framework to copy a builder.
• Deserialize lets Serde construct the value from TOML.
• Serialize and Convert support HeddleMD’s unit-conversion path.

#[serde(deny_unknown_fields)] rejects misspelled or unsupported configuration keys instead of silently ignoring them.

Modules and public re-exports

Rust source files belong to modules. Adding a file is not enough by itself: the parent module must declare it with pub mod my_thermostat;. A statement such as pub use my_thermostat::MyThermostatBuilder; makes the named type available through the parent module. The per-feature guides identify the required declarations and registration entry.

Data layout on the GPU

HeddleMD stores particle data as a structure of arrays (SoA). Instead of an array in which every element contains one particle’s position, velocity, and mass, the engine keeps separate arrays such as positions_x, positions_y, and positions_z.

This layout allows adjacent GPU threads to read adjacent memory locations. Such reads are called coalesced and are substantially more efficient than scattered reads. New device code should use the existing ParticleBuffers layout rather than constructing an array of particle structs.

Deterministic parallel calculations

GPU threads do not finish in a predictable order. This matters because floating-point addition is not associative: (a + b) + c may differ in its least significant bits from a + (b + c).

HeddleMD therefore avoids floating-point sums whose order depends on thread scheduling. Its principal strategies are:

• Sort neighbor lists so that interactions are visited in a fixed order.
• Use fixed-topology reductions when summing floating-point values.
• Use integer fixed-point accumulators where many threads must contribute to the same result. Integer addition is independent of arrival order.
• Avoid floating-point atomicAdd when multiple threads contribute to one value.
• Use seeded, counter-based random-number generation rather than wall-clock or thread-local state.

An atomic operation safely updates shared memory when threads run concurrently. It prevents lost updates, but a floating-point atomic addition does not fix the order of those updates. Thread safety alone is therefore not enough for bitwise reproducibility.

Timesteps, schedules, and CUDA graphs

An integrator does not directly run an entire timestep. It returns a StepPlan, an ordered list of operations such as kicks, drifts, force evaluations, and thermostat or barostat coupling points. The runner validates this schedule and executes it.

This separation makes the schedule the authoritative definition of the physics. It also permits two execution paths:

• The per-step launch path has the CPU launch each kernel separately.
• A CUDA graph records a compatible sequence of device operations and replays it with less CPU launch overhead.

A component is graph-compatible only when all of its per-step work can be recorded as device operations. Reading a device scalar back to Rust, doing host-side arithmetic with it, and then launching another kernel creates host work in the middle of the sequence and normally makes that component incompatible with graph capture. See Schedules and CUDA Graphs for the complete execution model.

Registries, builders, and slots

Most extensible HeddleMD components use three related concepts:

• A registry stores the available implementations.
• A builder validates configuration parameters and constructs one implementation.
• A slot is the constructed runtime component used by a simulation phase, such as its thermostat or barostat.

This design avoids a central conditional containing every supported method. Adding a built-in generally means implementing the appropriate traits and adding its builder to one registry roster. The Extending HeddleMD overview explains this process and introduces the configuration and testing conventions shared by all extension points.

Schedules and CUDA Graphs

The overview sketched a timestep as an ordered list of operations and noted that it runs on the GPU. This page fills in both ideas: what a timestep is — the schedule — and how that schedule is actually launched on the hardware — either a per-step launch loop or a replayed CUDA graph. The two are separate concerns. The schedule defines the physics; the launch mechanism only decides how efficiently that same physics reaches the GPU.

The schedule: a timestep as a list of operations

A textbook velocity-Verlet step is “half-kick, drift, recompute forces, half-kick.” A real step usually does more: it may rescale velocities for a thermostat, rescale the simulation box for a barostat, and project velocities and positions back onto rigid-bond constraints — and each of those slots into the kick/drift sequence at a method-specific point.

Rather than write a bespoke, hand-ordered routine for every combination, HeddleMD builds each timestep as a single ordered list of typed operations called a StepPlan: a half-kick, a drift, a force evaluation, a constraint projection, a barostat point, and so on. Every operation declares the state it reads and the state it writes — its footprint. One executor, run_step, walks the list in order, and before the run begins the list is checked against those footprints for dependency mistakes. An operation may not, for example, read a force that a preceding position update has already invalidated without a force re-evaluation in between; such a plan is rejected rather than silently producing wrong dynamics. Validating Timestep Dependencies explains this check in full and works through a complete example.

This list is the definition of the step. An integrator does not run a timestep by calling kernels itself; it hands back a StepPlan, and the runner executes it. The payoff is that thermostats, barostats, and constraints compose with any integrator without special-casing, and the ordering rules live in exactly one place.

One principle protects that: any optimization that merges adjacent operations into a single GPU launch (fusion) must preserve the schedule’s order and the state each operation observes. Fusion may make the schedule cheaper to run; it may never become a second, competing definition of the step.

Two ways to run the same schedule

Once the schedule for a step is fixed, the runner has two ways to get its kernels onto the GPU.

The per-step launch loop. The straightforward way: every timestep the host walks the schedule and issues each kernel launch itself — on the order of 15–20 individual launches per step (the pair-force kernel, the integrator kicks and drift, the neighbour displacement check, and any thermostat/barostat/constraint kernels). This is simple, and every kernel is timed individually, but each launch carries a small fixed cost paid by the CPU-side driver. For a modestly sized system the per-atom compute of a step can be fast enough that these launch costs become a real fraction of the wall clock — and an MD run issues them millions of times.

Replayed CUDA graphs. A CUDA graph is a feature of the CUDA driver itself (no extra library) that lets you record a sequence of GPU operations once and then replay the whole sequence with a single launch. HeddleMD records one physical step’s kernel sequence into a graph at the start of a phase (one stage of a run — see Phases: the stages of a run), then replays that graph once per step. The ~15–20 per-step launches collapse into a single cuGraphLaunch, and the per-launch driver overhead essentially disappears. Picture it as recording a macro of one step’s GPU work and replaying the macro, instead of re-issuing every command by hand each step.

Graph mode is the default for any phase eligible for it; there is no opt-in. Two knobs under [simulation] control it: graph_batch_size (how many steps replay back-to-back between host check-ins, default 50) and cuda_graphs_disable (a diagnostic escape hatch that forces the per-step loop).

Batched replay

The host still has periodic work that cannot live inside the graph: checking whether atoms have drifted far enough to require rebuilding the neighbour list, and writing trajectory frames and log rows. Graph mode handles this in batches. The graph replays graph_batch_size steps in a row with almost no host involvement; only at the batch boundary does the host download the small neighbour-rebuild flag and write any output that is due. Between boundaries the GPU runs nearly unsupervised — which is where the speed-up comes from. Output cadences shrink a batch when needed so it always ends exactly on a step that produces a log row or trajectory frame.

When graph mode cannot be used

Recording a step into a graph works only if the entire step is device work — GPU kernels one after another, with nothing the CPU must compute or read back in the middle. Most of the engine is deliberately built this way. A few methods are not: the Nosé–Hoover-chain thermostat integrates its extended-system variables with arithmetic on the host and reads the kinetic energy back to the CPU every step, and the MTK-NPT integrator does similar host-side chain arithmetic inside its step. A step containing that kind of host work cannot be captured, so a phase using one of those slots runs on the per-step launch loop instead.

This fallback is automatic and silent — you do not configure it — and it changes nothing about the physics or the results; it only forgoes the launch-overhead saving for that phase. Each slot advertises whether it is graph-compatible, and the runner captures a graph only when every slot in the phase is.

The graph is a recording, not a second timestep

The key conceptual point: a captured graph is not a different definition of the timestep. It is a recording of the same schedule’s kernel sequence — the identical kernels, in the identical order, that the per-step loop would have launched. This is the schedule-is-authoritative principle again: the graph is an execution optimization over the one schedule, never a competing execution model.

That has a precise consequence for reproducibility, worth stating carefully:

• Run-to-run within one path is bit-exact. Two runs of the same config on the same GPU — both in graph mode, or both with cuda_graphs_disable = true — produce byte-identical trajectories and logs. This is the load-bearing invariant, and it is what surfaces a scheduling-dependent bug: such a bug would break byte-identity the moment thread scheduling varied between runs.
• The two paths agree only to floating-point rounding. The per-step loop and the graph replay are two different summation orders, not two spellings of one order. They run the same kernels and match to f32 rounding — and are bit-identical for regularly-packed systems — but on a disordered system they can differ by a single-precision quantum (order 1e-6). That difference is a benign consequence of summation order, not a race; it does mean you should not treat one path as a bit-exact reference for the other.

Stochastic methods keep their random-number counter in a one-element device buffer that both paths share, so a thermostat or barostat draws the same random sequence whichever way the step is launched — the RNG is never a source of divergence between the two paths.

The Reproducibility chapter states the user-facing guarantee and its limits; the exhaustive treatment of schedule orchestration and graph capture and replay lives in docs/architecture.md in the repository.

Validating Timestep Dependencies

The Schedules and CUDA Graphs page introduced a timestep as a single ordered list of operations — a StepPlan — and mentioned in passing that, before a run begins, HeddleMD checks that list for “dependency mistakes.” This page explains that check in full: what it protects against, how it works, and — through two worked examples — how it catches both a within-step multiple-timestep (RESPA) bug and an across-step barostat bug before the simulation ever runs.

The check is worth understanding for anyone adding or modifying an integrator, because an integrator’s whole job is to produce a correct StepPlan. The check is your safety net: it turns an ordering mistake that would otherwise produce silently wrong dynamics into a clear error at startup, naming the exact operation at fault.

Two kinds of state: held and computed

A running simulation keeps several pieces of per-atom and global state. For the purpose of this check they fall into two groups.

Held state is stored directly and is always available to read:

• the atom positions,
• the atom velocities,
• the periodic image flags (which periodic copy of the box each atom is in),
• the simulation box (the lattice).

Computed state is not stored independently — it is calculated from the positions and the box by a force evaluation:

• the forces on the atoms,
• and, when the force field is split for multiple-timestep integration, the separate fast-force and slow-force accumulators (in the code these are the Fast and Slow force classes; more on this below).

The distinction is the whole point. The moment you move the atoms — or change the box — the stored forces no longer describe the arrangement the atoms are now in. They have gone stale. They stay stale until a force evaluation recomputes them at the new positions. Reading stale forces is the single mistake this check exists to catch: a kick that pushes the atoms using forces that belong to where they used to be is wrong, and wrong in a way that a short test run may not obviously reveal.

Every operation declares what it touches

Each operation in a StepPlan declares two things: the state it reads and the state it writes. Together these are the operation’s footprint. The footprints are fixed by the kind of operation, so the engine knows them without running anything:

Two rows carry the key idea. A force evaluation is the only operation that produces forces: it reads the current positions and box and writes fresh forces. A half-kick consumes forces: it reads them to update the velocities. Every dependency the check reasons about is some arrangement of producers and consumers of that computed force state.

One subtlety that matters for the worked example: a force evaluation for a single class — say, only the fast forces — refreshes only that class. It does not revive the slow forces. This is exactly what multiple-timestep integration needs, and it is exactly where the bug below hides.

The check: two passes over the plan

The heart of the check is a single walk over the operations from first to last, keeping track of which computed force state is currently valid. Think of it as reading a recipe top to bottom and, at each step, asking “does this step use an ingredient that is still fresh?”

The walk starts from some set of forces assumed valid, and then, for each operation in order:

• Check its reads. Every piece of state the operation reads must currently be valid. If it reads a force that has gone stale, the walk stops and reports the offending operation.
• Mark forces stale if the atoms moved. If the operation writes the positions or the box, every kind of computed force is now stale — mark them all invalid.
• Record what it produces. If the operation is a force evaluation, mark the forces it computed as valid again.

The reads are checked before the same operation marks anything stale, so a fused kick-and-drift — which reads the forces and then moves the atoms — correctly reads the still-valid forces first and only afterward invalidates them.

A timestep is not run once: the same plan is replayed on every step. So the walk is run twice, from two different starting points, to catch two different kinds of mistake.

• The within-step pass starts by assuming all forces are valid — which is true on the very first step, whose forces come from a warm-up evaluation before the loop begins. This pass catches a force consumer placed after a move inside the same step with no force evaluation between them. Its failure is a stale resource error.
• The across-step pass recognizes that from the second step onward a step does not start with all forces valid — it starts with only the forces the previous run of the plan actually left valid at its end. So this pass starts the walk from “the forces this plan leaves valid at its end” and checks it again. It catches a plan whose opening operations read forces that the plan’s own closing operations invalidated. Its failure is a stale cached forces error.

A plan must pass both. The two worked examples below show one bug each: the first fails the within-step pass, the second fails the across-step pass.

Worked example 1: a within-step RESPA bug

Here is the within-step pass earning its keep on a real class of mistake.

Background: splitting the forces

Multiple-timestep integration (r-RESPA) speeds up a simulation by treating different parts of the force differently. The fast part — bonded terms and short-range interactions — changes quickly and is cheap, so it is evaluated several times per step on a small inner timestep. The slow part — chiefly long-range electrostatics — changes gradually and is expensive, so it is evaluated only once per outer step. That is the entire economy of RESPA: pay for the expensive forces rarely, and keep accuracy by refreshing the cheap forces often.

A correct RESPA outer step, with two inner substeps, is this list of operations:

1.  half-kick from slow force
2.  kick-and-drift from fast force     (inner substep 1)
3.  evaluate fast force
4.  half-kick from fast force
5.  kick-and-drift from fast force     (inner substep 2)
6.  evaluate fast force
7.  half-kick from fast force
8.  evaluate slow force
9.  half-kick from slow force

The slow force brackets the step (operations 1 and 9), and the inner loop (operations 2–7) advances the atoms while refreshing only the fast force. Operation 8 — re-evaluating the slow force once, at the end, before the final slow kick — is easy to overlook. It is a single expensive line in the middle of cheap inner-loop bookkeeping, and the developer already computed a slow force at the top of the step. Forgetting it is a natural mistake.

The bug

Suppose an integrator author omits operation 8. The plan they hand back is:

1.  half-kick from slow force
2.  kick-and-drift from fast force
3.  evaluate fast force
4.  half-kick from fast force
5.  kick-and-drift from fast force
6.  evaluate fast force
7.  half-kick from fast force
8.  half-kick from slow force          <-- was preceded by a slow-force evaluation

Physically, the final slow half-kick now pushes the atoms with a slow force that was computed back at the start of the step — before the two inner drifts moved every atom. The dynamics are wrong. A brief run might look plausible; the energy drift it introduces can take many steps to become obvious, and by then it is hard to trace back to a missing force evaluation.

Walking the buggy plan

The check tracks which forces are valid as it reads the plan top to bottom. The last column shows the state after the operation runs.

At operation 2 the kick-and-drift moves the atoms, so every force — total, fast, and slow — is marked stale. Operations 3 and 6 bring the fast force back, exactly as RESPA intends, but nothing ever brings the slow force back. When operation 8 tries to read the slow force, it is no longer valid. The check stops and reports:

The operation at index 7 (KickHalf) reads a stale resource: the slow-force accumulator.

(The operations are numbered from zero internally, so the eighth operation is reported as index 7.) The message names the operation kind and the exact piece of state that was stale, which points straight at the cause: a slow-force consumer with no slow-force evaluation between the drifts and itself.

The fix

Restore the missing slow-force evaluation before the final slow kick:

7.  half-kick from fast force
8.  evaluate slow force                <-- restored
9.  half-kick from slow force

Re-walking the plan, operation 8 reads the positions and box (always available) and writes a fresh slow force, so the slow force is valid again. Operation 9 then reads a valid slow force, and the walk reaches the end with no error. The plan is accepted and the simulation runs.

The lesson generalizes beyond RESPA: any time an integrator moves the atoms and later consumes a force, there must be a force evaluation in between for the force class it consumes. The within-step pass enforces exactly that, and the separate tracking of fast and slow forces is what lets it distinguish a correct multiple-timestep schedule from one that kicks on a stale accumulator.

Worked example 2: an across-step barostat bug

The within-step pass sees each step as a self-contained recipe. But a step is run over and over, and a mistake can hide in the seam between one step and the next. That is what the across-step pass is for.

The plan

A constant-pressure (NPT) step ends by letting a barostat rescale the simulation box — and every atom’s position with it — to nudge the system toward its target pressure. A minimal velocity-Verlet NPT step looks like:

1.  kick-and-drift from total force   (first half-kick, then move the atoms)
2.  evaluate forces
3.  half-kick from total force        (second half-kick)
4.  barostat point                    (rescale the box and every position)

Read on its own, this step is fine. Operation 1 reads the leftover forces from the previous step, operation 2 refreshes them for the new positions, operation 3 uses those fresh forces, and operation 4 rescales at the very end. The within-step pass is satisfied — nothing after operation 4 reads a force.

(The barostat point only rescales when a barostat is actually configured. On a constant-volume NVE or NVT run the same operation is present but does nothing — it is a placeholder — so it moves no atoms and none of what follows applies.)

The bug across the boundary

But the step repeats. When it runs again, operation 1 reads forces — and those forces are the ones operation 2 computed last time, before operation 4 rescaled the box and shifted every atom. They describe the configuration the system was in before the last rescale, not the one it is in now. The leading half-kick pushes the atoms with slightly wrong forces, and it does so on every step.

The across-step pass catches this. It first asks: what forces does this plan leave valid at its end? Walking the plan, operation 2 makes the forces valid, but operation 4 moves the atoms afterward — so by the end of the plan, no forces are valid. It then re-walks the plan starting from that leftover set (no valid forces) and immediately hits operation 1 reading the total force, which is not valid. The check stops and reports, in spirit:

The operation at index 0 (KickDrift) reads stale cached forces (the total force): the schedule’s own terminal operations invalidated it, with no trailing force evaluation.

The two ways to fix it are named in the message itself.

Two ways to make it valid

1. Add a force evaluation after the barostat. Append “evaluate forces” as a fifth operation. Now the plan leaves fresh forces at its end, the leading half-kick of the next step reads valid forces, and the across-step pass is satisfied. This is exactly correct — but it pays for a second full force evaluation every step, a steep price for what is usually a tiny correction.

2. Declare the barostat weak-coupling tolerant. The built-in pressure- coupling barostats (berendsen, c-rescale) rescale by a very small factor each step, so letting the next step’s leading half-kick consume the slightly-stale forces is a bounded, standard approximation — the same one long-established MD codes make. Such a barostat reports that it tolerates stale cached forces, and the across-step pass then accepts the leftover forces and lets the plan through. The exception is now written down in the barostat’s own code, rather than left as an unstated assumption that a reader has to know about.

An integrator that instead makes a large configuration change at the end of a step and does not declare tolerance is exactly the case the across-step pass exists to reject. And a barostat handled the measure-preserving way — where the box update happens in the middle of the step, before the force evaluation, rather than after it (as the mtk-npt integrator does) — leaves its forces valid at the plan’s end and needs no tolerance at all.

When the check runs, and why it is cheap

The check runs once per simulation phase, at setup — after the integrator produces its plan and before the timestep loop starts. (A phase is one stage of a run, with one fixed choice of integrator, thermostat, and barostat; see Phases: the stages of a run.) If the plan fails, the phase is aborted with an error carrying the offending operation and resource, and the loop never runs. There is no per-step cost.

Along with the plan, the check is handed a small amount of phase context: is a per-step barostat configured, and if so does it tolerate stale cached forces. The first tells the across-step pass whether the barostat point actually rescales (so an NVE/NVT plan is not wrongly flagged); the second tells it whether to accept the leftover forces a weak-coupling rescale leaves behind.

It is inexpensive and safe to run at startup because it reads only the plan and that context — the list of operations, their declared footprints, and two true/false facts about the barostat. It launches no GPU work and touches no positions, velocities, or forces. For the same reason it is easy to test directly — the test suite feeds it hand-written plans, correct and incorrect, and checks that it accepts or rejects each one, all without a GPU.

What the check does and does not cover

A few boundaries are worth knowing when you extend the engine.

• Integrator-private operations must declare their own footprint. Most operations have footprints fixed by their kind, but an integrator may emit a custom operation for a step that does not fit the standard shapes. When it does, it states that operation’s reads and writes explicitly, and the check holds it to the same rule as any built-in: a custom operation that reads a force after a drift, with no evaluation in between, is rejected. Declaring a custom operation’s footprint honestly is the author’s responsibility.

• Kinetic energy, virial, and potential energy are not tracked here. The check reasons about positions, velocities, box, and forces. The scalar quantities that a thermostat or barostat reduces — the kinetic energy, the virial — are computed inside those operations and are not part of this dependency check. The ordering rule that a constraint must publish its virial before a barostat reads it, for instance, is enforced elsewhere in the runner, not by this walk.

• Weak-coupling tolerance covers only forces. When a barostat declares that it tolerates stale cached forces, that applies to the forces alone. Keeping the neighbour list valid across a box change is a separate concern, handled by the rebuild trigger described in the neighbour-list documentation, not by this check.

• This is a correctness check, not an optimization. The check never reorders, merges, or removes operations. Every operation still runs as its own GPU launch, in the order the plan lists them. The check only confirms, before the run, that this order is free of stale-force mistakes.

Why this supports reproducibility

HeddleMD’s central guarantee is that two runs of the same input on the same GPU produce byte-for-byte identical output. That guarantee is only meaningful if the dynamics being reproduced are the correct dynamics. A stale-force bug would reproduce perfectly — it would give the same wrong answer every time — and byte-identity alone would never flag it. By rejecting an ill-ordered plan at startup, this check closes that gap: it ensures the schedule the engine faithfully and reproducibly executes is one whose force dependencies are sound. Together with the reproducible reductions described in the overview, it is part of what lets an integrator author compose kicks, drifts, force evaluations, thermostats, and barostats and trust that a mistake in their ordering becomes a clear error rather than a subtle, repeatable physical wrongness.

Testing extension components

HeddleMD tests an extension at several levels. A focused test can establish that one equation or kernel is correct, while a shared harness establishes that every registered implementation satisfies the same physical contract. Both are necessary: a component can pass an isolated test and still interact incorrectly with constraints, force aggregation, or the timestep runner.

This page introduces the two principal shared harnesses:

• the potential-consistency harness for pair and bonded potentials;
• the slot-conformance harness for thermostats and barostats.

It also explains how these harnesses fit with the feature-specific tests found under tests/.

Test layers

A new extension should normally be tested at four layers:

Do not use a long simulation to replace a small analytical test. A trajectory may reveal that something is wrong without identifying whether the cause is a sign error, unit conversion, reduction, placement, or random-number sequence. Conversely, do not rely only on direct kernel tests: they bypass the registry, runner, and other components with which the extension must compose.

Running the tests

The shared positive-path harnesses execute the real GPU pipeline and therefore require a supported NVIDIA GPU and CUDA environment. Run one integration-test target with Cargo:

cargo test --test potential_consistency
cargo test --test slot_conformance
cargo test --test barostats_berendsen

To run one named test while developing, add its name as a filter:

cargo test --test potential_consistency correct_pair_passes_finite_difference

Some helper and negative-control tests are CPU-only, but Cargo still builds the test target and its GPU-facing code. Treat a missing GPU or CUDA installation as an environment failure, not as evidence about the physical implementation.

Potential-consistency harness

The implementation lives in tests/common/consistency.rs, and tests/potential_consistency.rs drives its positive and negative scenarios. Its detailed requirements are recorded in rqm/forces/potential-consistency-harness.md.

The harness evaluates one isolated potential at a series of geometries. For each geometry it obtains the total energy, force on every atom, and scalar virial from the real force-field pipeline. It then checks invariants that do not depend on the potential’s particular functional form:

• Force–energy consistency: a central finite difference of the energy agrees with each Cartesian force component (F = -grad U).
• Newton’s third law: the forces sum to zero for the isolated interaction.
• Virial consistency: the reported scalar virial agrees with the energy change under an infinitesimal isotropic coordinate scaling.
• Reference values: optional analytical energies or generalized forces are reproduced at selected coordinates.
• Cutoff behavior for pair potentials: the switched interaction joins smoothly and vanishes at its cutoff. This check does not apply to fragments explicitly declared unbounded.

These checks catch common CUDA-fragment errors such as reversing a force sign, omitting the division by distance, splitting an energy twice, producing asymmetric forces, or reporting an inconsistent virial.

The fixture and evaluator

A ConsistencyFixture describes what the harness should evaluate. Its most important fields are:

• label, used in failure messages and coverage checks;
• shape, one of Pair, Bond, Angle, or Dihedral;
• samples, the distances, angles, or dihedrals at which to test;
• a geometry function that turns each sample coordinate into atom positions;
• a builder for an isolated ForceField containing the potential;
• tolerances and finite-difference step sizes;
• optional analytical ReferencePoint values;
• switching and cutoff information for pair potentials.

The shape constructors—ConsistencyFixture::pair, bond, angle, and dihedral—supply the standard geometry and topology for their interaction type. pair_correction supports an unbounded correction-only pair fragment.

GpuEvaluator is the bridge to the production code. For every requested geometry it creates fresh ParticleBuffers, runs ForceField::step(..., AggregateLevel::ForcesAndScalars), copies the forces, per-particle energies, and virials back to the host, and returns their totals as an Eval.

The checking functions are generic over the Evaluator trait. Negative tests can therefore use a small CPU evaluator containing a deliberate defect. This proves that the check fails for the mistake it is intended to detect without adding a broken CUDA kernel to the project.

Adding a pair-potential fixture

Add the new fixture to builtin_consistency_fixtures() in tests/common/consistency.rs. A pair fixture has the following conceptual form; the exact configuration types should match the new potential:

#![allow(unused)]
fn main() {
ConsistencyFixture::pair(
    "buckingham",
    |registry| registry.register(Box::new(BuckinghamBuilder)),
    particle_types,
    pair_interactions,
    charges,
    None,                 // optional SPME configuration
    cutoff,
    r_switch,
    vec![r1, r2, r3],     // informative, non-singular sample distances
)
.with_reference_points(vec![
    ReferencePoint {
        coordinate: r_ref,
        energy: Some(expected_energy),
        coord_force: Some(expected_radial_force),
        tol: Tolerance::new(relative_tolerance, absolute_tolerance),
    },
])
}

Choose samples that exercise the repulsive and attractive regions and, when applicable, the switched region. Avoid singular points and configurations whose expected force is so close to zero that a sign or scale error becomes invisible. Reference points are especially useful because finite-difference consistency alone cannot detect an energy and force that are both wrong by the same constant factor.

The fixture should construct a registry containing only the builder under test. This isolates its contribution from other potentials. The fixture constructors already provide a sufficiently large box and an all-pairs neighbor list for the small test system.

Coverage is enforced

assert_fixture_coverage compares fixture labels with the fragment labels reported by the built-in potential registry. Adding a built-in fragment without a fixture therefore fails the coverage test. A fixture whose label no longer corresponds to a registered fragment also fails.

assert_all_builtin_potentials_consistent runs the complete invariant suite over every built-in fixture. During development, an individual test such as check_force_energy can provide a narrower failure, but the complete sweep is the final shared-contract check.

The consistency harness does not replace tests for parameter claims, unit conversion, argument-schema strings, fragment metadata, exclusions, mixed particle types, or an empty interaction list. Keep those focused tests in the feature file and in tests/potential_claims.rs or tests/forces_framework.rs as appropriate.

Barostat tests

Barostats have two complementary test styles:

1. Direct feature tests exercise construction, reductions, scale factors, box mutation, random counters, and empty-system behavior.
2. The shared slot-conformance harness checks whether every registered barostat maintains a physically meaningful density in a complete constrained-water simulation.

The direct tests are grouped by implementation:

• tests/barostats_berendsen.rs for deterministic per-step coupling;
• tests/barostats_c_rescale.rs for stochastic per-step coupling;
• tests/barostats_mc.rs for periodic Monte Carlo volume moves.

Direct per-step tests

A typical direct fixture constructs:

• a GpuContext from init_device();
• a SlotConfig from a kind and TOML parameter string;
• a barostat through BarostatRegistry::with_builtins().build_optional(...);
• a small ParticleState with prescribed positions, velocities, masses, and per-particle virials;
• ParticleBuffers, a SimulationBox, and Timings.

Call the public Barostat trait method when testing behavior shared by the framework. Tests may inspect a concrete built-in type when verifying method-specific diagnostics, but prefer safe public accessors where available. An unsafe cast from Box<dyn Barostat> to a concrete type is tightly coupled to the registry roster and should not be copied into ordinary downstream code.

For a per-step barostat, focused tests should cover:

• registry discovery and construction, including zero particles;
• parameter validation and conversion;
• the deterministic kinetic-energy and virial reductions;
• the analytical scale factor for controlled kinetic energy, virial, box volume, and timestep values;
• changes to the box, positions, and box-generation counter;
• no changes to velocities unless the method explicitly requires them;
• no kernel launch or box mutation for an empty particle set;
• diagnostic and conserved-quantity bookkeeping;
• identical results from identical seeds and device-counter states for a stochastic method.

Tests should populate the per-particle virials deliberately. Calling barostat.apply with an arbitrary zeroed buffer does not verify whether the barostat consumes the pressure information produced by the force pipeline. At least one integration test should run a real force evaluation with AggregateLevel::ForcesAndScalars before applying the barostat.

Direct periodic-move tests

A periodic barostat is host-orchestrated and uses apply_move rather than the per-step apply hook. Its fixture also needs a mutable ForceField, molecule tables, and a topology-compatible state.

In addition to configuration and construction, test:

• the declared BarostatPeriodicity and move frequency;
• the proposal counter advancing once per attempted move;
• acceptance and rejection paths;
• exact restoration of positions, forces, and lattice after rejection;
• no velocity changes;
• rigid intramolecular geometry under molecular center-of-mass scaling;
• cutoff/box-width guards;
• deterministic decisions for identical seeds and counters;
• graph-mode and ordinary-launch behavior at host batch boundaries.

Slot-conformance harness for barostats

The shared implementation is in tests/common/slot_conformance.rs, and tests/slot_conformance.rs contains the test sweep. The harness runs each built-in barostat on dense, rigid SPC/E water with SETTLE constraints and a CSVR thermostat. It then checks that the mean density over the equilibrated part of the phase lies within the case’s declared interval.

The system is intentionally more demanding than a minimal smoke test. It is dense enough for pressure and constraint-virial errors to be visible. The barostat run uses a thermostat that has its own conformance case, helping to attribute a density failure to the barostat rather than uncontrolled thermal dynamics.

Each built-in barostat has a SlotCase in builtin_barostat_cases():

#![allow(unused)]
fn main() {
SlotCase {
    kind: "my-barostat",
    expect: Expect::HoldsDensity { lo, hi },
    note: "brief explanation of the method and the chosen bounds",
}
}

To add a built-in barostat to the harness:

1. Add a SlotCase with a scientifically justified density interval and a note explaining it.
2. Add a named #[test] in tests/slot_conformance.rs that finds the case and passes it to run_barostat_case.
3. Update the barostat sweep count checked by every_case_in_the_tables_is_driven_by_a_test.
4. Run the new test and the complete slot-conformance target.

assert_barostat_coverage compares the case table with BarostatRegistry::with_builtins(). A registered barostat without a case, or a stale case without a registered builder, fails immediately. The separate sweep-count test ensures that adding a table row is not enough: a named test must actually execute it.

The density bounds are regression and conformance bounds for the harness’s finite system, cutoff, and run length. They are not a claim that every valid barostat must reproduce an experimental bulk density within the same narrow range under arbitrary settings. Record the reason for the interval in the case’s note so future changes do not widen it without understanding what defect it was designed to catch.

The harness also contains negative controls that feed historically incorrect temperature or density values directly to its checking functions and assert that they fail. These tests demonstrate that the tolerances can discriminate the targeted defect; they do not require a GPU simulation.

Choosing tolerances

A useful tolerance is wide enough to accommodate the known numerical noise of the tested precision and finite system, but narrow enough to fail for a plausible implementation error.

Use these principles:

• Base analytical-test tolerances on the precision and reduction path, not on the size of the observed error after the test fails.
• Include an absolute floor when the expected value can approach zero.
• Use several sample coordinates so a cancellation at one point cannot hide a defect.
• Give every broad statistical interval a physical and finite-size rationale.
• Add a negative control when a tolerance protects against a known class of error. A check that has never been shown to reject a representative defect may provide false confidence.

Checklist for a new extension

• Configuration parsing, domain validation, unknown fields, and unit conversion are tested.
• Registry construction and an empty particle or interaction set are tested.
• Small analytical tests isolate the method’s equations and kernel outputs.
• The appropriate shared harness contains a fixture or conformance case.
• Registry-to-harness coverage checks pass.
• Tolerances are justified and discriminate meaningful defects.
• Stochastic tests control both the seed and counter state.
• At least one end-to-end test exercises composition with the runner and relevant constraints, thermostat, or force field.
• Deterministic repeated runs are compared where the extension changes per-step numerical work.

Extending HeddleMD

This section explains how to add new physical models to HeddleMD. It is aimed at scientists who are comfortable describing a molecular-dynamics algorithm but may be new to Rust and CUDA. Before continuing, read the Developer Guide overview, especially its explanations of host and device code, Rust traits, device buffers, and deterministic reductions.

Choose the guide that matches the physical role of your new method:

• Adding a thermostat
• Adding an integrator
• Adding a barostat
• Adding a pair potential
• Adding a bonded potential

The guides use the same overall workflow: choose a nearby built-in, describe and validate the configuration, implement the runtime behavior, register the builder, and test both the physics and the framework contract.

Shared terminology

Decide how the component is selected

HeddleMD has seven registries, collected in Registries in src/registries.rs:

The first six use named selection. For example, the configuration [phase.thermostat] kind = "csvr" asks the thermostat registry for the builder whose kind_name() is "csvr". Only one thermostat and one integrator can occupy their respective slots in a phase.

Potentials use compositional activation because a force field normally contains several contributions. Every potential builder examines the configuration. It returns a component when its interaction type is present and returns no component otherwise. The force field is the sum of all active potential components.

Understand the builder and runtime component

Most extensions contain two conceptually different Rust types:

1. The builder is small and normally contains no fields. It owns the configuration schema, validates values, and constructs the component.
2. The runtime component holds parameters, device buffers, and any state that changes during a simulation. It implements Thermostat, Barostat, Integrator, or Potential.

For named components, the builder also implements KindedBuilder, which supplies the configuration name and unit conversion. A typical declaration is struct MyThermostatBuilder;: the trailing semicolon defines a zero-field struct. #[derive(Clone)] is sufficient for the registry’s generated cloning support.

Define and validate configuration parameters

Named slots are initially parsed into a kind string and an untyped toml::Value containing the remaining parameters. The selected builder owns the interpretation of that value. Its validate_params method should:

1. Deserialize the TOML value into a typed parameter struct.
2. Reject unknown fields with #[serde(deny_unknown_fields)].
3. Check physical and numerical restrictions, such as finite positive temperatures, timescales, and cutoffs.
4. Return a user-facing ConfigError that identifies the invalid field.

Potentials route parameter tables through a parameter claim. A claim is a category and kind pair, such as (PairInteraction, "buckingham"). It tells the loader which builder validates and converts each interaction entry.

Preserve physical units

Configuration values may be expressed in SI or Hartree atomic units. Values must be converted to atomic units before the runtime component is built. Use the dimension types in crate::units rather than representing every value as a plain f64:

#![allow(unused)]
fn main() {
#[derive(
    Debug,
    Clone,
    serde::Deserialize,
    serde::Serialize,
    crate::units::Convert,
)]
#[serde(deny_unknown_fields)]
struct MyParams {
    temperature: crate::units::Temperature,
    tau: crate::units::Time,
    seed: u64,
}
}

Temperature and Time are converted according to the user’s unit system. The dimensionless integer seed is left unchanged. For a named builder, convert_params can normally delegate to convert_params_in_place::<MyParams>(units, params).

Do not manually apply conversion factors in build. That duplicates policy and makes it easy for SI and atomic-unit inputs to behave differently.

Add GPU work only when necessary

First look for an existing launch wrapper in src/gpu/kernels.rs. Many new thermostats and integrators can reuse the kinetic-energy, rescaling, kick, and drift kernels already provided by HeddleMD.

If the physical method requires a new CUDA kernel, the wiring is:

1. Add a .cu file under kernels/. build.rs discovers it and compiles it to embedded PTX automatically.
2. Declare its kernel handles with gpu_kernels! in the appropriate Rust subsystem file.
3. Add that handle set to the define_kernels! manifest in src/gpu/device.rs.
4. Add a launch wrapper in src/gpu/kernels.rs and re-export it from src/gpu/mod.rs.

Use the types and mathematical functions from kernels/precision.cuh, such as Real, R(x), and Real_sqrt, instead of spelling the precision as float or double. Launch wrappers should use the project’s gpu_launch! conventions and should not contain their own timing calls; the caller brackets the launch with the appropriate timing stage.

Pair and bonded potentials are a special case: their physical calculation is usually supplied as a CUDA source fragment that the framework combines and compiles at run time. Their individual guides explain that mechanism.

Preserve reproducibility

Every new component must preserve HeddleMD’s fixed-GPU reproducibility guarantee. In particular:

• Do not accumulate floating-point values in an order determined by thread arrival. Avoid multi-writer floating-point atomicAdd.
• Use existing deterministic reductions for kinetic energy, virial, and force contributions.
• Make per-particle or per-interaction CUDA functions pure: their outputs should depend only on their inputs, not on mutable shared state.
• For stochastic per-step device work, use Philox with an explicit seed and a device-resident counter. A host counter does not advance correctly when a captured CUDA graph is replayed.
• Allocate device scratch during construction and reuse it.

Run the same simulation twice and compare outputs in an end-to-end reproducibility test when the new algorithm changes per-step numerical work.

Register the builder

To include a component with HeddleMD, add its builder to the appropriate Builtins::builtins() roster shown in the table above. Module declarations and public re-exports are separate Rust steps; the feature-specific guide lists each required edit.

Libraries embedding HeddleMD can also create Registries and call methods such as register_thermostat or register_potential before invoking run_simulation_with_registries. Registration appends to the registry. A named registry uses the first builder matching a kind, so a later builder with the same name does not replace an earlier built-in. Potential registries are different: every builder is consulted during construction, while the first matching parameter claim owns validation and conversion. A duplicate potential claim can therefore construct an additional component without owning its configuration processing. Use distinct names and parameter claims unless constructing the corresponding registry from empty components.

Document and test the extension

Requirements specifications live under rqm/. For a substantial physical method, add a neighboring specification that records its equations, parameters, execution sequence, limitations, determinism requirements, and behavior for an empty particle set.

The Testing Extension Components chapter explains the potential-consistency and thermostat/barostat conformance harnesses. Read it before adding a built-in potential, thermostat, or barostat: registry coverage tests require every built-in to declare an appropriate shared fixture or conformance case.

Tests normally cover three layers:

1. Configuration: valid values, invalid domains, unknown fields, and unit conversion.
2. Framework integration: registry discovery, successful construction, empty-system behavior, and compatibility checks.
3. Physics and numerics: a known analytical case, conservation or target behavior where appropriate, and deterministic repeated runs.

Each feature-specific guide ends with a concise file checklist. Use it after understanding the implementation path, rather than treating it as a substitute for the preceding explanation.

Adding a thermostat

A thermostat couples particle motion to a target temperature. This guide explains how to add a thermostat that occupies the optional [phase.thermostat] configuration slot.

Use this extension point when the thermostat is separate from the integrator. It can be combined with velocity-verlet or respa. Do not use it for an integrator such as langevin-baoab or mtk-npt, which owns its thermostat and rejects a separately configured thermostat.

Read the extension overview first for shared configuration, registration, unit-conversion, and determinism conventions.

Choose the closest example

Start from the simplest existing thermostat that resembles your algorithm:

• src/integrator/berendsen.rs is the best starting point for a deterministic velocity-rescaling method. Its per-step calculation remains on the GPU.
• src/integrator/csvr.rs demonstrates a stochastic rescaling method, a device-resident Philox counter, and conserved-quantity bookkeeping.
• src/integrator/andersen.rs demonstrates a stochastic method that requires its own CUDA kernel.
• src/integrator/nose_hoover_chain.rs is an advanced example with coupling both before and after integration and host-side chain arithmetic. It cannot be captured in a CUDA graph.

Copying a nearby implementation is preferable to beginning with the full skeleton below: the existing file includes complete error handling, timing, and tests for its particular execution path.

Understand when a thermostat runs

A runtime thermostat implements the Thermostat trait in src/integrator/mod.rs. Two methods can couple it to a timestep:

• apply_pre runs before the integrator walk. Its default implementation does nothing. Override it only when the mathematical splitting requires an initial coupling operation, as in a symmetric Nosé–Hoover scheme.
• apply_post runs after the trailing velocity update. Every separately configured thermostat must implement this method.

The runner calls these methods only on coupling steps. The dt argument is already the effective coupling time, coupling_interval * simulation_timestep. Do not multiply it by the coupling interval again. The runner removes coupling_interval before passing the remaining parameter table to your builder.

Both methods must return Ok(()) without launching a kernel when buffers.particle_count() == 0. Empty particle sets are used by framework tests and must remain valid.

Thermostats may also provide log_column_names and log_column_values for diagnostic or conserved-quantity columns. Follow the naming and unit conventions of the closest existing thermostat.

Implement a velocity-rescaling thermostat

Most rescaling thermostats perform three operations:

1. Calculate the kinetic energy after the integrator’s final velocity update.
2. Calculate a velocity scale factor, conventionally called lambda, from the kinetic energy, target temperature, coupling time, and method-specific parameters.
3. Multiply every particle velocity by that factor.

The kinetic-energy calculation is a global reduction and the velocity update is a per-particle operation. They must remain separate kernel launches; trying to combine the reduction with the velocity update would require the order-dependent synchronization that HeddleMD avoids.

Preferred path: keep intermediate values on the GPU

Use this path when the scale-factor calculation can be expressed as a CUDA kernel:

1. Allocate a one-element CudaSlice<Real> for the kinetic energy and another for the scale factor when constructing the thermostat.
2. Call compute_kinetic_energy_on_device. It writes the result into the kinetic-energy buffer without copying it to Rust.
3. Launch a small method-specific kernel that reads the kinetic energy and writes the scale factor.
4. Call rescale_velocities_device_factor with the scale-factor buffer.

Because all per-step operations remain on the device, this form can retain the default graph_compatible() == true. Berendsen and CSVR follow this pattern.

Alternative path: calculate the factor on the host

If the algorithm requires substantial CPU-side state or arithmetic, call compute_kinetic_energy, which performs the same deterministic reduction and then copies the scalar result to the host. Calculate lambda in Rust and pass it to rescale_velocities.

This device-to-host copy and intervening Rust calculation cannot be recorded as an uninterrupted CUDA graph. The builder must therefore override graph_compatible and return false. Nosé–Hoover chain is the relevant example.

Do not mix the two paths accidentally. In particular, calling the host-returning compute_kinetic_energy while leaving graph compatibility at its default is incorrect.

Add stochastic sampling safely

Stochastic thermostats use the counter-based Philox generator. Store the configured seed as a u64; it is dimensionless and is not changed by unit conversion.

For per-step CUDA work, allocate a one-element draw counter on the device. The sampling kernel reads and advances that counter. A captured CUDA graph then advances the random stream on every replay. A counter stored only in Rust would be captured with a fixed value and would repeat random draws.

Use src/integrator/csvr.rs or src/integrator/andersen.rs as the complete pattern. The host implementation in src/integrator/philox.rs and the device helper in kernels/philox.cuh produce matching draws for the same keys.

Define the Rust types

The implementation file normally contains four pieces:

1. A typed parameter struct with physical dimension types.
2. A runtime thermostat holding converted parameters and device buffers.
3. An implementation of Thermostat for that runtime type.
4. A builder implementing KindedBuilder and ThermostatBuilder.

The following is an abbreviated outline, not compilable code:

#![allow(unused)]
fn main() {
#[derive(
    Debug,
    Clone,
    serde::Deserialize,
    serde::Serialize,
    crate::units::Convert,
)]
#[serde(deny_unknown_fields)]
struct MyThermostatParams {
    temperature: crate::units::Temperature,
    tau: crate::units::Time,
}

struct MyThermostat {
    temperature: f64,
    tau: f64,
    kinetic_energy: cudarc::driver::CudaSlice<Real>,
    scale_factor: cudarc::driver::CudaSlice<Real>,
}

impl Thermostat for MyThermostat {
    fn apply_post(
        &mut self,
        buffers: &mut ParticleBuffers,
        dt: Real,
        timings: &mut Timings,
    ) -> Result<(), ThermostatError> {
        if buffers.particle_count() == 0 {
            return Ok(());
        }

        // Time and launch the deterministic kinetic-energy reduction.
        compute_kinetic_energy_on_device(buffers, &mut self.kinetic_energy)?;

        // Launch a method-specific kernel that computes self.scale_factor.

        // Time and apply the device-resident scale factor.
        rescale_velocities_device_factor(buffers, &self.scale_factor)?;
        Ok(())
    }
}

#[derive(Debug, Clone)]
struct MyThermostatBuilder;

impl KindedBuilder for MyThermostatBuilder {
    fn kind_name(&self) -> &'static str {
        "my-thermostat"
    }

    fn convert_params(
        &self,
        units: crate::units::UnitSystem,
        params: &mut toml::Value,
    ) -> Result<(), ConfigError> {
        convert_params_in_place::<MyThermostatParams>(units, params)
    }
}
}

In validate_params, deserialize MyThermostatParams and reject non-finite or non-positive temperatures and timescales. In build, deserialize the already converted values, compute the thermal degrees of freedom consistently with the existing thermostats, allocate device buffers once, and return Box<dyn Thermostat>.

Register the thermostat

In src/integrator/mod.rs:

1. Declare the module with pub mod my_thermostat;.
2. Re-export the runtime type and builder with the other thermostat exports.
3. Add Box::new(MyThermostatBuilder) to the built-in thermostat roster.

The roster entry makes kind = "my-thermostat" resolvable. No central configuration enum or dispatch statement needs to change.

If you add a CUDA file, also follow the kernel wiring described in the extension overview. A thermostat that reuses existing kernels does not need changes to build.rs or the aggregate device-kernel manifest.

Write the specification and tests

Every built-in thermostat is covered by the shared slot-conformance harness in tests/common/slot_conformance.rs. Add a SlotCase to builtin_thermostat_cases() that names the new kind, specifies an Expect::HoldsTemperature tolerance, and explains the choice of tolerance. Then add a named test in tests/slot_conformance.rs that passes the case to run_thermostat_case, and update the thermostat sweep count in that file.

The coverage check compares the case table with ThermostatRegistry::with_builtins(), so registering a built-in thermostat without a case fails the test suite. The harness runs dense SPC/E water with SETTLE constraints; it is intended to expose coupling to the wrong degrees of freedom. See Testing Extension Components for the system, negative controls, and tolerance guidance.

Add rqm/integration/my-thermostat.md. Use rqm/integration/csvr.md as a structural example and document:

• the physical equations and coupling location;
• parameters, units, and domain restrictions;
• the per-step kernel sequence;
• random-number keys and counter advancement, if stochastic;
• any conserved or diagnostic quantity;
• empty-system behavior and CUDA-graph compatibility;
• the determinism guarantee.

Add an integration test modeled on the closest existing thermostat. At a minimum, test:

• builder discovery and successful construction;
• valid, invalid, unknown, and unit-converted parameters;
• no work for an empty particle set;
• relaxation or sampling toward the expected temperature;
• the analytical scale factor for a controlled state;
• repeated-run identity for stochastic methods;
• the declared CUDA-graph compatibility.

Completion checklist

• The method belongs in a separate thermostat slot rather than an integrator-owned thermostat.
• dt is used as the effective coupling time without additional scaling.
• Device scratch is allocated during construction, not during a step.
• apply_post handles an empty particle set without launching work.
• Host/device transfers agree with the graph_compatible declaration.
• Stochastic device work uses a device-resident Philox counter.
• Parameter units and physical domains are validated.
• The builder is declared, re-exported, and added to the built-in roster.
• A built-in thermostat has a SlotCase and named slot-conformance test.
• Physics, configuration, empty-state, and determinism tests are present.

Adding an integrator

An integrator defines the ordered velocity updates, position updates, and force evaluations that make up one timestep. The user selects it by name with [phase.integrator] kind = "...".

Read the extension overview first. This guide assumes you understand its descriptions of builders, runtime components, host and device work, and CUDA-graph compatibility.

Unlike a conventional integrator routine, a HeddleMD integrator does not run the whole timestep itself. It describes the timestep by returning a StepPlan, which is an ordered list of typed SubStep values. The shared runner executes that list. This allows the runner to insert configured thermostat, barostat, and constraint operations, validate data dependencies, and capture compatible device work in a CUDA graph.

Choose the closest example

• src/integrator/velocity_verlet.rs — the canonical symplectic template. Start here for an integrator with ordinary kick and drift operations.
• src/integrator/langevin_baoab.rs — owns its own thermostat; a Custom sub-step (the OU step) with a device-resident RNG counter.
• src/integrator/mtk_npt.rs — owns thermostat and barostat; reads virial every step (ForcesAndScalars); graph_compatible = false; log columns.
• src/integrator/respa.rs — multiple force evaluations per step and KickSource::Class impulse splitting.

Before writing Rust, express the algorithm as an ordered list of operations and identify which forces each velocity update consumes. This makes it much easier to construct a valid StepPlan.

Describe the timestep with a StepPlan

plan(dt) -> StepPlan must be a pure function of dt and the integrator’s fixed configuration. In this context, pure means that it does not modify the integrator, launch kernels, allocate buffers, or depend on a per-step counter. It must return the same sequence of operation types every time. The runner relies on that stable sequence when validating the phase and recording CUDA graphs.

The canonical symplectic (velocity-Verlet family) shape reads F(t) before the force evaluation and F(t+dt) after it:

ConstraintPoint { BeforeDrift }      // snapshot; no-op without a constraint slot
KickDrift { Total }                  // v += (F(t)/m)·dt/2 ; x += v·dt   (fused)
ConstraintPoint { AfterDrift }       // project positions
ForceEval { class: None, level }     // runner recomputes F(t+dt)
KickHalf  { Total }                  // v += (F(t+dt)/m)·dt/2
BarostatPoint                        // terminal; no-op without a per-step barostat
ConstraintPoint { AfterKick }        // project velocities — RATTLE-last

The following rules explain how that sequence is executed:

• execute() never sees ForceEval. Only the runner holds the ForceField, so it dispatches force evaluation itself. Your execute handles the kicks, drifts, and any integrator-private Custom steps, and returns IntegratorError::UnexpectedSubStep for anything it shouldn’t receive.
• Forces at the start of a step are already available. At step entry the runner treats the cached forces as valid F(t) (left by the previous step’s trailing ForceEval, or the runner’s one warm-up evaluation before the loop). So the leading KickDrift legitimately reads F(t) — do not prepend a ForceEval to “prime” forces; that would double-evaluate, and the schedule validator already assumes carried-in forces are valid.
• Placement markers are harmless when no matching component is configured. Emit ConstraintPoint, BarostatPoint, and (if you own thermostat placement) ThermostatHalf markers unconditionally, so one plan drives constrained and unconstrained, NVE and NPT runs alike. Order matters: a terminal BarostatPoint goes before the final AfterKick velocity projection (RATTLE-last).
• AggregateLevel controls which force-derived quantities are calculated. On the ForceEval, it selects the cheap ForcesOnly path unless the integrator itself needs energy/virial every step (NPT integrators emit ForcesAndScalars). Use None to defer to the runner, which upgrades to ForcesAndScalars on steps that log or write a frame.
• KickSource identifies the force buffer used by a kick. Total reads the combined force buffer and is executed by the integrator. Class(Fast|Slow) reads one class accumulator and is dispatched by the runner — use it only for multiple-timestep (RESPA) impulse splitting.

The schedule is validated at phase setup: StepPlan::validate() rejects a plan that reads force-derived state a preceding position/box mutation invalidated without an intervening ForceEval. Each built-in SubStep has a fixed resource footprint; a Custom step declares its own reads/writes (see src/integrator/op_model.rs and rqm/integration/op-model.md). Model a Custom step’s label-and-footprint as a single enum, like MTK’s MtkStep, so execute re-derives the same footprint it declared.

Declare compatibility with other components

The builder provides several methods that tell the configuration loader which combinations are valid. These methods receive the parsed parameters, so an answer may depend on the selected mode:

velocity-verlet computes supports_constraints as !lossless (the lossless mode can’t project). Set graph_compatible = false if execute does host-side scalar arithmetic or a device→host copy between launches (MTK does). A separate runtime guard, BarostatPlacementMissing, fires if a per-step barostat is configured with a plan that carries no BarostatPoint — so a barostat-hosting integrator must emit one.

If the integrator supports constraints, also implement ConstraintCapableIntegrator. The builder predicate says that the configuration is supported; this runtime trait lets the constraint framework exercise the constructed integrator and reject a state-dependent case.

Reuse device kernels and preserve determinism

Reuse the existing integrate.cu kernels where you can — vv_kick, vv_kick_drift, and the class-kick kernels are re-exported from crate::gpu, so a standard symplectic integrator needs no new CUDA. If you add one, follow the overview CUDA section. Keep RNG state in a device-resident counter (like Langevin’s draw_counter_device) so it survives CUDA-graph replay deterministically; a host counter would freeze under capture and break reproducibility.

Implement the Rust types

The main implementation belongs in a new file such as src/integrator/my_integrator.rs. It contains the parameter type, runtime state, Integrator implementation, and builder. The following skeleton is an abbreviated outline: comments replace method-specific calculations and error mapping, so it is not directly compilable.

#![allow(unused)]
fn main() {
use serde::Deserialize;
use crate::gpu::{GpuContext, ParticleBuffers, vv_kick, vv_kick_drift};
use crate::io::config::ConfigError;
use crate::pbc::SimulationBox;
use crate::timings::{KernelStage, Timings};
use crate::precision::Real;
use super::{Integrator, IntegratorBuilder, IntegratorError,
            KickSource, StepPlan, SubStep, ConstraintPhase};

#[derive(Debug, Clone, Deserialize)]         // add serde::Serialize + Convert
#[serde(deny_unknown_fields)]                // if any field is unit-bearing
pub struct MyIntegratorParams {
    #[serde(default)] pub lossless: bool,
    // pub tau: crate::units::Time,          // unit-bearing → dimension newtype
}

#[derive(Debug)]
pub struct MyIntegratorState { /* device scratch, RNG counters, … */ }

impl Integrator for MyIntegratorState {
    fn plan(&self, dt: Real) -> StepPlan {           // MUST be pure
        StepPlan { steps: vec![
            SubStep::ConstraintPoint { phase: ConstraintPhase::BeforeDrift, dt },
            SubStep::KickDrift { dt, label: "my_kick_drift", source: KickSource::Total },
            SubStep::ConstraintPoint { phase: ConstraintPhase::AfterDrift, dt },
            SubStep::ForceEval { class: None,
                                 level: Some(crate::forces::AggregateLevel::ForcesOnly) },
            SubStep::KickHalf { dt, label: "my_kick", source: KickSource::Total },
            SubStep::BarostatPoint { dt },           // terminal, before AfterKick
            SubStep::ConstraintPoint { phase: ConstraintPhase::AfterKick, dt },
        ]}
    }
    fn execute(&mut self, substep: &SubStep, buffers: &mut ParticleBuffers,
               sim_box: &mut SimulationBox, timings: &mut Timings)
        -> Result<(), IntegratorError> {
        if buffers.particle_count() == 0 { return Ok(()); }        // mandatory
        match substep {
            SubStep::KickDrift { dt, .. } => { /* time + vv_kick_drift(buffers, sim_box, *dt) */ Ok(()) }
            SubStep::KickHalf  { dt, .. } => { /* time + vv_kick(buffers, *dt) */ Ok(()) }
            other => Err(IntegratorError::UnexpectedSubStep { variant: other.variant_name() }),
        }
    }
}

#[derive(Debug, Clone)]
pub struct MyIntegratorBuilder;

impl crate::registry::KindedBuilder for MyIntegratorBuilder {
    fn kind_name(&self) -> &'static str { "my-integrator" }
    // convert_params via convert_params_in_place::<MyIntegratorParams> — unit-bearing only
}

impl IntegratorBuilder for MyIntegratorBuilder {
    fn validate_params(&self, params: &toml::Value) -> Result<(), ConfigError> { /* … */ Ok(()) }
    // override owns_thermostat / owns_barostat / supports_barostat /
    // supports_constraints / graph_compatible where they differ from defaults
    fn build(&self, gpu: &GpuContext, particle_count: usize, n_constraints: usize,
             params: &toml::Value) -> Result<Box<dyn Integrator>, IntegratorError> {
        Ok(Box::new(MyIntegratorState { /* … */ }))
    }
}
}

Register and document the integrator

Add rqm/integration/my-integrator.md for the requirements specification; follow rqm/integration/velocity-verlet.md and the op-model conventions in rqm/integration/op-model.md.

Then edit the following existing files:

• src/integrator/mod.rs — declare pub mod my_integrator;, add the pub use my_integrator::{MyIntegratorBuilder, MyIntegratorState}; re-export, and Box::new(MyIntegratorBuilder) in the vec![...] of impl Builtins for dyn IntegratorBuilder (roster at ~line 958).

• src/gpu/device.rs — only if you added a .cu: one field in the define_kernels! manifest.

No change needed to src/io/config.rs (compatibility checks call your builder’s predicates; kind dispatch is registry-driven) or build.rs.

Test the integrator

There is currently no registry-wide physical conformance harness for integrators analogous to the potential-consistency or thermostat/barostat slot-conformance harnesses. Integrators instead share schedule validation and framework tests. Add pure, GPU-free plan assertions to tests/integrator_framework.rs when the assertion is an invariant that every built-in integrator should satisfy; keep method-specific plan and numerical tests in the new integrator’s test file.

The testing chapter explains how shared-contract tests complement focused analytical and end-to-end tests. For an integrator, the end-to-end portion is particularly important because it exercises carried forces and composition with thermostat, barostat, and constraint markers.

tests/integrators_my_integrator.rs (model on tests/integrator_framework.rs and tests/integrators_respa.rs):

• Registry roster / construction succeeds for the new kind, including via Registries::with_builtins().
• Build the plan and assert StepPlan::validate() is Ok (pure, no GPU); assert has_barostat_points() if applicable.
• The predicates return their intended values, and co-configuring an incompatible slot yields the matching ConfigError.
• Empty-state (particle_count == 0) construction and step are no-ops.
• If NVE-like, an energy-drift / determinism end-to-end test.

Completion checklist

• plan() is pure — no mutation of self, no branching on a per-step counter. The runner relies on a stable shape for graph capture.
• The warm-up force evaluation remains the runner’s job, before the loop — don’t prime forces in your plan.
• execute handles particle_count == 0; construction with zero particles must still succeed.
• graph_compatible reflects the execution path — any host-side scalar read between launches means false, or graph replay freezes stale values.
• The kind name is unique — a later same-name registration never shadows a built-in, so pick a fresh name.
• Physical parameters use dimension types, derive Convert, and implement convert_params; otherwise SI inputs may be interpreted as atomic units.

Adding a barostat

A barostat couples a system to a pressure bath by changing the simulation box and the particle or molecular positions within it. The user selects a separate barostat with [phase.barostat] kind = "..."; it can be combined with velocity Verlet and a thermostat to perform NPT dynamics.

Read the extension overview first for the shared builder, configuration, unit-conversion, and GPU terminology.

Choose the execution model

HeddleMD supports two kinds of barostat operation. Decide which model matches the physical algorithm before implementing the trait:

• Per-step (periodicity() == EveryStep, the default) — override apply, which runs every step at a BarostatPoint in the plan. It reads the current virial and kinetic energy, computes a scale factor, mutates the box, and rescales positions. Use berendsen_barostat.rs as the simpler equilibration-only example or c_rescale_barostat.rs as a canonical, stochastic example.
• Periodic (periodicity() == EveryNSteps(f)) — override apply_move, a host-orchestrated Monte-Carlo-style move run every f steps at a batch boundary. It receives mutable access to the force field because it must evaluate the energy of a trial box. Use mc_barostat.rs as the example.

Implement apply for a per-step barostat or apply_move for a periodic barostat. Do not put both algorithms into one implementation.

Per-step barostats

The runner calls apply(&mut self, buffers, sim_box, dt, timings) at the barostat point in every timestep. A typical implementation performs these operations:

1. Reduce the kinetic energy and the total virial into your own length-1 device scratch buffers (compute_kinetic_energy_on_device, compute_total_virial_on_device — both re-exported from crate::gpu). You do not read buffers.virials directly; the virial reduction sums the per-particle scalar-virial shares the force pipeline wrote.
2. Compute the scale factor µ and mutate the box lattice — in the existing barostats this is one scalar kernel that also writes diagnostics.
3. Rescale positions with the device-resident factor (rescale_positions_device_factor).

How the barostat receives a current virial

The per-particle virial shares are only current if the step’s ForceEval ran at AggregateLevel::ForcesAndScalars. You get this without the integrator having to know: when a per-step barostat is active, the runner sets runner_needs_scalars, which upgrades every ForceEval to ForcesAndScalars and disables the forces-only CUDA graph for the phase. So c-rescale/Berendsen work even with plain velocity-Verlet (which otherwise emits ForcesOnly).

Where the barostat appears in a timestep

apply runs only at a SubStep::BarostatPoint marker. velocity-verlet emits a terminal one, placed before the final AfterKick velocity projection so a RATTLE projection stays last. If a per-step barostat is configured with an integrator whose plan carries no BarostatPoint, phase setup fails with BarostatPlacementMissing.

Notify dependent calculations when the box changes

Mutating the lattice through sim_box.lattice_device_mut() bumps the box generation counter. The neighbor-list pipeline observes that counter and rebuilds on the next step. SPME then observes the neighbor-list rebuild generation and re-sorts its cached particle order. You get that for free.

Periodic barostats

The runner calls apply_move(&mut self, force_field, buffers, sim_box, constraint, dt, timings) at the requested host interval. A Monte Carlo volume move normally follows this sequence:

1. Read the current potential energy U_old (the runner guarantees the pre-move step evaluated ForcesAndScalars).
2. Snapshot positions, forces, and the lattice (device-to-device copies).
3. Draw the trial from Philox (host counter is fine here — the move is not captured in a CUDA graph), propose a volume change, guard the new box width against the cutoff, scale the box and the molecular centres of mass (mc_barostat_scale_molecule_com), and re-evaluate the energy with force_field.step(..., ForcesAndScalars) to get U_new.
4. Metropolis accept/reject; on reject, restore the snapshots.

Override init_run to upload the per-molecule tables once the molecule list is known. Because the move scales whole molecules’ centres of mass rather than individual atoms, it preserves rigid geometry and composes with constrained water. A periodic barostat needs no BarostatPoint marker and does no per-step work, so it stays graph_compatible (its move runs between graph replays).

Current limitation: the runner currently passes None for the constraint argument to apply_move. If your move needs to reproject constraints after scaling, that wiring does not exist yet — treat it as a limitation.

Preserve determinism in stochastic methods

For a stochastic barostat, draw from Philox keyed by a config seed. A per-step barostat must keep its counter on device (the c-rescale pattern: the µ kernel reads and increments a length-1 draw_counter_device buffer) so it survives CUDA-graph replay. A periodic barostat may use a host counter since its move isn’t captured. Same (seed, counter) → same draw → identical box trajectory.

Implement the Rust types

Create src/integrator/my_barostat.rs containing the typed parameters, runtime state, Barostat implementation, and builder. The following per-step skeleton is abbreviated and is not directly compilable; the comments stand in for the method-specific equations, launches, and error mapping.

#![allow(unused)]
fn main() {
use serde::Deserialize;
use crate::gpu::{GpuContext, GpuError, ParticleBuffers, c_rescale_compute_mu,
                 compute_kinetic_energy_on_device, compute_total_virial_on_device,
                 rescale_positions_device_factor};
use crate::io::config::ConfigError;
use crate::pbc::SimulationBox;
use crate::registry::{KindedBuilder, convert_params_in_place};
use crate::timings::{KernelStage, Timings};
use crate::precision::Real;
use super::{Barostat, BarostatBuilder, BarostatError};

#[derive(Debug, Clone, Deserialize, serde::Serialize, crate::units::Convert)]
#[serde(deny_unknown_fields)]
pub struct MyBarostatParams {
    pub pressure: crate::units::Pressure,
    pub tau: crate::units::Time,
    pub compressibility: crate::units::InversePressure,
    pub temperature: crate::units::Temperature, // stochastic only
    pub seed: u64,                              // stochastic only
}

#[derive(Debug)]
pub struct MyBarostat { /* params + length-1 device scratch (ke, virial, mu,
                           diagnostics, draw_counter) allocated in new() */ }

impl Barostat for MyBarostat {
    // periodicity() defaults to EveryStep — do NOT override for per-step.
    fn apply(&mut self, buffers: &mut ParticleBuffers, sim_box: &mut SimulationBox,
             dt: Real, timings: &mut Timings) -> Result<(), BarostatError> {
        if buffers.particle_count() == 0 { return Ok(()); }        // mandatory
        // reduce KE and virial → compute µ + mutate lattice + diagnostics
        //  → rescale_positions_device_factor(buffers, &self.mu_device)
        Ok(())
    }
    fn flush_pending_injection(&mut self, device: &std::sync::Arc<cudarc::driver::CudaDevice>)
        -> Result<(), BarostatError> { /* drain the injection accumulator for the log */ Ok(()) }
    fn log_column_names(&self) -> &'static [(&'static str, crate::units::Dimension)] {
        use crate::units::Dimension;
        &[("pressure", Dimension::Pressure), ("box_volume", Dimension::Dimensionless)]
    }
    fn log_column_values(&self, _ke: f64, _pe: f64) -> Vec<f64> { vec![/* pressure, volume */] }
}

#[derive(Debug, Clone)]
pub struct MyBarostatBuilder;
impl KindedBuilder for MyBarostatBuilder {
    fn kind_name(&self) -> &'static str { "my-barostat" }
    fn convert_params(&self, units: crate::units::UnitSystem, params: &mut toml::Value)
        -> Result<(), ConfigError> { convert_params_in_place::<MyBarostatParams>(units, params) }
}
impl BarostatBuilder for MyBarostatBuilder {
    fn validate_params(&self, params: &toml::Value) -> Result<(), ConfigError> { /* … */ Ok(()) }
    // graph_compatible defaults to true — keep it for a pure-launch per-step barostat.
    fn build(&self, gpu: &GpuContext, particle_count: usize, _n_constraints: usize,
             params: &toml::Value) -> Result<Box<dyn Barostat>, BarostatError> {
        Ok(Box::new(/* MyBarostat::new(...)? */))
    }
}
}

For a periodic barostat, instead override periodicity() → BarostatPeriodicity::EveryNSteps(freq), init_run, and apply_move (taking &mut ForceField); leave apply at the no-op default. See mc_barostat.rs.

Register and document the barostat

Add rqm/integration/my-barostat.md; follow rqm/integration/c-rescale-barostat.md (per-step) or rqm/integration/mc-barostat.md (periodic), and add a row to the barostat table in rqm/integration/framework.md.

Then edit the following existing files:

• src/integrator/mod.rs — declare pub mod my_barostat;, add the pub use my_barostat::{MyBarostat, MyBarostatBuilder}; re-export, and Box::new(MyBarostatBuilder) in the vec![...] of impl Builtins for dyn BarostatBuilder (roster at ~line 1230). This is the only registration point — dispatch is fully generic over dyn Barostat.

• src/gpu/device.rs — only if you add a new .cu: one field in the define_kernels! manifest.

No change needed to src/io/config.rs (validate/convert are called by registry lookup) or src/runner.rs (barostat dispatch is generic).

Test the barostat

Use both focused component tests and the shared slot-conformance harness. For a built-in barostat, add a SlotCase to builtin_barostat_cases() in tests/common/slot_conformance.rs. The case names the registered kind, uses Expect::HoldsDensity with scientifically justified bounds, and records why those bounds are appropriate for the finite constrained-water system.

Add a named test in tests/slot_conformance.rs that passes the case to run_barostat_case, and update the barostat sweep count in that file. Registry coverage fails if a built-in barostat has no case, while the sweep-count test fails if the case is not driven by a named test. The testing chapter explains the dense SPC/E system, negative controls, and the distinction between conformance bounds and bulk experimental accuracy.

tests/barostats_my.rs (model on tests/barostats_c_rescale.rs / tests/barostats_mc.rs):

• Registry exposes the kind (BarostatRegistry::with_builtins().builders() scan for kind_name()).
• validate_params accepts/rejects the expected params (deny_unknown_fields rejects extras).
• apply launches exactly its expected kernel set and none of the integrator kernels; µ matches the analytical formula for a known (K, W).
• Empty-state (particle_count == 0) apply is a no-op and does not bump the box generation.
• Determinism: same seed → identical box trajectory across builds.
• For periodic: the draw counter increments per attempted move, velocities are untouched by a move, and rigid-molecule geometry is invariant under a scale.

Completion checklist

• The execution model and method agree — apply (per-step) or apply_move (periodic), consistent with periodicity().
• A per-step barostat has a BarostatPoint in the integrator plan, or setup fails with BarostatPlacementMissing. velocity-verlet provides one.
• apply or apply_move handles particle_count == 0 and does not mutate the box in that case.
• The specification identifies the ensemble behavior. A canonical barostat should carry conserved-quantity bookkeeping in its log column; Berendsen is explicitly equilibration-only.
• A per-step stochastic barostat keeps its RNG counter on the device so graph replay stays deterministic; a host counter is fine for a periodic move.
• The kind name is unique. A later same-name registration never shadows a built-in.
• A built-in barostat has a SlotCase and named conformance test.

Adding a pair potential

A pair potential describes a short-range, non-bonded interaction evaluated for atom pairs in the neighbor list. Buckingham dispersion-repulsion is one example. Read the extension overview first for the shared Rust, configuration, unit-conversion, and determinism terminology.

Potentials differ from integrators and thermostats because a force field can contain several of them at once. HeddleMD therefore asks every registered PotentialBuilder whether the configuration contains interactions it owns. Each matching builder constructs one runtime component, and the force field combines all active components.

Choose the closest example

Use Lennard-Jones in src/forces/lj.rs as the primary example. A fast pair potential is JIT-composed: rather than launching its own kernel, it contributes a fragment of CUDA C++ source that the framework concatenates with every other active pair potential’s fragment into one nvrtc-compiled pair-force kernel. Your job is to (a) write a builder that activates on [[pair_interactions]] entries of your kind, and (b) write the per-pair functor source. JIT composition means that the framework combines active CUDA source fragments and compiles the resulting kernel at run time.

Connect configuration to the builder

• Activation is by presence. Your builder’s build(cx) scans cx.pair_interactions for entries whose kind matches yours; if there are none it returns Ok(None) and contributes nothing. No central match, no enum — adding a kind is entirely open.
• Parameters are routed by a claim. The builder returns a PotentialParamsClaim { category: PairInteraction, kind: "buckingham" }. The loader uses it to call your convert_params (unit conversion) and validate_params (domain checks) on each matching entry. An entry whose kind no builder claims is rejected at load with ConfigError::UnknownKind — so registering your builder is what turns "buckingham" from an error into a recognized kind.
• Common vs typed fields. [[pair_interactions]] entries carry only between, kind, and cutoff as centrally-parsed common fields; everything else (including r_switch) lives in the opaque params your typed struct deserializes.

Implement the one-pair calculation

The composition framework lives in src/forces/jit_composed.rs. A potential returns its CUDA C++ fragment as a Rust String; the composer inserts it into a shared CUDA translation unit. There is no separate kernels/pair_compute.cuh file to include. The example below shows the fixed functor interface; use src/forces/lj.rs for the complete Rust-side wiring.

pair_force_fragment() returns a PairForceFragment with a functor struct and three __device__ methods whose signatures are fixed by the composer’s call sites. Mirror LJ exactly:

struct BuckinghamPairFunctor {
    unsigned int n_types;
    const Real *type_a; const Real *type_rho; const Real *type_c6;
    const Real *type_cutoff; const Real *type_switch;
    const unsigned int *excl_offsets; const unsigned int *excl_partners;
    const Real *excl_scales;

    __device__ inline unsigned int slot(unsigned int ti, unsigned int tj) const {
        return ti * n_types + tj;
    }
    // factor = -(1/r) dU/dr   (so F_vec = factor * r_vec)
    // energy = U(r) ;  virial = factor * r2
    __device__ inline void evaluate(
        Real r2, Real inv_r, Real r, Real /*qi*/, Real /*qj*/,
        unsigned int i_type, unsigned int j_type,
        unsigned int /*i*/, unsigned int /*j*/,
        Real &factor, Real &energy, Real &virial) const
    {
        unsigned int p = slot(i_type, j_type);
        Real A = type_a[p], rho = type_rho[p], c6 = type_c6[p];
        Real inv_r6 = inv_r*inv_r * inv_r*inv_r * inv_r*inv_r;
        Real e_exp  = A * Real_exp(-r * (R(1.0) / rho));
        factor = e_exp * (R(1.0)/rho) * inv_r - R(6.0) * c6 * inv_r6 * inv_r*inv_r;
        energy = e_exp - c6 * inv_r6;
        virial = factor * r2;
    }
    // only needed for per-pair cutoffs (CutoffHandling::PerPair):
    __device__ inline Real cutoff_squared(unsigned int i_type, unsigned int j_type,
                                           unsigned int, unsigned int) const {
        Real c = type_cutoff[slot(i_type, j_type)]; return c * c;
    }
    __device__ inline Real exclusion_scale(unsigned int i, unsigned int j) const {
        return heddle_jit_exclusion_scale(i, j, excl_offsets, excl_partners, excl_scales);
    }
};

The generated kernel expects the following conventions:

• factor is already divided by r — it is -(1/r)·dU/dr, because the composer multiplies it by the raw displacement components, not the unit vector.
• The composer applies the cutoff mask and the exclusion scale outside evaluate — your exclusion_scale just does the table lookup; don’t pre-multiply it into your outputs.
• Use the precision shim from kernels/precision.cuh: Real, R(x) for literals, and the Real_* intrinsics (Real_exp, Real_sqrt, …) — never float/double/expf directly.
• Namespace everything. All fragments share one translation unit; give the functor struct and any free helpers a slot-unique prefix (BuckinghamPairFunctor, heddle_buck_*), and reuse shared preamble helpers like heddle_jit_exclusion_scale rather than redefining them.

Generate the argument interface from one schema

Do not hand-write entry_point_args or functor_init_source. Define one KernelArgSchema (via KernelArgSchema::pair_force(...)) and generate both from it — schema.entry_point_args() and schema.functor_init_source() — and route bind_pair_force_args through a KernelArgBinder built from the same schema. That makes the kernel parameter list, the functor-field initializers, and the launch-time argument binding share one source of truth, so they cannot drift. The functor struct’s field names must equal the schema’s functor_field strings; nvrtc catches mismatches.

Set consumes_type_index: true on the fragment so the composer supplies i_type/j_type (the per-atom type_indices buffer is a composer common argument — don’t bind it yourself), and cutoff: to CutoffHandling::Uniform when every pair shares one cutoff, else PerPair.

Rely on deterministic framework accumulation

You get bit-exact summation for free: the composer converts each per-pair (factor, energy, virial) to integer fixed-point and atomicAdds into per-atom u64 accumulators, and integer addition is order-independent. Your one obligation is that evaluate be a pure per-pair function — compute one triple from r and the per-type params, with no state carried across pairs and no float accumulation of your own. All summation belongs to the framework.

Implement the Rust types

Create src/forces/buckingham.rs. It should follow lj.rs and contain the kind and label constants, typed parameters, resolved state, argument schema, potential traits, builder, and CUDA fragment. The following code is an abbreviated outline rather than a directly compilable example.

#![allow(unused)]
fn main() {
pub const BUCKINGHAM_KIND: &str = "buckingham";
const LABEL: &str = "buckingham";

#[derive(Debug, Clone, serde::Deserialize, serde::Serialize, crate::units::Convert)]
#[serde(deny_unknown_fields)]
pub struct BuckinghamPairParams {
    pub a: crate::units::Energy,
    pub rho: crate::units::Length,
    pub c6: f64,   // energy·length⁶ has no dimension newtype — see the caveat below
    #[serde(default)] pub r_switch: Option<crate::units::Length>,
}

#[derive(Debug)]
pub struct BuckinghamState { /* device param table + Arc<DeviceExclusionList> clone + cutoffs */ }

impl Potential for BuckinghamState {
    fn label(&self) -> &'static str { LABEL }
    fn max_cutoff(&self) -> Option<Real> { Some(self.max_cutoff) }   // MUST be Some
    // frequency_class() defaults to Fast — leave it.
    fn compute(&mut self, /* … */) -> Result<(), ForceFieldError> {
        unreachable!("JIT-composed pair slot; compute() is never called")
    }
    fn jit_participant(&self) -> Option<JitParticipant<'_>> {
        Some(JitParticipant::PairForce(self))
    }
}

impl PairForcePotential for BuckinghamState {
    fn pair_force_fragment(&self) -> PairForceFragment { /* build from arg_schema() */ }
    fn bind_pair_force_args(&self, _ctx: &PairForceBindContext<'_>, builder: &mut ForceLaunchBuilder) {
        // KernelArgBinder::new(&buck_arg_schema(), LABEL, builder); push param + exclusion buffers
    }
}

#[derive(Debug, Clone)]
pub struct BuckinghamBuilder;
impl PotentialBuilder for BuckinghamBuilder {
    fn build(&self, cx: &PotentialBuildContext<'_>) -> Result<Option<Box<dyn Potential>>, ForceFieldError> {
        let pairs: Vec<_> = cx.pair_interactions.iter().filter_map(resolve_buckingham_pair).collect();
        if pairs.is_empty() { return Ok(None); }               // activation
        Ok(Some(Box::new(/* BuckinghamState::from_config(...) */)))
    }
    fn params_claim(&self) -> Option<PotentialParamsClaim> {
        Some(PotentialParamsClaim {
            category: PotentialParamsCategory::PairInteraction, kind: BUCKINGHAM_KIND })
    }
    fn validate_params(&self, entry: PotentialConfigEntry<'_>) -> Result<(), ConfigError> {
        let PotentialConfigEntry::PairInteraction(p) = entry else { unreachable!() };
        // A/ρ/C₆ finite, ρ > 0, cross-field r_switch <= p.cutoff …
        Ok(())
    }
    fn convert_params(&self, units: crate::units::UnitSystem, params: &mut toml::Value)
        -> Result<(), ConfigError> {
        crate::registry::convert_params_in_place::<BuckinghamPairParams>(units, params)
    }
}
}

The per-type parameter table is SoA (one CudaSlice<Real> per field, row-major ti*n_types+tj, uploaded with htod_or_empty) — copy LennardJonesParameterTable in src/gpu/kernels.rs, or keep it local to buckingham.rs.

Register and document the potential

Add rqm/forces/buckingham-pair-force.md; follow rqm/forces/lj-pair-force.md and cross-reference rqm/forces/jit-composed-pair-force.md and rqm/forces/packed-neighbour-pair-force.md rather than restating the shared kernel contract.

Then edit the following existing files:

• src/forces/mod.rs — three lines: pub mod buckingham;, the pub use buckingham::{...} re-export, and Box::new(BuckinghamBuilder) in the vec![...] of impl Builtins for dyn PotentialBuilder (roster at ~line 386). Registration order is the slot evaluation order; since your kind is unique in its category, position only affects the slot index — next to LennardJonesBuilder is natural.

• src/gpu/kernels.rs / src/gpu/mod.rs — only if you add the param table there rather than locally.

No change needed to src/io/config.rs (PairInteractionConfig is open-shaped and routing is claim-driven), build.rs, or any .cu file (the functor is a runtime-compiled string; there is no new kernel).

Test the potential

Add a ConsistencyFixture::pair(...) entry to builtin_consistency_fixtures() in tests/common/consistency.rs. Use the potential’s fragment label, register only the builder under test, and provide sample distances in the physically important and switched regions. Include at least one analytical ReferencePoint: finite-difference consistency alone cannot detect an energy and force that are both wrong by the same factor.

The shared harness checks force–energy, Newton, virial, reference-point, and cutoff invariants through the real GPU force pipeline. assert_fixture_coverage compares fixture labels with the built-in fragment labels, so a built-in pair fragment without a fixture fails the suite. See Testing Extension Components for the fixture API, sample selection, tolerance guidance, and negative controls.

• In-file #[cfg(test)] — pin the exact entry_point_args / functor_init_source strings the schema emits, and assert the fragment declares consumes_type_index and uses slot(i_type, j_type). These run without a GPU.
• tests/potential_claims.rs — add the builder’s claim to the claim-coverage test and a build/validate/convert round-trip.
• tests/forces_framework.rs — a Buckingham-only force field that steps and writes non-zero forces.
• tests/io_config.rs — a kind = "buckingham" parse test.

Completion checklist

• compute() is bypassed for JIT pair slots — make it unreachable!() and put all physics in the fragment.
• max_cutoff() returns Some because it feeds the shared neighbor-list cutoff and the JIT prune constant; None drops your interaction or panics.
• The shared neighbor and exclusion lists are used — clone cx.device_exclusions (an Arc); never allocate your own DeviceExclusionList. Bind only your per-type param buffers and the exclusion buffers.
• The units of every parameter are documented. A (Energy) and ρ (Length) convert cleanly, but C₆ carries energy·length⁶, which has no dimension newtype. Either keep it a plain f64 (a Convert no-op) and document that it must be supplied in atomic units regardless of the file’s units selector, or add a new Dimension variant and newtype in src/units/mod.rs.
• Tests use the implemented fixed-point scale of 2^48 (some spec prose says 2^32); the code is authoritative, but you never touch this directly.
• A built-in fragment has a consistency fixture with analytical reference data.

Adding a bonded potential

A bonded potential describes an intramolecular interaction whose atom indices come explicitly from the topology: two atoms for a bond, three for an angle, or four for a dihedral. It does not search the neighbor list.

Read the extension overview and Adding a pair potential first. Bonded and pair potentials use the same builder, parameter-claim, and run-time CUDA composition machinery. This guide concentrates on the different topology and reduction contracts.

Choose the closest example

This page uses a "cubic" bond as its running example. Begin with src/forces/harmonic_bond.rs for the simplest complete pattern or src/forces/morse.rs for another two-particle potential. Angle and dihedral implementations follow the same workflow but use different functor signatures; the table near the end identifies their examples.

Understand the bonded execution path

The following pieces are the same as for a pair potential: a cloneable builder implementing PotentialBuilder with a params_claim, validate_params, convert_params, and build; a State implementing Potential; a CUDA functor emitted as a Rust source string and compiled by nvrtc; registration by adding one line to the potential Builtins roster.

The important differences are:

• The claim category is BondType (or AngleType / DihedralType), and it matches the kind field of [[bond_types]] entries.
• max_cutoff() returns None — bonded potentials don’t consume the neighbor list; they iterate an explicit per-bond index table.
• compute() is NOT bypassed. For pair slots the framework skips compute (the JIT kernel does everything). For bonded slots, the composed kernel writes each bond’s contribution into a per-bond scratch buffer, and then compute() runs the reduction step (reduce_bond_forces) that scatters those per-bond contributions into per-atom forces. So your compute must run exactly that reduction and nothing else — doing more double-applies forces.
• No new .cu file is needed. The reduction kernel reduce_bond_forces already exists and is re-exported from crate::gpu; the per-bond physics is the JIT source string.

Map topology entries to the potential

The user declares [[bond_types]] name = "CT-CT" kind = "cubic" with a params inline table, and references it by name from the .topology file’s [bonds] section (atom_i atom_j <bond_type_name>). The loader resolves the name to a global bond_type_index, canonicalizes each bond i < j, and sorts the bond list — that fixed order is the determinism anchor.

Your build(cx) filters cx.bond_list to the bonds whose type uses your kind (BondList::filter_by_type_index), which rebuilds the per-atom reduction map (atom_bond_offsets / atom_bond_indices) over just that subset. Parameter tables stay indexed by the global bond_type_index (rows for other kinds hold zero placeholders your slot never reads), because the index the kernel reads out of the bond table is the global one.

Two things you get for free: 1-2 exclusions (every bond contributes an implicit non-bonded exclusion, derived kind-agnostically from the full bond list), and — for dihedrals — 1-4 scaling from the dihedral type’s common scale_lj_14 / scale_coul_14 fields.

Implement the one-bond calculation

The bond functor computes one bond’s contribution from the scalar distance. The composer’s generated entry point computes the minimum-image displacement (dx, dy, dz) = r_i − r_j, r2, and r = sqrt(r2), then calls your evaluate with the fixed signature:

struct CubicBondFunctor {
    const Real *bond_k; const Real *bond_r0; const Real *bond_kcub;
    __device__ inline void evaluate(
        Real r2, Real r,
        unsigned int bond_type_index,
        Real dx, Real dy, Real dz,
        Real &fmag, Real &u_k, Real &w_k) const
    {
        if (r < R(1.0e-7)) { fmag = R(0.0); u_k = R(0.0); w_k = R(0.0); return; }
        Real k = bond_k[bond_type_index], r0 = bond_r0[bond_type_index],
             kcub = bond_kcub[bond_type_index];
        Real dr = r - r0;
        // U = ½ k dr² + ⅓ k kcub dr³ ;  dU/dr = k dr (1 + kcub dr)
        fmag = -k * dr * (R(1.0) + kcub * dr) / r;   // per-component factor (÷ r)
        u_k  = R(0.5)*k*dr*dr + (R(1.0)/R(3.0))*k*kcub*dr*dr*dr;  // FULL bond energy
        w_k  = fmag * r2;                             // FULL scalar virial
    }
};

The generated entry point expects these conventions:

• fmag is −(dU/dr)/r (already divided by r). The entry point writes force fmag·(dx,dy,dz) onto atom i and its negation onto atom j (Newton’s third law), into per-bond scratch slots 2k and 2k+1.
• u_k and w_k are the FULL bond values — the entry point applies the ½-each split across the two atoms. Don’t pre-halve them.
• Use the precision shim (Real, R(x), Real_*) and a unique struct/helper prefix, exactly as for pair fragments.

Generate the argument interface from one schema

As with pair potentials, generate entry_point_args and functor_init_source from one KernelArgSchema (built with KernelArgSchema::intramolecular(...)) and route bind_bonded_force_args through a KernelArgBinder on the same schema.

Preserve the fixed reduction order

Each bond is processed by one thread that writes only its own two scratch slots — no cross-thread contention. The reduction then sums each atom’s bonds in the fixed index order set by the sorted bond list, so the per-atom float sum is bit-identical run to run. Bonded potentials therefore need no fixed-point accumulator (unlike the pair path). Keep evaluate a pure per-bond function.

Implement the Rust types

Create src/forces/cubic_bond.rs and follow harmonic_bond.rs for the complete implementation. The following code is an abbreviated outline rather than a directly compilable example.

#![allow(unused)]
fn main() {
pub const CUBIC_BOND_KIND: &str = "cubic";
const LABEL: &str = "cubic_bond";

#[derive(Debug, Clone, serde::Deserialize, serde::Serialize, crate::units::Convert)]
#[serde(deny_unknown_fields)]
pub struct CubicBondParams {
    pub k: crate::units::Stiffness,        // energy / length²
    pub r0: crate::units::Length,
    pub kcub: crate::units::InverseLength,
}

#[derive(Debug)]
pub struct CubicBondState {
    /* device: bonds [i,j,type_idx]; atom_bond_offsets/indices over the filtered
       subset; bond_k/r0/kcub param tables sized to bond_types.len();
       per-bond scratch bond_pair_{x,y,z,energy,virial} of length 2*bond_count */
}

impl Potential for CubicBondState {
    fn label(&self) -> &'static str { LABEL }
    fn max_cutoff(&self) -> Option<Real> { None }          // bonded → no neighbor cutoff
    // frequency_class() defaults to Fast.
    fn compute(&mut self, /* … */ level: AggregateLevel) -> Result<(), ForceFieldError> {
        // if particle_count == 0 || bond_count == 0 { return Ok(()); }
        // reduce_bond_forces(..., write_scalars = matches!(level, ForcesAndScalars))
        Ok(())
    }
    fn jit_participant(&self) -> Option<JitParticipant<'_>> {
        Some(JitParticipant::Bonded(self))
    }
}

impl BondedPotential for CubicBondState {
    fn bonded_force_fragment(&self) -> BondedForceFragment { cubic_bonded_force_fragment() }
    fn bonded_scratch(&self) -> BondedScratchView<'_> { /* bonds + bond_pair_* + bond_count */ }
    fn bind_bonded_force_args(&self, _ctx: &ForceLaunchContext<'_>, builder: &mut ForceLaunchBuilder) {
        // KernelArgBinder::new(&cubic_arg_schema(), LABEL, builder); push bond_k/r0/kcub
    }
}

#[derive(Debug, Clone)]
pub struct CubicBondBuilder;
impl PotentialBuilder for CubicBondBuilder {
    fn build(&self, cx: &PotentialBuildContext<'_>) -> Result<Option<Box<dyn Potential>>, ForceFieldError> {
        let has = cx.bond_list.bonds.iter()
            .any(|b| cx.bond_types.get(b.bond_type_index as usize)
                       .is_some_and(|t| t.kind == CUBIC_BOND_KIND));
        if !has { return Ok(None); }                       // activation
        Ok(Some(Box::new(/* CubicBondState::new(cx.gpu, cx.bond_list, cx.bond_types)? */)))
    }
    fn params_claim(&self) -> Option<PotentialParamsClaim> {
        Some(PotentialParamsClaim {
            category: PotentialParamsCategory::BondType, kind: CUBIC_BOND_KIND })
    }
    fn validate_params(&self, entry: PotentialConfigEntry<'_>) -> Result<(), ConfigError> {
        let PotentialConfigEntry::BondType(bt) = entry else { unreachable!() };
        // deserialize bt.params into CubicBondParams; k, r0 finite & positive …
        Ok(())
    }
    fn convert_params(&self, units: crate::units::UnitSystem, params: &mut toml::Value)
        -> Result<(), ConfigError> {
        crate::registry::convert_params_in_place::<CubicBondParams>(units, params)
    }
}
}

Register and document the potential

Add rqm/forces/cubic-bond.md; follow rqm/forces/harmonic-bond.md and cross-reference rqm/forces/jit-composed-intramolecular.md.

Then edit the following existing files:

• src/forces/mod.rs — three lines: pub mod cubic_bond;, the pub use cubic_bond::{...} re-export, and Box::new(CubicBondBuilder) in the vec![...] of impl Builtins for dyn PotentialBuilder (roster at ~line 386), conventionally next to the other bonded builders. Because your kind is unique within the BondType category, position only affects the slot index.

No change needed to src/io/config.rs (BondTypeConfig is open-shaped; routing is claim-driven — an unclaimed kind is UnknownKind), build.rs, src/gpu/device.rs, or any .cu file. PotentialParamsCategory::BondType already exists.

Test the potential

Add the potential to builtin_consistency_fixtures() in tests/common/consistency.rs using ConsistencyFixture::bond, angle, or dihedral as appropriate. Use the fragment label, register only the builder under test, and choose sample coordinates that exercise nonzero forces and the important regions of the potential. Add analytical ReferencePoint values so a consistently mis-scaled energy and force cannot pass unnoticed.

assert_fixture_coverage compares the fixture labels with all registered built-in fragments. Registering a built-in bonded fragment without the corresponding fixture therefore fails the suite. Testing Extension Components explains the shape-specific geometries, invariant checks, tolerances, and CPU negative controls.

• In-file #[cfg(test)] — pin the exact entry_point_args / functor_init_source strings from cubic_arg_schema(); GPU-free.
• tests/forces_cubic_bond.rs (model on tests/forces_harmonic_bond.rs) — a two-particle force field: assert per-atom force/energy/virial against the analytic U(r) and a finite-difference of the energy, check F_i = −F_j, and the empty-bond-list no-op.
• tests/io_config.rs / tests/potential_claims.rs — kind = "cubic" validates; an unknown kind yields ConfigError::UnknownKind { slot: "bond_types" }.

Adapt the pattern for angles and dihedrals

Same recipe; swap the capability trait, claim category, fragment type, scratch view, and reduction kernel:

The angle functor returns forces on atoms i and k directly (the central atom j is inferred as −(F_i + F_k)); the dihedral functor takes three bond displacements and returns all four forces. Dihedral types additionally carry the scale_lj_14 / scale_coul_14 common fields, consumed centrally by the topology loader to derive scaled 1-4 exclusions — a new dihedral kind gets that automatically.

Completion checklist

• compute() runs the reduction, and only the reduction. It is not bypassed like a pair slot. Anything extra double-applies forces.
• max_cutoff() returns None because a non-None value inflates the shared neighbor-list cutoff for no reason.
• fmag is −(dU/dr)/r (pre-divided); u_k / w_k are full values (the entry point splits them). These two are the most common mistakes.
• Parameter tables are sized to the full bond_types.len() and indexed by the global type index — don’t compact them to the filtered subset.
• An empty bond list is a no-op: build returns Ok(None) when no bond uses your kind, and compute early-returns on bond_count == 0.
• The kind name is unique. A later same-name registration never shadows a built-in.
• A built-in fragment has the matching bond, angle, or dihedral consistency fixture.