Analyze metadynamics and plot a 1D free-energy profile¶

Use mlmd.analyze_metadynamics() to reconstruct a one-dimensional free-energy surface (FES), compare basin populations, inspect sampling diagnostics, and export results. Plot.fes_1d() plots the resulting Data object. Analysis does not require a model, GPU, or ASE calculator.

The runnable examples below use synthetic data to demonstrate the algorithms. They do not represent a converged molecular simulation or a measured reaction rate. Each example creates its own small input files and runs independently in the browser. Reconstructing a profile from HILLS requires the PLUMED command-line executable on your computer; existing FES grids and normalized COLVAR analysis run in the browser.

Choose the input that contains the information you need¶

Input	`AnalysisConfig` fields	What is calculated
Existing one-dimensional FES	`fes_file`, `cv`, `energy_column`	Load the grid, convert energies, set its minimum to zero, and analyze basins
Regularly sampled COLVAR with normalized bias	`colvar`, `rbias_column`, or `bias_column` together with `rct_column`	A weighted histogram and relative `-RT ln(probability density)`
HILLS from a supported one-dimensional run	`hills`, optionally `plumed` and `grid_range`	A profile reconstructed by PLUMED `sum_hills`
Independent infrequent-metadynamics trajectories	`colvar=[...]`, `infrequent=True`, explicit basins, raw bias or acceleration column	First-passage rate estimates and diagnostics; a FES needs suitable additional information

Supply cv and temperature, or provide a plumed_input from which they can be read. Input energies default to "kJ/mol"; "kcal/mol" and "eV" are accepted. Output arrays and exported energy columns are always kJ/mol. cv_unit labels the coordinate; it does not rescale it. time_unit accepts "fs", "ps", "ns", or "s" and defaults to "ps".

Start with an existing FES grid¶

PLUMED-style tables start with a FIELDS header. The default energy column is file.free; use energy_column="your_column" for another name. A supplied FES needs at least three distinct, finite coordinate points and finite energies. Its points are sorted into increasing coordinate order.

from pathlib import Path
import numpy as np
from matmalib import mlmd, Plot

q = np.linspace(-1.5, 1.5, 101)
energy = 8 * (q**2 - 1)**2 + 0.4 * q
np.savetxt("example-fes.dat", np.column_stack((q, energy)),
           header="#! FIELDS q file.free", comments="")

simulation = mlmd()
fes = simulation.analyze_metadynamics(mlmd.AnalysisConfig(
    fes_file="example-fes.dat", cv="q", cv_unit="dimensionless",
    temperature=300, energy_unit="kJ/mol",
    reactant=(-1.4, -0.5), product=(0.5, 1.4),
))
print("Points:", len(fes.x), "Minimum:", fes.free_energy.min(), fes.energy_unit)
print("Product minus reactant:",
      round(fes.summary["delta_g_product_minus_reactant_kj_mol"], 3), "kJ/mol")
print("Projected forward barrier:",
      round(fes.summary["delta_g_dagger_forward_projected_kj_mol"], 3), "kJ/mol")
fig, ax = Plot().fes_1d(fes, smooth_window=0)
ax.set_title("Synthetic double-well profile")

The right-hand basin has a higher free energy, so the reported product-minus-reactant difference is positive. This difference integrates each basin's probability, including its width:

ΔG(product − reactant) = −RT ln(Pproduct / Preactant)

It is not the difference between the two minimum energies. At fixed volume, this quantity is a Helmholtz free-energy difference, ΔF. The analysis does not add standard-state, concentration, or radial Jacobian corrections.

Provide both reactant=(low, high) and product=(low, high), or neither. Ranges must increase, be disjoint, lie inside the profile, and contain sampled probability. Without them, the implementation tries to infer two geometric basins separated by a fully sampled barrier at least RT above both minima. Check summary["basin_assignment"] and review inferred ranges against the chemistry.

A projected barrier is the intervening profile maximum minus a basin minimum. Its location appears as a dashed line. eyring_proxy_forward and eyring_proxy_reverse, when available, apply the Eyring expression to those projected barriers using kappa (default 1; allowed range 0 < kappa <= 1). They are proxies, not dynamically recovered rates. A one-dimensional projection can hide important coordinates.

Reweight a normalized COLVAR¶

For a metadynamics trajectory, the log weight is:

log(w) = [V(s,t) − c(t) + other biases being removed] / (RT)

PLUMED's metad.rbias contains V − c(t). Alternatively, supply both metad.bias and metad.rct; the analysis subtracts the latter. Raw metadynamics bias alone is insufficient for population reweighting. Do not supply rbias_column and rct_column together.

