Coverage for cvpack/residue_coordination.py: 100%
45 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 00:10 +0000
1"""
2.. class:: ResidueCoordination
3 :platform: Linux, MacOS, Windows
4 :synopsis: The number of contacts between two groups of residues
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10from __future__ import annotations
12import typing as t
14import openmm
15from openmm import unit as mmunit
16from openmm.app.topology import Residue
18from .collective_variable import CollectiveVariable
19from .units import Quantity, ScalarQuantity, value_in_md_units
22class ResidueCoordination(CollectiveVariable, openmm.CustomCentroidBondForce):
23 r"""
24 The number of contacts between two disjoint groups of residues:
26 .. math::
27 N({\bf r}) = \sum_{i \in {\bf G}_1} \sum_{j \in {\bf G}_2} S\left(
28 \frac{\|{\bf R}_j({\bf r}) - {\bf R}_i({\bf r})\|}{r_0}
29 \right)
31 where :math:`{\bf G}_1` and :math:`{\bf G}_2` are the two groups of residues,
32 :math:`{\bf R}_i` is the centroid of the residue :math:`i`, :math:`r_0` is the
33 threshold distance for defining a contact and :math:`S(x)` is a step function equal
34 to :math:`1` if a contact is made or equal to :math:`0` otherwise. For trajectory
35 analysis, it is fine to set :math:`S(x) = H(1-x)`, where `H` is the `Heaviside step
36 function <https://en.wikipedia.org/wiki/Heaviside_step_function>`_. For molecular
37 dynamics, however, :math:`S(x)` should be a continuous approximation of
38 :math:`H(1-x)` for :math:`x \geq 0`. By default :cite:`Iannuzzi_2003`, the
39 following function is used:
41 .. math::
43 S(x) = \frac{1-x^6}{1-x^{12}} = \frac{1}{1+x^6}
45 Parameters
46 ----------
47 residueGroup1
48 The residues in the first group.
49 residueGroup2
50 The residues in the second group.
51 pbc
52 Whether to use periodic boundary conditions in distance calculations. This
53 should be ``True`` for periodic systems.
54 stepFunction
55 The function "step(1-x)" (for analysis only) or a continuous approximation
56 thereof.
57 thresholdDistance
58 The threshold distance (:math:`r_0`) for considering two residues as being
59 in contact.
60 normalize
61 Whether the number of contacts should be normalized by the total number of
62 possible contacts.
63 weighByMass
64 Whether the centroid of each residue should be weighted by the mass of the
65 atoms in the residue.
66 includeHydrogens
67 Whether hydrogen atoms should be included in the centroid calculations.
68 name
69 The name of the collective variable.
71 Raises
72 ------
73 ValueError
74 If the two groups of residues are not disjoint
76 Example
77 -------
78 >>> import cvpack
79 >>> import itertools as it
80 >>> import openmm
81 >>> from openmm import app, unit
82 >>> from openmmtools import testsystems
83 >>> model = testsystems.LysozymeImplicit()
84 >>> group1 = list(it.islice(model.topology.residues(), 125, 142))
85 >>> print(*[r.name for r in group1]) # doctest: +ELLIPSIS
86 TRP ASP GLU ... ASN GLN THR
87 >>> group2 = list(it.islice(model.topology.residues(), 142, 156))
88 >>> print(*[r.name for r in group2]) # doctest: +ELLIPSIS
89 PRO ASN ARG ... ARG THR GLY
90 >>> residue_coordination = cvpack.ResidueCoordination(
91 ... group1,
92 ... group2,
93 ... stepFunction="step(1-x)",
94 ... thresholdDistance=0.6*unit.nanometers,
95 ... )
96 >>> residue_coordination.addToSystem(model.system)
97 >>> platform = openmm.Platform.getPlatformByName('Reference')
98 >>> context = openmm.Context(
99 ... model.system, openmm.CustomIntegrator(0), platform
100 ... )
101 >>> context.setPositions(model.positions)
102 >>> residue_coordination.getValue(context)
103 7.0 dimensionless
104 >>> residue_coordination.setReferenceValue(7 * unit.dimensionless)
105 >>> context.reinitialize(preserveState=True)
106 >>> residue_coordination.getValue(context)
107 0.99999... dimensionless
108 """
110 def __init__(
111 self,
112 residueGroup1: t.Iterable[Residue],
113 residueGroup2: t.Iterable[Residue],
114 pbc: bool = False,
115 stepFunction: str = "1/(1+x^6)",
116 thresholdDistance: ScalarQuantity = Quantity(1.0 * mmunit.nanometers),
117 normalize: bool = False,
118 weighByMass: bool = True,
119 includeHydrogens: bool = True,
120 name: str = "residue_coordination",
121 ) -> None:
122 residueGroup1 = list(residueGroup1)
123 residueGroup2 = list(residueGroup2)
124 if {r.index for r in residueGroup1} & {r.index for r in residueGroup2}:
125 raise ValueError("The two groups of residues must be disjoint")
126 nr1 = len(residueGroup1)
127 nr2 = len(residueGroup2)
128 self._ref_val = nr1 * nr2 if normalize else 1.0
129 threshold = thresholdDistance
130 if mmunit.is_quantity(threshold):
131 threshold = threshold.value_in_unit(mmunit.nanometers)
132 expression = (
133 f"({stepFunction})/refval"
134 f"; x=distance(g1,g2)/{threshold}"
135 f"; refval={self._ref_val}"
136 )
137 super().__init__(2, expression)
138 self.setUsesPeriodicBoundaryConditions(pbc)
139 for res in residueGroup1 + residueGroup2:
140 atoms = [
141 atom.index
142 for atom in res.atoms()
143 if includeHydrogens
144 or atom.element is None
145 or atom.element.atomic_number != 1 # hydrogen or deuterium
146 ]
147 self.addGroup(atoms, *([] if weighByMass else [[1] * len(atoms)]))
148 for idx1 in range(nr1):
149 for idx2 in range(nr1, nr1 + nr2):
150 self.addBond([idx1, idx2], [])
151 self._registerCV(
152 name,
153 mmunit.dimensionless,
154 residueGroup1=residueGroup1,
155 residueGroup2=residueGroup2,
156 pbc=pbc,
157 stepFunction=stepFunction,
158 thresholdDistance=thresholdDistance,
159 normalize=normalize,
160 weighByMass=weighByMass,
161 includeHydrogens=includeHydrogens,
162 )
164 def __getstate__(self) -> dict[str, t.Any]:
165 return {**super().__getstate__(), "referenceValue": self._ref_val}
167 def __setstate__(self, keywords: dict[str, t.Any]) -> None:
168 keywords = dict(keywords)
169 reference_value = keywords.pop("referenceValue", None)
170 super().__setstate__(keywords)
171 if reference_value is not None:
172 self.setReferenceValue(reference_value)
174 def getReferenceValue(self) -> mmunit.Quantity:
175 """
176 Get the reference value used for normalizing the residue coordination.
178 Returns
179 -------
180 mmunit.Quantity
181 The reference value.
182 """
183 return self._ref_val * self.getUnit()
185 def setReferenceValue(self, value: ScalarQuantity) -> None:
186 """
187 Set the reference value used for normalizing the residue coordination.
189 Parameters
190 ----------
191 value
192 The reference value.
193 """
194 expression = self.getEnergyFunction()
195 value = float(value_in_md_units(value))
196 self.setEnergyFunction(
197 expression.replace(f"refval={self._ref_val}", f"refval={value}")
198 )
199 self._ref_val = value
202ResidueCoordination.registerTag("!cvpack.ResidueCoordination")