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