This small example assigns weights 1, 2, 4, and 8 to four equally populated bins. The most heavily weighted bin becomes the free-energy minimum. Both supported ways of specifying normalization give the same result.

import numpy as np
from matmalib import mlmd, Plot

rt = 0.00831446261815324 * 300
q = np.tile([0.5, 1.5, 2.5, 3.5], 10)
rbias = np.tile(rt * np.log([1, 2, 4, 8]), 10)
rct = np.full(len(q), 20.0)
rows = np.column_stack((np.arange(len(q)), q, rbias, rbias + rct, rct))
np.savetxt("weighted-COLVAR", rows,
           header="#! FIELDS time q metad.rbias metad.bias metad.rct", comments="")

common = dict(colvar="weighted-COLVAR", cv="q", cv_unit="dimensionless",
              temperature=300, bins=4, grid_range=(0, 4), smooth_window=0)
fes = mlmd().analyze_metadynamics(mlmd.AnalysisConfig(
    **common, rbias_column="metad.rbias"))
equivalent = mlmd().analyze_metadynamics(mlmd.AnalysisConfig(
    **common, bias_column="metad.bias", rct_column="metad.rct"))
np.testing.assert_allclose(fes.free_energy, equivalent.free_energy)
print("Relative energies:", np.round(fes.free_energy, 3), "kJ/mol")
print("Weight effective sample size:",
      round(fes.summary["sampling"]["weight_effective_sample_size"], 2))
fig, ax = Plot().fes_1d(fes, show_basins=False)
ax.set_title("Known weights: 1, 2, 4, 8")

COLVAR must contain a time column with strictly increasing times and regularly spaced selected samples. start_time and end_time select an inclusive interval in the file's time units. At least two samples must remain. bins is an integer of at least 3; grid_range=(low, high) fixes the histogram edges. Keep this range and the number of bins fixed when comparing profiles.

The parser accepts repeated identical headers. When reading a COLVAR being written by a live simulation, it omits an unfinished final line and records a warning. Completed malformed rows, changing headers, and duplicate or decreasing times are errors. Do not concatenate restarted trajectories without first resolving time resets.

Inspect a histogram and its unsampled regions¶

An unbiased trajectory has unit weights. The same analysis API can represent that case with an explicitly zero bias column. Use zeros only when the samples are actually unbiased; they cannot correct a metadynamics run with missing bias information.

import numpy as np
from matmalib import mlmd, Plot

rng = np.random.default_rng(14)
q = np.r_[rng.normal(-0.9, 0.12, 900), rng.normal(0.9, 0.12, 600)]
rng.shuffle(q)
rows = np.column_stack((np.arange(len(q)) * 0.1, q, np.zeros(len(q))))
np.savetxt("unbiased-COLVAR", rows,
           header="#! FIELDS time q unbiased_weight", comments="")

fes = mlmd().analyze_metadynamics(mlmd.AnalysisConfig(
    colvar="unbiased-COLVAR", cv="q", temperature=300,
    rbias_column="unbiased_weight", bins=60, grid_range=(-1.5, 1.5),
    reactant=(-1.5, -0.3), product=(0.3, 1.5), smooth_window=0,
))
for name, value in fes.summary["sampling"].items():
    print(name, ":", value)
print("Finite bins:", np.isfinite(fes.free_energy).sum())
print("Forward barrier available:", "eyring_proxy_forward" in fes.summary)
fig, ax = Plot().fes_1d(fes)
ax.set_title("An unsampled gap stays a gap")

Empty histogram bins have +inf free energy and appear as gaps. The library neither adds pseudocounts nor interpolates a barrier across these bins. Basin populations can still be reported, but a missing path between their minima prevents a resolved barrier and Eyring proxy.

Read fes.summary["sampling"] before interpreting a COLVAR profile:

Diagnostic	Interpretation
`samples`, `start_time`, `end_time`	Number and time range of selected observations
`empty_bins`	Bins with no selected observations
`samples_outside_grid`	Selected samples excluded by the histogram range
`weight_effective_sample_size`	`(sum w)^2 / sum(w^2)`; a small value means a few weighted frames dominate

The effective sample size is a weight diagnostic over selected frames, including frames outside the grid. It does not correct for temporal correlation and is not an uncertainty estimate. Inspect fes.warnings, input field ranges in summary["inputs"], and any barrier_unavailable or profile_unavailable explanation as well.

