Coverage for cvpack/centroid_function.py: 89%
28 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 00:10 +0000
1"""
2.. class:: CentroidFunction
3 :platform: Linux, MacOS, Windows
4 :synopsis: A generic function of the centroids of groups of atoms
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10from __future__ import annotations
12import typing as t
14import numpy as np
15import openmm
16from numpy.typing import ArrayLike
17from openmm import unit as mmunit
19from .base_custom_function import BaseCustomFunction
20from .units import ScalarQuantity, VectorQuantity
23class CentroidFunction(BaseCustomFunction, openmm.CustomCentroidBondForce):
24 r"""
25 A generic function of the centroids of :math:`m \times n` atom groups split
26 into `m` collections of `n` groups:
28 .. math::
30 f({\bf r}) = \sum_{i=1}^m F\Big(
31 {\bf g}^i_1({\bf r}),
32 {\bf g}^i_2({\bf r}),
33 \dots,
34 {\bf g}^i_n({\bf r})
35 \Big)
37 where :math:`F` is a user-defined function and :math:`{\bf g}^i_j({\bf r})` is the
38 centroid of the :math:`j`-th group of atoms of the :math:`i`-th collection of
39 groups.
41 The centroid of a group of atoms is defined as:
43 .. math::
45 {\bf g}_j({\bf r}) = \frac{1}{N_j} \sum_{k=1}^{N_j} {\bf r}_{k,j}
47 where :math:`N_j` is the number of atoms in group :math:`j` and
48 :math:`{\bf r}_{k,j}` is the coordinate of the :math:`k`-th atom of group
49 :math:`j`. Optionally (and by default), the centroid can be weighted by the mass of
50 each atom in the group. In this case, it is redefined as:
52 .. math::
54 {\bf g}_j({\bf r}) = \frac{1}{M_j} \sum_{k=1}^{N_j} m_{k,j} {\bf r}_{k,j}
56 where :math:`M_j` is the total mass of atom group :math:`j` and :math:`m_{k,j}` is
57 the mass of the :math:`k`-th atom in group :math:`j`.
59 The function :math:`F` is defined as a string and can be any expression supported
60 by :OpenMM:`CustomCentroidBondForce`. If the expression contains named parameters,
61 the value of each parameter can be passed in one of three ways:
63 #. By a semicolon-separated definition in the function string, such as described
64 in the :OpenMM:`CustomCentroidBondForce` documentation. In this case, the
65 parameter value will be the same for all collections of atom groups.
67 #. By a 1D array or list of length :math:`m` passed as a keyword argument to
68 the :class:`CentroidFunction` constructor. In this case, each value will be
69 assigned to a different collection of atom groups.
71 #. By a scalar passed as a keyword argument to the :class:`CentroidFunction`
72 constructor. In this case, the parameter will apply to all collections of atom
73 groups and will become available for on-the-fly modification during a simulation
74 via the ``setParameter`` method of an :OpenMM:`Context` object. **Warning**:
75 other collective variables or :OpenMM:`Force` objects in the same system will
76 share the same values of equal-named parameters.
78 Parameters
79 ----------
80 function
81 The function to be evaluated. It must be a valid
82 :OpenMM:`CustomCentroidBondForce` expression.
83 unit
84 The unit of measurement of the collective variable. It must be compatible
85 with the MD unit system (mass in `daltons`, distance in `nanometers`, time
86 in `picoseconds`, temperature in `kelvin`, energy in `kilojoules_per_mol`,
87 angle in `radians`). If the collective variable does not have a unit, use
88 `dimensionless`.
89 groups
90 The groups of atoms to be used in the function. Each group must be specified
91 as a list of atom indices with arbitrary length.
92 collections
93 The indices of the groups in each collection, passed as a 2D array-like object
94 of shape `(m, n)`, where `m` is the number of collections and `n` is the number
95 of groups per collection. If a 1D object is passed, it is assumed that `m` is 1
96 and `n` is the length of the object.
97 periodicBounds
98 The periodic bounds of the collective variable if it is periodic, or `None` if
99 it is not.
100 pbc
101 Whether to use periodic boundary conditions. This should be ``True`` for
102 periodic systems.
103 weighByMass
104 Whether to define the centroid as the center of mass of the group instead of
105 the geometric center.
106 name
107 The name of the collective variable.
109 Keyword Args
110 ------------
111 **parameters
112 The named parameters of the function. Each parameter can be a 1D array-like
113 object or a scalar. In the former case, the array must have length :math:`m`
114 and each entry will be assigned to a different collection of atom groups. In
115 the latter case, it will define a global :OpenMM:`Context` parameter.
117 Raises
118 ------
119 ValueError
120 If the collections are not specified as a 1D or 2D array-like object.
121 ValueError
122 If group indices are out of bounds.
123 ValueError
124 If the unit of the collective variable is not compatible with the MD unit
125 system.
127 Example
128 -------
129 >>> import cvpack
130 >>> import numpy as np
131 >>> import openmm
132 >>> from openmm import unit
133 >>> from openmmtools import testsystems
134 >>> model = testsystems.LysozymeImplicit()
135 >>> residues = list(model.topology.residues())
136 >>> atoms = [[a.index for a in r.atoms()] for r in residues]
138 Compute the residue coordination between two helices:
140 >>> res_coord = cvpack.ResidueCoordination(
141 ... residues[115:124],
142 ... residues[126:135],
143 ... stepFunction="step(1-x)",
144 ... )
145 >>> res_coord.addToSystem(model.system)
146 >>> integrator = openmm.VerletIntegrator(0)
147 >>> platform = openmm.Platform.getPlatformByName('Reference')
148 >>> context = openmm.Context(model.system, integrator, platform)
149 >>> context.setPositions(model.positions)
150 >>> res_coord.getValue(context)
151 28.0 dimensionless
153 Recompute the residue coordination using the centroid function:
155 >>> groups = [atoms[115:124], atoms[126:135]]
156 >>> colvar = cvpack.CentroidFunction(
157 ... "step(1 - distance(g1, g2))",
158 ... unit.dimensionless,
159 ... atoms[115:124] + atoms[126:135],
160 ... [[i, j] for i in range(9) for j in range(9, 18)],
161 ... )
162 >>> colvar.addToSystem(model.system)
163 >>> context.reinitialize(preserveState=True)
164 >>> colvar.getValue(context)
165 28.0 dimensionless
166 """
168 def __init__(
169 self,
170 function: str,
171 unit: mmunit.Unit,
172 groups: t.Iterable[t.Iterable[int]],
173 collections: ArrayLike | None = None,
174 periodicBounds: VectorQuantity | None = None,
175 pbc: bool = False,
176 weighByMass: bool = True,
177 name: str = "centroid_function",
178 **parameters: ScalarQuantity | VectorQuantity,
179 ) -> None:
180 groups = [[int(atom) for atom in group] for group in groups]
181 num_groups = len(groups)
182 collections = np.atleast_2d(
183 np.arange(num_groups) if collections is None else collections
184 )
185 num_collections, groups_per_collection, *others = collections.shape
186 if others:
187 raise ValueError("Array `collections` cannot have more than 2 dimensions")
188 if np.any(collections < 0) or np.any(collections >= num_groups):
189 raise ValueError("Group index out of bounds")
190 super().__init__(groups_per_collection, function)
191 for group in groups:
192 self.addGroup(group, *([] if weighByMass else [[1] * len(group)]))
193 overalls, perbonds = self._extractParameters(num_collections, **parameters)
194 self._addParameters(overalls, perbonds, collections, pbc)
195 collections = [[int(atom) for atom in collection] for collection in collections]
196 self._registerCV(
197 name,
198 unit,
199 function=function,
200 unit=unit,
201 groups=groups,
202 collections=collections,
203 periodicBounds=periodicBounds,
204 pbc=pbc,
205 weighByMass=weighByMass,
206 **overalls,
207 **perbonds,
208 )
209 if periodicBounds is not None:
210 self._registerPeriodicBounds(*periodicBounds)
213CentroidFunction.registerTag("!cvpack.CentroidFunction")