Coverage for cvpack/composite_rmsd.py: 92%
38 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 00:10 +0000
1"""
2.. class:: CompositeRMSD
3 :platform: Linux, MacOS, Windows
4 :synopsis: Deviation of multiple corotating bodies from their reference structures
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10from __future__ import annotations
12import typing as t
14from openmm import unit as mmunit
16from .base_rmsd import BaseRMSD
17from .units import MatrixQuantity, VectorQuantity
19# Earlier versions compute wrong forces for groups of different sizes
20MIN_OPENMM_CPP_FORCES_VERSION = (0, 2, 0)
23def _stub(message: str) -> type:
24 class Stub:
25 def __init__(self, *_):
26 raise ImportError(message)
28 return Stub
31def _get_force_class() -> type:
32 """Return CompositeRMSDForce, or a stub that raises an informative ImportError."""
33 requirement = "openmm-cpp-forces 0.2.0 or later (mdtools conda channel or PyPI)"
34 try:
35 import openmmcppforces # noqa: PLC0415
36 except ImportError:
37 return _stub(f"CompositeRMSD requires {requirement}")
38 version = getattr(openmmcppforces, "__version__", None)
39 numbers = tuple(
40 int(part) for part in (version or "").split(".")[:3] if part.isdigit()
41 )
42 if numbers < MIN_OPENMM_CPP_FORCES_VERSION:
43 return _stub(
44 f"CompositeRMSD requires {requirement}, because earlier versions compute "
45 f"wrong forces, but version {version or 'unknown'} is installed"
46 )
47 return openmmcppforces.CompositeRMSDForce
50CompositeRMSDForce = _get_force_class()
53class CompositeRMSD(BaseRMSD, CompositeRMSDForce):
54 r"""
55 The composite root-mean-square deviation (RMSD) between the current and reference
56 coordinates of :math:`m` groups of atoms:
58 .. math::
60 d_{\rm crms}({\bf r}) = \sqrt{
61 \frac{1}{n} \min_{
62 \bf q \in \mathbb{R}^4 \atop \|{\bf q}\| = 1
63 } \sum_{j=1}^m \sum_{i \in {\bf g}_j} \left\|
64 {\bf A}({\bf q})\left({\bf r}_i - {\bf c}_j\right) -
65 {\bf r}_i^{\rm ref} + {\bf c}_j^{\rm ref}
66 \right\|^2
67 }
69 where each group :math:`{\bf g}_j` is a set of :math:`n_j` atom indices,
70 :math:`n = \sum_{j=1}^m n_j` is the total number of atoms in these groups,
71 :math:`{\bf A}(\bf q)` is the rotation matrix corresponding to a unit quaternion
72 :math:`{\bf q}`, :math:`{\bf r}_i` and :math:`{\bf r}_i^{\rm ref}` are the
73 positions of atom :math:`i` in the current and reference structures, respectively,
74 :math:`{\bf c}_j` is the centroid of group :math:`j`, given by
76 .. math::
78 {\bf c}_j = \frac{1}{n_j} \sum_{i \in {\bf g}_j} {\bf r}_i
80 and :math:`{\bf c}_j^{\rm ref}` is the centroid of the reference structure for
81 group :math:`j`, defined analogously to :math:`{\bf c}_j`.
83 .. warning::
85 To use this class, you must install the `openmm-cpp-forces`_ package,
86 version 0.2.0 or later.
88 .. _openmm-cpp-forces: https://anaconda.org/mdtools/openmm-cpp-forces
90 Parameters
91 ----------
92 referencePositions
93 The reference coordinates, which can be either a coordinate matrix or a mapping
94 from atom indices to coordinate vectors. It must contain all atoms in
95 ``groups``, and does not need to contain all atoms in the system. See
96 ``numAtoms`` below.
97 groups
98 A sequence of disjoint atom groups. Each group is a sequence of atom indices.
99 numAtoms
100 The total number of atoms in the system, including those that are not in
101 ``groups``. This argument is necessary only if ``referencePositions`` does not
102 contain all atoms in the system.
103 name
104 The name of the collective variable.
106 Raises
107 ------
108 ImportError
109 If the `openmm-cpp-forces`_ package (version 0.2.0 or later) is not
110 installed.
111 ValueError
112 If ``groups`` is not a sequence of disjoint atom groups.
113 ValueError
114 If referencePositions is a dict and numAtoms is not given
116 Example
117 -------
118 >>> import cvpack
119 >>> import openmm as mm
120 >>> import pytest
121 >>> from openmmtools import testsystems
122 >>> from openmm import unit
123 >>> model = testsystems.HostGuestVacuum()
124 >>> host_atoms, guest_atoms = (
125 ... [a.index for a in r.atoms()]
126 ... for r in model.topology.residues()
127 ... )
128 >>> try:
129 ... composite_rmsd = cvpack.CompositeRMSD(
130 ... model.positions,
131 ... [host_atoms, guest_atoms],
132 ... )
133 ... except ImportError:
134 ... pytest.skip("openmm-cpp-forces is not installed")
135 >>> composite_rmsd.addToSystem(model.system)
136 >>> context = mm.Context(
137 ... model.system,
138 ... mm.VerletIntegrator(1.0 * unit.femtoseconds),
139 ... mm.Platform.getPlatformByName('Reference'),
140 ... )
141 >>> context.setPositions(model.positions)
142 >>> composite_rmsd.getValue(context)
143 0.0 nm
144 >>> model.positions[guest_atoms] += 1.0 * unit.nanometers
145 >>> context.setPositions(model.positions)
146 >>> composite_rmsd.getValue(context)
147 0.0 nm
148 """
150 def __init__(
151 self,
152 referencePositions: MatrixQuantity | dict[int, VectorQuantity],
153 groups: t.Iterable[t.Iterable[int]],
154 numAtoms: int | None = None,
155 name: str = "composite_rmsd",
156 ) -> None:
157 if numAtoms is None and isinstance(referencePositions, dict):
158 raise ValueError("numAtoms is required when referencePositions is a dict")
159 num_atoms = numAtoms or len(referencePositions)
160 groups = [[int(atom) for atom in group] for group in groups]
161 if not all(groups):
162 # CompositeRMSDForce would interpret an empty group as all atoms
163 raise ValueError("The groups of atoms must not be empty")
164 defined_coords = self._getDefinedCoords(referencePositions, sum(groups, []))
165 all_coords = self._getAllCoords(defined_coords, num_atoms)
166 super().__init__(all_coords)
167 for group in groups:
168 self.addGroup(group)
169 self._registerCV(
170 name,
171 mmunit.nanometers,
172 referencePositions=defined_coords,
173 groups=groups,
174 numAtoms=num_atoms,
175 )
178CompositeRMSD.registerTag("!cvpack.CompositeRMSD")