Coverage for cvpack/shortest_distance.py: 100%

1""" 

2.. class:: ShortestDistance 

3 :platform: Linux, MacOS, Windows 

4 :synopsis: The number of contacts between two atom groups 

5 

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

7 

8""" 

9 

10import typing as t 

11 

12import numpy as np 

13import openmm 

14from openmm import unit as mmunit 

15 

16from .collective_variable import CollectiveVariable 

17from .units import Quantity 

18 

19 

20class ShortestDistance(CollectiveVariable, openmm.CustomCVForce): 

21 r""" 

22 A smooth approximation for the shortest distance between two atom groups: 

23 

24 .. math:: 

25 r_{\rm min}({\bf r}) = \frac{S_{rw}({\bf r})}{S_w({\bf r})} 

26 

27 with 

28 

29 .. math:: 

30 

31 S_{rw}({\bf r}) = r_c e^{-\beta} + 

32 \sum_{i \in {\bf g}_1} \sum_{j \in {\bf g}_2 \atop r_{ij} < r_c} 

33 r_{ij} e^{-\beta \frac{r_{ij}}{r_c}} 

34 

35 and 

36 

37 .. math:: 

38 

39 S_w({\bf r}) = e^{-\beta} + 

40 \sum_{i \in {\bf g}_1} \sum_{j \in {\bf g}_2 \atop r_{ij} < r_c} 

41 r_{ij} e^{-\beta \frac{r_{ij}}{r_c}} 

42 

43 where :math:`r_{ij} = \|{\bf r}_j - {\bf r}_i\|` is the distance between atoms 

44 :math:`i` and :math:`j`, :math:`{\bf g}_1` and :math:`{\bf g}_2` are the sets of 

45 indices of the atoms in the first and second groups, respectively, :math:`r_c` is 

46 the cutoff distance, and :math:`\beta` is a parameter that controls the degree of 

47 approximation. 

48 

49 The larger the value of :math:`\beta`, the closer the approximation to the true 

50 shortest distance. However, an excessively large value may lead to numerical 

51 instability. 

52 

53 The terms outside the summations guarantee that shortest distance is well-defined 

54 even when there are no atom pairs within the cutoff distance. In this case, the 

55 collective variable evaluates to the cutoff distance. 

56 

57 .. note:: 

58 

59 Atoms are allowed to be in both groups. In this case, terms for which 

60 :math:`i = j` are ignored. 

61 

62 Parameters 

63 ---------- 

64 group1 

65 The indices of the atoms in the first group. 

66 group2 

67 The indices of the atoms in the second group. 

68 numAtoms 

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

70 of the groups. This is necessary for the correct initialization of the 

71 collective variable. 

72 beta 

73 The parameter that controls the degree of approximation. 

74 cutoffDistance 

75 The cutoff distance for evaluating atom pairs. 

76 pbc 

77 Whether to consider periodic boundary conditions in distance calculations. 

78 name 

79 The name of the collective variable. 

80 

81 Example 

82 ------- 

83 >>> import cvpack 

84 >>> import numpy as np 

85 >>> from openmm import unit 

86 >>> from openmmtools import testsystems 

87 >>> model = testsystems.HostGuestVacuum() 

88 >>> group1, group2 = [], [] 

89 >>> for residue in model.topology.residues(): 

90 ... group = group1 if residue.name == "B2" else group2 

91 ... group.extend(atom.index for atom in residue.atoms()) 

92 >>> coords = model.positions.value_in_unit(mmunit.nanometers) 

93 >>> np.linalg.norm( 

94 ... coords[group1, None, :] - coords[None, group2, :], 

95 ... axis=-1, 

96 ... ).min() 

97 0.1757... 

98 >>> num_atoms = model.system.getNumParticles() 

99 >>> min_dist = cvpack.ShortestDistance(group1, group2, num_atoms) 

100 >>> min_dist.addToSystem(model.system) 

101 >>> platform = openmm.Platform.getPlatformByName("Reference") 

102 >>> integrator = openmm.VerletIntegrator(1.0 * mmunit.femtoseconds) 

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

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

105 >>> min_dist.getValue(context) 

106 0.1785... nm 

107 """ 

108 

109 def __init__( 

110 self, 

111 group1: t.Sequence[int], 

112 group2: t.Sequence[int], 

113 numAtoms: int, 

114 beta: float = 100, 

115 cutoffDistance: mmunit.Quantity = Quantity(1 * mmunit.nanometers), 

116 pbc: bool = True, 

117 name: str = "shortest_distance", 

118 ) -> None: 

119 rc = cutoffDistance 

120 if mmunit.is_quantity(rc): 

121 rc = rc.value_in_unit(mmunit.nanometers) 

122 weight = f"exp(-{beta/rc}*r)" 

123 forces = { 

124 "wrsum": openmm.CustomNonbondedForce(f"{weight}*r"), 

125 "wsum": openmm.CustomNonbondedForce(weight), 

126 } 

127 super().__init__(f"({rc*np.exp(-beta)}+wrsum)/({np.exp(-beta)}+wsum)") 

128 for cv, force in forces.items(): 

129 force.setNonbondedMethod( 

130 force.CutoffPeriodic if pbc else force.CutoffNonPeriodic 

131 ) 

132 force.setCutoffDistance(cutoffDistance) 

133 force.setUseSwitchingFunction(False) 

134 force.setUseLongRangeCorrection(False) 

135 for _ in range(numAtoms): 

136 force.addParticle([]) 

137 force.addInteractionGroup(group1, group2) 

138 self.addCollectiveVariable(cv, force) 

139 self._registerCV( 

140 name, 

141 mmunit.nanometers, 

142 group1=group1, 

143 group2=group2, 

144 numAtoms=numAtoms, 

145 beta=beta, 

146 cutoffDistance=rc, 

147 pbc=pbc, 

148 ) 

149 

150 

151ShortestDistance.registerTag("!cvpack.ShortestDistance")