Run metadynamics with PLUMED¶

mlmd.run_md can add a PLUMED bias to an ASE-compatible potential. The built-in configuration generates well-tempered metadynamics for one distance, angle, or torsion; custom PLUMED input supports additional actions and multiple collective variables (CVs). This tutorial covers configuration, native execution, outputs, continuous equilibration, and restart.

Metadynamics deposits history-dependent hills along selected CVs. Choosing a CV, Gaussian width, deposition interval, and bias strength requires physical knowledge of the process being sampled; faster motion in a chosen CV alone does not establish a converged free-energy surface. See the PLUMED METAD documentation for the method and action options.

1. Install native PLUMED¶

Use the repository's environment from its root directory:

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

It installs the MLMD dependencies plus the native PLUMED kernel, Python bindings, command-line tools, and FFTW. pip install 'matmalib[metadynamics]' supplies ASE and Python bindings for a supplied calculator, but the bindings alone do not include the native kernel. For AIMNet, use the full MLMD setup in the simulation tutorial.

Browser cells on this page only create configurations and inspect units. Running a simulation needs the local environment; PLUMED and AIMNet are not loaded in the browser.

2. Understand indices, units, and timing¶

The two interfaces intentionally use different conventions:

Quantity	MD/geometry interface	Generated PLUMED input
Atom selections	`fixed_atoms` uses zero-based positions	`MetadynamicsConfig.atoms` uses one-based positions
Length	Å	nm
Energy	eV	kJ/mol
Time	fs	ps in HILLS/COLVAR
Angles	Geometry methods commonly return degrees	`ANGLE`, `TORSION`, and their `sigma` use radians
Temperature	K	K
Frequency	Integration steps	`pace` and `colvar_stride` are integration-step counts

PLUMED's default units apply to input parameters and output files; custom input can change them with an initial UNITS directive. See PLUMED UNITS.

from matmalib import mlmd

meta = mlmd.MetadynamicsConfig(atoms=(3, 11), sigma=0.02, pace=300,
                               colvar_stride=10)
md = mlmd.MDConfig(steps=100000, timestep_fs=0.5,
                   mode="metadynamics", metadynamics=meta)
print("Selected zero-based positions:", tuple(i - 1 for i in meta.atoms))
print("Distance width:", 10 * meta.sigma, "angstrom")
print("Hill interval:", meta.pace * md.timestep_fs, "fs")
print("COLVAR interval:", meta.colvar_stride * md.timestep_fs / 1000, "ps")
print("Run duration:", md.steps * md.timestep_fs / 1000, "ps")

The example reports a 0.2 Å width, 150 fs hill interval, 0.005 ps COLVAR interval, and 50 ps duration. These are illustrative settings, not a validated protocol for a particular molecule. Changing the timestep changes every frequency expressed in physical time unless you adjust the corresponding step counts.

Use positions in the input atom ordering, not Atom.id, PDB serial numbers, or residue identifiers. Atom selections must contain the required number of unique, valid one-based indices.

3. Run a small installation check¶

Save the following as check_metadynamics.py and run python check_metadynamics.py in a new directory. This exercises native PLUMED with ASE's inexpensive EMT calculator and a small copper cell. It is an integration check, not a machine-learning model or a physical free-energy calculation. It avoids downloading model weights.

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

out = Path("metad_check")
out.mkdir(exist_ok=True)
simulation = mlmd(calculator=EMT())
meta = mlmd.MetadynamicsConfig(
    cv="DISTANCE", atoms=(1, 2),
    sigma=0.02, height_kj_mol=0.1, pace=10,
    biasfactor=10.0, colvar_stride=1,
    calc_rct=True, grid_auto=True,
    hills_file=out / "HILLS",
    colvar_file=out / "COLVAR",
    save_input=out / "plumed.dat",
    log_file=out / "plumed.log",
)
trajectory = simulation.run_md(
    bulk("Cu", cubic=True),
    mlmd.MDConfig(
        steps=100, timestep_fs=0.2, temperature_K=300.0,
        ensemble="NVT", thermostat="langevin", seed=17,
        mode="metadynamics", metadynamics=meta,
        trajectory=out / "trajectory.xyz",
        trajectory_interval=10,
        log_file=out / "md.log",
        log_interval=10,
    ),
)
print("Trajectory:", trajectory)
print("Checkpoint:", simulation.data.checkpoint_file)
print("Final step:", simulation.data.last_step)
print("Resolved PLUMED:")
print("\n".join(simulation.data.plumed_input))

