1"""Effective medium theory potential.""" 

2 

3from math import sqrt, exp, log 

4 

5import numpy as np 

6 

7from ase.data import chemical_symbols, atomic_numbers 

8from ase.units import Bohr 

9from ase.neighborlist import NeighborList 

10from ase.calculators.calculator import (Calculator, all_changes, 

11 PropertyNotImplementedError) 

12 

13 

14parameters = { 

15 # E0 s0 V0 eta2 kappa lambda n0 

16 # eV bohr eV bohr^-1 bohr^-1 bohr^-1 bohr^-3 

17 'Al': (-3.28, 3.00, 1.493, 1.240, 2.000, 1.169, 0.00700), 

18 'Cu': (-3.51, 2.67, 2.476, 1.652, 2.740, 1.906, 0.00910), 

19 'Ag': (-2.96, 3.01, 2.132, 1.652, 2.790, 1.892, 0.00547), 

20 'Au': (-3.80, 3.00, 2.321, 1.674, 2.873, 2.182, 0.00703), 

21 'Ni': (-4.44, 2.60, 3.673, 1.669, 2.757, 1.948, 0.01030), 

22 'Pd': (-3.90, 2.87, 2.773, 1.818, 3.107, 2.155, 0.00688), 

23 'Pt': (-5.85, 2.90, 4.067, 1.812, 3.145, 2.192, 0.00802), 

24 # extra parameters - just for fun ... 

25 'H': (-3.21, 1.31, 0.132, 2.652, 2.790, 3.892, 0.00547), 

26 'C': (-3.50, 1.81, 0.332, 1.652, 2.790, 1.892, 0.01322), 

27 'N': (-5.10, 1.88, 0.132, 1.652, 2.790, 1.892, 0.01222), 

28 'O': (-4.60, 1.95, 0.332, 1.652, 2.790, 1.892, 0.00850)} 

29 

30beta = 1.809 # (16 * pi / 3)**(1.0 / 3) / 2**0.5, preserve historical rounding 

31 

32 

33class EMT(Calculator): 

34 """Python implementation of the Effective Medium Potential. 

35 

36 Supports the following standard EMT metals: 

37 Al, Cu, Ag, Au, Ni, Pd and Pt. 

38 

39 In addition, the following elements are supported. 

40 They are NOT well described by EMT, and the parameters 

41 are not for any serious use: 

42 H, C, N, O 

43 

44 The potential takes a single argument, ``asap_cutoff`` 

45 (default: False). If set to True, the cutoff mimics 

46 how Asap does it; most importantly the global cutoff 

47 is chosen from the largest atom present in the simulation, 

48 if False it is chosen from the largest atom in the parameter 

49 table. True gives the behaviour of the Asap code and 

50 older EMT implementations, although the results are not 

51 bitwise identical. 

52 """ 

53 implemented_properties = ['energy', 'energies', 'forces', 

54 'stress', 'magmom', 'magmoms'] 

55 

56 nolabel = True 

57 

58 default_parameters = {'asap_cutoff': False} 

59 

60 def __init__(self, **kwargs): 

61 Calculator.__init__(self, **kwargs) 

62 

63 def initialize(self, atoms): 

64 self.par = {} 

65 self.rc = 0.0 

66 self.numbers = atoms.get_atomic_numbers() 

67 if self.parameters.asap_cutoff: 

68 relevant_pars = {} 

69 for symb, p in parameters.items(): 

70 if atomic_numbers[symb] in self.numbers: 

71 relevant_pars[symb] = p 

72 else: 

73 relevant_pars = parameters 

74 maxseq = max(par[1] for par in relevant_pars.values()) * Bohr 

75 rc = self.rc = beta * maxseq * 0.5 * (sqrt(3) + sqrt(4)) 

76 rr = rc * 2 * sqrt(4) / (sqrt(3) + sqrt(4)) 

77 self.acut = np.log(9999.0) / (rr - rc) 

78 if self.parameters.asap_cutoff: 

79 self.rc_list = self.rc * 1.045 

80 else: 

81 self.rc_list = self.rc + 0.5 

82 for Z in self.numbers: 

83 if Z not in self.par: 

84 sym = chemical_symbols[Z] 

85 if sym not in parameters: 

86 raise NotImplementedError('No EMT-potential for {0}' 

87 .format(sym)) 

88 p = parameters[sym] 

89 s0 = p[1] * Bohr 

90 eta2 = p[3] / Bohr 

91 kappa = p[4] / Bohr 

92 x = eta2 * beta * s0 

93 gamma1 = 0.0 

94 gamma2 = 0.0 

95 for i, n in enumerate([12, 6, 24]): 

96 r = s0 * beta * sqrt(i + 1) 

97 x = n / (12 * (1.0 + exp(self.acut * (r - rc)))) 

98 gamma1 += x * exp(-eta2 * (r - beta * s0)) 

99 gamma2 += x * exp(-kappa / beta * (r - beta * s0)) 

100 

101 self.par[Z] = {'E0': p[0], 

102 's0': s0, 

103 'V0': p[2], 

104 'eta2': eta2, 

105 'kappa': kappa, 

106 'lambda': p[5] / Bohr, 

107 'n0': p[6] / Bohr**3, 

108 'rc': rc, 

109 'gamma1': gamma1, 

110 'gamma2': gamma2} 

111 

112 self.ksi = {} 

113 for s1, p1 in self.par.items(): 

114 self.ksi[s1] = {} 

115 for s2, p2 in self.par.items(): 

116 self.ksi[s1][s2] = p2['n0'] / p1['n0'] 

117 

118 self.energies = np.empty(len(atoms)) 

119 self.forces = np.empty((len(atoms), 3)) 

