Coverage for cvpack/meta_collective_variable.py: 99%
73 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 00:10 +0000
1"""
2.. class:: MetaCollectiveVariable
3 :platform: Linux, MacOS, Windows
4 :synopsis: A function of other collective variables
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10from __future__ import annotations
12import typing as t
13from copy import copy
15import openmm
16from openmm import unit as mmunit
18from .collective_variable import CollectiveVariable
19from .units import Quantity, ScalarQuantity, VectorQuantity, in_md_units
20from .utils import compute_effective_mass, get_single_force_state
23class MetaCollectiveVariable(CollectiveVariable, openmm.CustomCVForce):
24 r"""
25 A collective variable that is a function of other collective variables:
27 .. math::
29 f({\bf r}) = F\left(c_1({\bf r}), c_2({\bf r}), \ldots, c_n({\bf r})\right)
31 where :math:`F(c_1,c_2,\ldots,c_n)` is a user-defined function and
32 :math:`c_i({\bf r})` is the value of the :math:`i`-th collective variable at a
33 given configuration :math:`{\bf r}`.
35 The function :math:`F` is defined as a string and can be any expression supported
36 by the :OpenMM:`CustomCVForce` class. If the expression contains named parameters,
37 the value of each parameter can be passed in one of two ways:
39 #. By a semicolon-separated definition in the function string, such as described
40 in the :OpenMM:`CustomCVForce` documentation. In this case, the parameter value
41 is fixed.
43 #. By a scalar passed as a keyword argument to the
44 :class:`MetaCollectiveVariable` constructor. In this case, the parameter will
45 become available for on-the-fly modification during a simulation via the
46 ``setParameter`` method of an :OpenMM:`Context` object. **Warning**: other
47 collective variables or :OpenMM:`Force` objects in the same system will share
48 the same values of equal-named parameters.
50 Parameters
51 ----------
52 function
53 The function to be evaluated. It must be a valid :OpenMM:`CustomCVForce`
54 expression.
55 variables
56 A sequence of :class:`CollectiveVariable` instances that represent the
57 collective variables on which the meta-collective variable depends. The name
58 of each collective variable must be unique and match a symbol used in the
59 function.
60 unit
61 The unit of measurement of the collective variable. It must be compatible
62 with the MD unit system (mass in `daltons`, distance in `nanometers`, time
63 in `picoseconds`, temperature in `kelvin`, energy in `kilojoules_per_mol`,
64 angle in `radians`). If the collective variable does not have a unit, use
65 `unit.dimensionless`.
66 periodicBounds
67 The periodic bounds of the collective variable if it is periodic, or `None` if
68 it is not periodic.
70 Keyword Args
71 ------------
72 **parameters
73 The named parameters of the function, if any. They will become settable context
74 parameters when this meta-collective variable is added to an :OpenMM:`System`.
75 The passed objects must be scalar quantities. Their values will be converted to
76 OpenMM's MD unit system to serve as default values for the context parameters.
78 Example
79 -------
80 >>> import cvpack
81 >>> import math
82 >>> import openmm
83 >>> from openmm import unit
84 >>> from openmmtools import testsystems
85 >>> model = testsystems.AlanineDipeptideVacuum()
86 >>> phi = cvpack.Torsion(4, 6, 8, 14, name="phi")
87 >>> driving_force = cvpack.MetaCollectiveVariable(
88 ... f"0.5*kappa*min(delta,{2*math.pi}-delta)^2"
89 ... "; delta=abs(phi-phi0)",
90 ... [phi],
91 ... unit.kilojoules_per_mole,
92 ... kappa = 1e3 * unit.kilojoules_per_mole/unit.radian**2,
93 ... phi0 = 120 * unit.degrees
94 ... )
95 >>> driving_force.addToSystem(model.system)
96 >>> integrator = openmm.VerletIntegrator(0)
97 >>> platform = openmm.Platform.getPlatformByName('Reference')
98 >>> context = openmm.Context(model.system, integrator, platform)
99 >>> context.setPositions(model.positions)
100 >>> driving_force.getParameterDefaultValues()
101 {'kappa': 1000.0 kJ/(mol rad**2), 'phi0': 2.094... rad}
102 >>> driving_force.getParameterValues(context)
103 {'kappa': 1000.0 kJ/(mol rad**2), 'phi0': 2.094... rad}
104 >>> driving_force.getValue(context)
105 548.3... kJ/mol
106 >>> driving_force.getInnerValues(context)
107 {'phi': -3.14... rad}
108 >>> driving_force.getInnerEffectiveMasses(context)
109 {'phi': 0.04795... nm**2 Da/(rad**2)}
110 >>> driving_force.getParameterDerivatives(context)
111 {'kappa': 0.548... rad**2, 'phi0': -1047.20... kJ/(mol rad)}
112 """
114 def __init__(
115 self,
116 function: str,
117 variables: t.Iterable[CollectiveVariable],
118 unit: mmunit.Unit,
119 periodicBounds: VectorQuantity | None = None,
120 name: str = "meta_collective_variable",
121 **parameters: ScalarQuantity,
122 ) -> None:
123 super().__init__(function)
124 variables = list(variables)
125 names = [cv.getName() for cv in variables]
126 if len(set(names)) != len(names):
127 raise ValueError(f"The names of the variables must be unique: {names}")
128 self._cvs = tuple(map(copy, variables))
129 self._parameters = {k: in_md_units(v) for k, v in parameters.items()}
130 for cv in self._cvs:
131 self.addCollectiveVariable(cv.getName(), cv)
132 for parameter, value in self._parameters.items():
133 self.addGlobalParameter(parameter, value / value.unit)
134 self.addEnergyParameterDerivative(parameter)
135 self._registerCV(
136 name,
137 unit,
138 function=function,
139 variables=variables,
140 unit=unit,
141 periodicBounds=periodicBounds,
142 **self._parameters,
143 )
144 if periodicBounds is not None:
145 self._registerPeriodicBounds(*periodicBounds)
147 def getInnerVariables(self) -> tuple[CollectiveVariable]:
148 """
149 Get the collective variables on which the meta-collective variable depends.
151 Returns
152 -------
153 Tuple[CollectiveVariable]
154 A tuple with the collective variables.
155 """
156 return self._cvs
158 def getInnerValues(self, context: openmm.Context) -> dict[str, Quantity]:
159 """
160 Get the values of the collective variables on which the meta-collective variable
161 depends. The values are returned as a dictionary with the names of the
162 collective variables as keys.
164 Parameters
165 ----------
166 context
167 The context in which the collective variables will be evaluated.
169 Returns
170 -------
171 Dict[str, Quantity]
172 A dictionary with the names of the collective variables as keys and their
173 values as values.
174 """
175 values = self.getCollectiveVariableValues(context)
176 return {
177 cv.getName(): Quantity(value, cv.getUnit())
178 for cv, value in zip(self._cvs, values)
179 }
181 def getInnerEffectiveMasses(self, context: openmm.Context) -> dict[str, Quantity]:
182 """
183 Get the effective masses of the collective variables on which the
184 meta-collective variable depends. The effective masses are calculated from the
185 forces acting on the particles that represent the collective variables.
187 Parameters
188 ----------
189 context
190 The context in which the collective variables will be evaluated.
192 Returns
193 -------
194 Dict[str, Quantity]
195 A dictionary with the names of the collective variables as keys and their
196 effective masses as values.
197 """
198 self.getCollectiveVariableValues(context) # update the inner context positions
199 inner_context = self.getInnerContext(context)
200 forces = inner_context.getSystem().getForces()
201 system = context.getSystem()
202 if any(map(system.isVirtualSite, range(system.getNumParticles()))):
203 # The inner context has no virtual sites, so the forces on them would not
204 # be distributed to the atoms that define them
205 inner_context = self._getContextWithVirtualSites(forces, context)
206 forces = inner_context.getSystem().getForces()
207 masses = [compute_effective_mass(force, inner_context) for force in forces]
208 return {
209 cv.getName(): Quantity(mass, cv.getMassUnit())
210 for cv, mass in zip(self._cvs, masses)
211 }
213 def _getContextWithVirtualSites(
214 self, forces: t.Sequence[openmm.Force], context: openmm.Context
215 ) -> openmm.Context:
216 # Built once per outer context, then updated with its current state
217 cached = getattr(self, "_virtual_site_context", None)
218 if cached is None or cached[0] is not context:
219 cached = (context, self._createContextWithVirtualSites(forces, context))
220 self._virtual_site_context = cached
221 new_context = cached[1]
222 state = context.getState(getPositions=True, getParameters=True)
223 new_context.setPeriodicBoxVectors(*state.getPeriodicBoxVectors())
224 new_context.setPositions(state.getPositions())
225 parameters = state.getParameters()
226 for name in new_context.getParameters():
227 if name in parameters:
228 new_context.setParameter(name, parameters[name])
229 return new_context
231 @staticmethod
232 def _createContextWithVirtualSites(
233 forces: t.Sequence[openmm.Force], context: openmm.Context
234 ) -> openmm.Context:
235 system = openmm.XmlSerializer.deserialize(
236 openmm.XmlSerializer.serialize(context.getSystem())
237 )
238 while system.getNumForces() > 0:
239 system.removeForce(0)
240 for group, force in enumerate(forces):
241 force_copy = openmm.XmlSerializer.deserialize(
242 openmm.XmlSerializer.serialize(force)
243 )
244 force_copy.setForceGroup(group)
245 system.addForce(force_copy)
246 return openmm.Context(
247 system, openmm.VerletIntegrator(1.0), context.getPlatform()
248 )
250 def getParameterDefaultValues(self) -> dict[str, Quantity]:
251 """
252 Get the default values of the named parameters of this meta-collective variable.
254 Returns
255 -------
256 Dict[str, Quantity]
257 A dictionary with the names of the named parameters as keys and their
258 default values as values.
259 """
260 return self._parameters.copy()
262 def getParameterValues(self, context: openmm.Context) -> dict[str, Quantity]:
263 """
264 Get the values of the named parameters of this meta-collective variable. The
265 values are returned as a dictionary with the names of the parameters as keys.
267 Parameters
268 ----------
269 context
270 The context in which the named parameters will be evaluated.
272 Returns
273 -------
274 Dict[str, Quantity]
275 A dictionary with the names of the named parameters as keys and their values
276 as values.
277 """
278 return {
279 name: Quantity(context.getParameter(name), parameter.unit)
280 for name, parameter in self._parameters.items()
281 }
283 def getParameterDerivatives(
284 self,
285 context: openmm.Context,
286 allowReinitialization: bool = False,
287 ) -> dict[str, Quantity]:
288 """
289 Get the derivatives of the named parameters of this meta-collective variable.
291 Returns
292 -------
293 Dict[str, Quantity]
294 A dictionary with the names of the named parameters as keys and their
295 derivatives as values.
296 """
297 state = get_single_force_state(
298 self, context, allowReinitialization, getParameterDerivatives=True
299 )
300 derivatives = state.getEnergyParameterDerivatives()
301 return {
302 name: Quantity(derivatives[name], self._unit / parameter.unit)
303 for name, parameter in self._parameters.items()
304 }
307MetaCollectiveVariable.registerTag("!cvpack.MetaCollectiveVariable")