PathInCVSpace¶

class cvpack.PathInCVSpace(metric, variables, milestones, sigma, scales=None, name=None)[source]¶

A metric of the system’s progress (\(s\)) or deviation (\(z\)) with respect to a path defined by a sequence of \(n\) milestones positioned in a collective variable space [6]:

\[s({\bf r}) = \frac{ \dfrac{\sum_{i=1}^n i w_i({\bf r})}{\sum_{i=1}^n w_i({\bf r})} - 1 }{n-1} \quad \text{or} \quad z({\bf r}) = - 2 \sigma ^2 \ln \sum_{i=1}^n w_i({\bf r})\]

with \(w_i({\bf r})\) being a Gaussian kernel centered at the \(i\)-th milestone, i.e.,

\[w_i({\bf r}) = \exp\left(\ -\frac{\|{\bf c}({\bf r}) - \hat{\bf c}_i\|^2}{2 \sigma^2} \right)\]

where \({\bf c}({\bf r})\) is a vector of collective variables, \(\hat{\bf c}_i\) is the location of the \(i\)-th milestone, and \(\sigma\) sets the width of the kernels. The squared norm of a vector \({\bf x}\) in the collective variable space is defined as

\[\|{\bf x}\|^2 = {\bf x}^T {\bf D}^{-2} {\bf x}\]

where \({\bf D}\) is a diagonal matrix whose each diagonal element is the characteristic scale of the corresponding collective variable, which makes \(\|{\bf x}\|^2\) dimensionless. Appropriate boundary conditions are used for periodic CVs.

Note

The kernel width \(\sigma\) is related to the parameter \(\lambda\) of Ref. [6] by \(\sigma = \frac{1}{\sqrt{2\lambda}}\).

Parameters:

metric (Metric) – The path-related metric to compute. Use cvpack.path.progress for computing \(s({\bf r})\) or cvpack.path.deviation for computing \(z({\bf r})\).
variables (Iterable[CollectiveVariable]) – The collective variables that define the space.
milestones (Quantity | ndarray | Sequence[Vec3] | Sequence[ndarray]) – The milestones in the collective variable space. The number of rows must be equal to the number of milestones and the number of columns must be equal to the number of collective variables.
sigma (float) – The width of the Gaussian kernels.
scales (Iterable[Quantity | Real] | None) – The characteristic scales for the collective variables. If not provided, the scales are assumed to be 1 (in standard MD units) for each collective variable.
name (str | None) – The name of the collective variable. If not provided, it is set to “path_progress_in_cv_space” or “path_deviation_in_cv_space” depending on the metric.

Raises:

ValueError – If the number of rows in the milestones matrix is less than 2
ValueError – If the number of columns in the milestones matrix is different from the number of collective variables
ValueError – If the metric is not cvpack.path.progress or cvpack.path.deviation

Examples

>>> import cvpack
>>> import openmm
>>> from openmmtools import testsystems
>>> import numpy as np
>>> model = testsystems.AlanineDipeptideVacuum()
>>> phi_atoms = ["ACE-C", "ALA-N", "ALA-CA", "ALA-C"]
>>> psi_atoms = ["ALA-N", "ALA-CA", "ALA-C", "NME-N"]
>>> atoms = [f"{a.residue.name}-{a.name}" for a in model.topology.atoms()]
>>> milestones = np.array(
...     [[1.3, -0.2], [1.2, 3.1], [-2.7, 2.9], [-1.3, 2.7]]
... )
>>> phi = cvpack.Torsion(*[atoms.index(atom) for atom in phi_atoms])
>>> psi = cvpack.Torsion(*[atoms.index(atom) for atom in psi_atoms])
>>> path_vars = []
>>> for metric in (cvpack.path.progress, cvpack.path.deviation):
...     var = cvpack.PathInCVSpace(metric, [phi, psi], milestones, np.pi / 6)
...     var.addToSystem(model.system)
...     path_vars.append(var)
>>> context = openmm.Context(model.system, openmm.VerletIntegrator(1.0))
>>> context.setPositions(model.positions)
>>> path_vars[0].getValue(context)
0.6... dimensionless
>>> path_vars[1].getValue(context)
0.2... dimensionless

Methods

addToSystem(system, setUnusedForceGroup=True)¶

Add this collective variable to an openmm.System.

Parameters:

system (System) – The system to which this collective variable should be added
setUnusedForceGroup (bool) – If True, the force group of this collective variable will be set to the first available force group in the system