Expect a streamed XYZ trajectory through step 100, an MD log, a final checkpoint, and PLUMED input/log/HILLS/COLVAR files. The simulation lasts only 20 fs. Re-running with the same names raises instead of silently replacing files; select a new output directory for another check.

For a molecular MLMD calculation, use simulation = mlmd(...) with appropriate AIMNet weights or another validated ASE calculator, and replace the copper structure with your equilibrated State or extended XYZ path. Choose the CV atoms for that system. The solvation tutorial creates equilibrated.xyz; retain its cell and PBC metadata for a periodic run.

4. Choose a generated CV and bias¶

Setting	Meaning
`cv="DISTANCE"`, `atoms=(i, j)`	One scalar separation; `sigma` and wall positions in nm.
`cv="ANGLE"`, `atoms=(i, j, k)`	Angle at the middle atom; width and wall positions in radians.
`cv="TORSION"`, `atoms=(i, j, k, l)`	Signed dihedral using the supplied ordering; width and wall positions in radians.
`height_kj_mol`	Positive initial hill height.
`pace`	Positive integer number of integration steps between hills.
`biasfactor`	Well-tempered factor, strictly greater than one.
`temperature_K`	Bias temperature; `None` uses the enclosing MD temperature.
`cv_label`, `metad_label`, `wall_label`	Distinct PLUMED action labels; defaults `cv`, `metad`, `wall`.
`metad_arg`, `wall_arg`	Default to the CV label; references must correspond to available actions.
`colvar_stride`, `colvar_arg`	PRINT frequency and optional comma-separated output references.

For nonperiodic systems, generated CV input includes NOPBC. For periodic systems PLUMED uses the box supplied at each step. Partial periodicity is rejected: use either no periodic axes or full three-dimensional PBC. If your CV needs molecule reconstruction or a special treatment of periodic images, express it in custom PLUMED input rather than assuming wrapped coordinates are sufficient.

The generated form defines one scalar CV. Use custom input for a multidimensional bias, combinations of CVs, or other PLUMED actions. metad_arg changes references; it does not create additional CV definitions.

Add an optional wall¶

from matmalib import mlmd

meta = mlmd.MetadynamicsConfig(
    cv="DISTANCE", atoms=(1, 2),
    wall_type="UPPER_WALLS",
    wall_at=0.6, wall_kappa=150000.0,
    wall_exp=2, wall_eps=1.0, wall_offset=0.0,
)
print("Wall location:", meta.wall_at * 10, "angstrom")
print("Wall action:", meta.wall_type)
print("Energy scale:", meta.wall_kappa, "in generated PLUMED energy units")

UPPER_WALLS or LOWER_WALLS enables a wall; wall_type=None disables it. wall_at/wall_offset follow the CV units. wall_kappa, wall_exp, and wall_eps control the PLUMED wall potential. A strong wall can produce large forces and restrict accessible states; choose its location and stiffness deliberately. The resulting wall.bias is printed by default and matters when interpreting or removing bias in analysis.

Decide whether to use a grid and normalized bias¶

calc_rct=True requests PLUMED's normalization estimate and adds metad.rbias and metad.rct to the default COLVAR columns. It requires a bias grid. Those columns enable normalized population reweighting in the analysis tutorial; raw metad.bias alone is not used as a normalized population weight.

Situation	Configuration
Fully periodic, fixed-cell scalar DISTANCE	Defaults `grid_auto=None, calc_rct=None` select an automatic grid and normalized bias when no explicit grid is supplied.
Same system, explicit bounds	Set `grid_min`, `grid_max`, `grid_bins`, and `grid_auto=False`; use `calc_rct=True` for normalized bias.
Nonperiodic distance, angle, or torsion	No automatic grid. Supply an appropriate complete grid for normalized bias, or use `calc_rct=False`.
NPT distance	Automatic bounds are disabled because the cell evolves. Use explicit bounds with `grid_auto=False`, or disable normalization and analyze HILLS.

