tools_cv — collective variables for proton transfer¶

Collective variables for proton transfer, and the PLUMED input that biases them.

A proton transfer will not happen on its own inside an MD run: the barrier is several times kBT and the event is rare enough that unbiased sampling never sees it. Every builder here writes a PLUMED script that biases the system along a collective variable so that it does, and returns alongside it the shell command that turns the resulting bias back into a free-energy surface – plumed sum_hills for well-tempered metadynamics, or the bundled reactiontools.opes scripts when OPES_METAD is used instead. Hand the script to whatever drives PLUMED (openmmnqe.run_openmm_prod for an OpenMM run, plumed_calculator() for an ASE one) and the surface it produces to plot_plumed_fes().

The builders fall into two families:

Coordination difference – plumed_input_1pt(), plumed_input_2pt_1d() and plumed_input_2pt_2d() bias one or two proton transfers directly, the CV running from +1 at the donor to -1 at the acceptor.
Path – plumed_input_neb_path() biases progress along a reference path instead, which is the option that copes when the reaction is not well described by any one geometric coordinate. Build the reference with reactiontools.tools_path.

plumed_input_steered() stands apart: it drags a CV rather than biasing it, and exists to generate the reference path the last family needs.

Building your own¶

A CV that encodes one particular system belongs with that study rather than here – the wobble base-pair builders that used to sit alongside these moved out to the study that uses them for exactly that reason. The plumbing under these builders is worth sharing though, so it is public: plumed_bias_and_fes() returns the bias line together with the command that reads the surface back, which is the pairing that goes wrong when the two are written separately; plumed_one_based(), plumed_units_header(), plumed_temperature_pair() and plumed_angle_radians() handle the conversions every script needs. Write the CV lines yourself, hand the label you biased to plumed_bias_and_fes(), and the result is a script indistinguishable from the ones here.

Units¶

These scripts are written for an external PLUMED, so they carry no UNITS line by default and are therefore in PLUMED’s own nanometres and kJ/mol – which is what it uses when driven from OpenMM. Pass units="ase" for angstrom and eV instead, which prepends PLUMED_ASE_UNITS and is what a run driven through plumed_calculator() needs.

Note that this is the opposite default to plumed_metad_input(), which writes the UNITS line unless told otherwise because it only ever feeds ASE. The two are not interchangeable: a script moved between them without changing units is wrong by a factor of ten in every length.

Lengths that come from the geometry are converted for you. Lengths you supply yourself – wall in plumed_input_neb_path() – are in whichever unit you asked for, and nothing checks that.

as_positions(source: Any) → ndarray[source]¶

Coordinates as a plain (n_atoms, 3) array of angstrom.

Accepts whatever holds a geometry: an array of coordinates, an ase.Atoms, an openmm.app.Modeller or anything else whose positions carry OpenMM units, or a path to a structure file ASE can read. OpenMM works in nanometres and ASE in angstrom; the conversion happens here so that no builder has to know which world it was called from, and so that nothing in this package imports OpenMM.

Parameters:: source (ase.Atoms or openmm.app.Modeller or array_like or str or path-like) – The geometry, or somewhere to read one from.
Returns:: numpy.ndarray – Coordinates in angstrom, of shape (n_atoms, 3).
Raises:: ValueError – If the coordinates are not (n_atoms, 3).

plumed_angle_radians(angle_lim: float) → float[source]¶

Convert a wall angle from degrees to the radians PLUMED wants.

Parameters:: angle_lim (float) – Angle in degrees.
Returns:: float – The angle in radians, rounded to two decimals as the scripts write it.

Build the metadynamics bias line and the matching FES-reconstruction command.

Every builder in this module goes through here, so that the bias it emits and the command that reads the surface back cannot disagree about the grid or the temperature – which is exactly what happened while three of them hand-rolled their own.

Parameters:

f_opes (bool) – If True, bias with OPES_METAD (reconstructed with the bundled OPES FES_from_State.py); otherwise use well-tempered METAD with plumed sum_hills.
arg (str) – The PLUMED ARG= value the bias acts on.
pace (int or float or str) – Values interpolated into the bias line. For a multi-dimensional bias pass pre-joined strings, e.g. f'{sigma},{sigma}'.
height (int or float or str) – Values interpolated into the bias line. For a multi-dimensional bias pass pre-joined strings, e.g. f'{sigma},{sigma}'.
sigma (int or float or str) – Values interpolated into the bias line. For a multi-dimensional bias pass pre-joined strings, e.g. f'{sigma},{sigma}'.
bias (int or float or str) – Values interpolated into the bias line. For a multi-dimensional bias pass pre-joined strings, e.g. f'{sigma},{sigma}'.
temperature (float) – Temperature in kelvin.
kt (float) – kBT in the script’s energy unit, for the --kt of the FES command.
grid_bin (int or str) – Grid bin count, pre-joined for a multi-dimensional bias.
grid_min (float or str or None, optional) – Grid bounds. When None the grid is left out of both commands, and PLUMED sizes it itself. Default is None.
grid_max (float or str or None, optional) – Grid bounds. When None the grid is left out of both commands, and PLUMED sizes it itself. Default is None.
label (str, optional) – Prefix of the bias line, being the label plus enough spaces to keep the script’s columns aligned. Default is 'metad: '.

Returns:

metad_line (str) – The PLUMED bias line.
fes_command (str) – Shell command reconstructing the free-energy surface from it.

plumed_input_1pt(geometry: Any, idx: Sequence[int], temperature: Any, r_0: float = 1.1, wall: float = 1.5, angle_lim: float = 130.0, pace: int = 500, height: float | None = None, sigma: float = 0.05, bias: float = 20.0, grid_min: float | None = -1.1, grid_max: float | None = 1.1, grid_bin: int = 200, kappa: float = 500.0, f_opes: bool = False, units: str = 'plumed') → tuple[str, str][source]¶

Build a PLUMED input that biases a single proton transfer with metadynamics.

The collective variable is the same one plumed_input_steered_pt() drags: the difference between the donor-hydrogen and acceptor-hydrogen coordination numbers, running from +1 (bonded to the donor) to -1 (bonded to the acceptor). Distance and angle walls keep the hydrogen bond intact while the proton moves. Use this to bias the transfer with metadynamics; use plumed_input_steered_pt() to drag it instead.

Parameters:

geometry (ase.Atoms or openmm.app.Modeller or array_like or str) – The reactant geometry, used to size the switching function and the walls. See as_positions().
idx (list of int) – Three 0-based atom indices, ordered donor, hydrogen, acceptor.
temperature (float or openmm.unit.Quantity) – Simulation temperature in kelvin, used for the BIASFACTOR/OPES scaling and to report --kt to the FES reconstruction command.
r_0 (float, optional) – Multiplier on the shorter of the two donor/acceptor-hydrogen distances that sets R_0 of the coordination switching function. Default is 1.1.
wall (float, optional) – Multiplier on the donor-acceptor distance that sets the upper wall keeping the hydrogen bond intact. Default is 1.5.
angle_lim (float, optional) – Lower wall on the donor-hydrogen-acceptor angle, in degrees. Default is 130.0.
pace (int, optional) – PACE of the metadynamics bias, in steps. Default is 500.
height (float, optional) – Gaussian height (HEIGHT) for standard METAD, or the BARRIER for OPES_METAD, in the script’s energy unit. None, the default, deposits 15 kJ/mol – written as 15.0 for units='plumed' and as its equivalent 0.155 for units='ase'.
sigma (float, optional) – Gaussian width of the bias, in CV units. Default is 0.05.
bias (float, optional) – Well-tempered BIASFACTOR. Ignored when f_opes is True. Default is 20.0.
grid_min (float, optional) – Bounds of the bias/FES grid. Default is -1.1 and 1.1.
grid_max (float, optional) – Bounds of the bias/FES grid. Default is -1.1 and 1.1.
grid_bin (int, optional) – Number of grid bins. Default is 200.
kappa (float, optional) – Spring constant of the distance and angle walls. Default is 500.0.
f_opes (bool, optional) – If True, bias with OPES_METAD instead of well-tempered METAD. Default is False.
units (str, optional) – 'plumed' (the default) for nanometres and kJ/mol, or 'ase' for angstrom and eV. See the module docstring.