Read defaults and walls from the saved PLUMED input¶

plumed_input="plumed.dat" can infer the CV, temperature, units, PRINT output path, bias labels, and supported wall definitions. Relative input paths are resolved against the saved PLUMED file. Explicit analysis settings override inferred defaults. Input parsing supports the common single-line and ACTION ... block syntax, with exactly one METAD action on one explicit CV; it is not a general PLUMED preprocessor.

This example writes a small saved input and a matching synthetic COLVAR. Its recorded wall energy is removed automatically in the default analysis. keep_wall=True retains that restraint, so the two curves answer different thermodynamic questions.

from pathlib import Path
import numpy as np
from matmalib import mlmd, Plot

Path("wall-example.dat").write_text(
    "UNITS LENGTH=nm ENERGY=kj/mol TIME=ps\n"
    "q: DISTANCE ATOMS=1,2\n"
    "metad: METAD ARG=q SIGMA=0.02 HEIGHT=1 PACE=100 BIASFACTOR=10 TEMP=300\n"
    "wall: UPPER_WALLS ARG=q AT=0.8 KAPPA=100 EXP=2\n"
    "PRINT ARG=q,metad.rbias,wall.bias FILE=wall-COLVAR STRIDE=10\n"
)
q = np.tile(np.linspace(0.1, 1.1, 11), 20)
wall = 100 * np.maximum(q - 0.8, 0)**2
np.savetxt("wall-COLVAR", np.column_stack((np.arange(len(q)), q, q * 0, wall)),
           header="#! FIELDS time q metad.rbias wall.bias", comments="")

common = dict(plumed_input="wall-example.dat", bins=11,
              grid_range=(0.05, 1.15), smooth_window=0)
removed = mlmd().analyze_metadynamics(mlmd.AnalysisConfig(**common))
retained = mlmd().analyze_metadynamics(mlmd.AnalysisConfig(**common, keep_wall=True))
fig, ax = Plot().fes_1d(removed, show_basins=False, label="Wall removed")
Plot().fes_1d(retained, ax=ax, show_basins=False,
              label="Wall retained", color="#d48122")
print("Inferred units:", removed.cv_unit, removed.summary["input_energy_unit"])

Without saved input, set additional_bias=("wall.bias",) to add that recorded bias to the population log weights. You are choosing which restraints to remove; the code cannot infer the desired physical ensemble from a column name alone.

For an existing FES or a HILLS reconstruction, upper_wall=(at, kappa) subtracts the quadratic energy kappa * max(q − at, 0)^2 in the input energy units. With hills and plumed_input, a compatible single-CV quadratic upper wall can be inferred automatically, including its EPS/OFFSET settings. keep_wall=True disables inferred wall removal. Do not combine it with upper_wall. Lower or more complex walls require an appropriate explicit treatment. Removing a wall's energy cannot recover a region the trajectory never sampled.

Reconstruct from HILLS on your computer¶

Use this local workflow when a run has HILLS but lacks normalized COLVAR weights:

from matmalib import mlmd, Plot

simulation = mlmd()
fes = simulation.analyze_metadynamics(mlmd.AnalysisConfig(
    hills="run/HILLS", plumed_input="run/plumed.dat",
    bins=200, grid_range=(0.15, 0.8),
    # plumed="/path/to/plumed",  # optional explicit executable
    output_dir="hills-analysis", plot=True,
))
fig, ax = Plot().fes_1d(fes, energy_unit="kcal/mol")

The executable is found beside the current Python executable or on PATH, unless plumed is supplied. This calls plumed sum_hills and reads its output. PLUMED already applies the well-tempered scaling: do not multiply the returned energies by the bias-factor correction again. If no grid_range is given, the range is inferred from the deposited hill centers. Inspect that choice rather than assuming it includes every region of interest.

Choose hills or fes_file, not both. grid_range and bins configure reconstruction or histogramming; they do not resample an existing FES file. grid_range is rejected with fes_file. Time-window selection and population-reweighting column arguments apply to COLVAR, not HILLS or supplied FES grids.

For older COLVAR files containing only raw bias, plumed_input plus the matching HILLS and explicit METAD grid bounds can recover missing c(t) using PLUMED. This recovery supports a nonperiodic, one-dimensional, single-walker, fixed-width well-tempered run with a constant bias factor greater than one. The grid must cover all samples and hill centers. Reconstructed and recorded raw bias must agree within 0.01 kJ/mol; inconsistent histories are rejected. Read summary["bias_reconstruction"] and its warnings. This path also requires local PLUMED.

