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 

10import typing as t 

11 

12from openmm import unit as mmunit 

13 

14from .base_path_cv import BasePathCV 

15from .path import Metric, progress 

16from .rmsd import RMSD 

17from .units.units import VectorQuantity 

18 

19 

20class PathInRMSDSpace(BasePathCV): 

21 r""" 

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

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

24 reference structures :cite:`Branduardi_2007`: 

25 

26 .. math:: 

27 

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

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

30 }{n-1} 

31 \quad \text{or} \quad 

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

33 

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

35 milestone, i.e., 

36 

37 .. math:: 

38 

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

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

41 \right) 

42 

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

44 distance between the current system state and the :math:`i`-th reference structure 

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

46 

47 .. note:: 

48 

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

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

51 

52 Parameters 

53 ---------- 

54 metric 

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

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

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

58 milestones 

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

60 atom indices to coordinate vectors. 

61 numAtoms 

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

63 any of the reference structures. 

64 sigma 

65 The width of the Gaussian kernels in nanometers 

66 name 

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

68 "path_progress_in_rmsd_space" or "path_deviation_in_rmsd_space" depending 

69 on the metric. 

70 

71 Raises 

72 ------ 

73 ValueError 

74 The number of milestones is less than 2 

75 ValueError 

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

77 

78 Examples 

79 -------- 

80 >>> import cvpack 

81 >>> import networkx as nx 

82 >>> import numpy as np 

83 >>> import openmm 

84 >>> from openmm import unit 

85 >>> from openmmtools import testsystems 

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

87 >>> model = testsystems.AlanineDipeptideVacuum() 

88 >>> atom1, atom2 = 8, 14 

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

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

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

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

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

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

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

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

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

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

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

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

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

102 ... frames.append(x.copy()) 

103 >>> milestones = [ 

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

105 ... for frame in frames 

106 ... ] 

107 >>> s, z = [ 

108 ... cvpack.PathInRMSDSpace( 

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

110 ... ) 

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

112 ... ] 

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

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

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

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

117 >>> s.getValue(context) 

118 0.172... dimensionless 

119 >>> z.getValue(context) 

120 -0.004... nm**2 

121 """ 

122 

123 def __init__( # pylint: disable=too-many-branches 

124 self, 

125 metric: Metric, 

126 milestones: t.Sequence[t.Dict[int, VectorQuantity]], 

127 numAtoms: int, 

128 sigma: mmunit.Quantity, 

129 name: t.Optional[str] = None, 

130 ) -> None: 

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

132 if mmunit.is_quantity(sigma): 

133 sigma = sigma.value_in_unit(mmunit.nanometers) 

134 n = len(milestones) 

135 if n < 2: 

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

137 collective_variables = { 

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

139 for i, reference in enumerate(milestones) 

140 } 

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

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

143 self._registerCV( 

144 name, 

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

146 metric=metric, 

147 milestones=[{k: list(v) for k, v in m.items()} for m in milestones], 

148 numAtoms=numAtoms, 

149 sigma=sigma, 

150 ) 

151 

152 

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