Coverage for /builds/kinetik161/ase/ase/calculators/acn.py: 93.83%

1import numpy as np 

2 

3import ase.units as units 

4from ase.calculators.calculator import Calculator, all_changes 

5from ase.data import atomic_masses 

6from ase.geometry import find_mic 

7 

8# Electrostatic constant 

9k_c = units.Hartree * units.Bohr 

10 

11# Force field parameters 

12q_me = 0.206 

13q_c = 0.247 

14q_n = -0.453 

15sigma_me = 3.775 

16sigma_c = 3.650 

17sigma_n = 3.200 

18epsilon_me = 0.7824 * units.kJ / units.mol 

19epsilon_c = 0.544 * units.kJ / units.mol 

20epsilon_n = 0.6276 * units.kJ / units.mol 

21r_mec = 1.458 

22r_cn = 1.157 

23r_men = r_mec + r_cn 

24m_me = atomic_masses[6] + 3 * atomic_masses[1] 

25m_c = atomic_masses[6] 

26m_n = atomic_masses[7] 

27 

28 

29def combine_lj_lorenz_berthelot(sigma, epsilon): 

30 """Combine LJ parameters according to the Lorenz-Berthelot rule""" 

31 sigma_c = np.zeros((len(sigma), len(sigma))) 

32 epsilon_c = np.zeros_like(sigma_c) 

33 

34 for ii in range(len(sigma)): 

35 sigma_c[:, ii] = (sigma[ii] + sigma) / 2 

36 epsilon_c[:, ii] = (epsilon[ii] * epsilon) ** 0.5 

37 return sigma_c, epsilon_c 

38 

39 

40class ACN(Calculator): 

41 implemented_properties = ['energy', 'forces'] 

42 nolabel = True 

43 

44 def __init__(self, rc=5.0, width=1.0): 

45 """Three-site potential for acetonitrile. 

46 

47 Atom sequence must be: 

48 MeCNMeCN ... MeCN or NCMeNCMe ... NCMe 

49 

50 When performing molecular dynamics (MD), forces are redistributed 

51 and only Me and N sites propagated based on a scheme for MD of 

52 rigid triatomic molecules from Ciccotti et al. Molecular Physics 

53 1982 (https://doi.org/10.1080/00268978200100942). Apply constraints 

54 using the FixLinearTriatomic to fix the geometry of the acetonitrile 

55 molecules. 

56 

57 rc: float 

58 Cutoff radius for Coulomb interactions. 

59 width: float 

60 Width for cutoff function for Coulomb interactions. 

61 

62 References: 

63 

64 https://doi.org/10.1080/08927020108024509 

65 """ 

66 self.rc = rc 

67 self.width = width 

68 self.forces = None 

69 Calculator.__init__(self) 

70 self.sites_per_mol = 3 

71 self.pcpot = None 

72 

73 def calculate(self, atoms=None, 

74 properties=['energy'], 

75 system_changes=all_changes): 

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

77 

78 Z = atoms.numbers 

79 masses = atoms.get_masses() 

80 if Z[0] == 7: 

81 n = 0 

82 me = 2 

83 sigma = np.array([sigma_n, sigma_c, sigma_me]) 

84 epsilon = np.array([epsilon_n, epsilon_c, epsilon_me]) 

85 else: 

86 n = 2 

87 me = 0 

88 sigma = np.array([sigma_me, sigma_c, sigma_n]) 

89 epsilon = np.array([epsilon_me, epsilon_c, epsilon_n]) 

90 assert (Z[n::3] == 7).all(), 'incorrect atoms sequence' 

91 assert (Z[1::3] == 6).all(), 'incorrect atoms sequence' 

92 assert (masses[n::3] == m_n).all(), 'incorrect masses' 

93 assert (masses[1::3] == m_c).all(), 'incorrect masses' 

94 assert (masses[me::3] == m_me).all(), 'incorrect masses' 

95 

96 R = self.atoms.positions.reshape((-1, 3, 3)) 

97 pbc = self.atoms.pbc 

98 cd = self.atoms.cell.diagonal() 

99 nm = len(R) 

100 

101 assert (self.atoms.cell == np.diag(cd)).all(), 'not orthorhombic' 

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

103 

104 charges = self.get_virtual_charges(atoms[:3]) 

105 

106 # LJ parameters 

107 sigma_co, epsilon_co = combine_lj_lorenz_berthelot(sigma, epsilon) 

108 

109 energy = 0.0 

110 self.forces = np.zeros((3 * nm, 3)) 

111 

112 for m in range(nm - 1): 

113 Dmm = R[m + 1:, 1] - R[m, 1] 

114 # MIC PBCs 

115 Dmm_min, Dmm_min_len = find_mic(Dmm, atoms.cell, pbc) 

116 shift = Dmm_min - Dmm 

117 

118 # Smooth cutoff 

119 cut, dcut = self.cutoff(Dmm_min_len) 

120 

121 for j in range(3): 

