tools_path — reference paths from steered MD

Reference paths for path collective variables, estimated from steered MD.

A PATHMSD simulation needs a reference path: an ordered set of frames walking the system from reactant to product. reactiontools.tools_reaction builds that path with a nudged elastic band, which is accurate but needs both endpoints optimised up front. The functions here take the cheaper route of dragging a collective variable across with a moving harmonic restraint (see plumed_input_steered()) and picking frames out of the trajectory that leaves behind. The result is a path in the full simulation environment – solvent included – rather than one interpolated in vacuum.

The typical sequence is

  1. plumed_input_steered_pt() writes the pulling script,

  2. an MD engine runs it – openmmnqe.run_openmm_steered for an OpenMM run,

  3. path_from_steered_md() selects the frames and writes neb_path.pdb,

  4. plumed_input_neb_path() biases along it, with the LAMBDA estimate_path_lambda() recommends.

Note that the path atoms should belong to a single molecule. A trajectory is written with molecules wrapped into the periodic box, so a path spanning two molecules – the two bases of a pair, say – can have them wrapped to opposite sides of the box, which no amount of alignment repairs.

Trajectory readers use MDTraj, which is installed with the package and imported only when needed. The frame-selection functions (select_frames_by_cv(), select_frames_by_msd(), cv_from_colvar()) remain NumPy-only.

cv_from_colvar(colvar_file: str, n_frames: int, cv_name: str | None = None) → ndarray[source]

Read a CV from a PLUMED COLVAR file, one value per trajectory frame.

PLUMED and OpenMM disagree about when to write: PRINT fires at step 0 and every stride thereafter, while a reporter first fires one interval in. So a COLVAR written with the same stride as the reporter holds exactly one row more than the trajectory has frames, and its first row is dropped here. Any other length is resampled onto the frames by assuming both cover the same span of time at a constant rate.

Parameters:

  • colvar_file (str) – Path to the COLVAR file written by the steered run.

  • n_frames (int) – Number of frames in the trajectory the values are wanted for.

  • cv_name (str or None, optional) – Field name of the CV, e.g. 'pt_cv'. If None, the first column after time is used. Default is None.

Returns:

numpy.ndarray – CV value for each trajectory frame, of length n_frames.

estimate_path_lambda(pdb_path: str, length_unit: str = 'nm') → float[source]

Estimate the LAMBDA a reference path should be given in PATHMSD.

LAMBDA sets how sharply path.sss tells neighbouring frames apart. Too small and the whole path reads as one blurred position; too large and the CV jumps between frames instead of moving smoothly along them. The recommendation is the usual rule of thumb, 2.3 / avg_msd, which puts the crossover between consecutive frames at roughly one unit of the exponent. A summary of the path and any warnings are printed as it goes.

Parameters:

  • pdb_path (str) – Multi-frame PDB holding the path, as written by convert_xyz_to_plumed_ref() or path_from_steered_md().

  • length_unit (str, optional) – Length unit the PLUMED script this LAMBDA is going into works in: 'nm' (the default, and what PLUMED uses when driven from OpenMM) or 'A' (what it uses when driven from ASE, via PLUMED_ASE_UNITS). LAMBDA has units of inverse squared length, so the two answers differ by a factor of a hundred and using the wrong one is not a small error.

Returns:

float – Recommended LAMBDA, to pass to the plumed_input_neb_path* builders in reactiontools.tools_cv.

Raises:

  • ImportError – If mdtraj is not installed.

  • ValueError – If length_unit is not 'nm' or 'A', or the path holds fewer than two frames.

Notes

A path of 15 to 30 frames behaves best. Fewer and consecutive frames are too far apart for the CV to interpolate between; more and PATHMSD spends its time on frames that add nothing.

path_from_steered_md(traj_file: str, template_pdb: str = 'index_atoms.pdb', output_file: str = 'neb_path.pdb', colvar_file: str | None = 'COLVAR_SMD', cv_name: str | None = None, n_images: int = 15, atom_indices: Sequence[int] | ndarray | None = None, top: str | None = None, cv_start: float | None = None, cv_stop: float | None = None, smooth: int = 0, align: bool = True, atom_line: str | tuple[str, ...] = 'HETATM', length_unit: str = 'nm') → float[source]

