Improves NumberOfContacts API and DOC

cvpack/number_of_contacts.py

@@ -7,7 +7,7 @@
 
 """
 
-from typing import Sequence
+import typing as t
 
 import openmm
 
@@ -24,26 +24,27 @@ class NumberOfContacts(openmm.CustomNonbondedForce, AbstractCollectiveVariable):
         N({\\bf r}) = \\sum_{i \\in {\\bf g}_1} \\sum_{j \\in {\\bf g}_2}
                         S\\left(\\frac{\\|{\\bf r}_j - {\\bf r}_i\\|}{r_0}\\right)
 
-    where :math:`r_0` is the threshold distance for defining a contact and :math:`S(x)` is a step
-    function equal to 1 if a contact is made or equal to 0 otherwise. In analysis, it is fine to
-    make :math:`S(x) = H(1-x)`, where `H` is the `Heaviside step function
-    <https://en.wikipedia.org/wiki/Heaviside_step_function>`_. In a simulation, however,
-    :math:`S(x)` should continuously approximate :math:`H(1-x)` for :math:`x \\geq 0`. By default
-    :cite:`Iannuzzi_2003`,
+    where :math:`r_0` is the threshold distance for defining a contact and :math:`S(x)` is
+    a step function equal to :math:`1` if a contact is made or equal to :math:`0` otherwise.
+    For trajectory analysis, it is fine to make :math:`S(x) = H(1-x)`, where `H` is the
+    `Heaviside step function <https://en.wikipedia.org/wiki/Heaviside_step_function>`_. For
+    molecular dynamics, however, :math:`S(x)` should be a continuous approximation of
+    :math:`H(1-x)` for :math:`x \\geq 0`. By default :cite:`Iannuzzi_2003`, the following
+    function is used:
 
     .. math::
 
         S(x) = \\frac{1-x^6}{1-x^{12}} = \\frac{1}{1+x^6}
 
-    Atom pairs are ignored for distances beyond a cutoff :math:`r_c`. To avoid discontinuities,
-    a switching function is applied at :math:`r_s \\leq r \\leq r_c` to make :math:`S(r/r_0)`
-    smoothly decay to zero.
+    In fact, a cutoff distance :math:`r_c = x_c r_0` (typically, :math:`x_c = 2`) is applied
+    so that :math:`S(x) = 0` for :math:`x \\geq x_c`. To avoid discontinuities, there is also
+    the option to smoothly switch off :math:`S(x)` starting from :math:`r_s = x_s r_0`
+    (typically, :math:`x_s = 1.5`) instead of doing it abruptly at :math:`r_c`.
 
     .. note::
 
-        The two groups are allowed to overlap. In this case, terms with :math:`j = i`
-        (self-contacts) are ignored and each combination with :math:`j \\neq i` is counted
-        only once.
+        Atoms are allowed to be in both groups. In this case, self-contacts (:math:`i = j`)
+        are ignored and each pair of distinct atoms (:math:`i \\neq j`) is counted only once.
 
     Parameters
     ----------
@@ -59,11 +60,15 @@ class NumberOfContacts(openmm.CustomNonbondedForce, AbstractCollectiveVariable):
             The function "step(1-x)" (for analysis only) or a continuous approximation
             thereof
         thresholdDistance
-            The threshold distance for considering two atoms as being in contact
-        cutoffDistance
-            The distance beyond which an atom pair will be ignored
-        switchingDistance
-            The distance beyond which a switching function will be applied
+            The threshold distance (:math:`r_0`) for considering two atoms as being in
+            contact
+        cutoffFactor
+            The factor :math:`x_c` that multiplies the threshold distance to define
+            the cutoff distance
+        switchFactor
+            The factor :math:`x_s` that multiplies the threshold distance to define
+            the distance at which the step function starts switching off smoothly.
+            If None, it switches off abruptly at the cutoff distance.
 
     Example
     -------
@@ -89,27 +94,27 @@ class NumberOfContacts(openmm.CustomNonbondedForce, AbstractCollectiveVariable):
     @mmunit.convert_quantities
     def __init__(  # pylint: disable=too-many-arguments
         self,
-        group1: Sequence[int],
-        group2: Sequence[int],
+        group1: t.Sequence[int],
+        group2: t.Sequence[int],
         numAtoms: int,
-        pbc: bool = False,
+        pbc: bool,
         stepFunction: str = "1/(1+x^6)",
         thresholdDistance: mmunit.ScalarQuantity = mmunit.Quantity(
             0.3, mmunit.nanometers
         ),
-        cutoffDistance: mmunit.ScalarQuantity = mmunit.Quantity(0.6, mmunit.nanometers),
-        switchingDistance: mmunit.ScalarQuantity = mmunit.Quantity(
-            0.5, mmunit.nanometers
-        ),
+        cutoffFactor: float = 2.0,
+        switchFactor: t.Optional[float] = 1.5,
     ) -> None:
         super().__init__(stepFunction + f"; x=r/{thresholdDistance}")
         nonbonded_method = self.CutoffPeriodic if pbc else self.CutoffNonPeriodic
         self.setNonbondedMethod(nonbonded_method)
         for _ in range(numAtoms):
             self.addParticle([])
-        self.setUseSwitchingFunction(True)
-        self.setCutoffDistance(cutoffDistance)
-        self.setSwitchingDistance(switchingDistance)
+        self.setCutoffDistance(cutoffFactor * thresholdDistance)
+        use_switching_function = switchFactor is not None
+        self.setUseSwitchingFunction(use_switching_function)
+        if use_switching_function:
+            self.setSwitchingDistance(switchFactor * thresholdDistance)
         self.setUseLongRangeCorrection(False)
         self.addInteractionGroup(group1, group2)
         self._registerCV(
@@ -120,6 +125,6 @@ def __init__(  # pylint: disable=too-many-arguments
             pbc,
             stepFunction,
             thresholdDistance,
-            cutoffDistance,
-            switchingDistance,
+            cutoffFactor,
+            switchFactor,
         )

cvpack/tests/test_cvpack.py

@@ -262,7 +262,13 @@ def test_number_of_contacts():
     distances = np.array([np.linalg.norm(pos[i] - pos[j]) for i, j in pairs])
     contacts = np.where(distances <= 0.6, 1 / (1 + (distances / 0.3) ** 6), 0)
     num_atoms = model.topology.getNumAtoms()
-    number_of_contacts = cvpack.NumberOfContacts(group1, group2, num_atoms, pbc=False)
+    number_of_contacts = cvpack.NumberOfContacts(
+        group1,
+        group2,
+        num_atoms,
+        pbc=False,
+        switchFactor=None,
+    )
     model.system.addForce(number_of_contacts)
     integrator = openmm.CustomIntegrator(0)
     platform = openmm.Platform.getPlatformByName("Reference")