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