Estimate a path collective variable from a steered MD trajectory.

Frames evenly spaced along the CV are pulled out of the trajectory, aligned, and written as the multi-model PDB that PATHMSD reads, in the same format convert_xyz_to_plumed_ref() produces for a NEB path. An XYZ copy is written alongside it for viewing.

Parameters:

  • traj_file (str) – Trajectory from the steered run, e.g. the smd_steps.pdb written by openmmnqe.run_openmm_steered.

  • template_pdb (str, optional) – PDB holding the path atoms, as written by openmmnqe.save_only_index_atoms. Its atom records are the template for the output. Default is 'index_atoms.pdb'. Note that this file is renumbered in place so that it and the path agree on atom numbering, which is what PLUMED expects of them.

  • output_file (str, optional) – Multi-model PDB to write the path to. Default is 'neb_path.pdb', which is what the PATHMSD inputs in reactiontools.tools_cv reference.

  • colvar_file (str or None, optional) – COLVAR file from the steered run. If None, frames are spaced by RMSD along the trajectory instead of by CV. Default is 'COLVAR_SMD'.

  • cv_name (str or None, optional) – Field name of the CV in colvar_file. If None, the first column after time is used. Default is None.

  • n_images (int, optional) – Number of frames in the path. Default is 15; PATHMSD behaves best with 15 to 30.

  • atom_indices (list of int or None, optional) – Atoms of the full system that make up the path, normally the same indices the template was written from. If None, every atom in the trajectory is used. Default is None.

  • top (str or None, optional) – Topology file for trajectory formats that carry none, such as DCD. Default is None.

  • cv_start (float or None, optional) – Ends of the CV range the path should span. If None, the first and last values in the COLVAR are used, i.e. the whole pull. Default is None.

  • cv_stop (float or None, optional) – Ends of the CV range the path should span. If None, the first and last values in the COLVAR are used, i.e. the whole pull. Default is None.

  • smooth (int, optional) – Number of neighbouring frames either side to average each path frame with. Default is 0, which keeps the frames as they were sampled. Try 2 or 3 if the path comes out jagged.

  • align (bool, optional) – Whether to superpose every frame on the first before selecting. Default is True.

  • atom_line (str or tuple of str, optional) – Record type the template’s atoms are written under. Default is 'HETATM', which is what OpenMM writes for the ligand-like residues these paths normally cover. ASE writes 'ATOM'.

  • length_unit (str, optional) – Length unit of the PLUMED run this path is for, passed on to estimate_path_lambda() to scale the recommended LAMBDA. Default is 'nm'. The written files are always in angstrom, which is what PDB and XYZ use; this affects only the returned number.

Returns:

float – The LAMBDA value recommended for this path, as reported by estimate_path_lambda().

Raises:

  • ImportError – If mdtraj is not installed.

  • ValueError – If the path atoms and the template PDB do not match, or the trajectory holds fewer frames than the path needs.

select_frames_by_cv(cv: Sequence[float] | ndarray, n_images: int, cv_start: float | None = None, cv_stop: float | None = None) → ndarray[source]

Choose the frames that are evenly spaced along a collective variable.

Parameters:

  • cv (array_like) – CV value for each frame of the trajectory.

  • n_images (int) – Number of frames to select.

  • cv_start (float or None, optional) – Ends of the CV range to cover. If None, the first and last values in cv are used. Default is None.

  • cv_stop (float or None, optional) – Ends of the CV range to cover. If None, the first and last values in cv are used. Default is None.

Returns:

numpy.ndarray – Indices of the selected frames, in path order.

select_frames_by_msd(xyz: ndarray, n_images: int) → ndarray[source]

Choose the frames that are evenly spaced along the trajectory itself.

Distance is measured as the RMSD between consecutive frames accumulated along the trajectory, which is the spacing PATHMSD cares about. Use this when there is no COLVAR to select on, and bear in mind that thermal jitter inflates the arc length of a noisy trajectory.

Parameters:

  • xyz (numpy.ndarray) – Coordinates with shape (n_frames, n_atoms, 3), already aligned.

  • n_images (int) – Number of frames to select.

Returns:

numpy.ndarray – Indices of the selected frames, in path order.