Coverage for /builds/ase/ase/ase/calculators/emt.py: 98.90%
182 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""Effective medium theory potential."""
4from collections import defaultdict
5from math import log, sqrt
7import numpy as np
9from ase.calculators.calculator import (
10 Calculator,
11 PropertyNotImplementedError,
12 all_changes,
13)
14from ase.data import atomic_numbers, chemical_symbols
15from ase.neighborlist import NeighborList
16from ase.units import Bohr
18parameters = {
19 # E0 s0 V0 eta2 kappa lambda n0
20 # eV bohr eV bohr^-1 bohr^-1 bohr^-1 bohr^-3
21 'Al': (-3.28, 3.00, 1.493, 1.240, 2.000, 1.169, 0.00700),
22 'Cu': (-3.51, 2.67, 2.476, 1.652, 2.740, 1.906, 0.00910),
23 'Ag': (-2.96, 3.01, 2.132, 1.652, 2.790, 1.892, 0.00547),
24 'Au': (-3.80, 3.00, 2.321, 1.674, 2.873, 2.182, 0.00703),
25 'Ni': (-4.44, 2.60, 3.673, 1.669, 2.757, 1.948, 0.01030),
26 'Pd': (-3.90, 2.87, 2.773, 1.818, 3.107, 2.155, 0.00688),
27 'Pt': (-5.85, 2.90, 4.067, 1.812, 3.145, 2.192, 0.00802),
28 # extra parameters - just for fun ...
29 'H': (-3.21, 1.31, 0.132, 2.652, 2.790, 3.892, 0.00547),
30 'C': (-3.50, 1.81, 0.332, 1.652, 2.790, 1.892, 0.01322),
31 'N': (-5.10, 1.88, 0.132, 1.652, 2.790, 1.892, 0.01222),
32 'O': (-4.60, 1.95, 0.332, 1.652, 2.790, 1.892, 0.00850)}
34beta = 1.809 # (16 * pi / 3)**(1.0 / 3) / 2**0.5, preserve historical rounding
37class EMT(Calculator):
38 """Python implementation of the Effective Medium Potential.
40 Supports the following standard EMT metals:
41 Al, Cu, Ag, Au, Ni, Pd and Pt.
43 In addition, the following elements are supported.
44 They are NOT well described by EMT, and the parameters
45 are not for any serious use:
46 H, C, N, O
48 Parameters
49 ----------
50 asap_cutoff : bool, default: False
51 If True the cutoff mimics how ASAP does it; most importantly the global
52 cutoff is chosen from the largest atom present in the simulation.
53 If False it is chosen from the largest atom in the parameter table.
54 True gives the behaviour of ASAP and older EMT implementations,
55 although the results are not bitwise identical.
57 Notes
58 -----
59 Formulation mostly follows Jacobsen *et al*. [1]_
60 `Documentation in ASAP can also be referred to <https://gitlab.com/asap/
61 asap/blob/master/docs/manual/potentials/emt.pdf>`_.
63 .. [1] K. W. Jacobsen, P. Stoltze, and J. K. Nørskov,
64 Surf. Sci. 366, 394 (1996).
65 """
66 implemented_properties = ['energy', 'free_energy', 'energies', 'forces',
67 'stress', 'magmom', 'magmoms']
69 nolabel = True
71 default_parameters = {'asap_cutoff': False}
73 def __init__(self, **kwargs):
74 Calculator.__init__(self, **kwargs)
76 def initialize(self, atoms):
77 self.rc, self.rc_list, self.acut = self._calc_cutoff(atoms)
79 numbers = atoms.get_atomic_numbers()
81 # ia2iz : map from idx of atoms to idx of atomic numbers in self.par
82 unique_numbers, self.ia2iz = np.unique(numbers, return_inverse=True)
83 self.par = defaultdict(lambda: np.empty(len(unique_numbers)))
84 for i, Z in enumerate(unique_numbers):
85 sym = chemical_symbols[Z]
86 if sym not in parameters:
87 raise NotImplementedError(f'No EMT-potential for {sym}')
88 p = parameters[sym]
89 s0 = p[1] * Bohr
90 eta2 = p[3] / Bohr
91 kappa = p[4] / Bohr
92 gamma1, gamma2 = self._calc_gammas(s0, eta2, kappa)
93 self.par['Z'][i] = Z
94 self.par['E0'][i] = p[0]
95 self.par['s0'][i] = s0
96 self.par['V0'][i] = p[2]
97 self.par['eta2'][i] = eta2
98 self.par['kappa'][i] = kappa
99 self.par['lambda'][i] = p[5] / Bohr
100 self.par['n0'][i] = p[6] / Bohr**3
101 self.par['inv12gamma1'][i] = 1.0 / (12.0 * gamma1)
102 self.par['neghalfv0overgamma2'][i] = -0.5 * p[2] / gamma2
104 self.chi = self.par['n0'][None, :] / self.par['n0'][:, None]
106 self.energies = np.empty(len(atoms))
107 self.forces = np.empty((len(atoms), 3))
108 self.stress = np.empty((3, 3))
109 self.deds = np.empty(len(atoms))
111 self.nl = NeighborList([0.5 * self.rc_list] * len(atoms),
112 self_interaction=False, bothways=True)
114 def _calc_cutoff(self, atoms):
115 """Calculate parameters of the logistic smoothing function etc.
117 The logistic smoothing function is given by
119 .. math:
121 w(r) = \\frac{1}{1 + \\exp a (r - r_\\mathrm{c})}
123 Returns
124 -------
125 rc : float
126 "Midpoint" of the logistic smoothing function, set to be the mean
127 of the 3rd and the 4th nearest-neighbor distances in FCC.
128 rc_list : float
129 Cutoff radius for the neighbor search, set to be slightly larger
130 than ``rc`` depending on ``asap_cutoff``.
131 acut : float
132 "Slope" of the smoothing function, set for the smoothing function
133 value to be ``1e-4`` at the 4th nearest-neighbor distance in FCC.
135 Notes
136 -----
137 ``maxseq`` is the present FCC Wigner-Seitz radius. ``beta * maxseq``
138 (`r1nn`) is the corresponding 1st nearest-neighbor distance in FCC.
139 The 2nd, 3rd, 4th nearest-neighbor distances in FCC are given using
140 ``r1nn`` by ``sqrt(2) * r1nn``, ``sqrt(3) * r1nn``, ``sqrt(4) * r1nn``,
141 respectively.
142 """
143 numbers = atoms.get_atomic_numbers()
144 if self.parameters['asap_cutoff']:
145 relevant_pars = {
146 symb: p
147 for symb, p in parameters.items()
148 if atomic_numbers[symb] in numbers
149 }
150 else:
151 relevant_pars = parameters
152 maxseq = max(par[1] for par in relevant_pars.values()) * Bohr
153 r1nn = beta * maxseq # 1st NN distance in FCC
154 rc = r1nn * 0.5 * (sqrt(3.0) + 2.0) # mean of 3NN and 4NN dists.
155 r4nn = r1nn * 2.0 # 4NN distance in FCC
156 eps = 1e-4 # value at r4nn, should be small
158 # "slope" is set so that the function value becomes eps at r4nn
159 acut = log(1.0 / eps - 1.0) / (r4nn - rc)
161 rc_list = rc * 1.045 if self.parameters['asap_cutoff'] else rc + 0.5
163 return rc, rc_list, acut
165 def _calc_gammas(self, s0, eta2, kappa):
166 n = np.array([12, 6, 24]) # numbers of 1, 2, 3NN sites in fcc
167 r = beta * s0 * np.sqrt([1.0, 2.0, 3.0]) # distances of 1, 2, 3NNs
168 w = 1.0 / (1.0 + np.exp(self.acut * (r - self.rc)))
169 x = n * w / 12.0
170 gamma1 = x @ np.exp(-eta2 * (r - beta * s0))
171 gamma2 = x @ np.exp(-kappa / beta * (r - beta * s0))
172 return gamma1, gamma2
174 def calculate(self, atoms=None, properties=['energy'],
175 system_changes=all_changes):
176 Calculator.calculate(self, atoms, properties, system_changes)
178 if 'numbers' in system_changes:
179 self.initialize(self.atoms)
181 self.nl.update(self.atoms)
183 self.energies[:] = 0.0
184 self.forces[:] = 0.0
185 self.stress[:] = 0.0
186 self.deds[:] = 0.0
188 natoms = len(self.atoms)
190 # store nearest neighbor info for all the atoms
191 # suffixes 's' and 'o': contributions from self and the other atoms
192 ps = {}
193 for a1 in range(natoms):
194 a2, d, r = self._get_neighbors(a1)
195 if len(a2) == 0:
196 continue
197 w, dwdroverw = self._calc_theta(r)
198 dsigma1s, dsigma1o = self._calc_dsigma1(a1, a2, r, w)
199 dsigma2s, dsigma2o = self._calc_dsigma2(a1, a2, r, w)
200 ps[a1] = {
201 'a2': a2,
202 'd': d,
203 'r': r,
204 'invr': 1.0 / r,
205 'w': w,
206 'dwdroverw': dwdroverw,
207 'dsigma1s': dsigma1s,
208 'dsigma1o': dsigma1o,
209 'dsigma2s': dsigma2s,
210 'dsigma2o': dsigma2o,
211 }
213 # deds is computed in _calc_e_c_a2
214 # since deds for all the atoms are used later in _calc_f_c_a2,
215 # _calc_e_c_a2 must be called beforehand for all the atoms
216 for a1, p in ps.items():
217 a2 = p['a2']
218 dsigma1s = p['dsigma1s']
219 self._calc_e_c_a2(a1, dsigma1s)
221 for a1, p in ps.items():
222 a2 = p['a2']
223 d = p['d']
224 invr = p['invr']
225 dwdroverw = p['dwdroverw']
226 dsigma1s = p['dsigma1s']
227 dsigma1o = p['dsigma1o']
228 dsigma2s = p['dsigma2s']
229 dsigma2o = p['dsigma2o']
230 self._calc_fs_c_a2(a1, a2, d, invr, dwdroverw, dsigma1s, dsigma1o)
231 self._calc_efs_a1(a1, a2, d, invr, dwdroverw, dsigma2s, dsigma2o)
233 # subtract E0 (ASAP convention)
234 self.energies -= self.par['E0'][self.ia2iz]
236 energy = np.add.reduce(self.energies, axis=0)
237 self.results['energy'] = self.results['free_energy'] = energy
238 self.results['energies'] = self.energies
239 self.results['forces'] = self.forces
241 if self.atoms.cell.rank == 3:
242 self.stress = (self.stress + self.stress.T) * 0.5 # symmetrize
243 self.stress /= self.atoms.get_volume()
244 self.results['stress'] = self.stress.flat[[0, 4, 8, 5, 2, 1]]
245 elif 'stress' in properties:
246 raise PropertyNotImplementedError
248 def _get_neighbors(self, a1):
249 positions = self.atoms.positions
250 cell = self.atoms.cell
251 neighbors, offsets = self.nl.get_neighbors(a1)
252 offsets = np.dot(offsets, cell)
253 d = positions[neighbors] + offsets - positions[a1]
254 r = np.sqrt(np.add.reduce(d**2, axis=1))
255 mask = r < self.rc_list
256 return neighbors[mask], d[mask], r[mask]
258 def _calc_theta(self, r):
259 """Calculate cutoff function and its r derivative"""
260 w = 1.0 / (1.0 + np.exp(self.acut * (r - self.rc)))
261 dwdroverw = self.acut * (w - 1.0)
262 return w, dwdroverw
264 def _calc_dsigma1(self, a1, a2, r, w):
265 """Calculate contributions of neighbors to sigma1"""
266 s0s = self.par['s0'][self.ia2iz[a1]]
267 s0o = self.par['s0'][self.ia2iz[a2]]
268 eta2s = self.par['eta2'][self.ia2iz[a1]]
269 eta2o = self.par['eta2'][self.ia2iz[a2]]
270 chi = self.chi[self.ia2iz[a1], self.ia2iz[a2]]
272 dsigma1s = np.exp(-eta2o * (r - beta * s0o)) * chi * w
273 dsigma1o = np.exp(-eta2s * (r - beta * s0s)) / chi * w
275 return dsigma1s, dsigma1o
277 def _calc_dsigma2(self, a1, a2, r, w):
278 """Calculate contributions of neighbors to sigma2"""
279 s0s = self.par['s0'][self.ia2iz[a1]]
280 s0o = self.par['s0'][self.ia2iz[a2]]
281 kappas = self.par['kappa'][self.ia2iz[a1]]
282 kappao = self.par['kappa'][self.ia2iz[a2]]
283 chi = self.chi[self.ia2iz[a1], self.ia2iz[a2]]
285 dsigma2s = np.exp(-kappao * (r / beta - s0o)) * chi * w
286 dsigma2o = np.exp(-kappas * (r / beta - s0s)) / chi * w
288 return dsigma2s, dsigma2o
290 def _calc_e_c_a2(self, a1, dsigma1s):
291 """Calculate E_c and the second term of E_AS and their s derivatives"""
292 e0s = self.par['E0'][self.ia2iz[a1]]
293 v0s = self.par['V0'][self.ia2iz[a1]]
294 eta2s = self.par['eta2'][self.ia2iz[a1]]
295 lmds = self.par['lambda'][self.ia2iz[a1]]
296 kappas = self.par['kappa'][self.ia2iz[a1]]
297 inv12gamma1s = self.par['inv12gamma1'][self.ia2iz[a1]]
299 sigma1 = np.add.reduce(dsigma1s)
300 ds = -1.0 * np.log(sigma1 * inv12gamma1s) / (beta * eta2s)
302 lmdsds = lmds * ds
303 expneglmdds = np.exp(-1.0 * lmdsds)
304 self.energies[a1] += e0s * (1.0 + lmdsds) * expneglmdds
305 self.deds[a1] += -1.0 * e0s * lmds * lmdsds * expneglmdds
307 sixv0expnegkppds = 6.0 * v0s * np.exp(-1.0 * kappas * ds)
308 self.energies[a1] += sixv0expnegkppds
309 self.deds[a1] += -1.0 * kappas * sixv0expnegkppds
311 self.deds[a1] /= -1.0 * beta * eta2s * sigma1 # factor from ds/dr
313 def _calc_efs_a1(self, a1, a2, d, invr, dwdroverw, dsigma2s, dsigma2o):
314 """Calculate the first term of E_AS and derivatives"""
315 neghalfv0overgamma2s = self.par['neghalfv0overgamma2'][self.ia2iz[a1]]
316 neghalfv0overgamma2o = self.par['neghalfv0overgamma2'][self.ia2iz[a2]]
317 kappas = self.par['kappa'][self.ia2iz[a1]]
318 kappao = self.par['kappa'][self.ia2iz[a2]]
320 es = neghalfv0overgamma2s * dsigma2s
321 eo = neghalfv0overgamma2o * dsigma2o
322 self.energies[a1] += 0.5 * np.add.reduce(es + eo, axis=0)
324 dedrs = es * (dwdroverw - kappao / beta)
325 dedro = eo * (dwdroverw - kappas / beta)
326 f = ((dedrs + dedro) * invr)[:, None] * d
327 self.forces[a1] += np.add.reduce(f, axis=0)
328 self.stress += 0.5 * np.dot(d.T, f) # compensate double counting
330 def _calc_fs_c_a2(self, a1, a2, d, invr, dwdroverw, dsigma1s, dsigma1o):
331 """Calculate forces and stress from E_c and the second term of E_AS"""
332 eta2s = self.par['eta2'][self.ia2iz[a1]]
333 eta2o = self.par['eta2'][self.ia2iz[a2]]
335 ddsigma1sdr = dsigma1s * (dwdroverw - eta2o)
336 ddsigma1odr = dsigma1o * (dwdroverw - eta2s)
337 dedrs = self.deds[a1] * ddsigma1sdr
338 dedro = self.deds[a2] * ddsigma1odr
339 f = ((dedrs + dedro) * invr)[:, None] * d
340 self.forces[a1] += np.add.reduce(f, axis=0)
341 self.stress += 0.5 * np.dot(d.T, f) # compensate double counting