Coverage for /builds/kinetik161/ase/ase/calculators/emt.py: 90.68%
161 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Effective medium theory potential."""
3from math import exp, 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 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', 'free_energy', 'energies', 'forces',
54 'stress', 'magmom', 'magmoms']
56 nolabel = True
58 default_parameters = {'asap_cutoff': False}
60 def __init__(self, **kwargs):
61 Calculator.__init__(self, **kwargs)
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 {}'
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))
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}
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']
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))
124 self.nl = NeighborList([0.5 * self.rc_list] * len(atoms),
125 self_interaction=False)
127 def calculate(self, atoms=None, properties=['energy'],
128 system_changes=all_changes):
129 Calculator.calculate(self, atoms, properties, system_changes)
131 if 'numbers' in system_changes:
132 self.initialize(self.atoms)
134 positions = self.atoms.positions
135 numbers = self.atoms.numbers
136 cell = self.atoms.cell
138 self.nl.update(self.atoms)
140 self.energy = 0.0
141 self.energies[:] = 0
142 self.sigma1[:] = 0.0
143 self.forces[:] = 0.0
144 self.stress[:] = 0.0
146 natoms = len(self.atoms)
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])
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
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])
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
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
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'])
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)
242def main():
243 import sys
245 from ase.io import read, write
246 inputfile = sys.argv[1]
247 outputfile = sys.argv[2]
248 atoms = read(inputfile)
249 atoms.calc = EMT()
250 atoms.get_stress()
251 write(outputfile, atoms)
254if __name__ == '__main__':
255 main()