Added contextParameter argument to Reporter

cvpack/reporter.py

@@ -44,6 +44,9 @@ class Reporter(mmapp.StateDataReporter):
     metaCVs
         The meta-collective variables whose inner variables will be reported. These must
         be the same objects added to the simulation's :OpenMM:`System`.
+    contextParameters
+        Global :OpenMM:`Context` parameters to be reported along with their units of
+        measurement. Use `unit.dimensionless` for parameters without units.
     step
         Whether to report the current step index.
     time
@@ -63,6 +66,8 @@ class Reporter(mmapp.StateDataReporter):
         If `variables` and `metaCVs` are both empty.
     ValueError
         If `values` and `masses` are both `False`.
+    ValueError
+        If any value `contextParameters` is not a unit of measurement.
 
     Examples
     --------
@@ -75,18 +80,20 @@ class Reporter(mmapp.StateDataReporter):
     >>> phi = cvpack.Torsion(6, 8, 14, 16, name="phi")
     >>> psi = cvpack.Torsion(8, 14, 16, 18, name="psi")
     >>> umbrella = cvpack.MetaCollectiveVariable(
-    ...     "50*(min(dphi,2*pi-dphi)^2+min(dpsi,2*pi-dpsi)^2)"
+    ...     "0.5*kappa*(min(dphi,2*pi-dphi)^2+min(dpsi,2*pi-dpsi)^2)"
     ...     "; dphi=abs(phi-5*pi/6); dpsi=abs(psi+5*pi/6)",
     ...     [phi, psi],
     ...     unit.kilojoules_per_mole,
     ...     name="umbrella",
     ...     pi=3.141592653589793,
+    ...     kappa=100 * unit.kilojoules_per_mole/unit.radian**2,
     ... )
     >>> reporter = cvpack.Reporter(
     ...     stdout,
     ...     100,
     ...     variables=[umbrella],
     ...     metaCVs=[umbrella],
+    ...     contextParameters={"kappa": unit.kilojoules_per_mole/unit.radian**2},
     ...     step=True,
     ...     values=True,
     ...     masses=True
@@ -103,25 +110,26 @@ class Reporter(mmapp.StateDataReporter):
     >>> simulation.context.setVelocitiesToTemperature(300 * unit.kelvin, 5678)
     >>> simulation.reporters.append(reporter)
     >>> simulation.step(1000)  # doctest: +SKIP
-    #"Step","umbrella (kJ/mol)",...,"psi mass (nm**2 Da/(rad**2))"
-    100,11.2618...,1.8774...e-05,2.3684...,0.04388...,-3.0217...,0.05552...
-    200,7.46382...,4.8803...e-05,2.8851...,0.04968...,-2.8971...,0.05230...
-    300,2.55804...,8.9086...e-05,2.4311...,0.04637...,-2.4905...,0.03721...
-    400,6.19945...,4.6824...e-05,2.9678...,0.05449...,-2.6577...,0.04840...
-    500,8.82728...,3.0253...e-05,2.5838...,0.04584...,-3.0367...,0.05485...
-    600,3.76160...,7.9099...e-05,2.7248...,0.04772...,-2.8706...,0.05212...
-    700,3.38895...,6.3273...e-05,2.5583...,0.04999...,-2.8714...,0.04674...
-    800,1.07166...,0.00028746...,2.7104...,0.05321...,-2.7314...,0.05418...
-    900,8.58602...,2.4899...e-05,2.4391...,0.04096...,-2.9917...,0.04675...
-    1000,5.8404...,5.1011...e-05,2.7584...,0.04951...,-2.9295...,0.05030...
+    #"Step","umbrella (kJ/mol)",...(nm**2 Da/(rad**2))","kappa (kJ/(mol rad**2))"
+    100,11.2618...,1.8774...e-05,2.3684...,0.04388...,-3.0217...,0.05552...,100.0
+    200,7.46382...,4.8803...e-05,2.8851...,0.04968...,-2.8971...,0.05230...,100.0
+    300,2.55804...,8.9086...e-05,2.4311...,0.04637...,-2.4905...,0.03721...,100.0
+    400,6.19945...,4.6824...e-05,2.9678...,0.05449...,-2.6577...,0.04840...,100.0
+    500,8.82728...,3.0253...e-05,2.5838...,0.04584...,-3.0367...,0.05485...,100.0
+    600,3.76160...,7.9099...e-05,2.7248...,0.04772...,-2.8706...,0.05212...,100.0
+    700,3.38895...,6.3273...e-05,2.5583...,0.04999...,-2.8714...,0.04674...,100.0
+    800,1.07166...,0.00028746...,2.7104...,0.05321...,-2.7314...,0.05418...,100.0
+    900,8.58602...,2.4899...e-05,2.4391...,0.04096...,-2.9917...,0.04675...,100.0
+    1000,5.8404...,5.1011...e-05,2.7584...,0.04951...,-2.9295...,0.05030...,100.0
     """
 
-    def __init__(
+    def __init__(  # pylint: disable=dangerous-default-value
         self,
         file: t.Union[str, io.TextIOBase],
         reportInterval: int,
         variables: t.Sequence[CollectiveVariable] = (),
         metaCVs: t.Sequence[MetaCollectiveVariable] = (),
+        contextParameters: t.Mapping[str, mmunit.Unit] = {},
         step: bool = False,
         time: bool = False,
         values: bool = False,
@@ -133,6 +141,9 @@ def __init__(
             raise ValueError("Arguments 'variables' and 'metaCVs' cannot be both empty")
         if not (values or masses):
             raise ValueError("Arguments 'values' and 'masses' cannot be both False")
+        for parameter, unit in contextParameters.items():
+            if not mmunit.is_unit(unit):
+                raise ValueError(f"Invalid unit '{unit}' for parameter '{parameter}'")
         super().__init__(
             file,
             reportInterval,
@@ -145,6 +156,7 @@ def __init__(
         self._meta_cvs = metaCVs
         self._values = values
         self._masses = masses
+        self._parameters = contextParameters
 
     def _constructReportValues(  # pylint: disable=too-many-branches
         self, simulation: mmapp.Simulation, state: openmm.State
@@ -168,6 +180,8 @@ def _constructReportValues(  # pylint: disable=too-many-branches
                     values.append(cv_values[name] / cv_values[name].unit)
                 if name in cv_masses:
                     values.append(cv_masses[name] / cv_masses[name].unit)
+        for parameter in self._parameters:
+            values.append(context.getParameter(parameter))
         return values
 
     def _constructHeaders(self) -> t.List[str]:
@@ -188,4 +202,6 @@ def add_headers(cv: CollectiveVariable) -> None:
         for cv in self._meta_cvs:
             for inner_cv in cv.getInnerVariables():
                 add_headers(inner_cv)
+        for parameter, unit in self._parameters.items():
+            headers.append(f"{parameter} ({unit.get_symbol()})")
         return headers

cvpack/tests/test_cvpack.py

@@ -864,9 +864,10 @@ def test_reporter():
     model = testsystems.AlanineDipeptideVacuum()
     phi = cvpack.Torsion(6, 8, 14, 16, name="phi")
     umbrella = cvpack.MetaCollectiveVariable(
-        f"50*min(delta,2*pi-delta)^2" "; delta=abs(phi-5*pi/6)" f"; pi={np.pi}",
+        f"0.5*kappa*min(delta,2*pi-delta)^2" "; delta=abs(phi-5*pi/6)" f"; pi={np.pi}",
         [phi],
         unit.kilojoules_per_mole,
+        kappa=100,
         name="umbrella",
     )
     with tempfile.TemporaryDirectory() as dirpath:
@@ -884,6 +885,7 @@ def test_reporter():
                 1,
                 [umbrella],
                 [umbrella],
+                {"kappa": unit.kilojoules_per_mole / unit.radian**2},
                 step=True,
                 values=True,
                 masses=True,
@@ -897,6 +899,7 @@ def test_reporter():
                     '"umbrella (kJ/mol)"',
                     '"umbrella mass (nm**2 mol**2 Da/(kJ**2))"',
                     '"phi (rad)"',
-                    '"phi mass (nm**2 Da/(rad**2))"\n',
+                    '"phi mass (nm**2 Da/(rad**2))"',
+                    '"kappa (kJ/(mol rad**2))"\n',
                 ]
             )