Coverage for cvpack/path_in_cv_space.py: 94%
33 statements
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-09 16:14 +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"""
10import typing as t
11from copy import deepcopy
13from openmm import unit as mmunit
15from .base_path_cv import BasePathCV
16from .collective_variable import CollectiveVariable
17from .path import Metric
18from .units.units import MatrixQuantity, ScalarQuantity, value_in_md_units
19from .utils import convert_to_matrix
22class PathInCVSpace(BasePathCV):
23 r"""
24 A metric of the system's progress (:math:`s`) or deviation (:math:`z`) with
25 respect to a path defined by a sequence of :math:`n` milestones positioned in a
26 collective variable space :cite:`Branduardi_2007`:
28 .. math::
30 s({\bf r}) = \frac{
31 \dfrac{\sum_{i=1}^n i w_i({\bf r})}{\sum_{i=1}^n w_i({\bf r})} - 1
32 }{n-1}
33 \quad \text{or} \quad
34 z({\bf r}) = - 2 \sigma ^2 \ln \sum_{i=1}^n w_i({\bf r})
36 with :math:`w_i({\bf r})` being a Gaussian kernel centered at the :math:`i`-th
37 milestone, i.e.,
39 .. math::
41 w_i({\bf r}) = \exp\left(\
42 -\frac{\|{\bf c}({\bf r}) - \hat{\bf c}_i\|^2}{2 \sigma^2}
43 \right)
45 where :math:`{\bf c}({\bf r})` is a vector of collective variables,
46 :math:`\hat{\bf c}_i` is the location of the :math:`i`-th milestone, and
47 :math:`\sigma` sets the width of the kernels. The squared norm of a vector
48 :math:`{\bf x}` in the collective variable space is defined as
50 .. math::
52 \|{\bf x}\|^2 = {\bf x}^T {\bf D}^{-2} {\bf x}
54 where :math:`{\bf D}` is a diagonal matrix whose each diagonal element is the
55 characteristic scale of the corresponding collective variable, which makes
56 :math:`\|{\bf x}\|^2` dimensionless. Appropriate boundary conditions are used for
57 periodic CVs.
59 .. note::
61 The kernel width :math:`\sigma` is related to the parameter :math:`\lambda` of
62 Ref. :cite:`Branduardi_2007` by :math:`\sigma = \frac{1}{\sqrt{2\lambda}}`.
64 Parameters
65 ----------
66 metric
67 The path-related metric to compute. Use ``cvpack.path.progress`` for
68 computing :math:`s({\bf r})` or ``cvpack.path.deviation`` for computing
69 :math:`z({\bf r})`.
70 variables
71 The collective variables that define the space.
72 milestones
73 The milestones in the collective variable space. The number of rows must be
74 equal to the number of milestones and the number of columns must be equal to
75 the number of collective variables.
76 sigma
77 The width of the Gaussian kernels.
78 scales
79 The characteristic scales for the collective variables. If not provided, the
80 scales are assumed to be 1 (in standard MD units) for each collective variable.
81 name
82 The name of the collective variable. If not provided, it is set to
83 "path_progress_in_cv_space" or "path_deviation_in_cv_space" depending on the
84 metric.
86 Raises
87 ------
88 ValueError
89 If the number of rows in the milestones matrix is less than 2
90 ValueError
91 If the number of columns in the milestones matrix is different from the number
92 of collective variables
93 ValueError
94 If the metric is not `cvpack.path.progress` or `cvpack.path.deviation`
96 Examples
97 --------
98 >>> import cvpack
99 >>> import openmm
100 >>> from openmmtools import testsystems
101 >>> import numpy as np
102 >>> model = testsystems.AlanineDipeptideVacuum()
103 >>> phi_atoms = ["ACE-C", "ALA-N", "ALA-CA", "ALA-C"]
104 >>> psi_atoms = ["ALA-N", "ALA-CA", "ALA-C", "NME-N"]
105 >>> atoms = [f"{a.residue.name}-{a.name}" for a in model.topology.atoms()]
106 >>> milestones = np.array(
107 ... [[1.3, -0.2], [1.2, 3.1], [-2.7, 2.9], [-1.3, 2.7]]
108 ... )
109 >>> phi = cvpack.Torsion(*[atoms.index(atom) for atom in phi_atoms])
110 >>> psi = cvpack.Torsion(*[atoms.index(atom) for atom in psi_atoms])
111 >>> path_vars = []
112 >>> for metric in (cvpack.path.progress, cvpack.path.deviation):
113 ... var = cvpack.PathInCVSpace(metric, [phi, psi], milestones, np.pi / 6)
114 ... var.addToSystem(model.system)
115 ... path_vars.append(var)
116 >>> context = openmm.Context(model.system, openmm.VerletIntegrator(1.0))
117 >>> context.setPositions(model.positions)
118 >>> path_vars[0].getValue(context)
119 0.6... dimensionless
120 >>> path_vars[1].getValue(context)
121 0.2... dimensionless
122 """
124 def __init__( # pylint: disable=too-many-branches
125 self,
126 metric: Metric,
127 variables: t.Iterable[CollectiveVariable],
128 milestones: MatrixQuantity,
129 sigma: float,
130 scales: t.Optional[t.Iterable[ScalarQuantity]] = None,
131 name: t.Optional[str] = None,
132 ) -> None:
133 name = self._generateName(metric, name, "cv")
134 variables = list(variables)
135 cv_scales = [1.0] * len(variables) if scales is None else list(scales)
136 milestones, n, numvars = convert_to_matrix(milestones)
137 if numvars != len(variables):
138 raise ValueError("Wrong number of columns in the milestones matrix.")
139 if n < 2:
140 raise ValueError("At least two rows are required in the milestones matrix.")
141 squared_distances = []
142 periods = {}
143 for i, variable in enumerate(variables):
144 values = variable.getPeriodicBounds()
145 if values is not None:
146 periods[i] = value_in_md_units(values[1] - values[0])
147 for i, values in enumerate(milestones):
148 deltas = [f"({value}-cv{j})" for j, value in enumerate(values)]
149 for j, period in periods.items():
150 deltas[j] = f"min(abs{deltas[j]},{period}-abs{deltas[j]})"
151 squared_distances.append(
152 "+".join(
153 f"{delta}^2" if scale == 1.0 else f"({delta}/{scale})^2"
154 for delta, scale in zip(deltas, cv_scales)
155 )
156 )
157 collective_variables = {
158 f"cv{i}": deepcopy(variable) for i, variable in enumerate(variables)
159 }
160 super().__init__(metric, sigma, squared_distances, collective_variables)
161 self._registerCV(
162 name,
163 mmunit.dimensionless,
164 metric=metric,
165 variables=variables,
166 milestones=milestones.tolist(),
167 sigma=sigma,
168 scales=scales,
169 )
172PathInCVSpace.registerTag("!cvpack.PathInCVSpace")