Coverage for /builds/ase/ase/ase/calculators/eam.py: 87.50%

1# fmt: off 

2 

3# flake8: noqa 

4"""Calculator for the Embedded Atom Method Potential""" 

5 

6# eam.py 

7# Embedded Atom Method Potential 

8# These routines integrate with the ASE simulation environment 

9# Paul White (Oct 2012) 

10# UNCLASSIFIED 

11# License: See accompanying license files for details 

12 

13import os 

14 

15import numpy as np 

16from scipy.interpolate import InterpolatedUnivariateSpline as spline 

17 

18from ase.calculators.calculator import Calculator, all_changes 

19from ase.data import chemical_symbols 

20from ase.neighborlist import NeighborList 

21from ase.units import Bohr, Hartree 

22from ase.stress import full_3x3_to_voigt_6_stress 

23 

24 

25class EAM(Calculator): 

26 r""" 

27 

28 EAM Interface Documentation 

29 

30Introduction 

31============ 

32 

33The Embedded Atom Method (EAM) [1]_ is a classical potential which is 

34good for modelling metals, particularly fcc materials. Because it is 

35an equiaxial potential the EAM does not model directional bonds 

36well. However, the Angular Dependent Potential (ADP) [2]_ which is an 

37extended version of EAM is able to model directional bonds and is also 

38included in the EAM calculator. 

39 

40Generally all that is required to use this calculator is to supply a 

41potential file or as a set of functions that describe the potential. 

42The files containing the potentials for this calculator are not 

43included but many suitable potentials can be downloaded from The 

44Interatomic Potentials Repository Project at 

45https://www.ctcms.nist.gov/potentials/ 

46 

47Theory 

48====== 

49 

50A single element EAM potential is defined by three functions: the 

51embedded energy, electron density and the pair potential. A two 

52element alloy contains the individual three functions for each element 

53plus cross pair interactions. The ADP potential has two additional 

54sets of data to define the dipole and quadrupole directional terms for 

55each alloy and their cross interactions. 

56 

57The total energy `E_{\rm tot}` of an arbitrary arrangement of atoms is 

58given by the EAM potential as 

59 

60.. math:: 

61 E_\text{tot} = \sum_i F(\bar\rho_i) + \frac{1}{2}\sum_{i\ne j} \phi(r_{ij}) 

62 

63and 

64 

65.. math:: 

66 \bar\rho_i = \sum_j \rho(r_{ij}) 

67 

68where `F` is an embedding function, namely the energy to embed an atom `i` in 

69the combined electron density `\bar\rho_i` which is contributed from 

70each of its neighbouring atoms `j` by an amount `\rho(r_{ij})`, 

71`\phi(r_{ij})` is the pair potential function representing the energy 

72in bond `ij` which is due to the short-range electro-static 

73interaction between atoms, and `r_{ij}` is the distance between an 

74atom and its neighbour for that bond. 

75 

76The ADP potential is defined as 

77 

78.. math:: 

79 E_\text{tot} = \sum_i F(\bar\rho_i) + \frac{1}{2}\sum_{i\ne j} \phi(r_{ij}) 

80 + \frac{1}{2} \sum_{i,\alpha} (\mu_i^\alpha)^2 

81 + \frac{1}{2} \sum_{i,\alpha,\beta} (\lambda_i^{\alpha\beta})^2 

82 - \frac{1}{6} \sum_i \nu_i^2 

83 

84where `\mu_i^\alpha` is the dipole vector, `\lambda_i^{\alpha\beta}` 

85is the quadrupole tensor and `\nu_i` is the trace of 

86`\lambda_i^{\alpha\beta}`. 

87 

88The fs potential is defined as 

89 

90.. math:: 

91 E_i = F_\alpha (\sum_{j\neq i} \rho_{\alpha \beta}(r_{ij})) 

92 + \frac{1}{2}\sum_{j\neq i}\phi_{\alpha \beta}(r_{ij}) 

93 

94where `\alpha` and `\beta` are element types of atoms. This form is similar to 

95original EAM formula above, except that `\rho` and `\phi` are determined 

96by element types. 

97 

98Running the Calculator 

99====================== 

100 

101EAM calculates the cohesive atom energy, forces and stress. Internally the 

102potential functions are defined by splines which may be directly 

103supplied or created by reading the spline points from a data file from 

104which a spline function is created. The LAMMPS compatible ``.alloy``, ``.fs`` 

105and ``.adp`` formats are supported. The LAMMPS ``.eam`` format is 

106slightly different from the ``.alloy`` format and is currently not 

107supported. 

108 

109For example:: 

110 

111 from ase.calculators.eam import EAM 

112 

113 mishin = EAM(potential='Al99.eam.alloy') 

114 mishin.write_potential('new.eam.alloy') 

115 mishin.plot() 

116 

117 slab.calc = mishin 

118 slab.get_potential_energy() 

119 slab.get_forces() 

120 slab.get_stress() 

121 

122The breakdown of energy contribution from the indvidual components are 

123stored in the calculator instance ``.results['energy_components']`` 

124 

125Arguments 

126========= 

127 

128========================= ==================================================== 

129Keyword Description 

130========================= ==================================================== 

131``potential`` file of potential in ``.eam``, ``.alloy``, ``.adp`` or ``.fs`` 

132 format or file object 

133 (This is generally all you need to supply). 

134 For file object the ``form`` argument is required 

135 

136``elements[N]`` array of N element abbreviations 

137 

138``embedded_energy[N]`` arrays of embedded energy functions 

139 

140``electron_density[N]`` arrays of electron density functions 

141 

142``phi[N,N]`` arrays of pair potential functions 

143 

144``d_embedded_energy[N]`` arrays of derivative embedded energy functions 

145 

146``d_electron_density[N]`` arrays of derivative electron density functions 

147 

148``d_phi[N,N]`` arrays of derivative pair potentials functions 

149 

150``d[N,N], q[N,N]`` ADP dipole and quadrupole function 

151 

152``d_d[N,N], d_q[N,N]`` ADP dipole and quadrupole derivative functions 

153 

154``skin`` skin distance passed to NeighborList(). If no atom 

155 has moved more than the skin-distance since the last 

156 call to the ``update()`` method then the neighbor 

157 list can be reused. Defaults to 1.0. 

158 

159``form`` the form of the potential 

160 ``eam``, ``alloy``, ``adp`` or 

161 ``fs``. This will be determined from the file suffix 

162 or must be set if using equations or file object 

163 

164========================= ==================================================== 

165 

166 

167Additional parameters for writing potential files 

168================================================= 

169 

170The following parameters are only required for writing a potential in 

171``.alloy``, ``.adp`` or ``fs`` format file. 

172 

173========================= ==================================================== 

174Keyword Description 

175========================= ==================================================== 

176``header`` Three line text header. Default is standard message. 

177 

178``Z[N]`` Array of atomic number of each element 

179 

180``mass[N]`` Atomic mass of each element 

181 

182``a[N]`` Array of lattice parameters for each element 

183 

184``lattice[N]`` Lattice type 

185 

186``nrho`` No. of rho samples along embedded energy curve 

187 

188``drho`` Increment for sampling density 

189 

190``nr`` No. of radial points along density and pair 

191 potential curves 

192 

193``dr`` Increment for sampling radius 

194 

195========================= ==================================================== 

196 

197Special features 

198================ 

199 

200``.plot()`` 

201 Plots the individual functions. This may be called from multiple EAM 

202 potentials to compare the shape of the individual curves. This 

203 function requires the installation of the Matplotlib libraries. 

204 

205Notes/Issues 

206============= 

207 

208* Although currently not fast, this calculator can be good for trying 

209 small calculations or for creating new potentials by matching baseline 

210 data such as from DFT results. The format for these potentials is 

211 compatible with LAMMPS_ and so can be used either directly by LAMMPS or 

212 with the ASE LAMMPS calculator interface. 

213 

214* Supported formats are the LAMMPS_ ``.alloy`` and ``.adp``. The 

215 ``.eam`` format is currently not supported. The form of the 

216 potential will be determined from the file suffix. 

217 

218* Any supplied values will override values read from the file. 

219 

220* The derivative functions, if supplied, are only used to calculate 

221 forces. 

222 

223* There is a bug in early versions of scipy that will cause eam.py to 

224 crash when trying to evaluate splines of a potential with one 

225 neighbor such as caused by evaluating a dimer. 

226 

227.. _LAMMPS: http://lammps.sandia.gov/ 

228 

229.. [1] M.S. Daw and M.I. Baskes, Phys. Rev. Letters 50 (1983) 

230 1285. 

231 

232.. [2] Y. Mishin, M.J. Mehl, and D.A. Papaconstantopoulos, 

233 Acta Materialia 53 2005 4029--4041. 

234 

235 

236End EAM Interface Documentation 

237 """ 

