-
Notifications
You must be signed in to change notification settings - Fork 1
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Includes PathCV that uses reference structures as milestones (#93)
* Defined base class for path-related CVs * Fixed doctest * Initial implementation of PathInRMSDSpace * Included doctest * Improved doctest * Fixed base path cv usage * Added unit test * Included missing package * Fixed tests * Fixed dependencies * Using absolute error in assertion * Updated README.md
- Loading branch information
Showing
12 changed files
with
313 additions
and
48 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,71 @@ | ||
| """ | ||
| .. class:: BasePathCV | ||
| :platform: Linux, MacOS, Windows | ||
| :synopsis: A base class for path-related collective variables | ||
| .. classauthor:: Charlles Abreu <craabreu@gmail.com> | ||
| """ | ||
|
|
||
| import typing as t | ||
| from collections import OrderedDict | ||
|
|
||
| import openmm | ||
|
|
||
| from .collective_variable import CollectiveVariable | ||
| from .path import Metric, deviation, progress | ||
|
|
||
|
|
||
| class BasePathCV(openmm.CustomCVForce, CollectiveVariable): | ||
| """ | ||
| A base class for path-related collective variables | ||
| Parameters | ||
| ---------- | ||
| metric | ||
| A measure of progress or deviation with respect to a path in CV space | ||
| sigma | ||
| The width of the Gaussian kernels | ||
| squared_distances | ||
| Expressions for the squared distance to each milestone | ||
| variables | ||
| A dictionary of collective variables used in the expressions for the squared | ||
| distances | ||
| """ | ||
|
|
||
| def __init__( | ||
| self, | ||
| metric: Metric, | ||
| sigma: float, | ||
| squared_distances: t.Sequence[str], | ||
| variables: t.Dict[str, CollectiveVariable], | ||
| ) -> None: | ||
| n = len(squared_distances) | ||
| definitions = OrderedDict( | ||
| {f"x{i}": sqdist for i, sqdist in enumerate(squared_distances)} | ||
| ) | ||
| definitions["lambda"] = 1 / (2 * sigma**2) | ||
| definitions["xmin0"] = "min(x0,x1)" | ||
| for i in range(n - 2): | ||
| definitions[f"xmin{i+1}"] = f"min(xmin{i},x{i+2})" | ||
| for i in range(n): | ||
| definitions[f"w{i}"] = f"exp(lambda*(xmin{n - 2}-x{i}))" | ||
| definitions["wsum"] = "+".join(f"w{i}" for i in range(n)) | ||
| expressions = [f"{key}={value}" for key, value in definitions.items()] | ||
| if metric == progress: | ||
| numerator = "+".join(f"{i}*w{i}" for i in range(1, n)) | ||
| expressions.append(f"({numerator})/({n - 1}*wsum)") | ||
| else: | ||
| expressions.append(f"xmin{n - 2}-log(wsum)/lambda") | ||
| super().__init__("; ".join(reversed(expressions))) | ||
| for name, variable in variables.items(): | ||
| self.addCollectiveVariable(name, variable) | ||
|
|
||
| def _generateName(self, metric: Metric, name: str, kind: str) -> str: | ||
| if metric not in (progress, deviation): | ||
| raise ValueError( | ||
| "Invalid metric. Use 'cvpack.path.progress' or 'cvpack.path.deviation'." | ||
| ) | ||
| if name is None: | ||
| return f"path_{metric.name}_in_{kind}_space" | ||
| return name |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,153 @@ | ||
| """ | ||
| .. class:: PathInRMSDSpace | ||
| :platform: Linux, MacOS, Windows | ||
| :synopsis: A metric of progress or deviation with respect to a path in RMSD space | ||
| .. classauthor:: Charlles Abreu <craabreu@gmail.com> | ||
| """ | ||
|
|
||
| import typing as t | ||
|
|
||
| from openmm import unit as mmunit | ||
|
|
||
| from .base_path_cv import BasePathCV | ||
| from .path import Metric, progress | ||
| from .rmsd import RMSD | ||
| from .units.units import VectorQuantity | ||
|
|
||
|
|
||
| class PathInRMSDSpace(BasePathCV): | ||
| r""" | ||
| A metric of the system's progress (:math:`s`) or deviation (:math:`z`) with | ||
| respect to a path defined by a sequence of :math:`n` milestones defined as | ||
| reference structures :cite:`Branduardi_2007`: | ||
| .. math:: | ||
| s({\bf r}) = \frac{ | ||
| \dfrac{\sum_{i=1}^n i w_i({\bf r})}{\sum_{i=1}^n w_i({\bf r})} - 1 | ||
| }{n-1} | ||
| \quad \text{or} \quad | ||
| z({\bf r}) = - 2 \sigma ^2 \ln \sum_{i=1}^n w_i({\bf r}) | ||
| with :math:`w_i({\bf r})` being a Gaussian kernel centered at the :math:`i`-th | ||
| milestone, i.e., | ||
| .. math:: | ||
| w_i({\bf r}) = \exp\left(\ | ||
| -\frac{d^2_{\rm rms}({\bf r},{\bf r}^{\rm ref}_i)}{2 \sigma^2} | ||
| \right) | ||
| where :math:`d_{\rm rms}({\bf r},{\bf r}^{\rm ref}_i)` is the root-mean-square | ||
| distance between the current system state and the :math:`i`-th reference structure | ||
| and :math:`\sigma` sets the width of the kernels. | ||
| .. note:: | ||
| The kernel width :math:`\sigma` is related to the parameter :math:`\lambda` of | ||
| Ref. :cite:`Branduardi_2007` by :math:`\sigma = \frac{1}{\sqrt{2\lambda}}`. | ||
| Parameters | ||
| ---------- | ||
| metric | ||
| The path-related metric to compute. Use ``cvpack.path.progress`` for | ||
| computing :math:`s({\bf r})` or ``cvpack.path.deviation`` for computing | ||
| :math:`z({\bf r})`. | ||
| milestones | ||
| A sequence of reference structures, each represented as a dictionary mapping | ||
| atom indices to coordinate vectors. | ||
| numAtoms | ||
| The total number of atoms in the system, including those that are not in | ||
| any of the reference structures. | ||
| sigma | ||
| The width of the Gaussian kernels in nanometers | ||
| name | ||
| The name of the collective variable. If not provided, it is set to | ||
| "path_progress_in_rmsd_space" or "path_deviation_in_rmsd_space" depending | ||
| on the metric. | ||
| Raises | ||
| ------ | ||
| ValueError | ||
| The number of milestones is less than 2 | ||
| ValueError | ||
| If the metric is not `cvpack.path.progress` or `cvpack.path.deviation` | ||
| Examples | ||
| -------- | ||
| >>> import cvpack | ||
| >>> import networkx as nx | ||
| >>> import numpy as np | ||
| >>> import openmm | ||
| >>> from openmm import unit | ||
| >>> from openmmtools import testsystems | ||
| >>> from scipy.spatial.transform import Rotation | ||
| >>> model = testsystems.AlanineDipeptideVacuum() | ||
| >>> atom1, atom2 = 8, 14 | ||
| >>> graph = model.mdtraj_topology.to_bondgraph() | ||
| >>> nodes = list(graph.nodes) | ||
| >>> graph.remove_edge(nodes[atom1], nodes[atom2]) | ||
| >>> movable = list(nx.connected_components(graph))[1] | ||
| >>> x = model.positions / model.positions.unit | ||
| >>> x0 = x[atom1, :] | ||
| >>> vector = x[atom2, :] - x0 | ||
| >>> vector /= np.linalg.norm(vector) | ||
| >>> rotation = Rotation.from_rotvec((np.pi / 6) * vector) | ||
| >>> atoms = [nodes.index(atom) for atom in movable] | ||
| >>> frames = [x.copy()] | ||
| >>> for _ in range(6): | ||
| ... x[atoms, :] = x0 + rotation.apply(x[atoms, :] - x0) | ||
| ... frames.append(x.copy()) | ||
| >>> milestones = [ | ||
| ... {i: row for i, row in enumerate(frame)} | ||
| ... for frame in frames | ||
| ... ] | ||
| >>> s, z = [ | ||
| ... cvpack.PathInRMSDSpace( | ||
| ... metric, milestones, len(x), 0.5 * unit.angstrom | ||
| ... ) | ||
| ... for metric in (cvpack.path.progress, cvpack.path.deviation) | ||
| ... ] | ||
| >>> s.addToSystem(model.system) | ||
| >>> z.addToSystem(model.system) | ||
| >>> context = openmm.Context(model.system, openmm.VerletIntegrator(0.001)) | ||
| >>> context.setPositions(model.positions) | ||
| >>> s.getValue(context) | ||
| 0.172... dimensionless | ||
| >>> z.getValue(context) | ||
| -0.004... nm**2 | ||
| """ | ||
|
|
||
| def __init__( # pylint: disable=too-many-branches | ||
| self, | ||
| metric: Metric, | ||
| milestones: t.Sequence[t.Dict[int, VectorQuantity]], | ||
| numAtoms: int, | ||
| sigma: mmunit.Quantity, | ||
| name: t.Optional[str] = None, | ||
| ) -> None: | ||
| name = self._generateName(metric, name, "rmsd") | ||
| if mmunit.is_quantity(sigma): | ||
| sigma = sigma.value_in_unit(mmunit.nanometers) | ||
| n = len(milestones) | ||
| if n < 2: | ||
| raise ValueError("At least two reference structures are required.") | ||
| collective_variables = { | ||
| f"rmsd{i}": RMSD(reference, reference.keys(), numAtoms, name=f"rmsd{i}") | ||
| for i, reference in enumerate(milestones) | ||
| } | ||
| squared_distances = [f"rmsd{i}^2" for i in range(n)] | ||
| super().__init__(metric, sigma, squared_distances, collective_variables) | ||
| self._registerCV( | ||
| name, | ||
| mmunit.dimensionless if metric == progress else mmunit.nanometers**2, | ||
| metric=metric, | ||
| milestones=[{k: list(v) for k, v in m.items()} for m in milestones], | ||
| numAtoms=numAtoms, | ||
| sigma=sigma, | ||
| ) | ||
|
|
||
|
|
||
| PathInRMSDSpace.registerTag("!cvpack.PathInRMSDSpace") |
Oops, something went wrong.