1import numpy as np 

2 

3from ase import units 

4from ase.calculators.calculator import Calculator, all_changes 

5from ase.calculators.tip3p import rOH, angleHOH, TIP3P 

6 

7__all__ = ['rOH', 'angleHOH', 'TIP4P', 'sigma0', 'epsilon0'] 

8 

9# Electrostatic constant and parameters: 

10k_c = units.Hartree * units.Bohr 

11qH = 0.52 

12A = 600e3 * units.kcal / units.mol 

13B = 610 * units.kcal / units.mol 

14sigma0 = (A / B)**(1 / 6.) 

15epsilon0 = B**2 / (4 * A) 

16# https://doi.org/10.1063/1.445869 

17 

18 

19class TIP4P(TIP3P): 

20 def __init__(self, rc=7.0, width=1.0): 

21 """ TIP4P potential for water. 

22 

23 https://doi.org/10.1063/1.445869 

24 

25 Requires an atoms object of OHH,OHH, ... sequence 

26 Correct TIP4P charges and LJ parameters set automatically. 

27 

28 Virtual interaction sites implemented in the following scheme: 

29 Original atoms object has no virtual sites. 

30 When energy/forces are requested: 

31 

32 * virtual sites added to temporary xatoms object 

33 * energy / forces calculated 

34 * forces redistributed from virtual sites to actual atoms object 

35 

36 This means you do not get into trouble when propagating your system 

37 with MD while having to skip / account for massless virtual sites. 

38 

39 This also means that if using for QM/MM MD with GPAW, the EmbedTIP4P 

40 class must be used. 

41 """ 

42 

43 TIP3P.__init__(self, rc, width) 

44 self.atoms_per_mol = 3 

45 self.sites_per_mol = 4 

46 self.energy = None 

47 self.forces = None 

48 

49 def calculate(self, atoms=None, 

50 properties=['energy', 'forces'], 

51 system_changes=all_changes): 

52 Calculator.calculate(self, atoms, properties, system_changes) 

53 

54 assert (atoms.numbers[::3] == 8).all() 

55 assert (atoms.numbers[1::3] == 1).all() 

56 assert (atoms.numbers[2::3] == 1).all() 

57 

58 xpos = self.add_virtual_sites(atoms.positions) 

59 xcharges = self.get_virtual_charges(atoms) 

60 

61 cell = atoms.cell 

62 pbc = atoms.pbc 

63 

64 natoms = len(atoms) 

65 nmol = natoms // 3 

66 

67 self.energy = 0.0 

