Coverage for cvpack/helix_hbond_content.py: 94%
35 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 00:10 +0000
1"""
2.. class:: HelixHBondContent
3 :platform: Linux, MacOS, Windows
4 :synopsis: Alpha-helix hydrogen-bond content of a sequence of residues
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10from __future__ import annotations
12import re as regex
13import typing as t
15import openmm
16from openmm import app as mmapp
17from openmm import unit as mmunit
19from .collective_variable import CollectiveVariable
20from .units import Quantity, ScalarQuantity
21from .utils import check_contiguous_residues
24class HelixHBondContent(CollectiveVariable, openmm.CustomBondForce):
25 r"""
26 The alpha-helix hydrogen-bond content of a sequence of `n` residues:
28 .. math::
30 \alpha_{\rm HB}({\bf r}) = \sum_{i=5}^n B_m\left(
31 \frac{\| {\bf r}^{\rm H}_i - {\bf r}^{\rm O}_{i-4} \|}{d_{\rm HB}}
32 \right)
34 where :math:`{\bf r}^{\rm H}_k` and :math:`{\bf r}^{\rm O}_k` are the positions
35 of the hydrogen and oxygen atoms bonded, respectively, to the backbone nitrogen and
36 carbon atoms of residue :math:`k`. In addition, :math:`d_{\rm HB}` is the threshold
37 distance for a hydrogen bond and :math:`B_m(x)` is a smooth step function given by
39 .. math::
40 B_m(x) = \frac{1}{1 + x^{2m}}
42 where :math:`m` is an integer parameter that controls its steepness.
44 Optionally, this collective variable can be normalized to the range :math:`[0, 1]`.
46 .. note::
48 The residues must be a contiguous sequence from a single chain, ordered from the
49 N- to the C-terminus. Proline residues, which have no amide hydrogen, are
50 skipped as donors, and the normalization uses the number of remaining
51 donor-acceptor pairs.
53 Parameters
54 ----------
55 residues
56 The residues in the sequence.
57 pbc
58 Whether to use periodic boundary conditions.
59 thresholdDistance
60 The threshold distance for a hydrogen bond.
61 halfExponent
62 The parameter :math:`m` of the step function.
63 normalize
64 Whether to normalize the collective variable to the range :math:`[0, 1]`.
65 name
66 The name of the collective variable.
68 Example
69 -------
70 >>> import cvpack
71 >>> import openmm
72 >>> from openmm import app, unit
73 >>> from openmmtools import testsystems
74 >>> model = testsystems.LysozymeImplicit()
75 >>> residues = list(model.topology.residues())[59:80]
76 >>> print(*[r.name for r in residues])
77 LYS ASP GLU ... ILE LEU ARG
78 >>> helix_content = cvpack.HelixHBondContent(residues)
79 >>> helix_content.addToSystem(model.system)
80 >>> platform = openmm.Platform.getPlatformByName('Reference')
81 >>> integrator = openmm.VerletIntegrator(0)
82 >>> context = openmm.Context(model.system, integrator, platform)
83 >>> context.setPositions(model.positions)
84 >>> helix_content.getValue(context)
85 15.880... dimensionless
86 """
88 def __init__(
89 self,
90 residues: t.Sequence[mmapp.topology.Residue],
91 pbc: bool = False,
92 thresholdDistance: ScalarQuantity = Quantity(0.33 * mmunit.nanometers),
93 halfExponent: int = 3,
94 normalize: bool = False,
95 name: str = "helix_hbond_content",
96 ) -> None:
97 residues = list(residues)
98 check_contiguous_residues(residues)
99 if normalize and len(residues) < 5: # noqa: PLR2004
100 raise ValueError("Normalization requires at least 5 residues")
102 def find_atom(residue: mmapp.topology.Residue, pattern: t.Pattern) -> int:
103 for atom in residue.atoms():
104 if regex.match(pattern, atom.name):
105 return atom.index
106 raise ValueError(
107 f"Could not find atom matching regex "
108 f"'{pattern.pattern}'"
109 f" in residue {residue.name}{residue.id}"
110 )
112 pairs = [
113 (residues[i - 4], residues[i])
114 for i in range(4, len(residues))
115 if residues[i].name != "PRO" # proline has no amide hydrogen
116 ]
117 if normalize and not pairs:
118 raise ValueError("No hydrogen-bond donors found in the residue sequence")
119 numerator = 1 / len(pairs) if normalize else 1
120 threshold = thresholdDistance
121 if mmunit.is_quantity(threshold):
122 threshold = threshold.value_in_unit_system(mmunit.md_unit_system)
123 super().__init__(f"{numerator}/(1+x^{2 * halfExponent}); x=r/{threshold}")
124 hydrogen_pattern = regex.compile(r"\b(H|1H|HN1|HT1|H1|HN)\b")
125 oxygen_pattern = regex.compile(r"\b(O|OCT1|OC1|OT1|O1)\b")
126 for acceptor, donor in pairs:
127 self.addBond(
128 find_atom(acceptor, oxygen_pattern),
129 find_atom(donor, hydrogen_pattern),
130 [],
131 )
132 self.setUsesPeriodicBoundaryConditions(pbc)
133 self._registerCV(
134 name,
135 mmunit.dimensionless,
136 residues=residues,
137 pbc=pbc,
138 thresholdDistance=thresholdDistance,
139 halfExponent=halfExponent,
140 normalize=normalize,
141 )
144HelixHBondContent.registerTag("!cvpack.HelixHBondContent")