RMSD

class cvpack.RMSD(referencePositions, group, numAtoms=None, name='rmsd')[source]

The root-mean-square deviation (RMSD) between the current and reference coordinates of a group of n atoms:

\[d_{\rm rms}({\bf r}) = \sqrt{ \frac{1}{n} \min_{ \bf q \in \mathbb{R}^4 \atop \|{\bf q}\| = 1 } \sum_{i=1}^n \left\| {\bf A}({\bf q}) \hat{\bf r}_i - \hat{\bf r}_i^{\rm ref} \right\|^2 }\]

where \(\hat{\bf r}_i\) is the position of the \(i\)-th atom in the group relative to the group’s center of geometry (centroid), the superscript \(\rm ref\) denotes the reference structure, \({\bf q}\) is a unit quaternion, and \({\bf A}({\bf q})\) is the rotation matrix corresponding to \({\bf q}\).

Warning

Periodic boundary conditions are not supported. This is not a problem if all atoms in the group belong to the same molecule. If they belong to distinct molecules, it is possible to circumvent the issue by calling the method getNullBondForce() and adding the resulting force to the system.

Parameters:

  • referencePositions (Quantity | ndarray | Sequence[Vec3] | Sequence[ndarray] | Dict[int, Quantity | ndarray | Vec3]) – The reference coordinates, which can be either a coordinate matrix or a mapping from atom indices to coordinate vectors. It must contain all atoms in group, and does not need to contain all atoms in the system. See numAtoms below.

  • groups – A sequence of atom indices.

  • numAtoms (int | None) – The total number of atoms in the system, including those that are not in group. This argument is necessary only if referencePositions does not contain all atoms in the system.

  • name (str) – The name of the collective variable.

Example

>>> import cvpack
>>> import openmm
>>> from openmm import app, unit
>>> from openmmtools import testsystems
>>> model = testsystems.AlanineDipeptideImplicit()
>>> num_atoms = model.topology.getNumAtoms()
>>> group = list(range(num_atoms))
>>> rmsd = cvpack.RMSD(model.positions, group, num_atoms)
>>> rmsd.addToSystem(model.system)
>>> platform = openmm.Platform.getPlatformByName('Reference')
>>> integrator = openmm.VerletIntegrator(2*unit.femtoseconds)
>>> context = openmm.Context(model.system, integrator, platform)
>>> context.setPositions(model.positions)
>>> integrator.step(1000)
>>> rmsd.getValue(context)
0.123138... nm

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

getNullBondForce()[source]

Get a null bond force that creates a connected graph with all the atoms in the group.

Returns:

The null bond force.

Return type:

openmm.HarmonicBondForce

Example

>>> import cvpack
>>> import openmm
>>> from openmm import app, unit
>>> from openmmtools import testsystems
>>> model = testsystems.WaterBox(
...     box_edge=10*unit.angstroms,
...     cutoff=5*unit.angstroms,
... )
>>> group = [
...     atom.index
...     for atom in model.topology.atoms()
...     if atom.residue.index < 3
... ]
>>> rmsd = cvpack.RMSD(
...     model.positions, group, model.topology.getNumAtoms()
... )
>>> rmsd.addToSystem(model.system)
>>> _ = model.system.addForce(rmsd.getNullBondForce())
>>> integrator = openmm.VerletIntegrator(2*unit.femtoseconds)
>>> platform = openmm.Platform.getPlatformByName('Reference')
>>> context = openmm.Context(model.system, integrator, platform)
>>> context.setPositions(model.positions)
>>> integrator.step(100)
>>> rmsd.getValue(context)
0.10436... nm
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