Coverage for cvpack/number_of_contacts.py: 100%
35 statements
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-09 16:14 +0000
1"""
2.. class:: NumberOfContacts
3 :platform: Linux, MacOS, Windows
4 :synopsis: The number of contacts between two atom groups
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10import typing as t
11from numbers import Real
13import openmm
14from openmm import unit as mmunit
16from .collective_variable import CollectiveVariable
17from .units import Quantity, ScalarQuantity
18from .utils import evaluate_in_context
21class NumberOfContacts(CollectiveVariable, openmm.CustomNonbondedForce):
22 r"""
23 The number of contacts between two atom groups:
25 .. math::
26 N({\bf r}) = \sum_{i \in {\bf g}_1} \sum_{j \in {\bf g}_2}
27 S\left(\frac{\|{\bf r}_j - {\bf r}_i\|}{r_0}\right)
29 where :math:`r_0` is the threshold distance for defining a contact and :math:`S(x)`
30 is a step function equal to :math:`1` if a contact is made or equal to :math:`0`
31 otherwise. For trajectory analysis, it is fine to make :math:`S(x) = H(1-x)`, where
32 `H` is the `Heaviside step function
33 <https://en.wikipedia.org/wiki/Heaviside_step_function>`_. For molecular dynamics,
34 however, :math:`S(x)` should be a continuous approximation of :math:`H(1-x)` for
35 :math:`x \geq 0`. By default :cite:`Iannuzzi_2003`, the following function is used:
37 .. math::
39 S(x) = \frac{1-x^6}{1-x^{12}} = \frac{1}{1+x^6}
41 In fact, a cutoff distance :math:`r_c = x_c r_0` (typically, :math:`x_c = 2`) is
42 applied so that :math:`S(x) = 0` for :math:`x \geq x_c`. To avoid discontinuities,
43 there is also the option to smoothly switch off :math:`S(x)` starting from
44 :math:`r_s = x_s r_0` (typically, :math:`x_s = 1.5`) instead of doing it abruptly at
45 :math:`r_c`.
47 .. note::
49 Atoms are allowed to be in both groups. In this case, self-contacts
50 (:math:`i = j`) are ignored and each pair of distinct atoms (:math:`i \neq j`)
51 is counted only once.
53 .. note::
54 Any non-exclusion exceptions involving atoms in :math:`{\bf g}_1` and
55 :math:`{\bf g}_2` in the provided :class:`openmm.NonbondedForce` are turned
56 into exclusions in this collective variable.
58 Parameters
59 ----------
60 group1
61 The indices of the atoms in the first group.
62 group2
63 The indices of the atoms in the second group.
64 nonbondedForce
65 The :class:`openmm.NonbondedForce` object from which the total number of
66 atoms, the exclusions, and whether to use periodic boundary conditions are
67 taken.
68 reference
69 A dimensionless reference value to which the collective variable should be
70 normalized. One can also provide an :OpenMM:`Context` object from which to
71 obtain the reference number of contacts.
72 stepFunction
73 The function "step(1-x)" (for analysis only) or a continuous approximation
74 thereof.
75 thresholdDistance
76 The threshold distance (:math:`r_0`) for considering two atoms as being in
77 contact.
78 cutoffFactor
79 The factor :math:`x_c` that multiplies the threshold distance to define
80 the cutoff distance.
81 switchFactor
82 The factor :math:`x_s` that multiplies the threshold distance to define
83 the distance at which the step function starts switching off smoothly.
84 If None, it switches off abruptly at the cutoff distance.
85 name
86 The name of the collective variable.
88 Example
89 -------
90 >>> import cvpack
91 >>> from openmm import unit
92 >>> from openmmtools import testsystems
93 >>> model = testsystems.HostGuestExplicit()
94 >>> group1, group2 = [], []
95 >>> for residue in model.topology.residues():
96 ... if residue.name != "HOH":
97 ... group = group1 if residue.name == "B2" else group2
98 ... group.extend(atom.index for atom in residue.atoms())
99 >>> forces = {f.getName(): f for f in model.system.getForces()}
100 >>> nc = cvpack.NumberOfContacts(
101 ... group1,
102 ... group2,
103 ... forces["NonbondedForce"],
104 ... stepFunction="step(1-x)",
105 ... )
106 >>> nc.addToSystem(model.system)
107 >>> platform = openmm.Platform.getPlatformByName("Reference")
108 >>> integrator = openmm.VerletIntegrator(1.0 * mmunit.femtoseconds)
109 >>> context = openmm.Context(model.system, integrator, platform)
110 >>> context.setPositions(model.positions)
111 >>> nc.getValue(context)
112 30.0... dimensionless
113 >>> nc_normalized = cvpack.NumberOfContacts(
114 ... group1,
115 ... group2,
116 ... forces["NonbondedForce"],
117 ... stepFunction="step(1-x)",
118 ... reference=context,
119 ... )
120 >>> nc_normalized.addToSystem(model.system)
121 >>> context.reinitialize(preserveState=True)
122 >>> nc_normalized.getValue(context)
123 0.99999... dimensionless
124 """
126 def __init__(
127 self,
128 group1: t.Sequence[int],
129 group2: t.Sequence[int],
130 nonbondedForce: openmm.NonbondedForce,
131 reference: t.Union[Real, openmm.Context] = 1.0,
132 stepFunction: str = "1/(1+x^6)",
133 thresholdDistance: ScalarQuantity = Quantity(0.3 * mmunit.nanometers),
134 cutoffFactor: float = 2.0,
135 switchFactor: t.Optional[float] = 1.5,
136 name: str = "number_of_contacts",
137 ) -> None:
138 num_atoms = nonbondedForce.getNumParticles()
139 pbc = nonbondedForce.usesPeriodicBoundaryConditions()
140 threshold = thresholdDistance
141 if mmunit.is_quantity(threshold):
142 threshold = threshold.value_in_unit(mmunit.nanometers)
143 expression = f"({stepFunction})/1; x=r/{threshold}"
144 super().__init__(expression)
145 nonbonded_method = self.CutoffPeriodic if pbc else self.CutoffNonPeriodic
146 self.setNonbondedMethod(nonbonded_method)
147 for _ in range(num_atoms):
148 self.addParticle([])
149 for index in range(nonbondedForce.getNumExceptions()):
150 i, j, *_ = nonbondedForce.getExceptionParameters(index)
151 self.addExclusion(i, j)
152 self.setCutoffDistance(cutoffFactor * thresholdDistance)
153 use_switching_function = switchFactor is not None
154 self.setUseSwitchingFunction(use_switching_function)
155 if use_switching_function:
156 self.setSwitchingDistance(switchFactor * thresholdDistance)
157 self.setUseLongRangeCorrection(False)
158 self.addInteractionGroup(group1, group2)
159 if isinstance(reference, openmm.Context):
160 reference = evaluate_in_context(self, reference)
161 self.setEnergyFunction(expression.replace("/1;", f"/{reference};"))
162 self._registerCV(
163 name,
164 mmunit.dimensionless,
165 group1=group1,
166 group2=group2,
167 nonbondedForce=nonbondedForce,
168 reference=reference,
169 stepFunction=stepFunction,
170 thresholdDistance=thresholdDistance,
171 cutoffFactor=cutoffFactor,
172 switchFactor=switchFactor,
173 )
176NumberOfContacts.registerTag("!cvpack.NumberOfContacts")