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 

10import typing as t 

11 

12import openmm 

13from openmm import unit as mmunit 

14 

15from .base_rmsd import BaseRMSD 

16from .units import MatrixQuantity, VectorQuantity 

17 

18 

19class RMSD(BaseRMSD, openmm.RMSDForce): 

20 r""" 

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

22 coordinates of a group of `n` atoms: 

23 

24 .. math:: 

25 

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

27 \frac{1}{n} \min_{ 

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

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

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

31 \right\|^2 

32 } 

33 

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

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

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

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

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

39 

40 .. warning:: 

41 

42 Periodic boundary conditions are `not supported 

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

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

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

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

47 

48 Parameters 

49 ---------- 

50 referencePositions 

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

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

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

54 groups 

55 A sequence of atom indices. 

56 numAtoms 

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

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

59 contain all atoms in the system. 

60 name 

61 The name of the collective variable. 

62 

63 Example 

64 ------- 

65 >>> import cvpack 

66 >>> import openmm 

67 >>> from openmm import app, unit 

68 >>> from openmmtools import testsystems 

69 >>> model = testsystems.AlanineDipeptideImplicit() 

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

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

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

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

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

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

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

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

78 >>> integrator.step(1000) 

79 >>> rmsd.getValue(context) 

80 0.123138... nm 

81 

82 """ 

83 

84 def __init__( 

85 self, 

86 referencePositions: t.Union[MatrixQuantity, t.Dict[int, VectorQuantity]], 

87 group: t.Iterable[int], 

88 numAtoms: t.Optional[int] = None, 

89 name: str = "rmsd", 

90 ) -> None: 

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

92 num_atoms = numAtoms or len(referencePositions) 

93 defined_coords = self._getDefinedCoords(referencePositions, group) 

94 all_coords = self._getAllCoords(defined_coords, num_atoms) 

95 super().__init__(all_coords, group) 

96 self._registerCV( 

97 name, 

98 mmunit.nanometers, 

99 referencePositions=defined_coords, 

100 group=group, 

101 numAtoms=num_atoms, 

102 ) 

103 

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

105 """ 

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

107 group. 

108 

109 Returns 

110 ------- 

111 openmm.HarmonicBondForce 

112 The null bond force. 

113 

114 Example 

115 ------- 

116 >>> import cvpack 

117 >>> import openmm 

118 >>> from openmm import app, unit 

119 >>> from openmmtools import testsystems 

120 >>> model = testsystems.WaterBox( 

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

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

123 ... ) 

124 >>> group = [ 

125 ... atom.index 

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

127 ... if atom.residue.index < 3 

128 ... ] 

129 >>> rmsd = cvpack.RMSD( 

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

131 ... ) 

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

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

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

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

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

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

138 >>> integrator.step(100) 

139 >>> rmsd.getValue(context) 

140 0.10436... nm 

141 

142 """ 

143 force = openmm.HarmonicBondForce() 

144 group = self._args["group"] 

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

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

147 return force 

148 

149 

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