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

1# fmt: off 

2 

3import numpy as np 

4 

5import ase.units as units 

6from ase.calculators.calculator import Calculator, all_changes 

7from ase.data import atomic_masses 

8from ase.geometry import find_mic 

9 

10# Electrostatic constant 

11k_c = units.Hartree * units.Bohr 

12 

13# Force field parameters 

14q_me = 0.206 

15q_c = 0.247 

16q_n = -0.453 

17sigma_me = 3.775 

18sigma_c = 3.650 

19sigma_n = 3.200 

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

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

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

23r_mec = 1.458 

24r_cn = 1.157 

25r_men = r_mec + r_cn 

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

27m_c = atomic_masses[6] 

28m_n = atomic_masses[7] 

29 

30 

31def combine_lj_lorenz_berthelot(sigma, epsilon): 

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

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

34 epsilon_c = np.zeros_like(sigma_c) 

35 

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

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

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

39 return sigma_c, epsilon_c 

40 

41 

42class ACN(Calculator): 

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

44 nolabel = True 

45 

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

47 """Three-site potential for acetonitrile. 

48 

49 Atom sequence must be: 

50 MeCNMeCN ... MeCN or NCMeNCMe ... NCMe 

51 

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

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

54 rigid triatomic molecules from Ciccotti et al. Molecular Physics 

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

56 using the FixLinearTriatomic to fix the geometry of the acetonitrile 

57 molecules. 

58 

59 rc: float 

60 Cutoff radius for Coulomb interactions. 

61 width: float 

62 Width for cutoff function for Coulomb interactions. 

63 

64 References 

65 ---------- 

66 

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

68 """ 

69 self.rc = rc 

70 self.width = width 

71 self.forces = None 

72 Calculator.__init__(self) 

73 self.sites_per_mol = 3 

74 self.pcpot = None 

75 

76 def calculate(self, atoms=None, 

77 properties=['energy'], 

78 system_changes=all_changes): 

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

80 

81 Z = atoms.numbers 

82 masses = atoms.get_masses() 

83 if Z[0] == 7: 

84 n = 0 

85 me = 2 

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

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

88 else: 

89 n = 2 

90 me = 0 

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

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

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

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

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

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

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

98 

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

100 pbc = self.atoms.pbc 

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

102 nm = len(R) 

103 

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

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

106 

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

108 

109 # LJ parameters 

110 sigma_co, epsilon_co = combine_lj_lorenz_berthelot(sigma, epsilon) 

111 

112 energy = 0.0 

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

114 

115 for m in range(nm - 1): 

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

117 # MIC PBCs 

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

119 shift = Dmm_min - Dmm 

120 

121 # Smooth cutoff 

122 cut, dcut = self.cutoff(Dmm_min_len) 

123 

124 for j in range(3): 

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

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

127 D_len = D_len2**0.5 

128 # Coulomb interactions 

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

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

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

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

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

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

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

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

137 # LJ interactions 

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

139 c12 = c6**2 

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

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

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

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

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

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

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

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

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

149 

150 if self.pcpot: 

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

152 self.atoms.positions) 

153 energy += e 

154 self.forces += f 

155 

156 self.results['energy'] = energy 

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

158 

159 def redistribute_forces(self, forces): 

160 return forces 

161 

162 def get_molcoms(self, nm): 

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

164 for m in range(nm): 

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

166 return molcoms 

167 

168 def cutoff(self, d): 

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

170 x2 = d < self.rc 

171 x12 = np.logical_and(x1, x2) 

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

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

174 cut[x2] = 1.0 

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

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

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

178 return cut, dtdd 

179 

180 def embed(self, charges): 

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

182 self.pcpot = PointChargePotential(charges) 

183 return self.pcpot 

184 

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

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

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

188 system_changes.append('positions') 

189 return system_changes 

190 

191 def add_virtual_sites(self, positions): 

192 return positions # no virtual sites 

193 

194 def get_virtual_charges(self, atoms): 

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

196 Z = atoms.numbers 

197 if Z[0] == 7: 

198 n = 0 

199 me = 2 

200 else: 

201 n = 2 

202 me = 0 

203 charges[me::3] = q_me 

204 charges[1::3] = q_c 

205 charges[n::3] = q_n 

206 return charges 

207 

208 

209class PointChargePotential: 

210 def __init__(self, mmcharges): 

211 """Point-charge potential for ACN. 

212 

213 Only used for testing QMMM. 

214 """ 

215 self.mmcharges = mmcharges 

216 self.mmpositions = None 

217 self.mmforces = None 

218 

219 def set_positions(self, mmpositions): 

220 self.mmpositions = mmpositions 

221 

222 def calculate(self, qmcharges, qmpositions): 

223 energy = 0.0 

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

225 qmforces = np.zeros_like(qmpositions) 

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

227 d = qmpositions - R 

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

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

230 energy += e.sum() 

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

232 qmforces += f 

233 F -= f.sum(0) 

234 self.mmpositions = None 

235 return energy, qmforces 

236 

237 def get_forces(self, calc): 

238 return self.mmforces