68 self.forces = np.zeros((4 * natoms // 3, 3)) 

69 

70 C = cell.diagonal() 

71 assert (cell == np.diag(C)).all(), 'not orthorhombic' 

72 assert ((C >= 2 * self.rc) | ~pbc).all(), 'cutoff too large' 

73 

74 # Get dx,dy,dz from first atom of each mol to same atom of all other 

75 # and find min. distance. Everything moves according to this analysis. 

76 for a in range(nmol - 1): 

77 D = xpos[(a + 1) * 4::4] - xpos[a * 4] 

78 shift = np.zeros_like(D) 

79 for i, periodic in enumerate(pbc): 

80 if periodic: 

81 shift[:, i] = np.rint(D[:, i] / C[i]) * C[i] 

82 q_v = xcharges[(a + 1) * 4:] 

83 

84 # Min. img. position list as seen for molecule !a! 

85 position_list = np.zeros(((nmol - 1 - a) * 4, 3)) 

86 

87 for j in range(4): 

88 position_list[j::4] += xpos[(a + 1) * 4 + j::4] - shift 

89 

90 # Make the smooth cutoff: 

91 pbcRoo = position_list[::4] - xpos[a * 4] 

92 pbcDoo = np.sum(np.abs(pbcRoo)**2, axis=-1)**(1 / 2) 

93 x1 = pbcDoo > self.rc - self.width 

94 x2 = pbcDoo < self.rc 

95 x12 = np.logical_and(x1, x2) 

96 y = (pbcDoo[x12] - self.rc + self.width) / self.width 

97 t = np.zeros(len(pbcDoo)) 

98 t[x2] = 1.0 

99 t[x12] -= y**2 * (3.0 - 2.0 * y) 

100 dtdd = np.zeros(len(pbcDoo)) 

101 dtdd[x12] -= 6.0 / self.width * y * (1.0 - y) 

102 self.energy_and_forces(a, xpos, position_list, q_v, nmol, t, dtdd) 

103 

104 if self.pcpot: 

105 e, f = self.pcpot.calculate(xcharges, xpos) 

106 self.energy += e 

107 self.forces += f 

108 

109 f = self.redistribute_forces(self.forces) 

110 

111 self.results['energy'] = self.energy 

112 self.results['forces'] = f 

113 

114 def energy_and_forces(self, a, xpos, position_list, q_v, nmol, t, dtdd): 

115 """ energy and forces on molecule a from all other molecules. 

116 cutoff is based on O-O Distance. """ 

117 

118 # LJ part - only O-O interactions 

119 epsil = np.tile([epsilon0], nmol - 1 - a) 

120 sigma = np.tile([sigma0], nmol - 1 - a) 

121 DOO = position_list[::4] - xpos[a * 4] 

122 d2 = (DOO**2).sum(1) 

123 d = np.sqrt(d2) 

124 e_lj = 4 * epsil * (sigma**12 / d**12 - sigma**6 / d**6) 

125 f_lj = (4 * epsil * (12 * sigma**12 / d**13 - 

126 6 * sigma**6 / d**7) * t - 

127 e_lj * dtdd)[:, np.newaxis] * DOO / d[:, np.newaxis] 

128 

129 self.forces[a * 4] -= f_lj.sum(0) 

130 self.forces[(a + 1) * 4::4] += f_lj 

131 

132 # Electrostatics 

133 e_elec = 0 

134 all_cut = np.repeat(t, 4) 

135 for i in range(4): 

136 D = position_list - xpos[a * 4 + i] 

137 d2_all = (D**2).sum(axis=1) 

138 d_all = np.sqrt(d2_all) 

139 e = k_c * q_v[i] * q_v / d_all 

140 e_elec += np.dot(all_cut, e).sum() 

141 e_f = e.reshape(nmol - a - 1, 4).sum(1) 

142 F = (e / d_all * all_cut)[:, np.newaxis] * D / d_all[:, np.newaxis] 

143 FOO = -(e_f * dtdd)[:, np.newaxis] * DOO / d[:, np.newaxis] 

144 self.forces[(a + 1) * 4 + 0::4] += FOO 

145 self.forces[a * 4] -= FOO.sum(0) 

146 self.forces[(a + 1) * 4:] += F 

147 self.forces[a * 4 + i] -= F.sum(0) 

148 

149 self.energy += np.dot(e_lj, t) + e_elec 

150 

151 def add_virtual_sites(self, pos): 

152 # Order: OHHM,OHHM,... 

153 # DOI: 10.1002/(SICI)1096-987X(199906)20:8 

154 b = 0.15 

155 xatomspos = np.zeros((4 * len(pos) // 3, 3)) 

156 for w in range(0, len(pos), 3): 

157 r_i = pos[w] # O pos 

158 r_j = pos[w + 1] # H1 pos 

159 r_k = pos[w + 2] # H2 pos 

160 n = (r_j + r_k) / 2 - r_i 

161 n /= np.linalg.norm(n) 

162 r_d = r_i + b * n 

163 

164 x = 4 * w // 3 

165 xatomspos[x + 0] = r_i 

166 xatomspos[x + 1] = r_j 

167 xatomspos[x + 2] = r_k 

168 xatomspos[x + 3] = r_d 

169 

170 return xatomspos 

171 

172 def get_virtual_charges(self, atoms): 

173 charges = np.empty(len(atoms) * 4 // 3) 

174 charges[0::4] = 0.00 # O 

175 charges[1::4] = qH # H1 

176 charges[2::4] = qH # H2 

177 charges[3::4] = - 2 * qH # X1 

178 return charges 

179 

180 def redistribute_forces(self, forces): 

181 f = forces 

182 b = 0.15 

183 a = 0.5 

184 pos = self.atoms.positions 

185 for w in range(0, len(pos), 3): 

186 r_i = pos[w] # O pos 

187 r_j = pos[w + 1] # H1 pos 

188 r_k = pos[w + 2] # H2 pos 

189 r_ij = r_j - r_i 

190 r_jk = r_k - r_j 

191 r_d = r_i + b * (r_ij + a * r_jk) / np.linalg.norm(r_ij + a * r_jk) 

192 r_id = r_d - r_i 

193 gamma = b / np.linalg.norm(r_ij + a * r_jk) 

194 

195 x = w * 4 // 3 

196 Fd = f[x + 3] # force on M 

197 F1 = (np.dot(r_id, Fd) / np.dot(r_id, r_id)) * r_id 

198 Fi = Fd - gamma * (Fd - F1) # Force from M on O 

199 Fj = (1 - a) * gamma * (Fd - F1) # Force from M on H1 

200 Fk = a * gamma * (Fd - F1) # Force from M on H2 

201 

202 f[x] += Fi 

203 f[x + 1] += Fj 

204 f[x + 2] += Fk 

205 

206 # remove virtual sites from force array 

207 f = np.delete(f, list(range(3, f.shape[0], 4)), axis=0) 

208 return f