NumberOfContacts¶

class cvpack.NumberOfContacts(group1, group2, nonbondedForce, reference=1.0, stepFunction='1/(1+x^6)', thresholdDistance=0.3 nm, cutoffFactor=2.0, switchFactor=1.5, name='number_of_contacts')[source]¶

The number of contacts between two atom groups:

\[N({\bf r}) = \sum_{i \in {\bf g}_1} \sum_{j \in {\bf g}_2} S\left(\frac{\|{\bf r}_j - {\bf r}_i\|}{r_0}\right)\]

where \(r_0\) is the threshold distance for defining a contact and \(S(x)\) is a step function equal to \(1\) if a contact is made or equal to \(0\) otherwise. For trajectory analysis, it is fine to make \(S(x) = H(1-x)\), where H is the Heaviside step function. For molecular dynamics, however, \(S(x)\) should be a continuous approximation of \(H(1-x)\) for \(x \geq 0\). By default [5], the following function is used:

\[S(x) = \frac{1-x^6}{1-x^{12}} = \frac{1}{1+x^6}\]

In fact, a cutoff distance \(r_c = x_c r_0\) (typically, \(x_c = 2\)) is applied so that \(S(x) = 0\) for \(x \geq x_c\). To avoid discontinuities, there is also the option to smoothly switch off \(S(x)\) starting from \(r_s = x_s r_0\) (typically, \(x_s = 1.5\)) instead of doing it abruptly at \(r_c\).

Note

Atoms are allowed to be in both groups. In this case, self-contacts (\(i = j\)) are ignored and each pair of distinct atoms (\(i \neq j\)) is counted only once.

Note

Any non-exclusion exceptions involving atoms in \({\bf g}_1\) and \({\bf g}_2\) in the provided openmm.NonbondedForce are turned into exclusions in this collective variable.

Parameters:

group1 (Sequence[int]) – The indices of the atoms in the first group.
group2 (Sequence[int]) – The indices of the atoms in the second group.
nonbondedForce (NonbondedForce) – The openmm.NonbondedForce object from which the total number of atoms, the exclusions, and whether to use periodic boundary conditions are taken.
reference (Real | Context) – A dimensionless reference value to which the collective variable should be normalized. One can also provide an openmm.Context object from which to obtain the reference number of contacts.
stepFunction (str) – The function “step(1-x)” (for analysis only) or a continuous approximation thereof.
thresholdDistance (Quantity | Real) – The threshold distance (\(r_0\)) for considering two atoms as being in contact.
cutoffFactor (float) – The factor \(x_c\) that multiplies the threshold distance to define the cutoff distance.
switchFactor (float | None) – The factor \(x_s\) that multiplies the threshold distance to define the distance at which the step function starts switching off smoothly. If None, it switches off abruptly at the cutoff distance.
name (str) – The name of the collective variable.

Example

>>> import cvpack
>>> from openmm import unit
>>> from openmmtools import testsystems
>>> model = testsystems.HostGuestExplicit()
>>> group1, group2 = [], []
>>> for residue in model.topology.residues():
...     if residue.name != "HOH":
...         group = group1 if residue.name == "B2" else group2
...         group.extend(atom.index for atom in residue.atoms())
>>> forces = {f.getName(): f for f in model.system.getForces()}
>>> nc = cvpack.NumberOfContacts(
...     group1,
...     group2,
...     forces["NonbondedForce"],
...     stepFunction="step(1-x)",
... )
>>> nc.addToSystem(model.system)
>>> platform = openmm.Platform.getPlatformByName("Reference")
>>> integrator = openmm.VerletIntegrator(1.0 * mmunit.femtoseconds)
>>> context = openmm.Context(model.system, integrator, platform)
>>> context.setPositions(model.positions)
>>> nc.getValue(context)
30.0... dimensionless
>>> nc_normalized = cvpack.NumberOfContacts(
...     group1,
...     group2,
...     forces["NonbondedForce"],
...     stepFunction="step(1-x)",
...     reference=context,
... )
>>> nc_normalized.addToSystem(model.system)
>>> context.reinitialize(preserveState=True)
>>> nc_normalized.getValue(context)
0.99999... 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

NumberOfContacts¶

CVPack v0.12.2

Table of Contents