Returns:

plumed_input (str) – The PLUMED input script.
fes_command (str) – Shell command that reconstructs the free-energy surface from the bias written by plumed_input: plumed sum_hills, or the bundled OPES FES_from_State.py when f_opes is True.

plumed_input_2pt_1d(geometry: Any, idx1: Sequence[int], idx2: Sequence[int], temperature: Any, r_0: float = 1.1, wall: float = 1.5, angle_lim: float = 130.0, pace: int = 500, height: float | None = None, sigma: float = 0.05, bias: float = 20.0, grid_min: float | None = -1.1, grid_max: float | None = 1.1, grid_bin: int = 200, kappa: float = 500.0, f_opes: bool = False, units: str = 'plumed') → tuple[str, str][source]¶

Build a PLUMED input that biases two proton transfers along one coordinate.

Each transfer gets its own coordination-difference CV, as in plumed_input_1pt(), and the two are averaged into a single concerted coordinate. That is the right picture when the two protons are expected to move together; use plumed_input_2pt_2d() when whether they do is the question.

Parameters:

geometry (ase.Atoms or openmm.app.Modeller or array_like or str) – The reactant geometry, used to size the switching functions and the walls. See as_positions().
idx1 (list of int) – Three 0-based atom indices each, ordered donor, hydrogen, acceptor, for the first and second proton transfer.
idx2 (list of int) – Three 0-based atom indices each, ordered donor, hydrogen, acceptor, for the first and second proton transfer.
temperature (float or openmm.unit.Quantity) – Simulation temperature in kelvin.
r_0 (float, optional) – Multiplier on the shorter donor/acceptor-hydrogen distance in each pair that sets R_0 of its coordination switching function. Default is 1.1.
wall (float, optional) – Multiplier on the larger of the two donor-acceptor distances that sets the upper wall keeping both hydrogen bonds intact. Default is 1.5.
angle_lim (float, optional) – Lower wall on each donor-hydrogen-acceptor angle, in degrees. Default is 130.0.
pace (int, optional) – PACE of the metadynamics bias, in steps. Default is 500.
height (float, optional) – Gaussian height (HEIGHT) for standard METAD, or the BARRIER for OPES_METAD, in the script’s energy unit. None, the default, deposits 15 kJ/mol – written as 15.0 for units='plumed' and as its equivalent 0.155 for units='ase'.
sigma (float, optional) – Gaussian width of the bias, in CV units. Default is 0.05.
bias (float, optional) – Well-tempered BIASFACTOR. Ignored when f_opes is True. Default is 20.0.
grid_min (float, optional) – Bounds of the bias/FES grid. Default is -1.1 and 1.1.
grid_max (float, optional) – Bounds of the bias/FES grid. Default is -1.1 and 1.1.
grid_bin (int, optional) – Number of grid bins. Default is 200.
kappa (float, optional) – Spring constant of the distance and angle walls. Default is 500.0.
f_opes (bool, optional) – If True, bias with OPES_METAD instead of well-tempered METAD. Default is False.
units (str, optional) – 'plumed' (the default) for nanometres and kJ/mol, or 'ase' for angstrom and eV. See the module docstring.

Returns:

plumed_input (str) – The PLUMED input script.
fes_command (str) – Shell command that reconstructs the free-energy surface from the bias written by plumed_input.

plumed_input_2pt_2d(geometry: Any, idx1: Sequence[int], idx2: Sequence[int], temperature: Any, r_0: float = 1.1, wall: float = 1.5, angle_lim: float = 130.0, pace: int = 500, height: float | None = None, sigma: float = 0.05, bias: float = 20.0, grid_min: float = -1.1, grid_max: float = 1.1, grid_bin: int = 200, kappa: float = 500.0, f_opes: bool = False, units: str = 'plumed') → tuple[str, str][source]¶

Build a PLUMED input that biases two proton transfers on a 2-D surface.

Each proton transfer gets its own donor-hydrogen/acceptor-hydrogen coordination-difference CV, as in plumed_input_1pt(), and the two are kept as separate axes (ARG=cv_diff1,cv_diff2) so the free-energy surface resolves how the transfers are correlated. Use plumed_input_2pt_1d() instead when only the combined coordinate is needed.