Automatic distance bounds are derived from the fixed cell's half-cell corners, converted to nm, and extended by a margin related to sigma. grid_spacing defaults to 0.002 nm and controls the resulting bin count. Automatic grid_min must be finite and at most zero. For explicit grids, both grid_max and a positive integer grid_bins are required, with grid_min < grid_max; setting calc_rct=True without a complete grid raises.

For example, an explicit distance grid could be configured as follows after choosing bounds that contain your sampled distances and hill support:

meta = mlmd.MetadynamicsConfig(
    atoms=(3, 11), sigma=0.02,
    grid_auto=False, grid_min=0.0, grid_max=3.0, grid_bins=1500,
    calc_rct=True,
    hills_file="npt_HILLS", colvar_file="npt_COLVAR",
    save_input="npt_plumed.dat",
)

The 3 nm bound is an example, not an automatic guarantee for an expanding NPT cell. For a torsion, use a domain consistent with its angular periodicity and the native PLUMED action; do not reuse a distance grid. Matmalib's FES analysis is currently one-dimensional and does not support a basin crossing a periodic boundary.

NPT additionally requires full PBC, a finite-volume cell, a stress-capable potential, and no fixed-atom, SHAKE, or center-of-mass constraints. The PLUMED virial contributes to pressure/stress, and the bias is updated once per integration step even if the integrator requests several force evaluations. A time-dependent bias also means that choosing ensemble="NVE" does not yield conservation of the unbiased system's energy.

5. Supply custom PLUMED input¶

Use either plumed_file or input_lines. These replace generated CV/bias/wall/PRINT commands; fields such as atoms, sigma, hills_file, and colvar_file no longer define the custom actions or their destinations.

This configuration defines two distances and a two-dimensional bias. It assumes at least four atoms; pair choices and hill widths must suit your system:

from pathlib import Path
from matmalib import mlmd

out = Path("custom_bias").resolve()
out.mkdir(exist_ok=True)
meta = mlmd.MetadynamicsConfig(
    input_lines=(
        "UNITS LENGTH=nm ENERGY=kj/mol TIME=ps",
        "d1: DISTANCE ATOMS=1,2",
        "d2: DISTANCE ATOMS=3,4",
        f"m: METAD ARG=d1,d2 SIGMA=0.02,0.02 HEIGHT=1 PACE=300 "
        f"BIASFACTOR=10 TEMP=300 FILE={out / 'HILLS'}",
        f"PRINT ARG=d1,d2,m.bias STRIDE=10 FILE={out / 'COLVAR'}",
    ),
    save_input=out / "plumed.dat",
    log_file=out / "plumed.log",
)
# In your prepared workflow:
# simulation.run_md(ready, mlmd.MDConfig(
#     mode="metadynamics", metadynamics=meta,
#     trajectory=out / "trajectory.xyz", log_file=out / "md.log",
# ))

Alternatively, mlmd.MetadynamicsConfig(plumed_file="my_plumed.dat", save_input="resolved.dat") reads an existing file. Comments and blank lines are stripped, and PLUMED ... blocks are flattened; nested or unterminated blocks are rejected. save_input stores that resolved input, also available as simulation.data.plumed_input.

Custom commands control their own paths, units, periodic treatment, and restart directives. Relative output paths are interpreted by the native run, not automatically relative to plumed_file; use explicit paths for reproducibility. Matmalib's generated-output path protections do not inspect every custom FILE action, so use a fresh directory and avoid accidental reuse. The one-dimensional analysis API cannot directly reconstruct this two-dimensional example; use PLUMED tools or a separately prepared one-dimensional marginal grid.

6. Keep a live bias through equilibration¶

Ordinary preparation is unbiased unless its own configuration explicitly requests metadynamics. To build the bias during equilibration and continue it during production, use one shared MetadynamicsConfig and continue_equilibration_bias=True in a single run_md call:

from matmalib import mlmd

