SheetRMSDContent

class cvpack.SheetRMSDContent(residues, numAtoms, parallel=False, blockSizes=None, thresholdRMSD=0.08 nm, stepFunction='(1+x^4)/(1+x^4+x^8)', normalize=False, name='sheet_rmsd_content')[source]

The beta-sheet RMSD content of n residues [3]:

\[\beta_{\rm rmsd}({\bf r}) = \sum_{i=1}^{n-2-d_{\rm min}} \sum_{j=i+d_{\rm min}}^{n-2} S\left( \frac{r_{\rm rmsd}({\bf g}_{i,j}, {\bf g}_{\rm ref})}{r_0} \right)\]

where \({\bf g}_{i,j}\) represents the coordinates of a group of atoms selected from residues i to i+2 and residues j to j+2 of the sequence, \(d_{\rm min}\) is the minimum distance between integers i and j, \({\bf g}_{\rm ref}\) represents the coordinates of the same atoms in an ideal beta-sheet configuration, \(r_{\rm rmsd}({\bf g}, {\bf g}_{\rm ref})\) is the root-mean-square deviation (RMSD) between groups \({\bf g}\) and \({\bf g}_{\rm ref}\), \(r_0\) is a threshold RMSD value, and \(S(x)\) is a smooth step function whose default form is

\[S(x) = \frac{1 + x^4}{1 + x^4 + x^8}\]

The sequence of residues must be a contiguous subset of a protein chain, ordered from the N-terminus to the C-terminus. The beta-sheet RMSD content is defined for both antiparallel and parallel beta sheets, with different ideal configurations \({\bf g}_{\rm ref}\) and minimum distances \(d_{\rm min}\) in each case.

Every group \({\bf g}_{i,j}\) is formed by the N, \({\rm C}_\alpha\), \({\rm C}_\beta\), C, and O atoms of the six residues involved, thus comprising 30 atoms in total. In the case of glycine, the missing \({\rm C}_\beta\) atom is replaced by the hydrogen atom in its position (HA3 in PDB naming, HA1 in CHARMM naming), whereas PLUMED uses HA2. The RMSD is then defined as

\[r_{\rm rmsd}({\bf g}, {\bf g}_{\rm ref}) = \sqrt{\frac{1}{30} \sum_{k=1}^{30} \left\| \hat{\bf r}_k({\bf g}) - {\bf A}({\bf g})\hat{\bf r}_k({\bf g}_{\rm ref}) \right\|^2}\]

where \(\hat{\bf r}_k({\bf g})\) is the position of the \(k\)-th atom in a group \({\bf g}\) relative to the group’s center of geometry and \({\bf A}({\bf g})\) is the rotation matrix that minimizes the RMSD between \({\bf g}\) and \({\bf g}_{\rm ref}\).

In addition, one can choose to partition the sequence of residues into \(m\) blocks with \(n_1, n_2, \ldots, n_m\) residues. In this case, only groups \({\bf g}_{i,j}\) composed of three consecutive residues from a block and three consecutive residues from the next block will be considered. The definition of the RMSD content is then modified to

\[\beta_{\rm rmsd}({\bf r}) = \sum_{k=1}^{m-1} \sum_{i=l_{k-1}+1}^{l_{k}-2} \sum_{j=l_k+1}^{l_{k+1}-2} S\left( \frac{r_{\rm rmsd}({\bf g}_{i,j}, {\bf g}_{\rm ref})}{r_0} \right)\]

where \(l_k = \sum_{i=1}^k n_i\) is the index of the last residue in block \(k\). All blocks must be contiguous subsets of the same protein chain, ordered from the N-terminus to the C-terminus. However, each block can be separated from the next one by any number of residues.

Optionally, the beta-sheet RMSD content can be normalized to the range \([0, 1]\). This is done by dividing its value by the number of \({\bf g}_{i,j}\) groups. For a single, contiguous sequence of residues, this number is

\[N_{\rm groups} = \frac{(n - 2 - d_{\rm min})(n - 1 - d_{\rm min})}{2}\]

In the case of a partitioned sequence, the number of groups is given by

\[N_{\rm groups} = \sum_{k=1}^{m-1} (n_k - 2)(n_{k+1} - 2)\]

This collective variable was introduced in Ref. [3]. The default step function shown above is identical to the one in the original paper, but written in a numerically safe form. The ideal beta-sheet configurations are the same as those used for the collective variables ANTIBETARMSD and PARABETARMSD in PLUMED v2.8.1. As in PLUMED, there are two ideal parallel configurations, which differ mainly in which strand comes first, and the RMSD of a group in a parallel sheet is the smaller of its RMSDs to them.

Note

The present implementation is limited to \(1 \leq N_{\rm groups} \leq 1024\) for antiparallel sheets and \(1 \leq N_{\rm groups} \leq 512\) for parallel sheets.

Parameters:

  • residues (t.Sequence[mmapp.topology.Residue]) – The residue sequence or residue blocks to be used in the calculation.

  • numAtoms (int) – The total number of atoms in the system (required by OpenMM).

  • parallel (bool) – Whether to consider a parallel beta sheet instead of an antiparallel one.

  • blockSizes (t.Sequence[int] | None) – The number of residues in each block. If None, a single contiguous sequence of residues is assumed.

  • thresholdRMSD (ScalarQuantity) – The threshold RMSD value for considering a group of residues as a close match to an ideal beta sheet.

  • stepFunction (str) – The form of the step function \(S(x)\).

  • normalize (bool) – Whether to normalize the collective variable to the range \([0, 1]\).

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

Example

>>> import itertools as it
>>> import cvpack
>>> import openmm
>>> from openmm import app, unit
>>> from openmmtools import testsystems
>>> model = testsystems.SrcImplicit()
>>> residues = list(it.islice(model.topology.residues(), 68, 82))
>>> print(*[r.name for r in residues])
TYR ALA VAL ... VAL THR GLU
>>> sheet_content = cvpack.SheetRMSDContent(
...     residues, model.system.getNumParticles()
... )
>>> sheet_content.getNumResidueBlocks()
28
>>> sheet_content.addToSystem(model.system)
>>> platform = openmm.Platform.getPlatformByName('Reference')
>>> integrator = openmm.VerletIntegrator(0)
>>> context = openmm.Context(model.system, integrator, platform)
>>> context.setPositions(model.positions)
>>> sheet_content.getValue(context)
1.0465... dimensionless
>>> blockwise_sheet_content = cvpack.SheetRMSDContent(
...     residues[:5] + residues[-5:],
...     model.system.getNumParticles(),
...     blockSizes=[5, 5],
... )
>>> blockwise_sheet_content.getNumResidueBlocks()
9
>>> blockwise_sheet_content.addToSystem(model.system)
>>> context.reinitialize(preserveState=True)
>>> blockwise_sheet_content.getValue(context)
0.9859... 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 shares 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("not water and name N CA C")
>>> 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.04795... nm**2 Da/(rad**2)
>>> psi.getEffectiveMass(context)
0.05119... 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 shares 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("not water and name N CA C")
>>> 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