Coverage for cvpack/helix_hbond_content.py: 100%

1""" 

2.. class:: HelixHBondContent 

3 :platform: Linux, MacOS, Windows 

4 :synopsis: Alpha-helix hydrogen-bond content of a sequence of residues 

5 

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

7 

8""" 

9 

10import re as regex 

11import typing as t 

12 

13import openmm 

14from openmm import app as mmapp 

15from openmm import unit as mmunit 

16 

17from .collective_variable import CollectiveVariable 

18from .units import Quantity, ScalarQuantity 

19 

20 

21class HelixHBondContent(CollectiveVariable, openmm.CustomBondForce): 

22 r""" 

23 The alpha-helix hydrogen-bond content of a sequence of `n` residues: 

24 

25 .. math:: 

26 

27 \alpha_{\rm HB}({\bf r}) = \sum_{i=5}^n B_m\left( 

28 \frac{\| {\bf r}^{\rm H}_i - {\bf r}^{\rm O}_{i-4} \|}{d_{\rm HB}} 

29 \right) 

30 

31 where :math:`{\bf r}^{\rm H}_k` and :math:`{\bf r}^{\rm O}_k` are the positions 

32 of the hydrogen and oxygen atoms bonded, respectively, to the backbone nitrogen and 

33 carbon atoms of residue :math:`k`. In addition, :math:`d_{\rm HB}` is the threshold 

34 distance for a hydrogen bond and :math:`B_m(x)` is a smooth step function given by 

35 

36 .. math:: 

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

38 

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

40 

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

42 

43 .. note:: 

44 

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

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

47 of residues is 37. 

48 

49 Parameters 

50 ---------- 

51 residues 

52 The residues in the sequence. 

53 pbc 

54 Whether to use periodic boundary conditions. 

55 thresholdDistance 

56 The threshold distance for a hydrogen bond. 

57 halfExponent 

58 The parameter :math:`m` of the step function. 

59 name 

60 The name of the collective variable. 

61 

62 Example 

63 ------- 

64 >>> import cvpack 

65 >>> import openmm 

66 >>> from openmm import app, unit 

67 >>> from openmmtools import testsystems 

68 >>> model = testsystems.LysozymeImplicit() 

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

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

71 LYS ASP GLU ... ILE LEU ARG 

72 >>> helix_content = cvpack.HelixHBondContent(residues) 

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

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

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

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

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

78 >>> helix_content.getValue(context) 

79 15.880... dimensionless 

80 """ 

81 

82 def __init__( 

83 self, 

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

85 pbc: bool = False, 

86 thresholdDistance: ScalarQuantity = Quantity(0.33 * mmunit.nanometers), 

87 halfExponent: int = 3, 

88 normalize: bool = False, 

89 name: str = "helix_hbond_content", 

90 ) -> None: 

91 def find_atom(residue: mmapp.topology.Residue, pattern: t.Pattern) -> int: 

92 for atom in residue.atoms(): 

93 if regex.match(pattern, atom.name): 

94 return atom.index 

95 raise ValueError( 

96 f"Could not find atom matching regex " 

97 f"'{pattern.pattern}'" 

98 f" in residue {residue.name}{residue.id}" 

99 ) 

100 

101 numerator = 1 / (len(residues) - 4) if normalize else 1 

102 threshold = thresholdDistance 

103 if mmunit.is_quantity(threshold): 

104 threshold = threshold.value_in_unit_system(mmunit.md_unit_system) 

105 super().__init__(f"{numerator}/(1+x^{2*halfExponent}); x=r/{threshold}") 

106 hydrogen_pattern = regex.compile(r"\b(H|1H|HN1|HT1|H1|HN)\b") 

107 oxygen_pattern = regex.compile(r"\b(O|OCT1|OC1|OT1|O1)\b") 

108 for i in range(4, len(residues)): 

109 self.addBond( 

110 find_atom(residues[i - 4], oxygen_pattern), 

111 find_atom(residues[i], hydrogen_pattern), 

112 [], 

113 ) 

114 self.setUsesPeriodicBoundaryConditions(pbc) 

115 self._registerCV( 

116 name, 

117 mmunit.dimensionless, 

118 residues=residues, 

119 pbc=pbc, 

120 thresholdDistance=thresholdDistance, 

121 halfExponent=halfExponent, 

122 normalize=normalize, 

123 ) 

124 

125 

126HelixHBondContent.registerTag("!cvpack.HelixHBondContent")