238 

239 implemented_properties = ['energy', 'free_energy', 'forces', 'stress'] 

240 

241 default_parameters = dict( 

242 skin=1.0, 

243 potential=None, 

244 header=[b'EAM/ADP potential file\n', 

245 b'Generated from eam.py\n', 

246 b'blank\n']) 

247 

248 def __init__(self, restart=None, 

249 ignore_bad_restart_file=Calculator._deprecated, 

250 label=os.curdir, atoms=None, form=None, **kwargs): 

251 

252 self.form = form 

253 

254 if 'potential' in kwargs: 

255 self.read_potential(kwargs['potential']) 

256 

257 Calculator.__init__(self, restart, ignore_bad_restart_file, 

258 label, atoms, **kwargs) 

259 

260 valid_args = ('potential', 'elements', 'header', 'drho', 'dr', 

261 'cutoff', 'atomic_number', 'mass', 'a', 'lattice', 

262 'embedded_energy', 'electron_density', 'phi', 

263 # derivatives 

264 'd_embedded_energy', 'd_electron_density', 'd_phi', 

265 'd', 'q', 'd_d', 'd_q', # adp terms 

266 'skin', 'Z', 'nr', 'nrho', 'mass') 

267 

268 # set any additional keyword arguments 

269 for arg, val in self.parameters.items(): 

