Source code for openpathsampling.collectivevariables.plumed_wrapper

import openpathsampling as paths
from openpathsampling.netcdfplus import StorableNamedObject
import openpathsampling.engines as peng
try:
    import plumed
except ImportError:
    pass
import numpy as np
import warnings
import copy
import os
import sys


[docs]class PLUMEDCV(paths.collectivevariable.CoordinateFunctionCV): """Make `CollectiveVariable` computed by PLUMED [PLU2a]_ according to the command `name`: `definition`, where `name` is a PLUMED label and `definition` contains all PLUMED keywords. Takes an `openpathsampling.engines.trajectory.Trajectory` as input. References ---------- .. [PLU2a] G.A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, G. Bussi, PLUMED2: New feathers for an old bird, Comp. Phys. Comm. 185, 604 (2014); https://doi.org/10.1016/j.cpc.2013.09.018 Examples -------- >>> # To create a `CollectiveVariable` which calculates the dihedral psi >>> # formed by the atoms [7,9,15,17] in Ala dipeptide: >>> from openpathsampling import PLUMEDCV, PLUMEDInterface >>> plmd = PLUMEDInterface(top) >>> # top is an `openpathsampling.engines.topology.MDTrajTopology` >>> psi_plumed = PLUMEDCV("psi",plmd,"TORSION ATOMS=7,9,15,17") >>> print psi_plumed(traj) # returns psi values for the trajectory """
[docs] def __init__(self, name, plmd, definition, components=None, cv_requires_lists=True, cv_wrap_numpy_array=True, cv_scalarize_numpy_singletons=True, **kwargs ): """ Parameters ---------- name : string A descriptive name of the PLUMED collective variable, equivalent to a `label` in a PLUMED input file. plmd : :obj:`openpathsampling.collectivevariable.PLUMEDInterface` An interface to the Cython PLUMED wrapper. If the PLUMED collective variable is a function of previously defined ones, or if it is defined based on group/virtual atoms, the `plmd` interface must be the same one that was used for the preceding instantiantions. definition : string The PLUMED keywords that define the collective variable (see http://www.plumed.org/documentation). components : list of string The components (either default of customized) of the PLUMED collective variable (see http://www.plumed.org/documentation). cv_requires_lists cv_wrap_numpy_array cv_scalarize_numpy_singletons kwargs """ super(PLUMEDCV, self).__init__( name, f=PLUMEDCV.compute_cv, cv_requires_lists=cv_requires_lists, cv_wrap_numpy_array=cv_wrap_numpy_array, cv_scalarize_numpy_singletons=cv_scalarize_numpy_singletons, **kwargs ) self.plmd = plmd self.definition = definition if components is None: components = [] self.components = components self.topology = plmd.topology self.var = self.create_plumed_var(name, definition)
def create_plumed_var(self, name, definition): """Create a PLUMED collective variable. Parameters ---------- name : string A descriptive name of the PLUMED collective variable, equivalent to a `label` in a PLUMED input file. definition : string The PLUMED keywords that define the collective variable (see http://www.plumed.org/documentation). Returns ------- data : array Array to store the computed PLUMED collective variable. """ self.plmd.cmd("readInputLine", name + ": " + definition) if (len(self.components)): cdata = [] for c in self.components: shape = np.zeros(1, dtype=np.int_) self.plmd.cmd("getDataRank " + name + "." + c, shape) data = np.zeros((1)) self.plmd.cmd("setMemoryForData " + name + "." + c, data) cdata.append(data) return cdata else: shape = np.zeros(1, dtype=np.int_) self.plmd.cmd("getDataRank " + name, shape) data = np.zeros((1)) self.plmd.cmd("setMemoryForData " + name, data) return data def compute_cv(self, trajectory): """Compute a PLUMED collective variable. Parameters ---------- trajectory : :obj:`openpathsampling.engines.trajectory.Trajectory` The trajectory along which the collective variable is to be computed. Returns ------- cv : array Computed values of the PLUMED collective variable along the `openpathsampling.engines.trajectory.Trajectory` """ cv = [] try: masses = np.array([a.element.mass for a in trajectory.topology.mdtraj.atoms], dtype=np.float64) except AttributeError: # pragma: no cover masses = np.ones(self.topology.n_atoms, dtype=np.float64) warnings.warn("No masses found in topology. All masses set to one") charges = np.zeros(self.topology.n_atoms, dtype=np.float64) # non-essential # warnings.warn("All charges set to zero") forces = np.zeros((self.topology.n_atoms, 3)) virial = np.zeros((3, 3), dtype=np.float64) bias = np.zeros((1), dtype=np.float64) # non-essential for step, snapshot in enumerate(trajectory): self.plmd.cmd("setStep", step) if snapshot.box_vectors != None: box = np.array(snapshot.box_vectors, dtype=np.float64) self.plmd.cmd("setBox", box) positions = snapshot.xyz.astype(np.float64) self.plmd.cmd("setBox", box) self.plmd.cmd("setPositions", positions) self.plmd.cmd("setMasses", masses) self.plmd.cmd("setForces", forces) self.plmd.cmd("setVirial", virial) self.plmd.cmd("setCharges", charges) # non-essential self.plmd.cmd("getBias", bias) # non-essential self.plmd.cmd("calc") cv.append(copy.deepcopy(self.var)) cv = np.array(cv) return cv def _eval(self, trajectory): trajectory = peng.Trajectory(trajectory) return self.cv_callable(self, trajectory) def to_dict(self): return { 'name': self.name, 'plmd': self.plmd, 'definition': self.definition, 'components': self.components, 'kwargs': self.kwargs, 'cv_requires_lists': self.cv_requires_lists, 'cv_wrap_numpy_array': self.cv_wrap_numpy_array, 'cv_scalarize_numpy_singletons': self.cv_scalarize_numpy_singletons }
[docs]class PLUMEDInterface(StorableNamedObject): """Interfaces the Cython PLUMED wrapper [PLU2b]_ located at `/path/to/plumed2/python` and allows to set and get non-`PLUMEDCV` commands (i.e., non-outputting). This includes groups of atoms, centers of mass, include files, etc. Requires PLUMED development version (see https://github.com/plumed/plumed2) and sourcing `/path/to/plumed2/sourceme.sh`. References ---------- .. [PLU2b] G.A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, G. Bussi, PLUMED2: New feathers for an old bird, Comp. Phys. Comm. 185, 604 (2014); https://doi.org/10.1016/j.cpc.2013.09.018 Examples -------- >>> # To group the atoms [7,9,15,17] corresponding to the dihedral psi >>> # in Ala dipeptide: >>> from openpathsampling import PLUMEDCV, PLUMEDInterface >>> plmd = PLUMEDInterface(top) >>> # top is an `openpathsampling.engines.topology.MDTrajTopology` >>> plmd.set("group","GROUP ATOMS=7,9,15,17") >>> psi_plumed = PLUMEDCV("psi",plmd,"TORSION ATOMS=group") >>> print psi_plumed(traj) # returns psi values for the trajectory >>> pld.get() # returns (('group', 'GROUP ATOMS=7,9,15,17')) """
[docs] def __init__(self, topology, pathtoplumed="", timestep=1., kbt=1., molinfo="", logfile="plumed.log"): """ Parameters ---------- topology : :obj:`openpathsampling.engines.topology.MDTrajTopology` pathtoplumed : string path to the PLUMED installation timestep : double Time step size of the simulation in PLUMED default units (ps). kbt : double :math:`$k_BT$` in PLUMED default units (kJ/mol). molinfo : string A PDB file containing information about the molecule. (see https://plumed.github.io/doc-v2.4/user-doc/html/_m_o_l_i_n_f_o.html). logfile : string Name of the PLUMED log file. """ self.interface = plumed.Plumed() # 8 is default size of real self.pathtoplumed = pathtoplumed self.topology = topology self.timestep = timestep self.kbt = kbt self.molinfo = molinfo self.logfile = logfile self._commandlist = [] self._init_plumed()
def cmd(self, *args, **kwargs): self.interface.cmd(*args, **kwargs) def _init_plumed(self): if self.pathtoplumed != "": # pragma: no cover #os.system("source " + self.pathtoplumed + "/sourceme.sh") sys.path.append(self.pathtoplumed + "/python") warnings.warn("Sourced PLUMED from: " + self.pathtoplumed) else: warnings.warn("Using currently sourced PLUMED from: " + plumed.sys.prefix + "/lib/libplumedKernel.dylib") #os.environ["PLUMED_KERNEL"][:-30]) not set in conda-forge plumed self.cmd("setMDEngine", "python") self.cmd("setTimestep", self.timestep) self.cmd("setKbT", self.kbt) self.cmd("setNatoms", self.topology.n_atoms) self.cmd("setLogFile", self.logfile) self.cmd("init") if (self.molinfo != ""): self.cmd("readInputLine", "MOLINFO STRUCTURE=" + self.molinfo) def set(self, name, definition): """Set a non-outputting command in the `PLUMEDInterface`. Parameters ---------- name : string Equivalent to a `label` in a PLUMED input file. Not required for all commands (can be left empty). definition: string The PLUMED keywords that define the command (see http://www.plumed.org/documentation). """ if (name != ""): self.cmd("readInputLine", name + ": " + definition) else: self.cmd("readInputLine", definition) self._commandlist.append((name, definition)) def get(self): """Get commands set in the `PLUMEDInterface`. Returns ------- list of tuples list of tuples (`label` and `definition`) ran in the `PLUMEDInterface` using the `set` function. """ return tuple(self._commandlist) def to_dict(self): return { 'pathtoplumed': self.pathtoplumed, 'topology': self.topology, 'timestep': self.timestep, 'kbt': self.kbt, 'molinfo': self.molinfo, 'logfile': self.logfile, '_commandlist': self._commandlist } @classmethod def from_dict(cls, dct): pathtoplumed = dct['pathtoplumed'] topology = dct['topology'] timestep = dct['timestep'] kbt = dct['kbt'] molinfo = dct['molinfo'] logfile = dct['logfile'] obj = cls(pathtoplumed=pathtoplumed, topology=topology, timestep=timestep, kbt=kbt, molinfo=molinfo, logfile=logfile) _commandlist = dct['_commandlist'] for name, definition in _commandlist: obj.set(name, definition) return obj