Coverage for cvpack/reporting/meta_cv_writer.py: 100%
53 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 00:10 +0000
1"""
2.. class:: MetaCVWriter
3 :platform: Linux, MacOS, Windows
4 :synopsis: A custom writer for reporting meta-collective variable data
6.. moduleauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10import typing as t
12import openmm as mm
13from openmm import unit as mmunit
15from ..meta_collective_variable import MetaCollectiveVariable
16from .custom_writer import CustomWriter
19class MetaCVWriter(CustomWriter):
20 """
21 A custom writer for reporting meta-collective variable data.
23 Parameters
24 ----------
25 metaCV
26 The meta-collective variable whose associated values will be reported.
27 values
28 The names of the inner variables whose values will be reported.
29 emasses
30 The names of the inner variables whose effective masses will be reported.
31 parameters
32 The names of the parameters whose values will be reported.
33 derivatives
34 The names of the parameters with respect to which the derivatives of the
35 meta-collective variable will be reported.
36 allowReinitialization
37 Whether to allow context reinitialization when the meta-collective variable
38 shares its force group with other forces (needed for derivatives only).
40 Raises
41 ------
42 ValueError
43 If any name does not match an inner variable or a parameter
45 Examples
46 --------
47 >>> import cvpack
48 >>> import openmm
49 >>> from math import pi
50 >>> from cvpack.reporting import StateDataReporter, MetaCVWriter
51 >>> from openmm import app, unit
52 >>> from sys import stdout
53 >>> from openmmtools import testsystems
54 >>> model = testsystems.AlanineDipeptideVacuum()
55 >>> phi = cvpack.Torsion(4, 6, 8, 14, name="phi")
56 >>> psi = cvpack.Torsion(6, 8, 14, 16, name="psi")
57 >>> umbrella = cvpack.MetaCollectiveVariable(
58 ... f"0.5*kappa*(min(dphi,{2*pi}-dphi)^2+min(dpsi,{2*pi}-dpsi)^2)"
59 ... "; dphi=abs(phi-phi0); dpsi=abs(psi-psi0)",
60 ... [phi, psi],
61 ... unit.kilojoules_per_mole,
62 ... name="umbrella",
63 ... kappa=100 * unit.kilojoules_per_mole/unit.radian**2,
64 ... phi0=5*pi/6 * unit.radian,
65 ... psi0=-5*pi/6 * unit.radian,
66 ... )
67 >>> reporter = StateDataReporter(
68 ... stdout,
69 ... 100,
70 ... writers=[
71 ... MetaCVWriter(
72 ... umbrella,
73 ... values=["phi", "psi"],
74 ... emasses=["phi", "psi"],
75 ... parameters=["phi0", "psi0"],
76 ... derivatives=["phi0", "psi0"],
77 ... ),
78 ... ],
79 ... step=True,
80 ... )
81 >>> integrator = openmm.LangevinIntegrator(
82 ... 300 * unit.kelvin,
83 ... 1 / unit.picosecond,
84 ... 2 * unit.femtosecond,
85 ... )
86 >>> integrator.setRandomNumberSeed(1234)
87 >>> umbrella.addToSystem(model.system)
88 >>> simulation = app.Simulation(model.topology, model.system, integrator)
89 >>> simulation.context.setPositions(model.positions)
90 >>> simulation.context.setVelocitiesToTemperature(300 * unit.kelvin, 5678)
91 >>> simulation.reporters.append(reporter)
92 >>> simulation.step(1000) # doctest: +SKIP
93 #"Step","phi (rad)",...,"d[umbrella]/d[psi0] (kJ/(mol rad))"
94 100,2.6766...,-20.519...
95 200,2.3694...,0.7804...
96 300,3.1342...,49.272...
97 400,2.8179...,-10.012...
98 500,2.7565...,-8.9529...
99 600,2.6956...,23.049...
100 700,2.9331...,25.068...
101 800,2.6988...,-0.4481...
102 900,2.7477...,9.5528...
103 1000,2.7744...,21.986...
104 """
106 def __init__(
107 self,
108 metaCV: MetaCollectiveVariable,
109 values: t.Sequence[str] = (),
110 emasses: t.Sequence[str] = (),
111 parameters: t.Sequence[str] = (),
112 derivatives: t.Sequence[str] = (),
113 allowReinitialization: bool = False,
114 ) -> None:
115 inner_cvs = {cv.getName(): cv for cv in metaCV.getInnerVariables()}
116 all_parameters = metaCV.getParameterDefaultValues()
118 def select(names, available, kind):
119 unknown = [name for name in names if name not in available]
120 if unknown:
121 raise ValueError(f"Unknown {kind}: {', '.join(unknown)}")
122 return {name: available[name] for name in names}
124 self._meta_cv = metaCV
125 self._values = list(select(values, inner_cvs, "inner variables").values())
126 self._masses = list(select(emasses, inner_cvs, "inner variables").values())
127 self._parameters = select(parameters, all_parameters, "parameters")
128 self._derivatives = select(derivatives, all_parameters, "parameters")
129 self._allow_reinitialization = allowReinitialization
131 def getHeaders(self) -> list[str]:
132 headers = []
134 def add_header(name: str, unit: mmunit.Unit) -> None:
135 headers.append(f"{name} ({unit.get_symbol()})")
137 for cv in self._values:
138 add_header(cv.getName(), cv.getUnit())
139 for cv in self._masses:
140 add_header(f"emass[{cv.getName()}]", cv.getMassUnit())
141 for name, quantity in self._parameters.items():
142 add_header(name, quantity.unit)
143 for name, quantity in self._derivatives.items():
144 add_header(
145 f"d[{self._meta_cv.getName()}]/d[{name}]",
146 self._meta_cv.getUnit() / quantity.unit,
147 )
148 return headers
150 def getValues(self, simulation: mm.app.Simulation) -> list[float]:
151 context = simulation.context
152 values = []
153 if self._values:
154 inner_values = self._meta_cv.getInnerValues(context)
155 for cv in self._values:
156 values.append(inner_values[cv.getName()] / cv.getUnit())
157 if self._masses:
158 inner_masses = self._meta_cv.getInnerEffectiveMasses(context)
159 for cv in self._masses:
160 values.append(inner_masses[cv.getName()] / cv.getMassUnit())
161 if self._parameters:
162 parameters = self._meta_cv.getParameterValues(context)
163 for name, quantity in self._parameters.items():
164 values.append(parameters[name] / quantity.unit)
165 if self._derivatives:
166 derivatives = self._meta_cv.getParameterDerivatives(
167 context, self._allow_reinitialization
168 )
169 for name, quantity in self._derivatives.items():
170 values.append(
171 derivatives[name] / (self._meta_cv.getUnit() / quantity.unit)
172 )
173 return values