AttractionStrength

class cvpack.AttractionStrength(group1, group2, nonbondedForce, contrastGroup=None, reference=1.0 kJ/mol, contrastScaling=1.0, name='attraction_strength')

The strength of the attraction between two atom groups:

\[\begin{split}S_{\rm attr}({\bf r}) = S_{1,2}({\bf r}) = &-\frac{1}{E_{\rm ref}} \Bigg[ \sum_{i \in {\bf g}_1} \sum_{\substack{j \in {\bf g}_2 \\ j \neq i}} \epsilon_{ij} u_{\rm disp}\left( \frac{\|{\bf r}_i - {\bf r}_j\|}{\sigma_{ij}} \right) \\ &+\sum_{i \in {\bf g}_1} \sum_{\substack{j \in {\bf g}_2 \\ q_jq_i < 0}} \frac{q_i q_j}{4 \pi \varepsilon_0 r_{\rm c}} u_{\rm elec}\left( \frac{\|{\bf r}_i - {\bf r}_j\|}{r_{\rm c}} \right) \Bigg]\end{split}\]

where \({\bf g}_1\) and \({\bf g}_2\) are the two atom groups, \(r_{\rm c}\) is the cutoff distance, \(\varepsilon_0\) is the permittivity of empty space, \(q_i\) is the charge of atom \(i\), and \(E_{\rm ref}\) is a reference value (in energy units per mole). The Lennard-Jones parameters are given by the Lorentz-Berthelot mixing rule, i.e. \(\epsilon_{ij} = \sqrt{\epsilon_i \epsilon_j}\), and \(\sigma_{ij} = (\sigma_i + \sigma_j)/2\).

Optionally, one can provide a third atom group \({\bf g}_3\) to contrast the attraction strength between \({\bf g}_1\) and \({\bf g}_2\) with that between \({\bf g}_1\) and \({\bf g}_3\). One can also provide a scaling factor \(\alpha\) to balance the contributions of the two interactions. In this case, the collective variable becomes

\[S_{\rm attr}({\bf r}) = S_{1,2}({\bf r}) - \alpha S_{1,3}({\bf r})\]

Note

Groups \({\bf g}_1\) and \({\bf g}_2\) can overlap or even be the same, in which case the collective variable will measure the strength of the self-attraction of \({\bf g}_1\). On the other hand, the contrast group \({\bf g}_3\) cannot overlap with neither \({\bf g}_1\) nor \({\bf g}_2\).

The function \(u_{\rm disp}(x)\) is a Lennard-Jones-like reduced potential with a highly softened repulsion part, defined as

\[u_{\rm disp}(x) = 4\left(\frac{1}{y^2} - \frac{1}{y}\right), \quad \text{where} \quad y = |x^6 - 2| + 2\]

The function \(u_{\rm elec}(x)\) provides a Coulomb-like decay with reaction-field screening, defined as

\[u_{\rm elec}(x) = \frac{1}{x} + \frac{x^2 - 3}{2}\]

The screening considers a perfect conductor as the surrounding medium [2].

Note

Only attractive electrostatic interactions are considered (\(q_i q_i < 0\)), which gives \(S_{\rm attr}({\bf r})\) a lower bound of zero. The upper bound will depends on the system details, the chosen groups of atoms, and the adopted reference value.

The Lennard-Jones parameters, atomic charges, cutoff distance, boundary conditions, as well as whether to use a switching function and its corresponding switching distance, are taken from openmm.NonbondedForce object.

Note

Any non-exclusion exceptions involving an atom in \({\bf g}_1\) and an atom in either \({\bf g}_2\) or \({\bf g}_3\) are turned into exclusions in this collective variable.

Parameters:

• group1 (Iterable[int]) – The first atom group.

• group2 (Iterable[int]) – The second atom group.

• nonbondedForce (NonbondedForce) – The openmm.NonbondedForce object from which to collect the necessary parameters.

• contrastGroup (Iterable[int] | None) – An optional third atom group to contrast the attraction strength between \({\bf g}_1\) and \({\bf g}_2\) with that between \({\bf g}_1\) and \({\bf g}_3\).

• reference (Quantity | Real | Context) – A reference value (in energy units per mole) to which the collective variable should be normalized. One can also provide an openmm.Context object from which to extract a reference attraction strength. The extracted value will be \(S_{1,2}({\bf r})\) for \(E_{\rm ref} = 1\), regardless of whether contrastGroup is provided or not.

• contrastScaling (float) – A scaling factor \(\alpha\) to balance the contributions of the two interactions. The default is 1.0.

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

Examples

>>> import cvpack
>>> from openmm import unit
>>> from openmmtools import testsystems
>>> model = testsystems.HostGuestExplicit()
>>> guest = [a.index for a in model.topology.atoms() if a.residue.name == "B2"]
>>> host = [a.index for a in model.topology.atoms() if a.residue.name == "CUC"]
>>> forces = {f.getName(): f for f in model.system.getForces()}
>>> cv1 = cvpack.AttractionStrength(guest, host, forces["NonbondedForce"])
>>> cv1.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)
>>> cv1.getValue(context)
4912.5... dimensionless

>>> water = [a.index for a in model.topology.atoms() if a.residue.name == "HOH"]
>>> cv2 = cvpack.AttractionStrength(guest, water, forces["NonbondedForce"])
>>> cv2.addToSystem(model.system)
>>> context.reinitialize(preserveState=True)
>>> cv2.getValue(context)
2063.3... dimensionless

>>> cv3 = cvpack.AttractionStrength(guest, host, forces["NonbondedForce"], water)
>>> cv3.addToSystem(model.system)
>>> context.reinitialize(preserveState=True)
>>> cv3.getValue(context)
2849.17... dimensionless
>>> print(cv1.getValue(context) - cv2.getValue(context))
2849.17... dimensionless

>>> cv4 = cvpack.AttractionStrength(
...     guest, host, forces["NonbondedForce"], water, contrastScaling=0.5
... )
>>> cv4.addToSystem(model.system)
>>> context.reinitialize(preserveState=True)
>>> cv4.getValue(context)
3880.8... dimensionless
>>> 1 * cv1.getValue(context) - 0.5 * cv2.getValue(context)
3880.8...

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