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