Parameters:

geometry (ase.Atoms or openmm.app.Modeller or array_like or str) – The reactant geometry, used to size the switching functions and the walls. See as_positions().
idx1 (list of int) – Three 0-based atom indices each, ordered donor, hydrogen, acceptor, for the first and second proton transfer.
idx2 (list of int) – Three 0-based atom indices each, ordered donor, hydrogen, acceptor, for the first and second proton transfer.
temperature (float or openmm.unit.Quantity) – Simulation temperature in kelvin.
r_0 (float, optional) – Multiplier on the shorter donor/acceptor-hydrogen distance in each pair that sets R_0 of its coordination switching function. Default is 1.1.
wall (float, optional) – Multiplier on the larger of the two donor-acceptor distances that sets the upper wall keeping both hydrogen bonds intact. Default is 1.5.
angle_lim (float, optional) – Lower wall on each donor-hydrogen-acceptor angle, in degrees. Default is 130.0.
pace (int, optional) – PACE of the metadynamics bias, in steps. Default is 500.
height (float, optional) – Gaussian height (HEIGHT) for standard METAD, or the BARRIER for OPES_METAD, in the script’s energy unit. None, the default, deposits 15 kJ/mol – written as 15.0 for units='plumed' and as its equivalent 0.155 for units='ase'.
sigma (float, optional) – Gaussian width of the bias along each axis, in CV units. Default is 0.05.
bias (float, optional) – Well-tempered BIASFACTOR. Ignored when f_opes is True. Default is 20.0.
grid_min (float, optional) – Bounds of the bias/FES grid, applied to both axes. Default is -1.1 and 1.1.
grid_max (float, optional) – Bounds of the bias/FES grid, applied to both axes. Default is -1.1 and 1.1.
grid_bin (int, optional) – Number of grid bins per axis. Default is 200.
kappa (float, optional) – Spring constant of the distance and angle walls. Default is 500.0.
f_opes (bool, optional) – If True, bias with OPES_METAD instead of well-tempered METAD. Default is False.
units (str, optional) – 'plumed' (the default) for nanometres and kJ/mol, or 'ase' for angstrom and eV. See the module docstring.

Returns:

plumed_input (str) – The PLUMED input script.
fes_command (str) – Shell command that reconstructs the free-energy surface from the bias written by plumed_input.

plumed_input_neb_path(temperature: Any, wall: float = 0.1, pace: int = 500, height: float | None = None, sigma: float = 0.1, bias: float = 5.0, grid_min: float | None = 0.0, grid_max: float | None = 26.0, grid_bin: int = 500, kappa: float = 500.0, lambda_val: float = 250.0, neigh_size: int = 8, f_opes: bool = False, units: str = 'plumed') → tuple[str, str][source]¶

Build a PLUMED input that biases progress along a reference path.

When the reaction is not well described by any one geometric coordinate – a transfer that only happens after the surroundings rearrange, say – there may be no CV to bias. PATHMSD sidesteps that by measuring progress along a reference path instead: path.sss says how far along it the system is and path.zzz how far off it has strayed. The bias goes on the first and a wall on the second, so the sampling follows the path rather than wandering away from it.

The reference is read from neb_path.pdb and the alignment template from index_atoms.pdb, both by those names in the working directory. Build them with path_from_steered_md() or convert_xyz_to_plumed_ref().

Parameters:

temperature (float or openmm.unit.Quantity) – Simulation temperature in kelvin.
wall (float, optional) – Upper wall on path.zzz, the distance off the path, as an absolute squared length. Default is 0.1.
pace (int, optional) – PACE of the metadynamics bias, in steps. Default is 500.
height (float, optional) – Gaussian height (HEIGHT) for standard METAD, or the BARRIER for OPES_METAD, in the script’s energy unit. None, the default, deposits 15 kJ/mol – written as 15.0 for units='plumed' and as its equivalent 0.155 for units='ase'.
sigma (float, optional) – Gaussian width of the bias, in units of path node index. Default is 0.1.
bias (float, optional) – Well-tempered BIASFACTOR. Ignored when f_opes is True. Default is 5.0.
grid_min (float, optional) – Bounds of the path.sss bias/FES grid, in units of path node index. Default is 0.0 and 26.0.
grid_max (float, optional) – Bounds of the path.sss bias/FES grid, in units of path node index. Default is 0.0 and 26.0.
grid_bin (int, optional) – Number of grid bins. Default is 500.
kappa (float, optional) – Spring constant of the path.zzz wall. Default is 500.0.
lambda_val (float, optional) – LAMBDA parameter of PATHMSD, controlling how sharply path.sss distinguishes neighbouring frames; see estimate_path_lambda(), which recommends one for a given path. Default is 250.0.
neigh_size (int, optional) – NEIGH_SIZE of PATHMSD, the number of reference frames considered when computing the path distance. Default is 8.
f_opes (bool, optional) – If True, bias with OPES_METAD instead of well-tempered METAD. Default is False.
units (str, optional) – 'plumed' (the default) for nanometres and kJ/mol, or 'ase' for angstrom and eV. See the module docstring, and note that LAMBDA and wall both scale with the choice.

Returns:

plumed_input (str) – The PLUMED input script.
fes_command (str) – Shell command that reconstructs the free-energy surface from the bias written by plumed_input.

plumed_input_steered(cv_block: str, cv_start: float, cv_stop: float, steps: int, cv_name: str = 'cv', kappa: float = 2000.0, stride: int = 100, steps_equil: int = 0, steps_relax: int = 0, colvar_file: str = 'COLVAR_SMD', extra_lines: str | None = None, units: str = 'plumed') → tuple[str, int][source]¶

Build a PLUMED input that drags a collective variable (steered MD).

The centre of a harmonic restraint is moved linearly from cv_start to cv_stop over steps MD steps, optionally after holding at the starting value for steps_equil steps and before holding at the final value for steps_relax steps. The trajectory this produces is a first guess at the reaction path, which path_from_steered_md() turns into a reference for PATHMSD.

Parameters:

cv_block (str) – PLUMED lines defining the CV to steer, ending with an action labelled cv_name.
cv_start (float) – Value of the CV the restraint starts at, normally the value of the reactant.
cv_stop (float) – Value of the CV the restraint finishes at, normally the value of the product.
steps (int) – Number of MD steps spent pulling from cv_start to cv_stop. Pulling slowly costs more time but leaves a path that is closer to the free energy valley.
cv_name (str, optional) – Label of the CV defined in cv_block. Must be a plain label rather than a component such as path.sss, because it is used to name the restraint’s output components. Default is 'cv'.
kappa (float, optional) – Spring constant of the moving restraint in kJ/mol per CV unit squared. Default is 2000.0. Too soft and the system lags behind the restraint, too stiff and the pulling heats it.
stride (int, optional) – How often the CV is written to colvar_file, in steps. Default is 100. Match this to the reporter interval of the MD run so that every trajectory frame has a CV value.
steps_equil (int, optional) – Steps held at cv_start before pulling starts. Default is 0.
steps_relax (int, optional) – Steps held at cv_stop after pulling finishes. Default is 0.
colvar_file (str, optional) – File the CV, the restraint centre and the work are written to. Default is 'COLVAR_SMD'.
extra_lines (str or None, optional) – Further PLUMED lines (walls, restraints, extra prints) inserted after the CV definition. Default is None.
units (str, optional) – 'plumed' (the default) for nanometres and kJ/mol, or 'ase' for angstrom and eV. See the module docstring.

Returns:

plumed_input (str) – The PLUMED input script.
n_steps (int) – Total number of MD steps the pulling schedule covers, i.e. steps_equil + steps + steps_relax. Pass this to the MD run so the simulation does not stop mid-pull.

plumed_input_steered_pt(geometry: Any, idx: Sequence[int], steps: int, r_0: float = 1.1, wall: float = 1.5, angle_lim: float = 130.0, kappa: float = 2000.0, stride: int = 100, cv_start: float | None = None, cv_stop: float | None = None, steps_equil: int = 0, steps_relax: int = 0, colvar_file: str = 'COLVAR_SMD', wall_kappa: float = 500.0, units: str = 'plumed') → tuple[str, int][source]¶

Build a steered MD input that pulls a proton across a hydrogen bond.

