Coverage for cvpack/sheet_rmsd_content.py: 96%
27 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 00:10 +0000
1"""
2.. class:: SheetRMSDContent
3 :platform: Linux, MacOS, Windows
4 :synopsis: Beta-sheet RMSD content of a sequence of residues
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10from __future__ import annotations
12import typing as t
14import numpy as np
15from openmm import app as mmapp
16from openmm import unit as mmunit
18from .base_rmsd_content import BaseRMSDContent
19from .units import Quantity, ScalarQuantity
20from .utils import check_contiguous_residues
22PARABETA_POSITIONS = BaseRMSDContent._loadPositions("ideal_parallel_beta_sheet.csv")
23PARABETA_SWAPPED_POSITIONS = BaseRMSDContent._loadPositions(
24 "ideal_parallel_beta_sheet_swapped.csv"
25)
26ANTIBETA_POSITIONS = BaseRMSDContent._loadPositions("ideal_antiparallel_beta_sheet.csv")
29class SheetRMSDContent(BaseRMSDContent):
30 r"""
31 The beta-sheet RMSD content of `n` residues :cite:`Pietrucci_2009`:
33 .. math::
35 \beta_{\rm rmsd}({\bf r}) = \sum_{i=1}^{n-2-d_{\rm min}}
36 \sum_{j=i+d_{\rm min}}^{n-2} S\left(
37 \frac{r_{\rm rmsd}({\bf g}_{i,j}, {\bf g}_{\rm ref})}{r_0}
38 \right)
40 where :math:`{\bf g}_{i,j}` represents the coordinates of a group of atoms
41 selected from residues i to i+2 and residues j to j+2 of the sequence,
42 :math:`d_{\rm min}` is the minimum distance between integers i and j,
43 :math:`{\bf g}_{\rm ref}` represents the coordinates of the same atoms in an
44 ideal beta-sheet configuration, :math:`r_{\rm rmsd}({\bf g}, {\bf g}_{\rm ref})`
45 is the root-mean-square deviation (RMSD) between groups :math:`{\bf g}` and
46 :math:`{\bf g}_{\rm ref}`, :math:`r_0` is a threshold RMSD value, and
47 :math:`S(x)` is a smooth step function whose default form is
49 .. math::
50 S(x) = \frac{1 + x^4}{1 + x^4 + x^8}
52 The sequence of residues must be a contiguous subset of a protein chain, ordered
53 from the N-terminus to the C-terminus. The beta-sheet RMSD content is defined for
54 both antiparallel and parallel beta sheets, with different ideal configurations
55 :math:`{\bf g}_{\rm ref}` and minimum distances :math:`d_{\rm min}` in each case.
57 Every group :math:`{\bf g}_{i,j}` is formed by the N, :math:`{\rm C}_\alpha`,
58 :math:`{\rm C}_\beta`, C, and O atoms of the six residues involved, thus
59 comprising 30 atoms in total. In the case of glycine, the missing
60 :math:`{\rm C}_\beta` atom is replaced by the hydrogen atom in its position (HA3 in
61 PDB naming, HA1 in CHARMM naming), whereas PLUMED uses HA2. The RMSD is then
62 defined as
64 .. math::
66 r_{\rm rmsd}({\bf g}, {\bf g}_{\rm ref}) =
67 \sqrt{\frac{1}{30} \sum_{k=1}^{30} \left\|
68 \hat{\bf r}_k({\bf g}) -
69 {\bf A}({\bf g})\hat{\bf r}_k({\bf g}_{\rm ref})
70 \right\|^2}
72 where :math:`\hat{\bf r}_k({\bf g})` is the position of the :math:`k`-th atom in
73 a group :math:`{\bf g}` relative to the group's center of geometry and
74 :math:`{\bf A}({\bf g})` is the rotation matrix that minimizes the RMSD between
75 :math:`{\bf g}` and :math:`{\bf g}_{\rm ref}`.
77 In addition, one can choose to partition the sequence of residues into :math:`m`
78 blocks with :math:`n_1, n_2, \ldots, n_m` residues. In this case, only groups
79 :math:`{\bf g}_{i,j}` composed of three consecutive residues from a block and
80 three consecutive residues from the next block will be considered. The definition
81 of the RMSD content is then modified to
83 .. math::
85 \beta_{\rm rmsd}({\bf r}) = \sum_{k=1}^{m-1} \sum_{i=l_{k-1}+1}^{l_{k}-2}
86 \sum_{j=l_k+1}^{l_{k+1}-2} S\left(
87 \frac{r_{\rm rmsd}({\bf g}_{i,j}, {\bf g}_{\rm ref})}{r_0}
88 \right)
90 where :math:`l_k = \sum_{i=1}^k n_i` is the index of the last residue in block
91 :math:`k`. All blocks must be contiguous subsets of the same protein chain, ordered
92 from the N-terminus to the C-terminus. However, each block can be separated from the
93 next one by any number of residues.
95 Optionally, the beta-sheet RMSD content can be normalized to the range
96 :math:`[0, 1]`. This is done by dividing its value by the number of
97 :math:`{\bf g}_{i,j}` groups. For a single, contiguous sequence of residues, this
98 number is
100 .. math::
102 N_{\rm groups} = \frac{(n - 2 - d_{\rm min})(n - 1 - d_{\rm min})}{2}
104 In the case of a partitioned sequence, the number of groups is given by
106 .. math::
108 N_{\rm groups} = \sum_{k=1}^{m-1} (n_k - 2)(n_{k+1} - 2)
110 This collective variable was introduced in Ref. :cite:`Pietrucci_2009`. The default
111 step function shown above is identical to the one in the original paper, but written
112 in a numerically safe form. The ideal beta-sheet configurations are the same as
113 those used for the collective variables `ANTIBETARMSD`_ and `PARABETARMSD`_ in
114 `PLUMED v2.8.1 <https://github.com/plumed/plumed2>`_. As in PLUMED, there are two
115 ideal parallel configurations, which differ mainly in which strand comes first,
116 and the RMSD of a group in a parallel sheet is the smaller of its RMSDs to them.
118 .. _ANTIBETARMSD:
119 https://www.plumed.org/doc-v2.8/user-doc/html/_a_n_t_i_b_e_t_a_r_m_s_d.html
121 .. _PARABETARMSD:
122 https://www.plumed.org/doc-v2.8/user-doc/html/_p_a_r_a_b_e_t_a_r_m_s_d.html
124 .. note::
126 The present implementation is limited to :math:`1 \leq N_{\rm groups} \leq
127 1024` for antiparallel sheets and :math:`1 \leq N_{\rm groups} \leq 512` for
128 parallel sheets.
130 Parameters
131 ----------
132 residues
133 The residue sequence or residue blocks to be used in the calculation.
134 numAtoms
135 The total number of atoms in the system (required by OpenMM).
136 parallel
137 Whether to consider a parallel beta sheet instead of an antiparallel one.
138 blockSizes
139 The number of residues in each block. If ``None``, a single contiguous
140 sequence of residues is assumed.
141 thresholdRMSD
142 The threshold RMSD value for considering a group of residues as a close
143 match to an ideal beta sheet.
144 stepFunction
145 The form of the step function :math:`S(x)`.
146 normalize
147 Whether to normalize the collective variable to the range :math:`[0, 1]`.
148 name
149 The name of the collective variable.
151 Example
152 -------
153 >>> import itertools as it
154 >>> import cvpack
155 >>> import openmm
156 >>> from openmm import app, unit
157 >>> from openmmtools import testsystems
158 >>> model = testsystems.SrcImplicit()
159 >>> residues = list(it.islice(model.topology.residues(), 68, 82))
160 >>> print(*[r.name for r in residues])
161 TYR ALA VAL ... VAL THR GLU
162 >>> sheet_content = cvpack.SheetRMSDContent(
163 ... residues, model.system.getNumParticles()
164 ... )
165 >>> sheet_content.getNumResidueBlocks()
166 28
167 >>> sheet_content.addToSystem(model.system)
168 >>> platform = openmm.Platform.getPlatformByName('Reference')
169 >>> integrator = openmm.VerletIntegrator(0)
170 >>> context = openmm.Context(model.system, integrator, platform)
171 >>> context.setPositions(model.positions)
172 >>> sheet_content.getValue(context)
173 1.0465... dimensionless
174 >>> blockwise_sheet_content = cvpack.SheetRMSDContent(
175 ... residues[:5] + residues[-5:],
176 ... model.system.getNumParticles(),
177 ... blockSizes=[5, 5],
178 ... )
179 >>> blockwise_sheet_content.getNumResidueBlocks()
180 9
181 >>> blockwise_sheet_content.addToSystem(model.system)
182 >>> context.reinitialize(preserveState=True)
183 >>> blockwise_sheet_content.getValue(context)
184 0.9859... dimensionless
185 """
187 def __init__(
188 self,
189 residues: t.Sequence[mmapp.topology.Residue],
190 numAtoms: int,
191 parallel: bool = False,
192 blockSizes: t.Sequence[int] | None = None,
193 thresholdRMSD: ScalarQuantity = Quantity(0.08 * mmunit.nanometers),
194 stepFunction: str = "(1+x^4)/(1+x^4+x^8)",
195 normalize: bool = False,
196 name: str = "sheet_rmsd_content",
197 ) -> None:
198 residues = list(residues)
199 if blockSizes is None:
200 check_contiguous_residues(residues)
201 min_distance = 6 if parallel else 5
202 residue_groups = [
203 [i, i + 1, i + 2, j, j + 1, j + 2]
204 for i in range(len(residues) - 2 - min_distance)
205 for j in range(i + min_distance, len(residues) - 2)
206 ]
207 elif sum(blockSizes) == len(residues):
208 bounds = np.insert(np.cumsum(blockSizes), 0, 0)
209 for k in range(len(blockSizes)):
210 check_contiguous_residues(residues[bounds[k] : bounds[k + 1]])
211 residue_groups = [
212 [i, i + 1, i + 2, j, j + 1, j + 2]
213 for k in range(len(blockSizes) - 1)
214 for i in range(bounds[k], bounds[k + 1] - 2)
215 for j in range(bounds[k + 1], bounds[k + 2] - 2)
216 ]
217 else:
218 raise ValueError(
219 f"The sum of block sizes ({sum(blockSizes)}) and the "
220 f"number of residues ({len(residues)}) must be equal."
221 )
223 super().__init__(
224 residue_groups,
225 (
226 [PARABETA_POSITIONS, PARABETA_SWAPPED_POSITIONS]
227 if parallel
228 else [ANTIBETA_POSITIONS]
229 ),
230 residues,
231 numAtoms,
232 thresholdRMSD,
233 stepFunction,
234 normalize,
235 )
236 self._registerCV(
237 name,
238 mmunit.dimensionless,
239 residues=residues,
240 numAtoms=numAtoms,
241 parallel=parallel,
242 blockSizes=blockSizes,
243 thresholdRMSD=thresholdRMSD,
244 stepFunction=stepFunction,
245 normalize=normalize,
246 )
249SheetRMSDContent.registerTag("!cvpack.SheetRMSDContent")