r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np

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

8# Electrostatic constant

9k_c = units.Hartree * units.Bohr

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]

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)

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

40class ACN(Calculator):

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

42 nolabel = True

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

45 """Three-site potential for acetonitrile.

47 Atom sequence must be:

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

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.

57 rc: float

58 Cutoff radius for Coulomb interactions.

59 width: float

60 Width for cutoff function for Coulomb interactions.

62 References:

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

73 def calculate(self, atoms=None,

74 properties=['energy'],

75 system_changes=all_changes):

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

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'

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

97 pbc = self.atoms.pbc

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

99 nm = len(R)

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

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

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

106 # LJ parameters

107 sigma_co, epsilon_co = combine_lj_lorenz_berthelot(sigma, epsilon)

109 energy = 0.0

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

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

118 # Smooth cutoff

119 cut, dcut = self.cutoff(Dmm_min_len)

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)

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

153 self.results['energy'] = energy

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

156 def redistribute_forces(self, forces):

157 return forces

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

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

177 def embed(self, charges):

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

179 self.pcpot = PointChargePotential(charges)

180 return self.pcpot

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

189 return positions # no virtual sites

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

206class PointChargePotential:

207 def __init__(self, mmcharges):

208 """Point-charge potential for ACN.

210 Only used for testing QMMM.

211 """

212 self.mmcharges = mmcharges

213 self.mmpositions = None

214 self.mmforces = None

216 def set_positions(self, mmpositions):

217 self.mmpositions = mmpositions

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

234 def get_forces(self, calc):

235 return self.mmforces