Coverage for /builds/hweiske/ase/ase/calculators/emt.py: 98.90%
182 statements
« prev ^ index » next coverage.py v7.2.7, created at 2024-04-22 11:22 +0000
1"""Effective medium theory potential."""
2from collections import defaultdict
3from math import log, sqrt
5import numpy as np
7from ase.calculators.calculator import (Calculator,
8 PropertyNotImplementedError,
9 all_changes)
10from ase.data import atomic_numbers, chemical_symbols
11from ase.neighborlist import NeighborList
12from ase.units import Bohr
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)}
30beta = 1.809 # (16 * pi / 3)**(1.0 / 3) / 2**0.5, preserve historical rounding
33class EMT(Calculator):
34 """Python implementation of the Effective Medium Potential.
36 Supports the following standard EMT metals:
37 Al, Cu, Ag, Au, Ni, Pd and Pt.
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
44 Parameters
45 ----------
46 asap_cutoff : bool, default: False
47 If True the cutoff mimics how ASAP does it; most importantly the global
48 cutoff is chosen from the largest atom present in the simulation.
49 If False it is chosen from the largest atom in the parameter table.
50 True gives the behaviour of ASAP and older EMT implementations,
51 although the results are not bitwise identical.
53 Notes
54 -----
55 Formulation mostly follows Jacobsen *et al*. [1]_
56 `Documentation in ASAP can also be referred to <https://gitlab.com/asap/
57 asap/blob/master/docs/manual/potentials/emt.pdf>`_.
59 .. [1] K. W. Jacobsen, P. Stoltze, and J. K. Nørskov,
60 Surf. Sci. 366, 394 (1996).
61 """
62 implemented_properties = ['energy', 'free_energy', 'energies', 'forces',
63 'stress', 'magmom', 'magmoms']
65 nolabel = True
67 default_parameters = {'asap_cutoff': False}
69 def __init__(self, **kwargs):
70 Calculator.__init__(self, **kwargs)
72 def initialize(self, atoms):
73 self.rc, self.rc_list, self.acut = self._calc_cutoff(atoms)
75 numbers = atoms.get_atomic_numbers()
77 # ia2iz : map from idx of atoms to idx of atomic numbers in self.par
78 unique_numbers, self.ia2iz = np.unique(numbers, return_inverse=True)
79 self.par = defaultdict(lambda: np.empty(len(unique_numbers)))
80 for i, Z in enumerate(unique_numbers):
81 sym = chemical_symbols[Z]
82 if sym not in parameters:
83 raise NotImplementedError(f'No EMT-potential for {sym}')
84 p = parameters[sym]
85 s0 = p[1] * Bohr
86 eta2 = p[3] / Bohr
87 kappa = p[4] / Bohr
88 gamma1, gamma2 = self._calc_gammas(s0, eta2, kappa)
89 self.par['Z'][i] = Z
90 self.par['E0'][i] = p[0]
91 self.par['s0'][i] = s0
92 self.par['V0'][i] = p[2]
93 self.par['eta2'][i] = eta2
94 self.par['kappa'][i] = kappa
95 self.par['lambda'][i] = p[5] / Bohr
96 self.par['n0'][i] = p[6] / Bohr**3
97 self.par['inv12gamma1'][i] = 1.0 / (12.0 * gamma1)
98 self.par['neghalfv0overgamma2'][i] = -0.5 * p[2] / gamma2
100 self.chi = self.par['n0'][None, :] / self.par['n0'][:, None]
102 self.energies = np.empty(len(atoms))
103 self.forces = np.empty((len(atoms), 3))
104 self.stress = np.empty((3, 3))
105 self.deds = np.empty(len(atoms))
107 self.nl = NeighborList([0.5 * self.rc_list] * len(atoms),
108 self_interaction=False, bothways=True)
110 def _calc_cutoff(self, atoms):
111 """Calculate parameters of the logistic smoothing function etc.
113 The logistic smoothing function is given by
115 .. math:
117 w(r) = \\frac{1}{1 + \\exp a (r - r_\\mathrm{c})}
119 Returns
120 -------
121 rc : float
122 "Midpoint" of the logistic smoothing function, set to be the mean
123 of the 3rd and the 4th nearest-neighbor distances in FCC.
124 rc_list : float
125 Cutoff radius for the neighbor search, set to be slightly larger
126 than ``rc`` depending on ``asap_cutoff``.
127 acut : float
128 "Slope" of the smoothing function, set for the smoothing function
129 value to be ``1e-4`` at the 4th nearest-neighbor distance in FCC.
131 Notes
132 -----
133 ``maxseq`` is the present FCC Wigner-Seitz radius. ``beta * maxseq``
134 (`r1nn`) is the corresponding 1st nearest-neighbor distance in FCC.
135 The 2nd, 3rd, 4th nearest-neighbor distances in FCC are given using
136 ``r1nn`` by ``sqrt(2) * r1nn``, ``sqrt(3) * r1nn``, ``sqrt(4) * r1nn``,
137 respectively.
138 """
139 numbers = atoms.get_atomic_numbers()
140 if self.parameters['asap_cutoff']:
141 relevant_pars = {
142 symb: p
143 for symb, p in parameters.items()
144 if atomic_numbers[symb] in numbers
145 }
146 else:
147 relevant_pars = parameters
148 maxseq = max(par[1] for par in relevant_pars.values()) * Bohr
149 r1nn = beta * maxseq # 1st NN distance in FCC
150 rc = r1nn * 0.5 * (sqrt(3.0) + 2.0) # mean of 3NN and 4NN dists.
151 r4nn = r1nn * 2.0 # 4NN distance in FCC
152 eps = 1e-4 # value at r4nn, should be small
154 # "slope" is set so that the function value becomes eps at r4nn
155 acut = log(1.0 / eps - 1.0) / (r4nn - rc)
157 rc_list = rc * 1.045 if self.parameters['asap_cutoff'] else rc + 0.5
159 return rc, rc_list, acut
161 def _calc_gammas(self, s0, eta2, kappa):
162 n = np.array([12, 6, 24]) # numbers of 1, 2, 3NN sites in fcc
163 r = beta * s0 * np.sqrt([1.0, 2.0, 3.0]) # distances of 1, 2, 3NNs
164 w = 1.0 / (1.0 + np.exp(self.acut * (r - self.rc)))
165 x = n * w / 12.0
166 gamma1 = x @ np.exp(-eta2 * (r - beta * s0))
167 gamma2 = x @ np.exp(-kappa / beta * (r - beta * s0))
168 return gamma1, gamma2
170 def calculate(self, atoms=None, properties=['energy'],
171 system_changes=all_changes):
172 Calculator.calculate(self, atoms, properties, system_changes)
174 if 'numbers' in system_changes:
175 self.initialize(self.atoms)
177 self.nl.update(self.atoms)
179 self.energies[:] = 0.0
180 self.forces[:] = 0.0
181 self.stress[:] = 0.0
182 self.deds[:] = 0.0
184 natoms = len(self.atoms)
186 # store nearest neighbor info for all the atoms
187 # suffixes 's' and 'o': contributions from self and the other atoms
188 ps = {}
189 for a1 in range(natoms):
190 a2, d, r = self._get_neighbors(a1)
191 if len(a2) == 0:
192 continue
193 w, dwdroverw = self._calc_theta(r)
194 dsigma1s, dsigma1o = self._calc_dsigma1(a1, a2, r, w)
195 dsigma2s, dsigma2o = self._calc_dsigma2(a1, a2, r, w)
196 ps[a1] = {
197 'a2': a2,
198 'd': d,
199 'r': r,
200 'invr': 1.0 / r,
201 'w': w,
202 'dwdroverw': dwdroverw,
203 'dsigma1s': dsigma1s,
204 'dsigma1o': dsigma1o,
205 'dsigma2s': dsigma2s,
206 'dsigma2o': dsigma2o,
207 }
209 # deds is computed in _calc_e_c_a2
210 # since deds for all the atoms are used later in _calc_f_c_a2,
211 # _calc_e_c_a2 must be called beforehand for all the atoms
212 for a1, p in ps.items():
213 a2 = p['a2']
214 dsigma1s = p['dsigma1s']
215 self._calc_e_c_a2(a1, dsigma1s)
217 for a1, p in ps.items():
218 a2 = p['a2']
219 d = p['d']
220 invr = p['invr']
221 dwdroverw = p['dwdroverw']
222 dsigma1s = p['dsigma1s']
223 dsigma1o = p['dsigma1o']
224 dsigma2s = p['dsigma2s']
225 dsigma2o = p['dsigma2o']
226 self._calc_fs_c_a2(a1, a2, d, invr, dwdroverw, dsigma1s, dsigma1o)
227 self._calc_efs_a1(a1, a2, d, invr, dwdroverw, dsigma2s, dsigma2o)
229 # subtract E0 (ASAP convention)
230 self.energies -= self.par['E0'][self.ia2iz]
232 energy = np.add.reduce(self.energies, axis=0)
233 self.results['energy'] = self.results['free_energy'] = energy
234 self.results['energies'] = self.energies
235 self.results['forces'] = self.forces
237 if self.atoms.cell.rank == 3:
238 self.stress = (self.stress + self.stress.T) * 0.5 # symmetrize
239 self.stress /= self.atoms.get_volume()
240 self.results['stress'] = self.stress.flat[[0, 4, 8, 5, 2, 1]]
241 elif 'stress' in properties:
242 raise PropertyNotImplementedError
244 def _get_neighbors(self, a1):
245 positions = self.atoms.positions
246 cell = self.atoms.cell
247 neighbors, offsets = self.nl.get_neighbors(a1)
248 offsets = np.dot(offsets, cell)
249 d = positions[neighbors] + offsets - positions[a1]
250 r = np.sqrt(np.add.reduce(d**2, axis=1))
251 mask = r < self.rc_list
252 return neighbors[mask], d[mask], r[mask]
254 def _calc_theta(self, r):
255 """Calculate cutoff function and its r derivative"""
256 w = 1.0 / (1.0 + np.exp(self.acut * (r - self.rc)))
257 dwdroverw = self.acut * (w - 1.0)
258 return w, dwdroverw
260 def _calc_dsigma1(self, a1, a2, r, w):
261 """Calculate contributions of neighbors to sigma1"""
262 s0s = self.par['s0'][self.ia2iz[a1]]
263 s0o = self.par['s0'][self.ia2iz[a2]]
264 eta2s = self.par['eta2'][self.ia2iz[a1]]
265 eta2o = self.par['eta2'][self.ia2iz[a2]]
266 chi = self.chi[self.ia2iz[a1], self.ia2iz[a2]]
268 dsigma1s = np.exp(-eta2o * (r - beta * s0o)) * chi * w
269 dsigma1o = np.exp(-eta2s * (r - beta * s0s)) / chi * w
271 return dsigma1s, dsigma1o
273 def _calc_dsigma2(self, a1, a2, r, w):
274 """Calculate contributions of neighbors to sigma2"""
275 s0s = self.par['s0'][self.ia2iz[a1]]
276 s0o = self.par['s0'][self.ia2iz[a2]]
277 kappas = self.par['kappa'][self.ia2iz[a1]]
278 kappao = self.par['kappa'][self.ia2iz[a2]]
279 chi = self.chi[self.ia2iz[a1], self.ia2iz[a2]]
281 dsigma2s = np.exp(-kappao * (r / beta - s0o)) * chi * w
282 dsigma2o = np.exp(-kappas * (r / beta - s0s)) / chi * w
284 return dsigma2s, dsigma2o
286 def _calc_e_c_a2(self, a1, dsigma1s):
287 """Calculate E_c and the second term of E_AS and their s derivatives"""
288 e0s = self.par['E0'][self.ia2iz[a1]]
289 v0s = self.par['V0'][self.ia2iz[a1]]
290 eta2s = self.par['eta2'][self.ia2iz[a1]]
291 lmds = self.par['lambda'][self.ia2iz[a1]]
292 kappas = self.par['kappa'][self.ia2iz[a1]]
293 inv12gamma1s = self.par['inv12gamma1'][self.ia2iz[a1]]
295 sigma1 = np.add.reduce(dsigma1s)
296 ds = -1.0 * np.log(sigma1 * inv12gamma1s) / (beta * eta2s)
298 lmdsds = lmds * ds
299 expneglmdds = np.exp(-1.0 * lmdsds)
300 self.energies[a1] += e0s * (1.0 + lmdsds) * expneglmdds
301 self.deds[a1] += -1.0 * e0s * lmds * lmdsds * expneglmdds
303 sixv0expnegkppds = 6.0 * v0s * np.exp(-1.0 * kappas * ds)
304 self.energies[a1] += sixv0expnegkppds
305 self.deds[a1] += -1.0 * kappas * sixv0expnegkppds
307 self.deds[a1] /= -1.0 * beta * eta2s * sigma1 # factor from ds/dr
309 def _calc_efs_a1(self, a1, a2, d, invr, dwdroverw, dsigma2s, dsigma2o):
310 """Calculate the first term of E_AS and derivatives"""
311 neghalfv0overgamma2s = self.par['neghalfv0overgamma2'][self.ia2iz[a1]]
312 neghalfv0overgamma2o = self.par['neghalfv0overgamma2'][self.ia2iz[a2]]
313 kappas = self.par['kappa'][self.ia2iz[a1]]
314 kappao = self.par['kappa'][self.ia2iz[a2]]
316 es = neghalfv0overgamma2s * dsigma2s
317 eo = neghalfv0overgamma2o * dsigma2o
318 self.energies[a1] += 0.5 * np.add.reduce(es + eo, axis=0)
320 dedrs = es * (dwdroverw - kappao / beta)
321 dedro = eo * (dwdroverw - kappas / beta)
322 f = ((dedrs + dedro) * invr)[:, None] * d
323 self.forces[a1] += np.add.reduce(f, axis=0)
324 self.stress += 0.5 * np.dot(d.T, f) # compensate double counting
326 def _calc_fs_c_a2(self, a1, a2, d, invr, dwdroverw, dsigma1s, dsigma1o):
327 """Calculate forces and stress from E_c and the second term of E_AS"""
328 eta2s = self.par['eta2'][self.ia2iz[a1]]
329 eta2o = self.par['eta2'][self.ia2iz[a2]]
331 ddsigma1sdr = dsigma1s * (dwdroverw - eta2o)
332 ddsigma1odr = dsigma1o * (dwdroverw - eta2s)
333 dedrs = self.deds[a1] * ddsigma1sdr
334 dedro = self.deds[a2] * ddsigma1odr
335 f = ((dedrs + dedro) * invr)[:, None] * d
336 self.forces[a1] += np.add.reduce(f, axis=0)
337 self.stress += 0.5 * np.dot(d.T, f) # compensate double counting