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