getEffectiveMass(context, allowReinitialization=False)¶

Compute the effective mass of this collective variable at a given openmm.Context.

The effective mass of a collective variable \(q({\bf r})\) is defined as [1]:

\[m_\mathrm{eff}({\bf r}) = \left( \sum_{i=1}^N \frac{1}{m_i} \left\| \frac{dq}{d{\bf r}_i} \right\|^2 \right)^{-1}\]

If this collective variable share the force group with other forces, then evaluating its effective mass requires reinitializing the openmm.Context twice at each call. This is inefficient and should be avoided. To allow this behavior, the allowReinitialization parameter must be set to True.

Note

By adding this collective variable to the system using the addToSystem() method, the force group of this collective variable is set to an available force group in the system by default.

Parameters:

context (Context) – The context at which this collective variable’s effective mass should be evaluated.
allowReinitialization (bool) – If True, the context will be reinitialized if necessary.

Returns:

The effective mass of this collective variable at the given context.

Return type:

Quantity

Example

In this example, we compute the effective masses of the backbone dihedral angles and the radius of gyration of an alanine dipeptide molecule in water:

>>> import cvpack
>>> import openmm
>>> from openmmtools import testsystems
>>> model = testsystems.AlanineDipeptideExplicit()
>>> top = model.mdtraj_topology
>>> backbone_atoms = top.select("name N C CA and resid 1 2")
>>> phi = cvpack.Torsion(*backbone_atoms[0:4])
>>> psi = cvpack.Torsion(*backbone_atoms[1:5])
>>> radius_of_gyration = cvpack.RadiusOfGyration(
...     top.select("not water")
... )
>>> for cv in [phi, psi, radius_of_gyration]:
...     cv.addToSystem(model.system)
>>> context = openmm.Context(
...     model.system, openmm.VerletIntegrator(0)
... )
>>> context.setPositions(model.positions)
>>> phi.getEffectiveMass(context)
0.05119... nm**2 Da/(rad**2)
>>> psi.getEffectiveMass(context)
0.05186... nm**2 Da/(rad**2)
>>> radius_of_gyration.getEffectiveMass(context)
30.946... Da

getMassUnit()¶

Get the unit of measurement of the effective mass of this collective variable.

Return type:: Unit

getName()¶

Get the name of this collective variable.

Return type:: str

getPeriodicBounds()¶

Get the periodic bounds of this collective variable.

Returns:: The periodic bounds of this collective variable or None if it is not periodic
Return type:: Union[Quantity, None]

getUnit()¶

Get the unit of measurement of this collective variable.

Return type:: Unit

getValue(context, allowReinitialization=False)¶

Evaluate this collective variable at a given openmm.Context.

If this collective variable share the force group with other forces, then its evaluation requires reinitializing the openmm.Context twice at each call. This is inefficient and should be avoided. To allow this behavior, the allowReinitialization parameter must be set to True.

Note

By adding this collective variable to the system using the addToSystem() method, the force group of this collective variable is set to an available force group in the system by default.

Parameters:

context (Context) – The context at which this collective variable should be evaluated.
allowReinitialization (bool) – If True, the context will be reinitialized if necessary.

Returns:

The value of this collective variable at the given context.

Return type:

Quantity

Example

In this example, we compute the values of the backbone dihedral angles and the radius of gyration of an alanine dipeptide molecule in water:

>>> import cvpack
>>> import openmm
>>> from openmmtools import testsystems
>>> model = testsystems.AlanineDipeptideExplicit()
>>> top = model.mdtraj_topology
>>> backbone_atoms = top.select("name N C CA and resid 1 2")
>>> phi = cvpack.Torsion(*backbone_atoms[0:4])
>>> psi = cvpack.Torsion(*backbone_atoms[1:5])
>>> radius_of_gyration = cvpack.RadiusOfGyration(
...     top.select("not water")
... )
>>> for cv in [phi, psi, radius_of_gyration]:
...     cv.addToSystem(model.system)
>>> context = openmm.Context(
...     model.system, openmm.VerletIntegrator(0)
... )
>>> context.setPositions(model.positions)
>>> phi.getValue(context)
3.1415... rad
>>> psi.getValue(context)
3.1415... rad
>>> radius_of_gyration.getValue(context)
0.29514... nm

PathInCVSpace¶

CVPack v0.12.2

Table of Contents