Coverage for cvpack/helix_angle_content.py: 100%

1""" 

2.. class:: HelixAngleContent 

3 :platform: Linux, MacOS, Windows 

4 :synopsis: Alpha-helix angle content of a sequence of residues 

5 

6.. classauthor:: Charlles Abreu <craabreu@gmail.com> 

7 

8""" 

9 

10from __future__ import annotations 

11 

12import typing as t 

13 

14import openmm 

15from openmm import app as mmapp 

16from openmm import unit as mmunit 

17 

18from .collective_variable import CollectiveVariable 

19from .units import Quantity, ScalarQuantity 

20from .utils import check_contiguous_residues 

21 

22 

23class HelixAngleContent(CollectiveVariable, openmm.CustomAngleForce): 

24 r""" 

25 The alpha-helix angle content of a sequence of `n` residues: 

26 

27 .. math:: 

28 

29 \alpha_{\theta}({\bf r}) = \sum_{i=2}^{n-1} B_m\left( 

30 \frac{\theta^\alpha_i({\bf r}) - \theta_{\rm ref}}{\theta_{\rm tol}} 

31 \right) 

32 

33 where :math:`\theta^\alpha_i` is the angle formed by the alpha-carbon atoms of 

34 residues :math:`i-1`, :math:`i`, and :math:`i+1`, :math:`\theta_{\rm ref}` is its 

35 reference value in an alpha helix, and :math:`\theta_{\rm tol}` is the threshold 

36 tolerance around this reference. :math:`B_m(x)` is a smooth boxcar function given by 

37 

38 .. math:: 

39 B_m(x) = \frac{1}{1 + x^{2m}} 

40 

41 where :math:`m` is an integer parameter that controls its steepness. Note that 

42 :math:`x` needs to be raised to an even power for :math:`B_m(x)` to be an even 

43 function. 

44 

45 Optionally, this collective variable can be normalized to the range :math:`[0, 1]`. 

46 

47 .. note:: 

48 

49 The residues must be a contiguous sequence from a single chain, ordered from the 

50 N- to the C-terminus. 

51 

52 Parameters 

53 ---------- 

54 residues 

55 The residues in the sequence. 

56 pbc 

57 Whether to use periodic boundary conditions. 

58 thetaReference 

59 The reference value of the 

60 :math:`{\rm C}_\alpha{\rm -C}_\alpha{\rm -C}_\alpha` angle in an alpha 

61 helix. 

62 tolerance 

63 The threshold tolerance around the reference values. 

64 halfExponent 

65 The parameter :math:`m` of the boxcar function. 

66 normalize 

67 Whether to normalize the collective variable to the range :math:`[0, 1]`. 

68 name 

69 The name of the collective variable. 

70 

71 Raises 

72 ------ 

73 ValueError 

74 If some residue does not contain a :math:`{\rm C}_\alpha` atom 

75 

76 Example 

77 ------- 

78 >>> import cvpack 

79 >>> import openmm 

80 >>> from openmm import app, unit 

81 >>> from openmmtools import testsystems 

82 >>> model = testsystems.LysozymeImplicit() 

83 >>> residues = list(model.topology.residues())[59:80] 

84 >>> print(*[r.name for r in residues]) # doctest: +ELLIPSIS 

85 LYS ASP GLU ... ILE LEU ARG 

86 >>> helix_content = cvpack.HelixAngleContent(residues) 

87 >>> helix_content.addToSystem(model.system) 

88 >>> platform = openmm.Platform.getPlatformByName('Reference') 

89 >>> integrator = openmm.VerletIntegrator(0) 

90 >>> context = openmm.Context(model.system, integrator, platform) 

91 >>> context.setPositions(model.positions) 

92 >>> helix_content.getValue(context) 

93 18.7605... dimensionless 

94 """ 

95 

96 def __init__( 

97 self, 

98 residues: t.Sequence[mmapp.topology.Residue], 

99 pbc: bool = False, 

100 thetaReference: ScalarQuantity = Quantity(88 * mmunit.degrees), 

101 tolerance: ScalarQuantity = Quantity(15 * mmunit.degrees), 

102 halfExponent: int = 3, 

103 normalize: bool = False, 

104 name: str = "helix_angle_content", 

105 ) -> None: 

106 residues = list(residues) 

107 check_contiguous_residues(residues) 

108 if normalize and len(residues) < 3: # noqa: PLR2004 

109 raise ValueError("Normalization requires at least 3 residues") 

110 

111 def find_alpha_carbon(residue: mmapp.topology.Residue) -> int: 

112 for atom in residue.atoms(): 

113 if atom.name == "CA": 

114 return atom.index 

115 raise ValueError( 

116 f"Could not find atom CA in residue {residue.name}{residue.id}" 

117 ) 

118 

119 num_angles = len(residues) - 2 

120 numerator = 1 / num_angles if normalize else 1 

121 theta0, tol = thetaReference, tolerance 

122 if mmunit.is_quantity(theta0): 

123 theta0 = theta0.value_in_unit(mmunit.radians) 

124 if mmunit.is_quantity(tol): 

125 tol = tol.value_in_unit(mmunit.radians) 

126 super().__init__( 

127 f"{numerator}/(1+x^{2 * halfExponent}); x=(theta-{theta0})/{tol}" 

128 ) 

129 for i in range(1, len(residues) - 1): 

130 self.addAngle( 

131 find_alpha_carbon(residues[i - 1]), 

132 find_alpha_carbon(residues[i]), 

133 find_alpha_carbon(residues[i + 1]), 

134 [], 

135 ) 

136 self.setUsesPeriodicBoundaryConditions(pbc) 

137 self._registerCV( 

138 name, 

139 mmunit.dimensionless, 

140 residues=residues, 

141 pbc=pbc, 

142 thetaReference=thetaReference, 

143 tolerance=tolerance, 

144 halfExponent=halfExponent, 

145 normalize=normalize, 

146 ) 

147 

148 

149HelixAngleContent.registerTag("!cvpack.HelixAngleContent")