270 if arg in valid_args: 

271 setattr(self, arg, val) 

272 else: 

273 raise RuntimeError( 

274 f'unknown keyword arg "{arg}" : not in {valid_args}') 

275 

276 def set_form(self, name): 

277 """set the form variable based on the file name suffix""" 

278 extension = os.path.splitext(name)[1] 

279 

280 if extension == '.eam': 

281 self.form = 'eam' 

282 elif extension == '.alloy': 

283 self.form = 'alloy' 

284 elif extension == '.adp': 

285 self.form = 'adp' 

286 elif extension == '.fs': 

287 self.form = 'fs' 

288 else: 

289 raise RuntimeError(f'unknown file extension type: {extension}') 

290 

291 def read_potential(self, filename): 

292 """Reads a LAMMPS EAM file in alloy or adp format 

293 and creates the interpolation functions from the data 

294 """ 

295 

296 if isinstance(filename, str): 

297 with open(filename) as fd: 

298 self._read_potential(fd) 

299 else: 

300 fd = filename 

301 self._read_potential(fd) 

302 

303 def _read_potential(self, fd): 

304 if self.form is None: 

305 self.set_form(fd.name) 

306 

307 lines = fd.readlines() 

308 

309 def lines_to_list(lines): 

310 """Make the data one long line so as not to care how its formatted 

311 """ 

312 data = [] 

313 for line in lines: 

314 data.extend(line.split()) 

315 return data 

316 

317 if self.form == 'eam': # single element eam file (aka funcfl) 

318 self.header = lines[:1] 

319 

320 data = lines_to_list(lines[1:]) 

321 

322 # eam form is just like an alloy form for one element 

323 

324 self.Nelements = 1 

325 self.elements = [chemical_symbols[int(data[0])]] 

326 self.Z = np.array([data[0]], dtype=int) 

327 self.mass = np.array([data[1]]) 

328 self.a = np.array([data[2]]) 

329 self.lattice = [data[3]] 

330 

331 self.nrho = int(data[4]) 

332 self.drho = float(data[5]) 

333 self.nr = int(data[6]) 

334 self.dr = float(data[7]) 

335 self.cutoff = float(data[8]) 

336 

337 n = 9 + self.nrho 

338 self.embedded_data = np.array([np.float64(data[9:n])]) 

339 

340 self.rphi_data = np.zeros([self.Nelements, self.Nelements, 

341 self.nr]) 

342 

343 effective_charge = np.float64(data[n:n + self.nr]) 

344 # convert effective charges to rphi according to 

345 # http://lammps.sandia.gov/doc/pair_eam.html 

346 self.rphi_data[0, 0] = Bohr * Hartree * (effective_charge**2) 

