AtomicFunction¶

class cvpack.AtomicFunction(function, unit, groups, periodicBounds=None, pbc=True, name='atomic_function', **parameters)[source]¶

A generic function of the coordinates of \(m\) groups of \(n\) atoms:

\[f({\bf r}) = \sum_{i=1}^m F\left( {\bf r}_{i,1}, {\bf r}_{i,2}, \dots, {\bf r}_{i,n} \right)\]

where \(F\) is a user-defined function and \({\bf r}_{i,j}\) is the coordinate of the \(j\)-th atom of the \(i\)-th group.

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:

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 groups of atoms.
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 group of atoms.
By a scalar passed as a keyword argument to the AtomicFunction constructor. In this case, the parameter will apply to all 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.CustomCompoundBondForce expression.
groups (ArrayLike) – The indices of the atoms in each group, passed as a 2D array-like object of shape (m, n), where m is the number of groups and n is the number of atoms per group. If a 1D object is passed, it is assumed that m is 1 and n is the length of the object.
unit (mmunit.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 unit.dimensionless.
periodicBounds (t.Optional[VectorQuantity]) – The lower and upper bounds of the collective variable if it is periodic, or None if it is not.
pbc (bool) – Whether to use periodic boundary conditions when computing atomic distances.
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 group of atoms. In the latter case, it will define a global openmm.Context parameter.

Raises:

ValueError – If the groups are not specified as a 1D or 2D array-like object.
ValueError – If the unit of the collective variable is not compatible with the MD unit system.

Example

>>> import cvpack
>>> import openmm
>>> import numpy as np
>>> from openmm import unit
>>> from openmmtools import testsystems
>>> model = testsystems.AlanineDipeptideVacuum()
>>> angle1 = cvpack.Angle(0, 5, 10)
>>> angle2 = cvpack.Angle(10, 15, 20)
>>> colvar = cvpack.AtomicFunction(
...     "(kappa/2)*(angle(p1, p2, p3) - theta0)^2",
...     unit.kilojoules_per_mole,
...     [[0, 5, 10], [10, 15, 20]],
...     kappa = 1000 * unit.kilojoules_per_mole/unit.radian**2,
...     theta0 = [np.pi/2, np.pi/3] * unit.radian,
... )
>>> for cv in [angle1, angle2, colvar]:
...     cv.addToSystem(model.system)
>>> integrator = openmm.VerletIntegrator(0)
>>> platform = openmm.Platform.getPlatformByName('Reference')
>>> context = openmm.Context(model.system, integrator, platform)
>>> context.setPositions(model.positions)
>>> theta1 = angle1.getValue(context) / openmm.unit.radian
>>> theta2 = angle2.getValue(context) / openmm.unit.radian
>>> 500*((theta1 - np.pi/2)**2 + (theta2 - np.pi/3)**2)
429.479...
>>> colvar.getValue(context)
429.479... kJ/mol
>>> context.setParameter(
...     "kappa",
...     2000 * unit.kilojoules_per_mole/unit.radian**2
... )
>>> colvar.getValue(context)
858.958... kJ/mol

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

classmethod fromOpenMMForce(force, unit, periodicBounds=None, pbc=False)[source]¶

Create an AtomicFunction from an openmm.Force.

Parameters:

force (Force) – The force to be converted
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 unit.dimensionless.
periodicBounds (Quantity | ndarray | Vec3 | None) – The lower and upper bounds of the collective variable if it is periodic, or None if it is not. This parameter is considered only if force is a custom openmm.Force.
pbc (bool) – Whether to use periodic boundary conditions when computing atomic distances.

Raises:

TypeError – If the force is not convertible to an AtomicFunction

Return type:

AtomicFunction

Examples

>>> import cvpack
>>> import numpy as np
>>> import openmm
>>> from openmm import unit
>>> from openmm import app
>>> from openmmtools import testsystems
>>> model = testsystems.LysozymeImplicit()
>>> residues = [r for r in model.topology.residues() if 59 <= r.index <= 79]
>>> helix_content = cvpack.HelixTorsionContent(residues)
>>> model.system.addForce(helix_content)
6
>>> num_atoms = model.system.getNumParticles()
>>> mean_x = openmm.CustomExternalForce("x/num_atoms")
>>> mean_x.addGlobalParameter("num_atoms", num_atoms)
0
>>> for i in range(num_atoms):
...     _ = mean_x.addParticle(i, [])
>>> model.system.addForce(mean_x)
7
>>> forces = {force.getName(): force for force in model.system.getForces()}
>>> copies = {
...     name: cvpack.AtomicFunction.fromOpenMMForce(
...         force, unit.kilojoules_per_mole
...     )
...     for name, force in forces.items()
...     if name in [
...         "HarmonicBondForce",
...         "HarmonicAngleForce",
...         "PeriodicTorsionForce",
...         "CustomExternalForce"
...     ]
... }
>>> copies["HelixTorsionContent"] = cvpack.AtomicFunction.fromOpenMMForce(
...     helix_content, unit.dimensionless
... )
>>> indices = {}
>>> for index, (name, force) in enumerate(copies.items(), start=1):
...     _ = model.system.addForce(force)
...     force.setForceGroup(index)
...     indices[name] = index
>>> platform = openmm.Platform.getPlatformByName('Reference')
>>> integrator = openmm.VerletIntegrator(0)
>>> context = openmm.Context(model.system, integrator, platform)
>>> context.setPositions(model.positions)
>>> for name in copies:
...    state = context.getState(getEnergy=True, groups={indices[name]})
...    value = state.getPotentialEnergy() / unit.kilojoules_per_mole
...    copy_value = copies[name].getValue(context)
...    print(f"{name}: original={value:.6f}, copy={copy_value}")
HarmonicBondForce: original=2094.312..., copy=2094.312... kJ/mol
HarmonicAngleForce: original=3239.795..., copy=3239.795... kJ/mol
PeriodicTorsionForce: original=4226.05..., copy=4226.05... kJ/mol
CustomExternalForce: original=5.02155..., copy=5.02155... kJ/mol
HelixTorsionContent: original=17.4528..., copy=17.4528... dimensionless

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

AtomicFunction¶

CVPack v0.12.2

Table of Contents