122 D = R[m + 1:] - R[m, j] + shift[:, np.newaxis] 

123 D_len2 = (D**2).sum(axis=2) 

124 D_len = D_len2**0.5 

125 # Coulomb interactions 

126 e = charges[j] * charges / D_len * k_c 

127 energy += np.dot(cut, e).sum() 

128 F = (e / D_len2 * cut[:, np.newaxis])[:, :, np.newaxis] * D 

129 Fmm = -(e.sum(1) * dcut / Dmm_min_len)[:, np.newaxis] * Dmm_min 

130 self.forces[(m + 1) * 3:] += F.reshape((-1, 3)) 

131 self.forces[m * 3 + j] -= F.sum(axis=0).sum(axis=0) 

132 self.forces[(m + 1) * 3 + 1::3] += Fmm 

133 self.forces[m * 3 + 1] -= Fmm.sum(0) 

134 # LJ interactions 

135 c6 = (sigma_co[:, j]**2 / D_len2)**3 

136 c12 = c6**2 

137 e = 4 * epsilon_co[:, j] * (c12 - c6) 

138 energy += np.dot(cut, e).sum() 

139 F = (24 * epsilon_co[:, j] * (2 * c12 - 

140 c6) / D_len2 * cut[:, np.newaxis])[:, :, np.newaxis] * D 

141 Fmm = -(e.sum(1) * dcut / Dmm_min_len)[:, np.newaxis] * Dmm_min 

142 self.forces[(m + 1) * 3:] += F.reshape((-1, 3)) 

143 self.forces[m * 3 + j] -= F.sum(axis=0).sum(axis=0) 

144 self.forces[(m + 1) * 3 + 1::3] += Fmm 

145 self.forces[m * 3 + 1] -= Fmm.sum(0) 

146 

147 if self.pcpot: 

148 e, f = self.pcpot.calculate(np.tile(charges, nm), 

149 self.atoms.positions) 

150 energy += e 

151 self.forces += f 

152 

153 self.results['energy'] = energy 

154 self.results['forces'] = self.forces 

155 

156 def redistribute_forces(self, forces): 

157 return forces 

158 

159 def get_molcoms(self, nm): 

160 molcoms = np.zeros((nm, 3)) 

161 for m in range(nm): 

162 molcoms[m] = self.atoms[m * 3:(m + 1) * 3].get_center_of_mass() 

163 return molcoms 

164 

165 def cutoff(self, d): 

166 x1 = d > self.rc - self.width 

167 x2 = d < self.rc 

168 x12 = np.logical_and(x1, x2) 

169 y = (d[x12] - self.rc + self.width) / self.width 

170 cut = np.zeros(len(d)) # cutoff function 

171 cut[x2] = 1.0 

172 cut[x12] -= y**2 * (3.0 - 2.0 * y) 

173 dtdd = np.zeros(len(d)) 

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

175 return cut, dtdd 

176 

177 def embed(self, charges): 

178 """Embed atoms in point-charges.""" 

179 self.pcpot = PointChargePotential(charges) 

180 return self.pcpot 

181 

182 def check_state(self, atoms, tol=1e-15): 

183 system_changes = Calculator.check_state(self, atoms, tol) 

184 if self.pcpot and self.pcpot.mmpositions is not None: 

185 system_changes.append('positions') 

186 return system_changes 

187 

188 def add_virtual_sites(self, positions): 

189 return positions # no virtual sites 

190 

191 def get_virtual_charges(self, atoms): 

192 charges = np.empty(len(atoms)) 

193 Z = atoms.numbers 

194 if Z[0] == 7: 

195 n = 0 

196 me = 2 

197 else: 

198 n = 2 

199 me = 0 

200 charges[me::3] = q_me 

201 charges[1::3] = q_c 

202 charges[n::3] = q_n 

203 return charges 

204 

205 

206class PointChargePotential: 

207 def __init__(self, mmcharges): 

208 """Point-charge potential for ACN. 

209 

210 Only used for testing QMMM. 

211 """ 

212 self.mmcharges = mmcharges 

213 self.mmpositions = None 

214 self.mmforces = None 

215 

216 def set_positions(self, mmpositions): 

217 self.mmpositions = mmpositions 

218 

219 def calculate(self, qmcharges, qmpositions): 

220 energy = 0.0 

221 self.mmforces = np.zeros_like(self.mmpositions) 

222 qmforces = np.zeros_like(qmpositions) 

223 for C, R, F in zip(self.mmcharges, self.mmpositions, self.mmforces): 

224 d = qmpositions - R 

225 r2 = (d**2).sum(1) 

226 e = units.Hartree * units.Bohr * C * r2**-0.5 * qmcharges 

227 energy += e.sum() 

228 f = (e / r2)[:, np.newaxis] * d 

229 qmforces += f 

230 F -= f.sum(0) 

231 self.mmpositions = None 

232 return energy, qmforces 

233 

234 def get_forces(self, calc): 

235 return self.mmforces