347 

348 self.density_data = np.array( 

349 [np.float64(data[n + self.nr:n + 2 * self.nr])]) 

350 

351 elif self.form in ['alloy', 'adp']: 

352 self.header = lines[:3] 

353 i = 3 

354 

355 data = lines_to_list(lines[i:]) 

356 

357 self.Nelements = int(data[0]) 

358 d = 1 

359 self.elements = data[d:d + self.Nelements] 

360 d += self.Nelements 

361 

362 self.nrho = int(data[d]) 

363 self.drho = float(data[d + 1]) 

364 self.nr = int(data[d + 2]) 

365 self.dr = float(data[d + 3]) 

366 self.cutoff = float(data[d + 4]) 

367 

368 self.embedded_data = np.zeros([self.Nelements, self.nrho]) 

369 self.density_data = np.zeros([self.Nelements, self.nr]) 

370 self.Z = np.zeros([self.Nelements], dtype=int) 

371 self.mass = np.zeros([self.Nelements]) 

372 self.a = np.zeros([self.Nelements]) 

373 self.lattice = [] 

374 d += 5 

375 

376 # reads in the part of the eam file for each element 

377 for elem in range(self.Nelements): 

378 self.Z[elem] = int(data[d]) 

379 self.mass[elem] = float(data[d + 1]) 

380 self.a[elem] = float(data[d + 2]) 

381 self.lattice.append(data[d + 3]) 

382 d += 4 

383 

384 self.embedded_data[elem] = np.float64( 

385 data[d:(d + self.nrho)]) 

386 d += self.nrho 

387 self.density_data[elem] = np.float64(data[d:(d + self.nr)]) 

388 d += self.nr 

389 

390 # reads in the r*phi data for each interaction between elements 

391 self.rphi_data = np.zeros([self.Nelements, self.Nelements, 

392 self.nr]) 

393 

394 for i in range(self.Nelements): 

395 for j in range(i + 1): 

396 self.rphi_data[j, i] = np.float64(data[d:(d + self.nr)]) 

397 d += self.nr 

398 

399 elif self.form == 'fs': 

400 self.header = lines[:3] 

401 i = 3 

402 

403 data = lines_to_list(lines[i:]) 

404 

405 self.Nelements = int(data[0]) 

406 d = 1 

407 self.elements = data[d:d + self.Nelements] 

408 d += self.Nelements 

409 

410 self.nrho = int(data[d]) 

411 self.drho = float(data[d + 1]) 

412 self.nr = int(data[d + 2]) 

413 self.dr = float(data[d + 3]) 

414 self.cutoff = float(data[d + 4]) 

415 

416 self.embedded_data = np.zeros([self.Nelements, self.nrho]) 

417 self.density_data = np.zeros([self.Nelements, self.Nelements, 

418 self.nr]) 

419 self.Z = np.zeros([self.Nelements], dtype=int) 

420 self.mass = np.zeros([self.Nelements]) 

421 self.a = np.zeros([self.Nelements]) 

422 self.lattice = [] 

423 d += 5 

424 

425 # reads in the part of the eam file for each element 

426 for elem in range(self.Nelements): 

427 self.Z[elem] = int(data[d]) 

428 self.mass[elem] = float(data[d + 1]) 

429 self.a[elem] = float(data[d + 2]) 

430 self.lattice.append(data[d + 3]) 

431 d += 4 

432 

433 self.embedded_data[elem] = np.float64( 

434 data[d:(d + self.nrho)]) 

435 d += self.nrho 

436 self.density_data[elem, :, :] = np.float64( 

437 data[d:(d + self.nr * self.Nelements)]).reshape([ 

438 self.Nelements, self.nr]) 

439 d += self.nr * self.Nelements 

440 

441 # reads in the r*phi data for each interaction between elements 

442 self.rphi_data = np.zeros([self.Nelements, self.Nelements, 

443 self.nr]) 

444 

445 for i in range(self.Nelements): 

446 for j in range(i + 1): 

447 self.rphi_data[j, i] = np.float64(data[d:(d + self.nr)]) 

448 d += self.nr 

449 

450 self.r = np.arange(0, self.nr) * self.dr 

451 self.rho = np.arange(0, self.nrho) * self.drho 

452 

453 # choose the set_splines method according to the type 

