Coverage for cvpack/path_in_rmsd_space.py: 95%

1""" 

2.. class:: PathInRMSDSpace 

3 :platform: Linux, MacOS, Windows 

4 :synopsis: A metric of progress or deviation with respect to a path in RMSD space 

5 

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

7 

8""" 

9 

10from __future__ import annotations 

11 

12import typing as t 

13 

14from openmm import unit as mmunit 

15 

16from .base_path_cv import BasePathCV 

17from .path import Metric, progress 

18from .rmsd import RMSD 

19from .units.units import VectorQuantity, value_in_md_units 

20 

21 

22class PathInRMSDSpace(BasePathCV): 

23 r""" 

24 A metric of the system's progress (:math:`s`) or deviation (:math:`z`) with 

25 respect to a path defined by a sequence of :math:`n` milestones defined as 

26 reference structures :cite:`Branduardi_2007`: 

27 

28 .. math:: 

29 

30 s({\bf r}) = \frac{ 

31 \dfrac{\sum_{i=1}^n i w_i({\bf r})}{\sum_{i=1}^n w_i({\bf r})} - 1 

32 }{n-1} 

33 \quad \text{or} \quad 

34 z({\bf r}) = - 2 \sigma ^2 \ln \sum_{i=1}^n w_i({\bf r}) 

35 

36 with :math:`w_i({\bf r})` being a Gaussian kernel centered at the :math:`i`-th 

37 milestone, i.e., 

38 

39 .. math:: 

40 

41 w_i({\bf r}) = \exp\left(\ 

42 -\frac{d^2_{\rm rms}({\bf r},{\bf r}^{\rm ref}_i)}{2 \sigma^2} 

43 \right) 

44 

45 where :math:`d_{\rm rms}({\bf r},{\bf r}^{\rm ref}_i)` is the root-mean-square 

46 deviation between the current system state and the :math:`i`-th reference structure 

47 and :math:`\sigma` sets the width of the kernels. 

48 

49 .. note:: 

50 

51 The kernel width :math:`\sigma` is related to the parameter :math:`\lambda` of 

52 Ref. :cite:`Branduardi_2007` by :math:`\sigma = \frac{1}{\sqrt{2\lambda}}`. 

53 

54 Parameters 

55 ---------- 

56 metric 

57 The path-related metric to compute. Use ``cvpack.path.progress`` for 

58 computing :math:`s({\bf r})` or ``cvpack.path.deviation`` for computing 

59 :math:`z({\bf r})`. 

60 milestones 

61 A sequence of reference structures, each represented as a dictionary mapping 

62 atom indices to coordinate vectors. 

63 numAtoms 

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

65 any of the reference structures. 

66 sigma 

67 The width of the Gaussian kernels in nanometers 

68 name 

69 The name of the collective variable. If not provided, it is set to 

70 "path_progress_in_rmsd_space" or "path_deviation_in_rmsd_space" depending 

71 on the metric. 

72 

73 Raises 

74 ------ 

75 ValueError 

76 If the number of milestones is less than 2 

77 ValueError 

78 If the metric is not `cvpack.path.progress` or `cvpack.path.deviation` 

79 

80 Examples 

81 -------- 

82 >>> import cvpack 

83 >>> import networkx as nx 

84 >>> import numpy as np 

85 >>> import openmm 

86 >>> from openmm import unit 

87 >>> from openmmtools import testsystems 

88 >>> from scipy.spatial.transform import Rotation 

89 >>> model = testsystems.AlanineDipeptideVacuum() 

90 >>> atom1, atom2 = 8, 14 

91 >>> graph = model.mdtraj_topology.to_bondgraph() 

92 >>> nodes = list(graph.nodes) 

93 >>> graph.remove_edge(nodes[atom1], nodes[atom2]) 

94 >>> movable = list(nx.connected_components(graph))[1] 

95 >>> x = model.positions / model.positions.unit 

96 >>> x0 = x[atom1, :] 

97 >>> vector = x[atom2, :] - x0 

98 >>> vector /= np.linalg.norm(vector) 

99 >>> rotation = Rotation.from_rotvec((np.pi / 6) * vector) 

100 >>> atoms = [nodes.index(atom) for atom in movable] 

101 >>> frames = [x.copy()] 

102 >>> for _ in range(6): 

103 ... x[atoms, :] = x0 + rotation.apply(x[atoms, :] - x0) 

104 ... frames.append(x.copy()) 

105 >>> milestones = [ 

106 ... {i: row for i, row in enumerate(frame)} 

107 ... for frame in frames 

108 ... ] 

109 >>> s, z = [ 

110 ... cvpack.PathInRMSDSpace( 

111 ... metric, milestones, len(x), 0.5 * unit.angstrom 

112 ... ) 

113 ... for metric in (cvpack.path.progress, cvpack.path.deviation) 

114 ... ] 

115 >>> s.addToSystem(model.system) 

116 >>> z.addToSystem(model.system) 

117 >>> context = openmm.Context(model.system, openmm.VerletIntegrator(0.001)) 

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

119 >>> s.getValue(context) 

120 0.172... dimensionless 

121 >>> z.getValue(context) 

122 -0.004... nm**2 

123 """ 

124 

125 def __init__( 

126 self, 

127 metric: Metric, 

128 milestones: t.Sequence[dict[int, VectorQuantity]], 

129 numAtoms: int, 

130 sigma: mmunit.Quantity, 

131 name: str | None = None, 

132 ) -> None: 

133 name = self._generateName(metric, name, "rmsd") 

134 if mmunit.is_quantity(sigma): 

135 sigma = sigma.value_in_unit(mmunit.nanometers) 

136 n = len(milestones) 

137 if n < 2: # noqa: PLR2004 

138 raise ValueError("At least two reference structures are required.") 

139 collective_variables = { 

140 f"rmsd{i}": RMSD(reference, reference.keys(), numAtoms, name=f"rmsd{i}") 

141 for i, reference in enumerate(milestones) 

142 } 

143 squared_distances = [f"rmsd{i}^2" for i in range(n)] 

144 super().__init__(metric, sigma, squared_distances, collective_variables) 

145 self._registerCV( 

146 name, 

147 mmunit.dimensionless if metric == progress else mmunit.nanometers**2, 

148 metric=metric, 

149 milestones=[ 

150 {int(k): [float(x) for x in value_in_md_units(v)] for k, v in m.items()} 

151 for m in milestones 

152 ], 

153 numAtoms=numAtoms, 

154 sigma=sigma, 

155 ) 

156 

157 

158PathInRMSDSpace.registerTag("!cvpack.PathInRMSDSpace")