Coverage for cvpack/composite_rmsd.py: 86%
22 statements
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-09 16:14 +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"""
10import typing as t
12from openmm import unit as mmunit
14from .base_rmsd import BaseRMSD
15from .units import MatrixQuantity, VectorQuantity
18class _Stub: # pylint: disable=too-few-public-methods
19 def __init__(self, _):
20 raise ImportError(
21 "CompositeRMSD requires the mdtools::openmm-cpp-forces conda package"
22 )
25try:
26 from openmmcppforces import CompositeRMSDForce
27except ImportError:
28 CompositeRMSDForce = _Stub
31class CompositeRMSD(BaseRMSD, CompositeRMSDForce):
32 r"""
33 The composite root-mean-square deviation (RMSD) between the current and reference
34 coordinates of :math:`m` groups of atoms:
36 .. math::
38 d_{\rm crms}({\bf r}) = \sqrt{
39 \frac{1}{n} \min_{
40 \bf q \in \mathbb{R}^4 \atop \|{\bf q}\| = 1
41 } \sum_{j=1}^m \sum_{i \in {\bf g}_j} \left\|
42 {\bf A}({\bf q})\left({\bf r}_i - {\bf c}_j\right) -
43 {\bf r}_i^{\rm ref} + {\bf c}_j^{\rm ref}
44 \right\|^2
45 }
47 where each group :math:`{\bf g}_j` is a set of :math:`n_j` atom indices,
48 :math:`n = \sum_{j=1}^m n_j` is the total number of atoms in these groups,
49 :math:`{\bf A}(\bf q)` is the rotation matrix corresponding to a unit quaternion
50 :math:`{\bf q}`, :math:`{\bf r}_i` and :math:`{\bf r}_i^{\rm ref}` are the
51 positions of atom :math:`i` in the current and reference structures, respectively,
52 :math:`{\bf c}_j` is the position of atom :math:`i`, given by
54 .. math::
56 {\bf c}_j = \frac{1}{n_j} \sum_{i \in {\bf g}_j} {\bf r}_i
58 and :math:`{\bf c}_j^{\rm ref}` is the centroid of the reference structure for
59 group :math:`j`, defined analogously to :math:`{\bf c}_j`.
61 .. warning::
63 To use this class, you must install the `openmm-cpp-forces`_ conda package.
65 .. _openmm-cpp-forces: https://anaconda.org/mdtools/openmm-cpp-forces
67 Parameters
68 ----------
69 referencePositions
70 The reference coordinates, which can be either a coordinate matrix or a mapping
71 from atom indices to coordinate vectors. It must contain all atoms in
72 ``groups``, and does not need to contain all atoms in the system. See
73 ``numAtoms`` below.
74 groups
75 A sequence of disjoint atom groups. Each group is a sequence of atom indices.
76 numAtoms
77 The total number of atoms in the system, including those that are not in
78 ``groups``. This argument is necessary only if ``referencePositions`` does not
79 contain all atoms in the system.
80 name
81 The name of the collective variable.
83 Raises
84 ------
85 ImportError
86 If the `openmm-cpp-forces`_ conda package is not installed.
87 ValueError
88 If ``groups`` is not a sequence of disjoint atom groups.
90 Example
91 -------
92 >>> import cvpack
93 >>> import openmm as mm
94 >>> import pytest
95 >>> from openmmtools import testsystems
96 >>> from openmm import unit
97 >>> model = testsystems.HostGuestVacuum()
98 >>> host_atoms, guest_atoms = (
99 ... [a.index for a in r.atoms()]
100 ... for r in model.topology.residues()
101 ... )
102 >>> try:
103 ... composite_rmsd = cvpack.CompositeRMSD(
104 ... model.positions,
105 ... [host_atoms, guest_atoms],
106 ... )
107 ... except ImportError:
108 ... pytest.skip("openmm-cpp-forces is not installed")
109 >>> composite_rmsd.addToSystem(model.system)
110 >>> context = mm.Context(
111 ... model.system,
112 ... mm.VerletIntegrator(1.0 * unit.femtoseconds),
113 ... mm.Platform.getPlatformByName('Reference'),
114 ... )
115 >>> context.setPositions(model.positions)
116 >>> composite_rmsd.getValue(context)
117 0.0 nm
118 >>> model.positions[guest_atoms] += 1.0 * unit.nanometers
119 >>> context.setPositions(model.positions)
120 >>> composite_rmsd.getValue(context)
121 0.0 nm
122 """
124 def __init__(
125 self,
126 referencePositions: t.Union[MatrixQuantity, t.Dict[int, VectorQuantity]],
127 groups: t.Iterable[t.Iterable[int]],
128 numAtoms: t.Optional[int] = None,
129 name: str = "composite_rmsd",
130 ) -> None:
131 num_atoms = numAtoms or len(referencePositions)
132 groups = [[int(atom) for atom in group] for group in groups]
133 defined_coords = self._getDefinedCoords(referencePositions, sum(groups, []))
134 all_coords = self._getAllCoords(defined_coords, num_atoms)
135 super().__init__(all_coords)
136 for group in groups:
137 self.addGroup(group)
138 self._registerCV(
139 name,
140 mmunit.nanometers,
141 referencePositions=defined_coords,
142 groups=groups,
143 numAtoms=num_atoms,
144 )
147CompositeRMSD.registerTag("!cvpack.CompositeRMSD")