454 if self.form == 'fs': 

455 self.set_fs_splines() 

456 else: 

457 self.set_splines() 

458 

459 if self.form == 'adp': 

460 self.read_adp_data(data, d) 

461 self.set_adp_splines() 

462 

463 def set_splines(self): 

464 # this section turns the file data into three functions (and 

465 # derivative functions) that define the potential 

466 self.embedded_energy = np.empty(self.Nelements, object) 

467 self.electron_density = np.empty(self.Nelements, object) 

468 self.d_embedded_energy = np.empty(self.Nelements, object) 

469 self.d_electron_density = np.empty(self.Nelements, object) 

470 

471 for i in range(self.Nelements): 

472 self.embedded_energy[i] = spline(self.rho, 

473 self.embedded_data[i], k=3) 

474 self.electron_density[i] = spline(self.r, 

475 self.density_data[i], k=3) 

476 self.d_embedded_energy[i] = self.deriv(self.embedded_energy[i]) 

477 self.d_electron_density[i] = self.deriv(self.electron_density[i]) 

478 

479 self.phi = np.empty([self.Nelements, self.Nelements], object) 

480 self.d_phi = np.empty([self.Nelements, self.Nelements], object) 

481 

482 # ignore the first point of the phi data because it is forced 

483 # to go through zero due to the r*phi format in alloy and adp 

484 for i in range(self.Nelements): 

485 for j in range(i, self.Nelements): 

486 self.phi[i, j] = spline( 

487 self.r[1:], 

488 self.rphi_data[i, j][1:] / self.r[1:], k=3) 

489 

490 self.d_phi[i, j] = self.deriv(self.phi[i, j]) 

491 

492 if j != i: 

493 self.phi[j, i] = self.phi[i, j] 

494 self.d_phi[j, i] = self.d_phi[i, j] 

495 

496 def set_fs_splines(self): 

497 self.embedded_energy = np.empty(self.Nelements, object) 

498 self.electron_density = np.empty( 

499 [self.Nelements, self.Nelements], object) 

500 self.d_embedded_energy = np.empty(self.Nelements, object) 

501 self.d_electron_density = np.empty( 

502 [self.Nelements, self.Nelements], object) 

503 

504 for i in range(self.Nelements): 

505 self.embedded_energy[i] = spline(self.rho, 

506 self.embedded_data[i], k=3) 

507 self.d_embedded_energy[i] = self.deriv(self.embedded_energy[i]) 

508 for j in range(self.Nelements): 

509 self.electron_density[i, j] = spline( 

510 self.r, self.density_data[i, j], k=3) 

511 self.d_electron_density[i, j] = self.deriv( 

512 self.electron_density[i, j]) 

513 

514 self.phi = np.empty([self.Nelements, self.Nelements], object) 

515 self.d_phi = np.empty([self.Nelements, self.Nelements], object) 

516 

517 for i in range(self.Nelements): 

518 for j in range(i, self.Nelements): 

519 self.phi[i, j] = spline( 

520 self.r[1:], 

521 self.rphi_data[i, j][1:] / self.r[1:], k=3) 

522 

523 self.d_phi[i, j] = self.deriv(self.phi[i, j]) 

524 

525 if j != i: 

526 self.phi[j, i] = self.phi[i, j] 

527 self.d_phi[j, i] = self.d_phi[i, j] 

528 

529 def set_adp_splines(self): 

530 self.d = np.empty([self.Nelements, self.Nelements], object) 

531 self.d_d = np.empty([self.Nelements, self.Nelements], object) 

532 self.q = np.empty([self.Nelements, self.Nelements], object) 

533 self.d_q = np.empty([self.Nelements, self.Nelements], object) 

534 

535 for i in range(self.Nelements): 

536 for j in range(i, self.Nelements): 

537 self.d[i, j] = spline(self.r[1:], self.d_data[i, j][1:], k=3) 

538 self.d_d[i, j] = self.deriv(self.d[i, j]) 

539 self.q[i, j] = spline(self.r[1:], self.q_data[i, j][1:], k=3) 

540 self.d_q[i, j] = self.deriv(self.q[i, j]) 

541 

542 # make symmetrical 

543 if j != i: 

544 self.d[j, i] = self.d[i, j] 

545 self.d_d[j, i] = self.d_d[i, j] 

