Coverage for cvpack/path_in_cv_space.py: 95%
44 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 00:10 +0000
1"""
2.. class:: PathInCVSpace
3 :platform: Linux, MacOS, Windows
4 :synopsis: A metric of progress or deviation with respect to a path in CV space
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10from __future__ import annotations
12import typing as t
13from copy import deepcopy
15from openmm import unit as mmunit
17from .base_path_cv import BasePathCV
18from .collective_variable import CollectiveVariable
19from .path import Metric
20from .units.units import MatrixQuantity, ScalarQuantity, value_in_md_units
21from .utils import convert_to_matrix
24def _milestone_values(milestones: MatrixQuantity) -> t.Any:
25 """Strip units from milestones given as a Quantity or as rows of Quantities."""
26 if mmunit.is_quantity(milestones):
27 return value_in_md_units(milestones)
28 return [
29 [value_in_md_units(value) for value in row]
30 if isinstance(row, t.Iterable)
31 else row
32 for row in milestones
33 ]
36class PathInCVSpace(BasePathCV):
37 r"""
38 A metric of the system's progress (:math:`s`) or deviation (:math:`z`) with
39 respect to a path defined by a sequence of :math:`n` milestones positioned in a
40 collective variable space :cite:`Branduardi_2007`:
42 .. math::
44 s({\bf r}) = \frac{
45 \dfrac{\sum_{i=1}^n i w_i({\bf r})}{\sum_{i=1}^n w_i({\bf r})} - 1
46 }{n-1}
47 \quad \text{or} \quad
48 z({\bf r}) = - 2 \sigma ^2 \ln \sum_{i=1}^n w_i({\bf r})
50 with :math:`w_i({\bf r})` being a Gaussian kernel centered at the :math:`i`-th
51 milestone, i.e.,
53 .. math::
55 w_i({\bf r}) = \exp\left(\
56 -\frac{\|{\bf c}({\bf r}) - \hat{\bf c}_i\|^2}{2 \sigma^2}
57 \right)
59 where :math:`{\bf c}({\bf r})` is a vector of collective variables,
60 :math:`\hat{\bf c}_i` is the location of the :math:`i`-th milestone, and
61 :math:`\sigma` sets the width of the kernels. The squared norm of a vector
62 :math:`{\bf x}` in the collective variable space is defined as
64 .. math::
66 \|{\bf x}\|^2 = {\bf x}^T {\bf D}^{-2} {\bf x}
68 where :math:`{\bf D}` is a diagonal matrix whose diagonal elements are the
69 characteristic scales of the corresponding collective variables, which makes
70 :math:`\|{\bf x}\|^2` dimensionless. Appropriate boundary conditions are used for
71 periodic CVs.
73 .. note::
75 The kernel width :math:`\sigma` is related to the parameter :math:`\lambda` of
76 Ref. :cite:`Branduardi_2007` by :math:`\sigma = \frac{1}{\sqrt{2\lambda}}`.
78 Parameters
79 ----------
80 metric
81 The path-related metric to compute. Use ``cvpack.path.progress`` for
82 computing :math:`s({\bf r})` or ``cvpack.path.deviation`` for computing
83 :math:`z({\bf r})`.
84 variables
85 The collective variables that define the space.
86 milestones
87 The milestones in the collective variable space. The number of rows must be
88 equal to the number of milestones and the number of columns must be equal to
89 the number of collective variables.
90 sigma
91 The width of the Gaussian kernels.
92 scales
93 The characteristic scales for the collective variables. If not provided, the
94 scales are assumed to be 1 (in standard MD units) for each collective variable.
95 name
96 The name of the collective variable. If not provided, it is set to
97 "path_progress_in_cv_space" or "path_deviation_in_cv_space" depending on the
98 metric.
100 Raises
101 ------
102 ValueError
103 If the number of rows in the milestones matrix is less than 2
104 ValueError
105 If the number of columns in the milestones matrix is different from the number
106 of collective variables
107 ValueError
108 If the metric is not `cvpack.path.progress` or `cvpack.path.deviation`
110 Examples
111 --------
112 >>> import cvpack
113 >>> import openmm
114 >>> from openmmtools import testsystems
115 >>> import numpy as np
116 >>> model = testsystems.AlanineDipeptideVacuum()
117 >>> phi_atoms = ["ACE-C", "ALA-N", "ALA-CA", "ALA-C"]
118 >>> psi_atoms = ["ALA-N", "ALA-CA", "ALA-C", "NME-N"]
119 >>> atoms = [f"{a.residue.name}-{a.name}" for a in model.topology.atoms()]
120 >>> milestones = np.array(
121 ... [[1.3, -0.2], [1.2, 3.1], [-2.7, 2.9], [-1.3, 2.7]]
122 ... )
123 >>> phi = cvpack.Torsion(*[atoms.index(atom) for atom in phi_atoms])
124 >>> psi = cvpack.Torsion(*[atoms.index(atom) for atom in psi_atoms])
125 >>> path_vars = []
126 >>> for metric in (cvpack.path.progress, cvpack.path.deviation):
127 ... var = cvpack.PathInCVSpace(metric, [phi, psi], milestones, np.pi / 6)
128 ... var.addToSystem(model.system)
129 ... path_vars.append(var)
130 >>> context = openmm.Context(model.system, openmm.VerletIntegrator(1.0))
131 >>> context.setPositions(model.positions)
132 >>> path_vars[0].getValue(context)
133 0.6... dimensionless
134 >>> path_vars[1].getValue(context)
135 0.2... dimensionless
136 """
138 def __init__(
139 self,
140 metric: Metric,
141 variables: t.Iterable[CollectiveVariable],
142 milestones: MatrixQuantity,
143 sigma: ScalarQuantity,
144 scales: t.Iterable[ScalarQuantity] | None = None,
145 name: str | None = None,
146 ) -> None:
147 name = self._generateName(metric, name, "cv")
148 variables = list(variables)
149 scales = None if scales is None else [value_in_md_units(s) for s in scales]
150 cv_scales = [1.0] * len(variables) if scales is None else scales
151 sigma = value_in_md_units(sigma)
152 milestones, n, numvars = convert_to_matrix(_milestone_values(milestones))
153 if numvars != len(variables):
154 raise ValueError("Wrong number of columns in the milestones matrix.")
155 if n < 2: # noqa: PLR2004
156 raise ValueError("At least two rows are required in the milestones matrix.")
157 squared_distances = []
158 periods = {}
159 for i, variable in enumerate(variables):
160 values = variable.getPeriodicBounds()
161 if values is not None:
162 periods[i] = value_in_md_units(values[1] - values[0])
163 lower_bounds = {
164 i: value_in_md_units(variables[i].getPeriodicBounds()[0]) for i in periods
165 }
166 for milestone in milestones:
167 values = list(milestone)
168 for j, period in periods.items():
169 # With the milestone in the CV's periodic range, |delta| <= period
170 values[j] = lower_bounds[j] + (values[j] - lower_bounds[j]) % period
171 deltas = [f"({value}-cv{j})" for j, value in enumerate(values)]
172 for j, period in periods.items():
173 deltas[j] = f"min(abs{deltas[j]},{period}-abs{deltas[j]})"
174 squared_distances.append(
175 "+".join(
176 f"{delta}^2" if scale == 1.0 else f"({delta}/{scale})^2"
177 for delta, scale in zip(deltas, cv_scales)
178 )
179 )
180 collective_variables = {
181 f"cv{i}": deepcopy(variable) for i, variable in enumerate(variables)
182 }
183 super().__init__(metric, sigma, squared_distances, collective_variables)
184 self._registerCV(
185 name,
186 mmunit.dimensionless,
187 metric=metric,
188 variables=variables,
189 milestones=milestones.tolist(),
190 sigma=sigma,
191 scales=scales,
192 )
195PathInCVSpace.registerTag("!cvpack.PathInCVSpace")