Periodic CVs and multidimensional profiles¶

Simulation setup can define a torsion CV, but this analysis currently has no periodic-boundary-crossing basin or barrier treatment. FES grids marked periodic are rejected. COLVAR histogramming treats its numeric coordinate as a line; it does not join -pi to pi. A basin split across that seam will therefore be misrepresented unless you first choose and validate a suitable unwrapped reaction path. Simply deleting periodic metadata does not make the analysis valid.

Multidimensional HILLS and repeated-coordinate FES tables are also unsupported. Prepare a physically appropriate one-dimensional marginal grid externally before using fes_file; do not flatten a multidimensional surface into a one-dimensional list. These restrictions are distinct from periodic simulation cells: a nonperiodic distance CV measured in a periodic cell can still be suitable for this analysis.

Compare time blocks before interpreting the profile¶

One profile does not establish convergence, and the API does not automatically calculate block uncertainties. Compare time windows on the same grid, inspect repeated transitions, and check stability against bin resolution and independent runs. The following creates three independent synthetic blocks with the same underlying distribution, then compares their basin free energies.

import numpy as np
from matmalib import mlmd, Plot

rng = np.random.default_rng(27)
centers = np.linspace(-1.45, 1.45, 30)
target_energy = 6 * (centers**2 - 1)**2 + 0.5 * centers
probability = np.exp(-target_energy / (0.00831446261815324 * 300))
probability /= probability.sum()
block_size = 2000
q = rng.choice(centers, size=3 * block_size, p=probability)
np.savetxt("blocks-COLVAR", np.column_stack((np.arange(len(q)), q, q * 0)),
           header="#! FIELDS time q unbiased_weight", comments="")

differences = []
ax = None
for block, color in enumerate(("#235789", "#d48122", "#38995b")):
    fes = mlmd().analyze_metadynamics(mlmd.AnalysisConfig(
        colvar="blocks-COLVAR", cv="q", temperature=300,
        rbias_column="unbiased_weight", bins=30, grid_range=(-1.5, 1.5),
        start_time=block * block_size, end_time=(block + 1) * block_size - 1,
        reactant=(-1.5, -0.4), product=(0.4, 1.5), smooth_window=0,
    ))
    differences.append(fes.summary["delta_g_product_minus_reactant_kj_mol"])
    fig, ax = Plot().fes_1d(fes, ax=ax, show_basins=False,
                          label=f"Block {block + 1}", color=color)
print("Block basin differences:", np.round(differences, 3), "kJ/mol")
print("Sample standard deviation:", round(float(np.std(differences, ddof=1)), 3))
ax.set_title("Compare profiles with identical histogram edges")

The displayed standard deviation describes these three block estimates. It is not automatically a confidence interval or a standard error for a correlated molecular trajectory. Real blocks must be long enough to address correlation and explore the relevant states; biased trajectories also need valid weights in each selected interval. Agreement among undersampled windows can be misleading.

Analyze separate equilibrium runs separately. A list of COLVAR files is reserved for independent infrequent first-passage runs with infrequent=True; it is not an equilibrium-pooling option.

Smooth the display and export results¶

smooth_window defaults to 11 points. Set it to 0 to disable smoothing, or to an odd integer of at least 5. The implementation fits a local cubic polynomial, not a Gaussian kernel. It smooths each finite sampled segment independently, reduces the window for short segments, and leaves segments shorter than five points unchanged. Gaps remain gaps. The smoothed curve is shifted to a zero minimum.

Smoothing affects the plotted curve and optional smoothed CSV. Basin populations, barriers, and rate estimates use the raw profile. A smooth curve is not evidence of convergence. There is no Gaussian-smoothing configuration option.

from pathlib import Path
import numpy as np
from matmalib import mlmd, Plot

q = np.linspace(-1.5, 1.5, 101)
energy_kcal = (8 * (q**2 - 1)**2 + 0.4 * np.sin(35 * q)) / 4.184
np.savetxt("noisy-fes.dat", np.column_stack((q, energy_kcal)),
           header="#! FIELDS q energy", comments="")

simulation = mlmd()
directory = simulation.analyze_metadynamics(mlmd.AnalysisConfig(
    fes_file="noisy-fes.dat", cv="q", temperature=300,
    energy_column="energy", energy_unit="kcal/mol", cv_unit="dimensionless",
    smooth_window=11, output_dir="analysis-export", plot=True, overwrite=True,
), output_format="files")
fes = simulation.data
print("Saved:", ", ".join(path.name for path in sorted(directory.iterdir())))
print("Stored energy units:", fes.energy_unit)
fig, ax = Plot().fes_1d(fes, energy_unit="kcal/mol", show_raw=True,
                       show_basins=False, label="Local cubic smoothing")
