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