546 self.q[j, i] = self.q[i, j] 

547 self.d_q[j, i] = self.d_q[i, j] 

548 

549 def read_adp_data(self, data, d): 

550 """read in the extra adp data from the potential file""" 

551 

552 self.d_data = np.zeros([self.Nelements, self.Nelements, self.nr]) 

553 # should be non symmetrical combinations of 2 

554 for i in range(self.Nelements): 

555 for j in range(i + 1): 

556 self.d_data[j, i] = data[d:d + self.nr] 

557 d += self.nr 

558 

559 self.q_data = np.zeros([self.Nelements, self.Nelements, self.nr]) 

560 # should be non symmetrical combinations of 2 

561 for i in range(self.Nelements): 

562 for j in range(i + 1): 

563 self.q_data[j, i] = data[d:d + self.nr] 

564 d += self.nr 

565 

566 def write_potential(self, filename, nc=1, numformat='%.8e'): 

567 """Writes out the potential in the format given by the form 

568 variable to 'filename' with a data format that is nc columns 

569 wide. Note: array lengths need to be an exact multiple of nc 

570 """ 

571 

572 with open(filename, 'wb') as fd: 

573 self._write_potential(fd, nc=nc, numformat=numformat) 

574 

575 def _write_potential(self, fd, nc, numformat): 

576 assert self.nr % nc == 0 

577 assert self.nrho % nc == 0 

578 

579 for line in self.header: 

580 fd.write(line) 

581 

582 fd.write(f'{self.Nelements} '.encode()) 

583 fd.write(' '.join(self.elements).encode() + b'\n') 

584 

585 fd.write(('%d %f %d %f %f \n' % 

586 (self.nrho, self.drho, self.nr, 

587 self.dr, self.cutoff)).encode()) 

588 

589 # start of each section for each element 

590# rs = np.linspace(0, self.nr * self.dr, self.nr) 

591# rhos = np.linspace(0, self.nrho * self.drho, self.nrho) 

592 

593 rs = np.arange(0, self.nr) * self.dr 

594 rhos = np.arange(0, self.nrho) * self.drho 

595 

596 for i in range(self.Nelements): 

597 fd.write(('%d %f %f %s\n' % 

598 (self.Z[i], self.mass[i], 

599 self.a[i], str(self.lattice[i]))).encode()) 

600 np.savetxt(fd, 

601 self.embedded_energy[i](rhos).reshape(self.nrho // nc, 

602 nc), 

603 fmt=nc * [numformat]) 

604 if self.form == 'fs': 

605 for j in range(self.Nelements): 

