Coverage for cvpack/helix_torsion_content.py: 100%

1""" 

2.. class:: HelixTorsionContent 

3 :platform: Linux, MacOS, Windows 

4 :synopsis: Alpha-helix torsion 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 HelixTorsionContent(CollectiveVariable, openmm.CustomTorsionForce): 

21 r""" 

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

23 

24 .. math:: 

25 

26 \alpha_{\phi,\psi}({\bf r}) = \frac{1}{2} \sum_{i=2}^{n-1} \left[ 

27 B_m\left( 

28 \frac{\phi_i({\bf r}) - \phi_{\rm ref}}{\theta_{\rm tol}} 

29 \right) + 

30 B_m\left( 

31 \frac{\psi_i({\bf r}) - \psi_{\rm ref}}{\theta_{\rm tol}} 

32 \right) 

33 \right] 

34 

35 where :math:`\phi_i({\bf r})` and :math:`\psi_i({\bf r})` are the Ramachandran 

36 dihedral angles of residue :math:`i`, :math:`\phi_{\rm ref}` and 

37 :math:`\psi_{\rm ref}` are their reference values in an alpha helix 

38 :cite:`Hovmoller_2002`, and :math:`\theta_{\rm tol}` is the threshold tolerance 

39 around these refenrences. :math:`B_m(x)` is a smooth boxcar function given by 

40 

41 .. math:: 

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

43 

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

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

46 function. 

47 

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

49 

50 .. note:: 

51 

52 The :math:`\phi` and :math:`\psi` angles of the first and last residues are 

53 not considered. They are used to compute the dihedral angles of the second and 

54 penultimate residues, respectively. 

55 

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

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

58 of residues is 37. 

59 

60 Parameters 

61 ---------- 

62 residues 

63 The residues in the sequence. 

64 pbc 

65 Whether to use periodic boundary conditions. 

66 phiReference 

67 The reference value of the phi dihedral angle in an alpha helix. 

68 psiReference 

69 The reference value of the psi dihedral angle in an alpha helix. 

70 tolerance 

71 The threshold tolerance around the reference values. 

72 halfExponent 

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

74 normalize 

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

76 name 

77 The name of the collective variable. 

78 

79 Raises 

80 ------ 

81 ValueError 

82 If some residue does not contain a :math:`\phi` or :math:`\psi` angle. 

83 

84 Example 

85 ------- 

86 >>> import cvpack 

87 >>> import openmm 

88 >>> from openmm import app, unit 

89 >>> from openmmtools import testsystems 

90 >>> model = testsystems.LysozymeImplicit() 

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

92 >>> print(*[r.name for r in residues]) 

93 LYS ASP GLU ... ILE LEU ARG 

94 >>> helix_content = cvpack.HelixTorsionContent(residues) 

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

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

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

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

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

100 >>> helix_content.getValue(context) 

101 17.452... dimensionless 

102 """ 

103 

104 def __init__( 

105 self, 

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

107 pbc: bool = False, 

108 phiReference: ScalarQuantity = Quantity(-63.8 * mmunit.degrees), 

109 psiReference: ScalarQuantity = Quantity(-41.1 * mmunit.degrees), 

110 tolerance: ScalarQuantity = Quantity(25 * mmunit.degrees), 

111 halfExponent: int = 3, 

112 normalize: bool = False, 

113 name: str = "helix_torsion_content", 

114 ) -> None: 

115 def find_atom(residue: mmapp.topology.Residue, name: str) -> int: 

116 for atom in residue.atoms(): 

117 if atom.name == name: 

118 return atom.index 

119 raise ValueError( 

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

121 ) 

122 

123 numerator = 1 / (2 * (len(residues) - 2)) if normalize else 1 / 2 

124 tol = tolerance 

125 if mmunit.is_quantity(tol): 

126 tol = tol.value_in_unit_system(mmunit.md_unit_system) 

127 super().__init__( 

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

129 ) 

130 self.addPerTorsionParameter("theta_ref") 

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

132 self.addTorsion( 

133 find_atom(residues[i - 1], "C"), 

134 find_atom(residues[i], "N"), 

135 find_atom(residues[i], "CA"), 

136 find_atom(residues[i], "C"), 

137 [phiReference], 

138 ) 

139 self.addTorsion( 

140 find_atom(residues[i], "N"), 

141 find_atom(residues[i], "CA"), 

142 find_atom(residues[i], "C"), 

143 find_atom(residues[i + 1], "N"), 

144 [psiReference], 

145 ) 

146 self.setUsesPeriodicBoundaryConditions(pbc) 

147 self._registerCV( 

148 name, 

149 mmunit.dimensionless, 

150 residues=residues, 

151 pbc=pbc, 

152 phiReference=phiReference, 

153 psiReference=psiReference, 

154 tolerance=tolerance, 

155 halfExponent=halfExponent, 

156 ) 

157 

158 

159HelixTorsionContent.registerTag("!cvpack.HelixTorsionContent")