The collective variable is the one plumed_input_1pt() biases, the difference between the donor-hydrogen and acceptor-hydrogen coordination numbers, and the same walls keep the donor and acceptor from drifting apart while the proton is dragged over. Unless they are given, the start and end values of the CV are read off geometry: the current value for the start, and its mirror image for the end, which is where the proton sits once transferred.

Parameters:

geometry (ase.Atoms or openmm.app.Modeller or array_like or str) – The reactant geometry, used to size the switching function and to work out where the pull should start. See as_positions().
idx (list of int) – Three 0-based atom indices, ordered donor, hydrogen, acceptor.
steps (int) – Number of MD steps spent pulling the proton across.
r_0 (float, optional) – Multiplier on the shorter of the two donor/acceptor-hydrogen distances that sets R_0 of the switching function. Default is 1.1.
wall (float, optional) – Multiplier on the donor-acceptor distance that sets the upper wall keeping the hydrogen bond intact. Default is 1.5.
angle_lim (float, optional) – Lower wall on the donor-hydrogen-acceptor angle, in degrees. Default is 130.0.
kappa (float, optional) – Spring constant of the moving restraint. Default is 2000.0.
stride (int, optional) – How often the CV is written, in steps. Default is 100.
cv_start (float or None, optional) – Explicit start and end values for the CV. If None, the start is computed from geometry and the end is its negative. Default is None.
cv_stop (float or None, optional) – Explicit start and end values for the CV. If None, the start is computed from geometry and the end is its negative. Default is None.
steps_equil (int, optional) – Steps held at the start and end values. Default is 0.
steps_relax (int, optional) – Steps held at the start and end values. Default is 0.
colvar_file (str, optional) – File the CV is written to. Default is 'COLVAR_SMD'.
wall_kappa (float, optional) – Spring constant of the distance and angle walls. Default is 500.0.
units (str, optional) – 'plumed' (the default) for nanometres and kJ/mol, or 'ase' for angstrom and eV. See the module docstring.

Returns:

plumed_input (str) – The PLUMED input script.
n_steps (int) – Total number of MD steps the pulling schedule covers.

plumed_one_based(indices: Iterable[int]) → list[int][source]¶

Convert 0-based atom indices to PLUMED’s 1-based convention.

Order is preserved, which is what a CV needs: the builders here read their indices positionally, as donor, hydrogen, acceptor. This is not interchangeable with plumed_selection(), which sorts and dedupes its indices to collapse them into ranges – right for the atom group of a COORDINATION, wrong for anything positional.

Parameters:: indices (iterable of int) – 0-based indices, as OpenMM and ASE number them.
Returns:: list of int – 1-based indices, as PLUMED numbers them, in the order given.

plumed_temperature_pair(temperature: Any, units: str) → tuple[float, float][source]¶

Temperature in kelvin and kBT in the script’s energy unit.

Parameters:

temperature (float or openmm.unit.Quantity) – Simulation temperature.
units (str) – 'plumed' or 'ase'.

Returns:

kelvin (float) – For the TEMP of the bias action.
kt (float) – For the --kt of the reconstruction command.

plumed_units_header(units: str) → str[source]¶

The UNITS line the script needs, if any.

Parameters:: units (str) – 'plumed' or 'ase'.
Returns:: str – The line plus a newline, or an empty string for PLUMED’s own units.

switching_value(r: float, r_0: float, nn: int = 6, mm: int | None = None) → float[source]¶

Evaluate PLUMED’s default rational switching function.

This is the same function COORDINATION applies to every pair distance, s(r) = (1 - (r/r_0)^nn) / (1 - (r/r_0)^mm), so it can be used to work out what a coordination-based CV is worth for a given geometry without running PLUMED – which is how plumed_input_steered_pt() decides where a pull should start and stop.

Parameters:

r (float) – Distance between the two atoms, in the same units as r_0.
r_0 (float) – The R_0 parameter of the switching function.
nn (int, optional) – Numerator exponent. Default is 6, as in PLUMED.
mm (int or None, optional) – Denominator exponent. If None, 2 * nn is used, again as in PLUMED.

Returns:

float – The value of the switching function, between 0 and 1.