Coverage for cvpack/helix_torsion_content.py: 100%
32 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 00:10 +0000
1"""
2.. class:: HelixTorsionContent
3 :platform: Linux, MacOS, Windows
4 :synopsis: Alpha-helix torsion 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
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 HelixTorsionContent(CollectiveVariable, openmm.CustomTorsionForce):
25 r"""
26 The alpha-helix Ramachandran content of a sequence of `n` residues:
28 .. math::
30 \alpha_{\phi,\psi}({\bf r}) = \frac{1}{2} \sum_{i=2}^{n-1} \left[
31 B_m\left(
32 \frac{\phi_i({\bf r}) - \phi_{\rm ref}}{\theta_{\rm tol}}
33 \right) +
34 B_m\left(
35 \frac{\psi_i({\bf r}) - \psi_{\rm ref}}{\theta_{\rm tol}}
36 \right)
37 \right]
39 where :math:`\phi_i({\bf r})` and :math:`\psi_i({\bf r})` are the Ramachandran
40 dihedral angles of residue :math:`i`, :math:`\phi_{\rm ref}` and
41 :math:`\psi_{\rm ref}` are their reference values in an alpha helix
42 :cite:`Hovmoller_2002`, and :math:`\theta_{\rm tol}` is the threshold tolerance
43 around these references. :math:`B_m(x)` is a smooth boxcar function given by
45 .. math::
46 B_m(x) = \frac{1}{1 + x^{2m}}
48 where :math:`m` is an integer parameter that controls its steepness. Note that
49 :math:`x` needs to be raised to an even power for :math:`B_m(x)` to be an even
50 function.
52 Optionally, this collective variable can be normalized to the range :math:`[0, 1]`.
54 .. note::
56 The :math:`\phi` and :math:`\psi` angles of the first and last residues are
57 not considered. They are used to compute the dihedral angles of the second and
58 penultimate residues, respectively.
60 The residues must be a contiguous sequence from a single chain, ordered from the
61 N- to the C-terminus.
63 Parameters
64 ----------
65 residues
66 The residues in the sequence.
67 pbc
68 Whether to use periodic boundary conditions.
69 phiReference
70 The reference value of the phi dihedral angle in an alpha helix.
71 psiReference
72 The reference value of the psi dihedral angle in an alpha helix.
73 tolerance
74 The threshold tolerance around the reference values.
75 halfExponent
76 The parameter :math:`m` of the boxcar function.
77 normalize
78 Whether to normalize the collective variable to the range :math:`[0, 1]`.
79 name
80 The name of the collective variable.
82 Raises
83 ------
84 ValueError
85 If some residue does not contain a :math:`\phi` or :math:`\psi` angle.
87 Example
88 -------
89 >>> import cvpack
90 >>> import openmm
91 >>> from openmm import app, unit
92 >>> from openmmtools import testsystems
93 >>> model = testsystems.LysozymeImplicit()
94 >>> residues = list(model.topology.residues())[59:80]
95 >>> print(*[r.name for r in residues])
96 LYS ASP GLU ... ILE LEU ARG
97 >>> helix_content = cvpack.HelixTorsionContent(residues)
98 >>> helix_content.addToSystem(model.system)
99 >>> platform = openmm.Platform.getPlatformByName('Reference')
100 >>> integrator = openmm.VerletIntegrator(0)
101 >>> context = openmm.Context(model.system, integrator, platform)
102 >>> context.setPositions(model.positions)
103 >>> helix_content.getValue(context)
104 17.452... dimensionless
105 """
107 def __init__(
108 self,
109 residues: t.Sequence[mmapp.topology.Residue],
110 pbc: bool = False,
111 phiReference: ScalarQuantity = Quantity(-63.8 * mmunit.degrees),
112 psiReference: ScalarQuantity = Quantity(-41.1 * mmunit.degrees),
113 tolerance: ScalarQuantity = Quantity(25 * mmunit.degrees),
114 halfExponent: int = 3,
115 normalize: bool = False,
116 name: str = "helix_torsion_content",
117 ) -> None:
118 residues = list(residues)
119 check_contiguous_residues(residues)
120 if normalize and len(residues) < 3: # noqa: PLR2004
121 raise ValueError("Normalization requires at least 3 residues")
123 def find_atom(residue: mmapp.topology.Residue, name: str) -> int:
124 for atom in residue.atoms():
125 if atom.name == name:
126 return atom.index
127 raise ValueError(
128 f"Could not find atom {name} in residue {residue.name}{residue.id}"
129 )
131 numerator = 1 / (2 * (len(residues) - 2)) if normalize else 1 / 2
132 tol = tolerance
133 if mmunit.is_quantity(tol):
134 tol = tol.value_in_unit_system(mmunit.md_unit_system)
135 super().__init__(
136 f"{numerator}/(1+x^{2 * halfExponent})"
137 f"; x=(delta-{2 * np.pi}*floor(delta/{2 * np.pi}+0.5))/{tol}"
138 "; delta=theta-theta_ref"
139 )
140 self.addPerTorsionParameter("theta_ref")
141 for i in range(1, len(residues) - 1):
142 self.addTorsion(
143 find_atom(residues[i - 1], "C"),
144 find_atom(residues[i], "N"),
145 find_atom(residues[i], "CA"),
146 find_atom(residues[i], "C"),
147 [phiReference],
148 )
149 self.addTorsion(
150 find_atom(residues[i], "N"),
151 find_atom(residues[i], "CA"),
152 find_atom(residues[i], "C"),
153 find_atom(residues[i + 1], "N"),
154 [psiReference],
155 )
156 self.setUsesPeriodicBoundaryConditions(pbc)
157 self._registerCV(
158 name,
159 mmunit.dimensionless,
160 residues=residues,
161 pbc=pbc,
162 phiReference=phiReference,
163 psiReference=psiReference,
164 tolerance=tolerance,
165 halfExponent=halfExponent,
166 normalize=normalize,
167 )
170HelixTorsionContent.registerTag("!cvpack.HelixTorsionContent")