tools_plumed — metadynamics support¶

PLUMED support: build the input, drive a biased run, sum the hills.

The three stages of a metadynamics calculation, in the order they happen. plumed_selection() and find_molecules() pick out the atoms a collective variable is built from, plumed_metad_input() turns a CV and a few METAD settings into the input PLUMED reads, plumed_calculator() wraps an ASE calculator in that bias so an ordinary ASE molecular-dynamics run becomes a biased one, and run_sum_hills() turns the hills it deposited into a free-energy surface for reactiontools.tools_fes to plot.

run_opes_fes() is the OPES_METAD counterpart of run_sum_hills(). OPES deposits no hills to add up, writing a running estimate of the bias to a STATE file instead, so the surface is read back out of that by one of the scripts bundled in reactiontools.opes. run_opes_reweighting() rebuilds the surface from the COLVAR samples and their bias column instead, which is how several independent walkers are combined into one estimate – combine_colvar_files() does the merge.

Only plumed_calculator() needs the plumed Python module; only run_sum_hills() needs the plumed executable. The rest is string handling and works without either.

PLUMED_ASE_UNITS = 'UNITS ENERGY=eV LENGTH=A TIME=fs'¶: PLUMED reads and writes its own units unless the input says otherwise, and they are not ASE’s: without this line a run driven from ASE reports lengths in nm and energies in kJ/mol. Declaring it makes COLVAR, HILLS and fes.dat come out in the Å and eV that the rest of this package, and plot_plumed() in particular, assume.

plumed_selection(indices: Iterable[int]) → str[source]¶

Format atom indices as a PLUMED ATOMS= selection string.

Parameters:: indices (iterable of int) – Zero-based atom indices.
Returns:: str – Comma-separated PLUMED selection using one-based indexing and compact ranges.
Raises:: ValueError – If no indices are given.

find_molecules(atoms: Atoms) → list[ndarray][source]¶

Return connected atom groups identified as molecules.

Parameters:: atoms (ase.Atoms) – Structure for which the bonded graph should be analysed.
Returns:: list of numpy.ndarray – Atom-index arrays, one per connected component.

plumed_metad_input(cvs: Sequence[str], sigma: float | Sequence[float], height: float, pace: int, biasfactor: float | None = None, temperature: float | None = None, hills: str = 'HILLS', colvar: str | None = 'COLVAR', stride: int = 10, units: bool = True, metad_extra: str | None = None, extra: Sequence[str] | None = None) → list[str][source]¶

Build the PLUMED input for a metadynamics run.

Assembles the lines plumed_calculator() takes: the units, the collective variables as given, a METAD action biasing all of them, and a PRINT collecting them into a COLVAR. The METAD file defaults to the name run_sum_hills() looks for, so the two ends of the workflow line up without being told to.

Parameters:

cvs (sequence of str) – Collective-variable lines in PLUMED syntax, each labelled, as in "d1: DISTANCE ATOMS=1,2". plumed_selection() formats the atom lists. Passed through untouched, so any PLUMED action can be used.
sigma (float or sequence of float) – Gaussian width per collective variable, in the units of that variable. A single value is used for all of them.
height (float) – Gaussian height in eV, or in kJ/mol if units is False.
pace (int) – Number of steps between deposited Gaussians.
biasfactor (float or None, optional) – Bias factor for well-tempered metadynamics, which must be greater than one. None, the default, deposits Gaussians of fixed height instead, which never stops filling and so has no converged surface to read off.
temperature (float or None, optional) – Simulation temperature in kelvin, written as TEMP. Required with biasfactor, and must match the temperature the dynamics actually runs at — nothing checks that for you.
hills (str, optional) – File the Gaussians are written to.
colvar (str or None, optional) – File the collective variables are printed to. None omits the PRINT.
stride (int, optional) – Number of steps between PRINT lines.
units (bool, optional) – Prepend PLUMED_ASE_UNITS, so the run reads and writes Å and eV. Leave it on unless you mean to work in PLUMED’s own nm and kJ/mol, in which case sigma, height and the output files are all in those units and plot_plumed() will read them wrongly.
metad_extra (str or None, optional) – Extra keywords appended to the METAD line, for the settings not given their own argument here — GRID_MIN, GRID_MAX and GRID_BIN above all, which a long run wants.
extra (sequence of str or None, optional) – Extra whole lines appended at the end, for walls, further prints, or anything else.

Returns:

list of str – The input lines, ready for plumed_calculator().

Raises:

ValueError – If no collective variable is given, if a line carries no label, if the number of sigmas does not match the number of variables, or if biasfactor is given without a temperature or is not above one.

Examples

Well-tempered metadynamics along one distance:

lines = plumed_metad_input(
    cvs=["d1: DISTANCE ATOMS=1,2"],
    sigma=0.05, height=0.02, pace=100,
    biasfactor=10, temperature=300)

plumed_calculator(atoms: Atoms, calc: Any, input_lines: Sequence[str], timestep: float, temperature: float | None = None, log: str = '', restart: bool = False) → Iterator[Plumed][source]¶

Bias an ASE calculator with PLUMED, for the length of the block.

