CentroidFunction

class cvpack.CentroidFunction(function, unit, groups, collections=None, periodicBounds=None, pbc=True, weighByMass=True, name='centroid_function', **parameters)[source]

A generic function of the centroids of \(m \times n\) atoms groups split into m collections of n groups:

\[f({\bf r}) = \sum_{i=1}^m F\Big( {\bf g}^i_1({\bf r}), {\bf g}^i_2({\bf r}), \dots, {\bf g}^i_n({\bf r}) \Big)\]

where \(F\) is a user-defined function and \({\bf g}^i_1({\bf r})\) is the centroid of the \(j\)-th group of atoms of the \(i\)-th collection of groups.

The centroid of a group of atoms is defined as:

\[{\bf g}_j({\bf r}) = \frac{1}{N_j} \sum_{k=1}^{N_j} {\bf r}_{k,j}\]

where \(N_j\) is the number of atoms in group \(j\) and \({\bf r}_{k,j}\) is the coordinate of the \(k\)-th atom of group \(j\). Optionally, the centroid can be weighted by the mass of each atom in the group. In this case, it is redefined as:

\[{\bf g}_j({\bf r}) = \frac{1}{M_j} \sum_{k=1}^{N_j} m_{k,j} {\bf r}_{k,j}\]

where \(M_j\) is the total mass of atom group \(j\) and \(m_{k,j}\) is the mass of the \(k\)-th atom in group \(j\).

The function \(F\) is defined as a string and can be any expression supported by openmm.CustomCompoundBondForce. If the expression contains named parameters, the value of each parameter can be passed in one of three ways:

  1. By a semicolon-separated definition in the function string, such as described in the openmm.CustomCompoundBondForce documentation. In this case, the parameter value will be the same for all collections of atom groups.

  2. By a 1D array or list of length \(m\) passed as a keyword argument to the AtomicFunction constructor. In this case, each value will be assigned to a different collection of atom groups.

  3. By a scalar passed as a keyword argument to the AtomicFunction constructor. In this case, the parameter will apply to all collections of atom groups and will become available for on-the-fly modification during a simulation via the setParameter method of an openmm.Context object. Warning: other collective variables or openmm.Force objects in the same system will share the same values of equal-named parameters.

Parameters:

  • function (str) – The function to be evaluated. It must be a valid openmm.CustomCentroidBondForce expression.

  • unit (Unit) – The unit of measurement of the collective variable. It must be compatible with the MD unit system (mass in daltons, distance in nanometers, time in picoseconds, temperature in kelvin, energy in kilojoules_per_mol, angle in radians). If the collective variables does not have a unit, use dimensionless.

  • groups (Iterable[Iterable[int]]) – The groups of atoms to be used in the function. Each group must be specified as a list of atom indices with arbitrary length.

  • collections (_SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None) – The indices of the groups in each collection, passed as a 2D array-like object of shape (m, n), where m is the number of collections and n is the number groups per collection. If a 1D object is passed, it is assumed that m is 1 and n is the length of the object.

  • periodicBounds (Quantity | ndarray | Vec3 | None) – The periodic bounds of the collective variable if it is periodic, or None if it is not.

  • pbc (bool) – Whether to use periodic boundary conditions.

  • weighByMass (bool) – Whether to define the centroid as the center of mass of the group instead of the geometric center.

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

Keyword Arguments:

**parameters – The named parameters of the function. Each parameter can be a 1D array-like object or a scalar. In the former case, the array must have length \(m\) and each entry will be assigned to a different collection of atom groups. In the latter case, it will define a global openmm.Context parameter.

Raises:

  • ValueError – If the collections are not specified as a 1D or 2D array-like object.

  • ValueError – If group indices are out of bounds.

  • ValueError – If the unit of the collective variable is not compatible with the MD unit system.

Example

>>> import cvpack
>>> import numpy as np
>>> import openmm
>>> from openmm import unit
>>> from openmmtools import testsystems
>>> model = testsystems.LysozymeImplicit()
>>> residues = list(model.topology.residues())
>>> atoms = [[a.index for a in r.atoms()] for r in residues]

Compute the residue coordination between two helices:

>>> res_coord = cvpack.ResidueCoordination(
...     residues[115:124],
...     residues[126:135],
...     stepFunction="step(1-x)",
... )
>>> res_coord.addToSystem(model.system)
>>> integrator = openmm.VerletIntegrator(0)
>>> platform = openmm.Platform.getPlatformByName('Reference')
>>> context = openmm.Context(model.system, integrator, platform)
>>> context.setPositions(model.positions)
>>> res_coord.getValue(context)
33.0 dimensionless

Recompute the residue coordination using the centroid function:

>>> groups = [atoms[115:124], atoms[126:135]]
>>> colvar = cvpack.CentroidFunction(
...    "step(1 - distance(g1, g2))",
...    unit.dimensionless,
...    atoms[115:124] + atoms[126:135],
...    [[i, j] for i in range(9) for j in range(9, 18)],
... )
>>> colvar.addToSystem(model.system)
>>> context.reinitialize(preserveState=True)
>>> colvar.getValue(context)
33.0 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