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 

10import typing as t 

11 

12import openmm 

13from openmm import app as mmapp 

14from openmm import unit as mmunit 

15 

16from .collective_variable import CollectiveVariable 

17from .units import Quantity, ScalarQuantity 

18 

19 

20class HelixAngleContent(CollectiveVariable, openmm.CustomAngleForce): 

21 r""" 

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

23 

24 .. math:: 

25 

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

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

28 \right) 

29 

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

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

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

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

34 

35 .. math:: 

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

37 

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

39 :math:`x` needs to be elevated to an even power for :math:`B_m(x)` to be an even 

40 function. 

41 

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

43 

44 .. note:: 

45 

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

47 N- to the C-terminus. Due to an OpenMM limitation, the maximum supported number 

48 of residues is 37. 

49 

50 Parameters 

51 ---------- 

52 residues 

53 The residues in the sequence. 

54 pbc 

55 Whether to use periodic boundary conditions. 

56 thetaReference 

57 The reference value of the 

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

59 helix. 

60 tolerance 

61 The threshold tolerance around the reference values. 

62 halfExponent 

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

64 normalize 

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

66 name 

67 The name of the collective variable. 

68 

69 Raises 

70 ------ 

71 ValueError 

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

73 

74 Example 

75 ------- 

76 >>> import cvpack 

77 >>> import openmm 

78 >>> from openmm import app, unit 

79 >>> from openmmtools import testsystems 

80 >>> model = testsystems.LysozymeImplicit() 

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

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

83 LYS ASP GLU ... ILE LEU ARG 

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

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

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

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

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

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

90 >>> helix_content.getValue(context) 

91 18.7605... dimensionless 

92 """ 

93 

94 def __init__( 

95 self, 

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

97 pbc: bool = False, 

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

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

100 halfExponent: int = 3, 

101 normalize: bool = False, 

102 name: str = "helix_angle_content", 

103 ) -> None: 

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

105 for atom in residue.atoms(): 

106 if atom.name == "CA": 

107 return atom.index 

108 raise ValueError( 

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

110 ) 

111 

112 num_angles = len(residues) - 2 

113 numerator = 1 / num_angles if normalize else 1 

114 theta0, tol = thetaReference, tolerance 

115 if mmunit.is_quantity(theta0): 

116 theta0 = theta0.value_in_unit(mmunit.radians) 

117 if mmunit.is_quantity(tol): 

118 tol = tol.value_in_unit(mmunit.radians) 

119 super().__init__( 

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

121 ) 

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

123 self.addAngle( 

124 find_alpha_carbon(residues[i - 1]), 

125 find_alpha_carbon(residues[i]), 

126 find_alpha_carbon(residues[i + 1]), 

127 [], 

128 ) 

129 self.setUsesPeriodicBoundaryConditions(pbc) 

130 self._registerCV( 

131 name, 

132 mmunit.dimensionless, 

133 residues=residues, 

134 pbc=pbc, 

135 thetaReference=thetaReference, 

136 tolerance=tolerance, 

137 halfExponent=halfExponent, 

138 normalize=normalize, 

139 ) 

140 

141 

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