Coverage for cvpack/rmsd.py: 100%

1""" 

2.. class:: RMSD 

3 :platform: Linux, MacOS, Windows 

4 :synopsis: Root-mean-square deviation of a group of atoms 

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 unit as mmunit 

16 

17from .base_rmsd import BaseRMSD 

18from .units import MatrixQuantity, VectorQuantity 

19 

20 

21class RMSD(BaseRMSD, openmm.RMSDForce): 

22 r""" 

23 The root-mean-square deviation (RMSD) between the current and reference 

24 coordinates of a group of `n` atoms: 

25 

26 .. math:: 

27 

28 d_{\rm rms}({\bf r}) = \sqrt{ 

29 \frac{1}{n} \min_{ 

30 \bf q \in \mathbb{R}^4 \atop \|{\bf q}\| = 1 

31 } \sum_{i=1}^n \left\| 

32 {\bf A}({\bf q}) \hat{\bf r}_i - \hat{\bf r}_i^{\rm ref} 

33 \right\|^2 

34 } 

35 

36 where :math:`\hat{\bf r}_i` is the position of the :math:`i`-th atom in the group 

37 relative to the group's center of geometry (centroid), the superscript 

38 :math:`\rm ref` denotes the reference structure, :math:`{\bf q}` is a unit 

39 quaternion, and :math:`{\bf A}({\bf q})` is the rotation matrix corresponding to 

40 :math:`{\bf q}`. 

41 

42 .. warning:: 

43 

44 Periodic boundary conditions are `not supported 

45 <https://github.com/openmm/openmm/issues/2913>`_. This is not a problem if all 

46 atoms in the group belong to the same molecule. If they belong to distinct 

47 molecules, it is possible to circumvent the issue by calling the method 

48 :func:`getNullBondForce` and adding the resulting force to the system. 

49 

50 Parameters 

51 ---------- 

52 referencePositions 

53 The reference coordinates, which can be either a coordinate matrix or a mapping 

54 from atom indices to coordinate vectors. It must contain all atoms in ``group``, 

55 and does not need to contain all atoms in the system. See ``numAtoms`` below. 

56 group 

57 A sequence of atom indices. 

58 numAtoms 

59 The total number of atoms in the system, including those that are not in 

60 ``group``. This argument is necessary only if ``referencePositions`` does not 

61 contain all atoms in the system. 

62 name 

63 The name of the collective variable. 

64 

65 Raises 

66 ------ 

67 ValueError 

68 If referencePositions is a dict and numAtoms is not given 

69 

70 Example 

71 ------- 

72 >>> import cvpack 

73 >>> import openmm 

74 >>> from openmm import app, unit 

75 >>> from openmmtools import testsystems 

76 >>> model = testsystems.AlanineDipeptideImplicit() 

77 >>> num_atoms = model.topology.getNumAtoms() 

78 >>> group = list(range(num_atoms)) 

79 >>> rmsd = cvpack.RMSD(model.positions, group, num_atoms) 

80 >>> rmsd.addToSystem(model.system) 

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

82 >>> integrator = openmm.VerletIntegrator(2*unit.femtoseconds) 

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

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

85 >>> integrator.step(1000) 

86 >>> rmsd.getValue(context) 

87 0.123138... nm 

88 

89 """ 

90 

91 def __init__( 

92 self, 

93 referencePositions: MatrixQuantity | dict[int, VectorQuantity], 

94 group: t.Iterable[int], 

95 numAtoms: int | None = None, 

96 name: str = "rmsd", 

97 ) -> None: 

98 group = list(map(int, group)) 

99 if not group: 

100 # OpenMM's RMSDForce would interpret an empty group as all atoms 

101 raise ValueError("The group of atoms must not be empty") 

102 if numAtoms is None and isinstance(referencePositions, dict): 

103 raise ValueError("numAtoms is required when referencePositions is a dict") 

104 num_atoms = numAtoms or len(referencePositions) 

105 defined_coords = self._getDefinedCoords(referencePositions, group) 

106 all_coords = self._getAllCoords(defined_coords, num_atoms) 

107 super().__init__(all_coords, group) 

108 self._registerCV( 

109 name, 

110 mmunit.nanometers, 

111 referencePositions=defined_coords, 

112 group=group, 

113 numAtoms=num_atoms, 

114 ) 

115 

116 def getNullBondForce(self) -> openmm.HarmonicBondForce: 

117 """ 

118 Get a null bond force that creates a connected graph with all the atoms in the 

119 group. 

120 

121 Returns 

122 ------- 

123 openmm.HarmonicBondForce 

124 The null bond force. 

125 

126 Example 

127 ------- 

128 >>> import cvpack 

129 >>> import openmm 

130 >>> from openmm import app, unit 

131 >>> from openmmtools import testsystems 

132 >>> model = testsystems.WaterBox( 

133 ... box_edge=10*unit.angstroms, 

134 ... cutoff=5*unit.angstroms, 

135 ... ) 

136 >>> group = [ 

137 ... atom.index 

138 ... for atom in model.topology.atoms() 

139 ... if atom.residue.index < 3 

140 ... ] 

141 >>> rmsd = cvpack.RMSD( 

142 ... model.positions, group, model.topology.getNumAtoms() 

143 ... ) 

144 >>> rmsd.addToSystem(model.system) 

145 >>> _ = model.system.addForce(rmsd.getNullBondForce()) 

146 >>> integrator = openmm.VerletIntegrator(2*unit.femtoseconds) 

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

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

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

150 >>> integrator.step(100) 

151 >>> rmsd.getValue(context) 

152 0.10436... nm 

153 

154 """ 

155 force = openmm.HarmonicBondForce() 

156 group = self._args["group"] 

157 for i, j in zip(group[:-1], group[1:]): 

158 force.addBond(i, j, 0.0, 0.0) 

159 return force 

160 

161 

162RMSD.registerTag("!cvpack.RMSD")