606 np.savetxt(fd, 

607 self.electron_density[i, j](rs).reshape( 

608 self.nr // nc, nc), 

609 fmt=nc * [numformat]) 

610 else: 

611 np.savetxt(fd, 

612 self.electron_density[i](rs).reshape( 

613 self.nr // nc, nc), 

614 fmt=nc * [numformat]) 

615 

616 # write out the pair potentials in Lammps DYNAMO setfl format 

617 # as r*phi for alloy format 

618 for i in range(self.Nelements): 

619 for j in range(i, self.Nelements): 

620 np.savetxt(fd, 

621 (rs * self.phi[i, j](rs)).reshape(self.nr // nc, 

622 nc), 

623 fmt=nc * [numformat]) 

624 

625 if self.form == 'adp': 

626 # these are the u(r) or dipole values 

627 for i in range(self.Nelements): 

628 for j in range(i + 1): 

629 np.savetxt(fd, self.d_data[i, j]) 

630 

631 # these are the w(r) or quadrupole values 

632 for i in range(self.Nelements): 

633 for j in range(i + 1): 

634 np.savetxt(fd, self.q_data[i, j]) 

635 

636 def update(self, atoms): 

637 # check all the elements are available in the potential 

638 self.Nelements = len(self.elements) 

639 elements = np.unique(atoms.get_chemical_symbols()) 

640 unavailable = np.logical_not( 

641 np.array([item in self.elements for item in elements])) 

642 

643 if np.any(unavailable): 

644 raise RuntimeError( 

645 f'These elements are not in the potential: {elements[unavailable]}') 

646 

647 # cutoffs need to be a vector for NeighborList 

648 cutoffs = self.cutoff * np.ones(len(atoms)) 

649 

650 # convert the elements to an index of the position 

651 # in the eam format 

652 self.index = np.array([self.elements.index(el) 

653 for el in atoms.get_chemical_symbols()]) 

654 self.pbc = atoms.get_pbc() 

655 

656 # since we need the contribution of all neighbors to the 

657 # local electron density we cannot just calculate and use 

658 # one way neighbors 

659 self.neighbors = NeighborList(cutoffs, 

660 skin=self.parameters.skin, 

661 self_interaction=False, 

662 bothways=True) 

663 self.neighbors.update(atoms) 

664 

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

666 system_changes=all_changes): 

667 """EAM Calculator 

668 

669 atoms: Atoms object 

670 Contains positions, unit-cell, ... 

671 properties: list of str 

672 List of what needs to be calculated. Can be any combination 

673 of 'energy', 'forces' 

674 system_changes: list of str 

675 List of what has changed since last calculation. Can be 

676 any combination of these five: 'positions', 'numbers', 'cell', 

677 'pbc', 'initial_charges' and 'initial_magmoms'. 

678 """ 

679 

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

681 

682 # we shouldn't really recalc if charges or magmos change 

683 if len(system_changes) > 0: # something wrong with this way 

684 self.update(self.atoms) 

685 self.calculate_energy(self.atoms) 

686 

687 if 'forces' in properties or 'stress' in properties: 

688 self.calculate_forces(self.atoms) 

689 

690 # check we have all the properties requested 

691 for property in properties: 

692 if property not in self.results: 

693 if property == 'energy': 

694 self.calculate_energy(self.atoms) 

695 

696 if property in ['forces', 'stress']: 

697 self.calculate_forces(self.atoms) 

698 

699 # we need to remember the previous state of parameters 

700# if 'potential' in parameter_changes and potential != None: 

701# self.read_potential(potential) 

702 

703 def calculate_energy(self, atoms): 

704 """Calculate the energy 

705 the energy is made up of the ionic or pair interaction and 

706 the embedding energy of each atom into the electron cloud 

707 generated by its neighbors 

708 """ 

709 

710 pair_energy = 0.0 

711 embedding_energy = 0.0 

712 mu_energy = 0.0 

713 lam_energy = 0.0 

714 trace_energy = 0.0 

715 

716 self.total_density = np.zeros(len(atoms)) 

717 if self.form == 'adp': 

718 self.mu = np.zeros([len(atoms), 3]) 

719 self.lam = np.zeros([len(atoms), 3, 3]) 

720 

721 for i in range(len(atoms)): # this is the atom to be embedded 

722 neighbors, offsets = self.neighbors.get_neighbors(i) 

723 offset = np.dot(offsets, atoms.get_cell()) 

724 

725 rvec = (atoms.positions[neighbors] + offset - 

726 atoms.positions[i]) 

727 

728 # calculate the distance to the nearest neighbors 

729 r = np.sqrt(np.sum(np.square(rvec), axis=1)) # fast 

730# r = np.apply_along_axis(np.linalg.norm, 1, rvec) # sloow 

731 

732 nearest = np.arange(len(r))[r <= self.cutoff] 

733 for j_index in range(self.Nelements): 

734 use = self.index[neighbors[nearest]] == j_index 

735 if not use.any(): 

736 continue 

737 pair_energy += np.sum(self.phi[self.index[i], j_index]( 

738 r[nearest][use])) / 2. 

739 

740 if self.form == 'fs': 

741 density = np.sum( 

742 self.electron_density[j_index, 

743 self.index[i]](r[nearest][use])) 

744 else: 

745 density = np.sum( 

746 self.electron_density[j_index](r[nearest][use])) 

747 self.total_density[i] += density 

748 

749 if self.form == 'adp': 

750 self.mu[i] += self.adp_dipole( 

751 r[nearest][use], 

752 rvec[nearest][use], 

753 self.d[self.index[i], j_index]) 

754 

755 self.lam[i] += self.adp_quadrupole( 

756 r[nearest][use], 

757 rvec[nearest][use], 

758 self.q[self.index[i], j_index]) 

759 

760 # add in the electron embedding energy 

761 embedding_energy += self.embedded_energy[self.index[i]]( 

762 self.total_density[i]) 

763 

764 components = dict(pair=pair_energy, embedding=embedding_energy) 

765 

766 if self.form == 'adp':