120 self.stress = np.empty((3, 3)) 

121 self.sigma1 = np.empty(len(atoms)) 

122 self.deds = np.empty(len(atoms)) 

123 

124 self.nl = NeighborList([0.5 * self.rc_list] * len(atoms), 

125 self_interaction=False) 

126 

127 def calculate(self, atoms=None, properties=['energy'], 

128 system_changes=all_changes): 

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

130 

131 if 'numbers' in system_changes: 

132 self.initialize(self.atoms) 

133 

134 positions = self.atoms.positions 

135 numbers = self.atoms.numbers 

136 cell = self.atoms.cell 

137 

138 self.nl.update(self.atoms) 

139 

140 self.energy = 0.0 

141 self.energies[:] = 0 

142 self.sigma1[:] = 0.0 

143 self.forces[:] = 0.0 

144 self.stress[:] = 0.0 

145 

146 natoms = len(self.atoms) 

147 

148 for a1 in range(natoms): 

149 Z1 = numbers[a1] 

150 p1 = self.par[Z1] 

151 ksi = self.ksi[Z1] 

152 neighbors, offsets = self.nl.get_neighbors(a1) 

153 offsets = np.dot(offsets, cell) 

154 for a2, offset in zip(neighbors, offsets): 

155 d = positions[a2] + offset - positions[a1] 

156 r = sqrt(np.dot(d, d)) 

157 if r < self.rc_list: 

158 Z2 = numbers[a2] 

159 p2 = self.par[Z2] 

160 self.interact1(a1, a2, d, r, p1, p2, ksi[Z2]) 

161 

162 for a in range(natoms): 

163 Z = numbers[a] 

164 p = self.par[Z] 

165 try: 

166 ds = -log(self.sigma1[a] / 12) / (beta * p['eta2']) 

167 except (OverflowError, ValueError): 

168 self.deds[a] = 0.0 

169 self.energy -= p['E0'] 

170 self.energies[a] -= p['E0'] 

171 continue 

172 x = p['lambda'] * ds 

173 y = exp(-x) 

174 z = 6 * p['V0'] * exp(-p['kappa'] * ds) 

175 self.deds[a] = ((x * y * p['E0'] * p['lambda'] + p['kappa'] * z) / 

176 (self.sigma1[a] * beta * p['eta2'])) 

177 E = p['E0'] * ((1 + x) * y - 1) + z 

178 self.energy += E 

179 self.energies[a] += E 

180 

181 for a1 in range(natoms): 

182 Z1 = numbers[a1] 

183 p1 = self.par[Z1] 

184 ksi = self.ksi[Z1] 

185 neighbors, offsets = self.nl.get_neighbors(a1) 

186 offsets = np.dot(offsets, cell) 

187 for a2, offset in zip(neighbors, offsets): 

188 d = positions[a2] + offset - positions[a1] 

189 r = sqrt(np.dot(d, d)) 

190 if r < self.rc_list: 

191 Z2 = numbers[a2] 

192 p2 = self.par[Z2] 

193 self.interact2(a1, a2, d, r, p1, p2, ksi[Z2]) 

194 

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

196 self.results['energies'] = self.energies 

197 self.results['free_energy'] = self.energy 

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

199 

200 if 'stress' in properties: 

201 if self.atoms.cell.rank == 3: 

202 self.stress += self.stress.T.copy() 

203 self.stress *= -0.5 / self.atoms.get_volume() 

204 self.results['stress'] = self.stress.flat[[0, 4, 8, 5, 2, 1]] 

205 else: 

206 raise PropertyNotImplementedError 

207 

208 def interact1(self, a1, a2, d, r, p1, p2, ksi): 

209 x = exp(self.acut * (r - self.rc)) 

210 theta = 1.0 / (1.0 + x) 

211 y1 = (0.5 * p1['V0'] * exp(-p2['kappa'] * (r / beta - p2['s0'])) * 

212 ksi / p1['gamma2'] * theta) 

213 y2 = (0.5 * p2['V0'] * exp(-p1['kappa'] * (r / beta - p1['s0'])) / 

214 ksi / p2['gamma2'] * theta) 

215 self.energy -= y1 + y2 

216 self.energies[a1] -= (y1 + y2) / 2 

217 self.energies[a2] -= (y1 + y2) / 2 

218 f = ((y1 * p2['kappa'] + y2 * p1['kappa']) / beta + 

219 (y1 + y2) * self.acut * theta * x) * d / r 

220 self.forces[a1] += f 

221 self.forces[a2] -= f 

222 self.stress -= np.outer(f, d) 

223 self.sigma1[a1] += (exp(-p2['eta2'] * (r - beta * p2['s0'])) * 

224 ksi * theta / p1['gamma1']) 

225 self.sigma1[a2] += (exp(-p1['eta2'] * (r - beta * p1['s0'])) / 

226 ksi * theta / p2['gamma1']) 

227 

228 def interact2(self, a1, a2, d, r, p1, p2, ksi): 

229 x = exp(self.acut * (r - self.rc)) 

230 theta = 1.0 / (1.0 + x) 

231 y1 = (exp(-p2['eta2'] * (r - beta * p2['s0'])) * 

232 ksi / p1['gamma1'] * theta * self.deds[a1]) 

233 y2 = (exp(-p1['eta2'] * (r - beta * p1['s0'])) / 

234 ksi / p2['gamma1'] * theta * self.deds[a2]) 

235 f = ((y1 * p2['eta2'] + y2 * p1['eta2']) + 

236 (y1 + y2) * self.acut * theta * x) * d / r 

237 self.forces[a1] -= f 

238 self.forces[a2] += f 

239 self.stress += np.outer(f, d)