Coverage for /builds/ase/ase/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 https://doi.org/10.1080/08927020108024509 

67 """ 

68 self.rc = rc 

69 self.width = width 

70 self.forces = None 

71 Calculator.__init__(self) 

72 self.sites_per_mol = 3 

73 self.pcpot = None 

74 

75 def calculate(self, atoms=None, 

76 properties=['energy'], 

77 system_changes=all_changes): 

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

79 

80 Z = atoms.numbers 

81 masses = atoms.get_masses() 

82 if Z[0] == 7: 

83 n = 0 

84 me = 2 

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

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

87 else: 

88 n = 2 

89 me = 0 

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

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

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

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

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

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

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

97 

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

99 pbc = self.atoms.pbc 

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

101 nm = len(R) 

102 

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

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

105 

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

107 

108 # LJ parameters 

109 sigma_co, epsilon_co = combine_lj_lorenz_berthelot(sigma, epsilon) 

110 

111 energy = 0.0 

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

113 

114 for m in range(nm - 1): 

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

116 # MIC PBCs 

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

118 shift = Dmm_min - Dmm 

119 

120 # Smooth cutoff 

121 cut, dcut = self.cutoff(Dmm_min_len) 

122 

123 for j in range(3): 

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

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

126 D_len = D_len2**0.5 

127 # Coulomb interactions 

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

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

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

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

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

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

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

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

136 # LJ interactions 

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

138 c12 = c6**2 

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

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

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

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

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

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

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

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

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

148 

149 if self.pcpot: 

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

151 self.atoms.positions) 

152 energy += e 

153 self.forces += f 

154 

155 self.results['energy'] = energy 

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

157 

158 def redistribute_forces(self, forces): 

159 return forces 

160 

161 def get_molcoms(self, nm): 

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

163 for m in range(nm): 

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

165 return molcoms 

166 

167 def cutoff(self, d): 

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

169 x2 = d < self.rc 

170 x12 = np.logical_and(x1, x2) 

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

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

173 cut[x2] = 1.0 

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

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

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

177 return cut, dtdd 

178 

179 def embed(self, charges): 

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

181 self.pcpot = PointChargePotential(charges) 

182 return self.pcpot 

183 

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

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

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

187 system_changes.append('positions') 

188 return system_changes 

189 

190 def add_virtual_sites(self, positions): 

191 return positions # no virtual sites 

192 

193 def get_virtual_charges(self, atoms): 

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

195 Z = atoms.numbers 

196 if Z[0] == 7: 

197 n = 0 

198 me = 2 

199 else: 

200 n = 2 

201 me = 0 

202 charges[me::3] = q_me 

203 charges[1::3] = q_c 

204 charges[n::3] = q_n 

205 return charges 

206 

207 

208class PointChargePotential: 

209 def __init__(self, mmcharges): 

210 """Point-charge potential for ACN. 

211 

212 Only used for testing QMMM. 

213 """ 

214 self.mmcharges = mmcharges 

215 self.mmpositions = None 

216 self.mmforces = None 

217 

218 def set_positions(self, mmpositions): 

219 self.mmpositions = mmpositions 

220 

221 def calculate(self, qmcharges, qmpositions): 

222 energy = 0.0 

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

224 qmforces = np.zeros_like(qmpositions) 

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

226 d = qmpositions - R 

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

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

229 energy += e.sum() 

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

231 qmforces += f 

232 F -= f.sum(0) 

233 self.mmpositions = None 

234 return energy, qmforces 

235 

236 def get_forces(self, calc): 

237 return self.mmforces