simulation = mlmd()
meta = mlmd.MetadynamicsConfig(
    atoms=(3, 11),
    hills_file="continuous_HILLS", colvar_file="continuous_COLVAR",
    save_input="continuous_plumed.dat", log_file="continuous_plumed.log",
)
equilibration = mlmd.EquilibrationConfig(
    steps=2000, timestep_fs=0.5, temperature_K=300,
    ensemble="NVT", thermostat="langevin", thermostat_time_fs=10,
    mode="metadynamics", metadynamics=meta,
    trajectory="biased_equilibration.xyz", log_file="biased_equilibration.log",
    seed=31,
)
simulation.run_md("equilibrated.xyz", mlmd.MDConfig(
    steps=100000, timestep_fs=0.5, temperature_K=300,
    ensemble="NVT", thermostat="langevin", thermostat_time_fs=100,
    mode="metadynamics", metadynamics=meta,
    equilibration=equilibration, continue_equilibration_bias=True,
    trajectory="biased_production.xyz", log_file="biased_production.log",
    seed=31,
))
print(simulation.data.start_step, simulation.data.last_step)  # 2000, 102000

The live PLUMED instance, hill history, velocities, RNG, and clock continue across the boundary. HILLS/COLVAR remain shared; MD trajectories and logs have separate paths. Thermostat damping may change. Both phases must use identical timestep, target temperature, ensemble, metadynamics settings, and potential; cell/PBC/charge/multiplicity overrides must agree. Seeds must match when both are set. NPT retains the cell reached during equilibration.

This also supports custom PLUMED input while the process remains alive. It is different from ending equilibrate and then starting an independent run_md: returning a State preserves molecular state but does not preserve a live PLUMED bias. Independent biased phases should use separate bias files unless you implement an explicit native restart workflow.

7. Restart generated metadynamics safely¶

By default, checkpoints are written every 1,000 steps and at completion, including biased runs. checkpoint_interval changes the interval and checkpoint_file=None disables them. Save simulation.data.checkpoint_file before performing another operation that replaces .data.

After the installation check above, run:

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

simulation = mlmd(calculator=EMT())  # same potential as the original check
continued = simulation.restart_md(
    "metad_check/trajectory.checkpoint.npz",
    steps=50,
    trajectory="metad_check/part2.xyz",
    log_file="metad_check/part2.log",
)
print(continued)
print(simulation.data.last_step)  # 150

steps means additional integration steps. Restart restores molecular state, constraints, random-generator state, thermostat/barostat variables, and the original step count. Use the same ASE version and model/calculator settings; model weights are not embedded in a checkpoint.

The new XYZ/log files are separate segments; existing files are protected. For generated PLUMED input, saved HILLS/COLVAR snapshots are copied to part2.plumed/, the saved input is updated to those paths, and native restart continues the bias and output history there. Original files are untouched. HILLS text precision can produce small numerical differences from an uninterrupted biased run.

overwrite=True is permission to replace output paths during a new run; it is not a continuation or append switch. Use restart_md for checkpoint continuation or continue_equilibration_bias for live phase continuation. Custom PLUMED runs still write molecular-state checkpoints, but restart_md rejects them because arbitrary actions may need additional state beyond HILLS/COLVAR. Handle custom native restart requirements explicitly rather than treating their checkpoint as a complete automatic bias restart.

8. Inspect the result before analysis¶

Output	Contents
Extended XYZ trajectory	Saved geometry, cell/PBC, masses, momenta, `step`, and `time_fs`; streamed to disk.
MD log	Potential/kinetic/total energy, temperature, maximum force, volume, and timing; unavailable model-disagreement fields are `nan`.
HILLS	Deposited-hill history, used by native bias restart and reconstruction.
COLVAR	CV and available bias columns at `colvar_stride`; generated times are in ps.
Saved PLUMED input	Resolved action definitions and paths for provenance and analysis.
`simulation.data`	Status, bounded running statistics, model/system/settings metadata, resolved input, and checkpoint paths.

In biased runs, potential and total energies in MD diagnostics include the applied bias; they are not an unbiased energy-conservation trace. Inspect CV coverage, numerical stability, hill history, and sensitivity to discarded equilibration and sampling duration before interpreting a FES. Continue to metadynamics analysis for normalized COLVAR reweighting, HILLS reconstruction, and one-dimensional free-energy plots.