Wraps calc in an ase.calculators.plumed.Plumed and hangs it on atoms, so an ordinary ASE dynamics run becomes a biased one: the integrator asks for forces as usual and PLUMED adds the bias on top.

A context manager because PLUMED buffers what it writes and only flushes on finalize. Run the dynamics inside the block; a run that returns without it leaves HILLS short of the hills it deposited last, and a free-energy surface summed from that is quietly wrong rather than obviously missing.

Parameters:

atoms (ase.Atoms) – Structure to simulate. Its calculator is swapped for the biased one for the length of the block and put back on the way out, so the atoms come back as they were, holding a calculator that still works.
calc (ase.calculators.Calculator) – Calculator computing the unbiased forces.
input_lines (sequence of str) – PLUMED input, one action per line, from plumed_metad_input() or written by hand. Hand-written input wants PLUMED_ASE_UNITS first.
timestep (float) – Timestep in ASE time units — the same value the integrator is given, for example 0.5 * ase.units.fs. PLUMED counts its own steps from it, so a mismatch silently misdates every hill.
temperature (float or None, optional) – Simulation temperature in kelvin, converted to the thermal energy PLUMED is told about. Needed by the methods that scale a bias by it, well-tempered metadynamics above all; None leaves ASE’s default, which those methods must not be run with. Match the temperature given to plumed_metad_input() and to the thermostat.
log (str, optional) – PLUMED’s own log file. Empty, the default, writes to stdout.
restart (bool, optional) – Continue a previous run, appending to its hills rather than starting the bias from nothing.

Yields:

ase.calculators.plumed.Plumed – The biased calculator, already attached to atoms.

Raises:

ImportError – If the plumed Python module is not installed.

Examples

Well-tempered metadynamics along one distance, at 300 K:

from ase import units
from ase.md.langevin import Langevin

lines = plumed_metad_input(cvs=["d1: DISTANCE ATOMS=1,2"],
                           sigma=0.05, height=0.02, pace=100,
                           biasfactor=10, temperature=300)

with plumed_calculator(atoms, calc, lines, timestep=0.5 * units.fs,
                       temperature=300):
    Langevin(atoms, 0.5 * units.fs, temperature_K=300,
             friction=0.01).run(10000)

run_sum_hills()

sum_hills_files(outfile: str | Path = 'fes.dat') → list[Path][source]¶

List the surfaces a strided run_sum_hills() wrote, in order.

--stride does not number the file it was given: it writes f"{outfile}{n}.dat", so the default outfile yields fes.dat0.dat, fes.dat1.dat and so on. Two things follow, and this exists because both are easy to get wrong. The obvious glob, fes*.dat, is right only by accident; and sorting the names puts fes.dat10.dat before fes.dat2.dat, which for a convergence series — where the order is the entire point — silently scrambles the answer.

Parameters:: outfile (str or path-like, optional) – The outfile given to run_sum_hills().
Returns:: list of pathlib.Path – The surfaces, ordered by the index PLUMED gave them, ready to hand to plot_fes_1d(). Empty if the run was not strided.

Examples

A convergence series, labelled by the simulated time each surface covers:

run_sum_hills(stride=100, outfile="fes")
surfaces = sum_hills_files("fes")
plot_fes_1d(surfaces,
            labels=[(i + 1) * 100 for i in range(len(surfaces))],
            label_template="{:g} hills")

Run plumed sum_hills to build a free-energy surface from the hills.

The paths are resolved by the plumed executable, so this acts on the current working directory unless absolute paths are given.

Parameters:

hills (str or path-like, optional) – Hills file written by the METAD action.
outfile (str or path-like, optional) – Free-energy surface file to write, as read by plot_plumed(). With stride it is a stem rather than a filename; see there.
mintozero (bool, optional) – Pass --mintozero so the surface minimum sits at zero.
stride (int or None, optional) – Write a surface every stride hills instead of one at the end, which is how a convergence series is made: a run is converged when the last few surfaces lie on top of each other. The files are named f"{outfile}{n}.dat", so the default outfile gives the unlovely fes.dat0.dat — pass outfile="fes" for fes0.dat. sum_hills_files() collects them in the right order either way.
nohistory (bool, optional) – With stride, make each surface from only the hills in its own interval rather than from everything up to it. Useful to watch where the bias is being deposited; not what a convergence series wants.
grid_min (float or sequence of float, optional) – Bounds of the output grid, one per collective variable. Worth setting for a series, since PLUMED otherwise picks bounds per surface from the hills it has so far and the surfaces come back on grids that do not line up.
grid_max (float or sequence of float, optional) – Bounds of the output grid, one per collective variable. Worth setting for a series, since PLUMED otherwise picks bounds per surface from the hills it has so far and the surfaces come back on grids that do not line up.
grid_bin (int or sequence of int, optional) – Number of bins per collective variable.
idw (str or sequence of str, optional) – Collective variables to keep, by label; the rest are integrated out, which needs kt. This is how a two-dimensional surface is projected onto one of its variables.
kt (float or None, optional) – Thermal energy for that integration, in the energy units of the hills file — eV for a run built by plumed_metad_input(), where ase.units.kB * 300 is 300 K. Only used with idw.
negbias (bool, optional) – Print the negative bias rather than the free energy.
extra (sequence of str, optional) – Further arguments appended to the command line, for the options without their own keyword here — --spacing, --fmt, --histo and the rest.
verbose (bool, optional) – Print the command being run.