ax.set_title("Input, display, and export units are explicit")

The generated CSV, JSON, and PNG files appear as downloads under the browser example. In a local run they are written into output_dir. Set overwrite=True only when replacing earlier analysis outputs is intended; it is enabled here to make repeated clicks work. Input files are protected from being overwritten by output filenames.

Output	Contents
`free_energy.csv`	Coordinate and raw relative free energy in kJ/mol
`free_energy_smoothed.csv`	Display-smoothed profile, only when smoothing is enabled
`summary.json`	Source files, method, units, diagnostics, basin results, and warnings
`free_energy.png`	Plot when `plot=True` and a profile is available

The default output_format="data" returns the Data object. It still writes files when output_dir is supplied. output_format="files" requires output_dir and returns a Path; the Data object remains available as simulation.data. plot=True also requires an output directory.

The returned object's main attributes are x, free_energy, free_energy_smoothed, two-column data, energy_unit, cv, cv_unit, temperature_K, summary, and warnings. Plot.fes_1d() returns (figure, axes) and retains them as plot.fig and plot.ax. Use energy_unit to change display units, smooth_window to override display smoothing, show_raw to show the original curve under a smoothed one, and show_basins to control basin shading and the barrier marker. ax, color, and label support overlays as shown above. Display changes do not alter the stored raw energy array.

Analyze independent infrequent-metadynamics passages¶

First-passage analysis has a different data contract from equilibrium population reweighting. It needs independent runs starting at time zero in the explicit reactant range, with strictly increasing times. Each run ends at its first recorded entry into the product range, or contributes right-censored exposure if no entry is observed.

Use the cumulative PLUMED acceleration column when available (acceleration_column="metad.acc"), or the raw nonnegative bias (bias_column="metad.bias"). A normalized rbias is not a kinetic acceleration bias. The estimator is events divided by total accelerated exposure. With no events, the result is a one-sided exponential rate bound instead of a measured zero rate.

import numpy as np
from matmalib import mlmd

rt = 0.00831446261815324 * 300
paths = []
for label, endpoint in (("event", 0.30), ("censored", 0.22)):
    path = f"passage-{label}.COLVAR"
    rows = np.column_stack(([0, 1, 2], [0.20, 0.23, endpoint],
                            np.full(3, rt * np.log(10))))
    np.savetxt(path, rows, header="#! FIELDS time q metad.bias", comments="")
    paths.append(path)

result = mlmd().analyze_metadynamics(mlmd.AnalysisConfig(
    colvar=paths, cv="q", cv_unit="nm", temperature=300, time_unit="ps",
    infrequent=True, bias_column="metad.bias",
    reactant=(0.18, 0.25), product=(0.28, 0.35),
))
for name in ("independent_runs", "events", "right_censored_runs", "rate_s_inv"):
    print(name, ":", result.summary[name])
print("FES available:", len(result.x) > 0)
print(result.summary["profile_unavailable"])

The constant synthetic bias gives an acceleration factor of 10: each two-picosecond record contributes 20 ps of accelerated exposure. One event over 40 ps gives 2.5e10 s⁻¹. This arithmetic example does not validate the physical assumptions of infrequent metadynamics.

The method assumes sufficiently infrequent deposition, negligible transition-state bias, independent equilibrated reactant starts, and exponential escape statistics. Active additional restraint-bias columns before first passage are rejected; adding a wall's energy to an acceleration factor is not a general kinetic correction. Output stride limits the event time resolution. Raw-bias integration should be checked against finer output; a recorded cumulative acceleration is preferable when available.

At least five uncensored independent runs enable a reproducible parametric-bootstrap exponential-fit diagnostic in the summary. Its p-value is a fit diagnostic, not proof of unbiased kinetics. kappa enters the reported activation-free-energy conversion and Eyring proxies, not the event/exposure rate estimator.

A kinetics-only result can have empty x and free_energy arrays. Do not pass it to Plot.fes_1d(). You can supply a suitable hills or fes_file alongside the independent COLVAR runs to analyze a profile as well. summary.json is still exportable when no FES is available.

Continue with the MLMD API and configuration reference and the script feature map. For the existing two-coordinate count-based landscape, see free-energy and population maps; Plot.fes() does not perform this metadynamics reweighting.