Coverage for cvpack/helix_rmsd_content.py: 100%
12 statements
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-09 16:14 +0000
1"""
2.. class:: HelixRMSDContent
3 :platform: Linux, MacOS, Windows
4 :synopsis: Alpha-helix RMSD content of a sequence of residues
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10import typing as t
12from openmm import app as mmapp
13from openmm import unit as mmunit
15from .base_rmsd_content import BaseRMSDContent
16from .units import Quantity, ScalarQuantity
18# pylint: disable=protected-access
19ALPHA_POSITIONS = BaseRMSDContent._loadPositions("ideal_alpha_helix.csv")
20# pylint: enable=protected-access
23class HelixRMSDContent(BaseRMSDContent):
24 r"""
25 The alpha-helix RMSD content of a sequence of `n` residues :cite:`Pietrucci_2009`:
27 .. math::
29 \alpha_{\rm rmsd}({\bf r}) = \sum_{i=1}^{n-5} S\left(
30 \frac{r_{\rm rmsd}({\bf g}_i, {\bf g}_{\rm ref})}{r_0}
31 \right)
33 where :math:`{\bf g}_i` represents a group of atoms selected from residues
34 :math:`i` to :math:`i+5` of the sequence, :math:`{\bf g}_{\rm ref}` represents
35 the same atoms in an ideal alpha-helix configuration,
36 :math:`r_{\rm rmsd}({\bf g}, {\bf g}_{\rm ref})` is the root-mean-square
37 distance (RMSD) between groups :math:`{\bf g}` and :math:`{\bf g}_{\rm ref}`,
38 :math:`r_0` is a threshold RMSD value, and :math:`S(x)` is a smooth step function
39 whose default form is
41 .. math::
42 S(x) = \frac{1 + x^4}{1 + x^4 + x^8}
44 The residues must be a contiguous sequence from a single chain, ordered from
45 the N-terminus to the C-terminus.
47 Every group :math:`{\bf g}_{i,j}` is formed by the N, :math:`{\rm C}_\alpha`,
48 :math:`{\rm C}_\beta`, C, and O atoms of the six residues involvend, thus
49 comprising 30 atoms in total. In the case of glycine, the missing
50 :math:`{\rm C}_\beta` atom is replaced by the corresponding H atom. The RMSD is
51 then defined as
53 .. math::
55 r_{\rm rmsd}({\bf g}, {\bf g}_{\rm ref}) =
56 \sqrt{\frac{1}{30} \sum_{j=1}^{30} \left\|
57 \hat{\bf r}_j({\bf g}) -
58 {\bf A}({\bf g})\hat{\bf r}_j({\bf g}_{\rm ref})
59 \right\|^2}
61 where :math:`\hat{\bf r}_k({\bf g})` is the position of the :math:`k`-th atom in
62 a group :math:`{\bf g}` relative to the group's center of geometry and
63 :math:`{\bf A}({\bf g})` is the rotation matrix that minimizes the RMSD between
64 :math:`{\bf g}` and :math:`{\bf g}_{\rm ref}`.
66 Optionally, the alpha-helix RMSD content can be normalized to the range
67 :math:`[0, 1]`. This is done by dividing its value by :math:`N_{\rm groups} =
68 n - 5`.
70 This collective variable was introduced in Ref. :cite:`Pietrucci_2009`. The default
71 step function shown above is identical to the one in the original paper, but written
72 in a numerically safe form. The ideal alpha-helix configuration is the same used for
73 the collective variable `ALPHARMSD`_ in `PLUMED v2.8.1
74 <https://github.com/plumed/plumed2>`_ .
76 .. note::
78 The present implementation is limited to :math:`1 \leq N_{\rm groups} \leq
79 1024`.
81 .. _ALPHARMSD: https://www.plumed.org/doc-v2.8/user-doc/html/_a_l_p_h_a_r_m_s_d.html
83 Parameters
84 ----------
85 residues
86 The residue sequence to be used in the calculation.
87 numAtoms
88 The total number of atoms in the system (required by OpenMM).
89 thresholdRMSD
90 The threshold RMSD value for considering a group of residues as a close
91 match to an ideal alpha-helix.
92 stepFunction
93 The form of the step function :math:`S(x)`.
94 normalize
95 Whether to normalize the collective variable to the range :math:`[0, 1]`.
96 name
97 The name of the collective variable.
99 Example
100 -------
101 >>> import itertools as it
102 >>> import cvpack
103 >>> import openmm
104 >>> from openmm import app, unit
105 >>> from openmmtools import testsystems
106 >>> model = testsystems.LysozymeImplicit()
107 >>> residues = list(it.islice(model.topology.residues(), 59, 80))
108 >>> print(*[r.name for r in residues])
109 LYS ASP GLU ... ILE LEU ARG
110 >>> helix_content = cvpack.HelixRMSDContent(
111 ... residues, model.system.getNumParticles()
112 ... )
113 >>> helix_content.getNumResidueBlocks()
114 16
115 >>> helix_content.addToSystem(model.system)
116 >>> normalized_helix_content = cvpack.HelixRMSDContent(
117 ... residues, model.system.getNumParticles(), normalize=True
118 ... )
119 >>> normalized_helix_content.addToSystem(model.system)
120 >>> platform = openmm.Platform.getPlatformByName('Reference')
121 >>> integrator = openmm.VerletIntegrator(0)
122 >>> context = openmm.Context(model.system, integrator, platform)
123 >>> context.setPositions(model.positions)
124 >>> helix_content.getValue(context)
125 15.98... dimensionless
126 >>> normalized_helix_content.getValue(context)
127 0.998... dimensionless
128 """
130 def __init__(
131 self,
132 residues: t.Sequence[mmapp.topology.Residue],
133 numAtoms: int,
134 thresholdRMSD: ScalarQuantity = Quantity(0.08 * mmunit.nanometers),
135 stepFunction: str = "(1+x^4)/(1+x^4+x^8)",
136 normalize: bool = False,
137 name: str = "helix_rmsd_content",
138 ) -> None:
139 residue_blocks = [
140 list(range(index, index + 6)) for index in range(len(residues) - 5)
141 ]
143 super().__init__(
144 residue_blocks,
145 ALPHA_POSITIONS,
146 residues,
147 numAtoms,
148 thresholdRMSD,
149 stepFunction,
150 normalize,
151 )
152 self._registerCV(
153 name,
154 mmunit.dimensionless,
155 residues=residues,
156 numAtoms=numAtoms,
157 thresholdRMSD=thresholdRMSD,
158 stepFunction=stepFunction,
159 normalize=normalize,
160 )
163HelixRMSDContent.registerTag("!cvpack.HelixRMSDContent")