Returns:

str – The command line that was run. With stride, the surfaces themselves are gathered by sum_hills_files().

Raises:

ValueError – If kt is given without idw, which would silently do nothing.
subprocess.CalledProcessError – If plumed exits non-zero.

combine_colvar_files(colvar_files: Sequence[str | Path], outfile: str | Path = 'COLVAR', sort_by_time: bool = True) → Path[source]¶

Merge several COLVAR files into one, keeping a single header.

The merge run_opes_reweighting() wants for independent walkers, done the two ways the bundled script’s own notes describe. Sorted by time (sort -gs COLVAR.*), the walkers interleave into one series, which is what stride and skiprows expect – a stride then samples across walkers, and skipped rows discard every walker’s transient together. Unsorted (cat COLVAR.*), the walkers stay contiguous, so blocks set to the number of walkers gives one block per walker and the error bars come from their scatter.

Parameters:

colvar_files (sequence of path-like) – The COLVAR files, one per walker. Their header blocks – the leading #! FIELDS line and any #! SET lines after it – must be identical, since one copy speaks for all of them.
outfile (str or path-like, optional) – The merged file to write. Must not be one of the inputs.
sort_by_time (bool, optional) – Stable-sort the rows on their first column. Default True; pass False for the contiguous per-walker layout blocks wants.

Returns:

pathlib.Path – The merged file.

Raises:

ValueError – If no files are given, outfile is one of the inputs, a file does not open with a #! FIELDS header, the headers disagree, or a row’s time is not a number.

Rebuild a free-energy surface from an OPES state file.

The OPES_METAD counterpart of run_sum_hills(), and the other half of what the f_opes switch on the reactiontools.tools_cv builders selects: they emit the bias action, this reads the surface back out of what it wrote.

Paths are resolved by the script, so this acts on the current working directory unless absolute paths are given.

Parameters:

state (str or path-like, optional) – State file written by the OPES_METAD action’s STATE_WFILE.
outfile (str or path-like, optional) – Free-energy surface file to write, as read by as_fes().
grid_min (float or sequence of float, optional) – Bounds of the output grid, one per collective variable.
grid_max (float or sequence of float, optional) – Bounds of the output grid, one per collective variable.
grid_bin (int or sequence of int, optional) – Number of bins per collective variable.
kt (float or None, optional) – Thermal energy in the energy units of the state file. See thermal_energy().
extra (sequence of str, optional) – Further arguments appended to the command line, for the options without their own keyword here – --deltaFat, --all_stored, --der and the rest.
verbose (bool, optional) – Print the command being run.

Returns:

str – The command line that was run.

Raises:

ValueError – If only one of grid_min and grid_max is given.
subprocess.CalledProcessError – If the script exits non-zero.

Rebuild a free-energy surface by reweighting COLVAR samples.

The reweighting counterpart of run_opes_fes(): that reads the bias’s own running estimate out of the STATE file, this rebuilds the surface from the samples and the bias they were collected under. For a single run the two should agree, and their disagreement is a convergence check. For several independent walkers this is the one that combines them – each sample carries its own walker’s bias, so the merged file combine_colvar_files() writes is reweighted as one: sorted by time for a single series (stride, skiprows), concatenated per walker for cross-walker error bars (blocks).

Paths are resolved by the script, so this acts on the current working directory unless absolute paths are given.

Parameters:

sigma (float or sequence of float) – Kernel bandwidth of the weighted density estimate, one value per collective variable. A fraction of the OPES SIGMA is the usual starting point.
kt (float) – Thermal energy in the energy units of the COLVAR file. See thermal_energy().
colvar (str or path-like, optional) – The COLVAR file, or the merged one combine_colvar_files() wrote.
outfile (str or path-like, optional) – Free-energy surface file to write, as read by as_fes().
cv (str, int, or sequence, optional) – Collective variable to bin, by column name or 1-based column number; two entries make a two-dimensional surface.
bias (str or None, optional) – Bias column to weight by; the script’s default matches any column containing .bias.
grid_min (float or sequence of float, optional) – Bounds of the output grid, one per collective variable.
grid_max (float or sequence of float, optional) – Bounds of the output grid, one per collective variable.
grid_bin (int or sequence of int, optional) – Number of bins per collective variable.
blocks (int or None, optional) – Number of blocks for the uncertainty column – one per walker on an unsorted merge.
stride (int or None, optional) – Also write the surface every stride samples, as a convergence series.
skiprows (int or None, optional) – Initial rows to discard, the transient before the bias settled.
extra (sequence of str, optional) – Further arguments appended to the command line.
verbose (bool, optional) – Print the command being run.

Returns:

str – The command line that was run.

Raises:

ValueError – If only one of grid_min and grid_max is given, or both blocks and stride are.
subprocess.CalledProcessError – If the script exits non-zero.