Coverage for /builds/kinetik161/ase/ase/calculators/eam.py: 82.56%

1# flake8: noqa 

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

3 

4# eam.py 

5# Embedded Atom Method Potential 

6# These routines integrate with the ASE simulation environment 

7# Paul White (Oct 2012) 

8# UNCLASSIFIED 

9# License: See accompanying license files for details 

10 

11import os 

12 

13import numpy as np 

14from scipy.interpolate import InterpolatedUnivariateSpline as spline 

15 

16from ase.calculators.calculator import Calculator, all_changes 

17from ase.neighborlist import NeighborList 

18from ase.units import Bohr, Hartree 

19 

20 

21class EAM(Calculator): 

22 r""" 

23 

24 EAM Interface Documentation 

25 

26Introduction 

27============ 

28 

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

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

31an equiaxial potential the EAM does not model directional bonds 

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

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

34included in the EAM calculator. 

35 

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

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

38The files containing the potentials for this calculator are not 

39included but many suitable potentials can be downloaded from The 

40Interatomic Potentials Repository Project at 

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

42 

43Theory 

44====== 

45 

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

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

48element alloy contains the individual three functions for each element 

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

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

51each alloy and their cross interactions. 

52 

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

54given by the EAM potential as 

55 

56.. math:: 

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

58 

59and 

60 

61.. math:: 

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

63 

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

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

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

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

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

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

70atom and its neighbour for that bond. 

71 

72The ADP potential is defined as 

73 

74.. math:: 

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

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

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

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

79 

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

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

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

83 

84The fs potential is defined as 

85 

86.. math:: 

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

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

89 

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

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

92by element types. 

93 

94Running the Calculator 

95====================== 

96 

97EAM calculates the cohesive atom energy and forces. Internally the 

98potential functions are defined by splines which may be directly 

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

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

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

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

103supported. 

104 

105For example:: 

106 

107 from ase.calculators.eam import EAM 

108 

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

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

111 mishin.plot() 

112 

113 slab.calc = mishin 

114 slab.get_potential_energy() 

115 slab.get_forces() 

116 

117The breakdown of energy contribution from the indvidual components are 

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

119 

120Arguments 

121========= 

122 

123========================= ==================================================== 

124Keyword Description 

125========================= ==================================================== 

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

127 format or file object 

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

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

130 

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

132 

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

134 

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

136 

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

138 

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

140 

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

142 

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

144 

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

146 

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

148 

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

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

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

152 list can be reused. Defaults to 1.0. 

153 

154``form`` the form of the potential 

155 ``eam``, ``alloy``, ``adp`` or 

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

157 or must be set if using equations or file object 

158 

159========================= ==================================================== 

160 

161 

162Additional parameters for writing potential files 

163================================================= 

164 

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

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

167 

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

169Keyword Description 

170========================= ==================================================== 

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

172 

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

174 

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

176 

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

178 

179``lattice[N]`` Lattice type 

180 

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

182 

183``drho`` Increment for sampling density 

184 

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

186 potential curves 

187 

188``dr`` Increment for sampling radius 

189 

190========================= ==================================================== 

191 

192Special features 

193================ 

194 

195``.plot()`` 

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

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

198 function requires the installation of the Matplotlib libraries. 

199 

200Notes/Issues 

201============= 

202 

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

204 small calculations or for creating new potentials by matching baseline 

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

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

207 with the ASE LAMMPS calculator interface. 

208 

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

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

211 potential will be determined from the file suffix. 

212 

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

214 

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

216 forces. 

217 

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

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

220 neighbor such as caused by evaluating a dimer. 

221 

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

223 

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

225 1285. 

226 

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

228 Acta Materialia 53 2005 4029--4041. 

229 

230 

231End EAM Interface Documentation 

232 """ 

233 

234 implemented_properties = ['energy', 'forces'] 

235 

236 default_parameters = dict( 

237 skin=1.0, 

238 potential=None, 

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

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

241 b'blank\n']) 

242 

243 def __init__(self, restart=None, 

244 ignore_bad_restart_file=Calculator._deprecated, 

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

246 

247 self.form = form 

248 

249 if 'potential' in kwargs: 

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

251 

252 Calculator.__init__(self, restart, ignore_bad_restart_file, 

253 label, atoms, **kwargs) 

254 

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

256 'cutoff', 'atomic_number', 'mass', 'a', 'lattice', 

257 'embedded_energy', 'electron_density', 'phi', 

258 # derivatives 

259 'd_embedded_energy', 'd_electron_density', 'd_phi', 

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

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

262 

263 # set any additional keyword arguments 

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

265 if arg in valid_args: 

266 setattr(self, arg, val) 

267 else: 

268 raise RuntimeError( 

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

270 

271 def set_form(self, name): 

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

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

274 

275 if extension == '.eam': 

276 self.form = 'eam' 

277 elif extension == '.alloy': 

278 self.form = 'alloy' 

279 elif extension == '.adp': 

280 self.form = 'adp' 

281 elif extension == '.fs': 

282 self.form = 'fs' 

283 else: 

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

285 

286 def read_potential(self, filename): 

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

288 and creates the interpolation functions from the data 

289 """ 

290 

291 if isinstance(filename, str): 

292 with open(filename) as fd: 

293 self._read_potential(fd) 

294 else: 

295 fd = filename 

296 self._read_potential(fd) 

297 

298 def _read_potential(self, fd): 

299 if self.form is None: 

300 self.set_form(fd.name) 

301 

302 lines = fd.readlines() 

303 

304 def lines_to_list(lines): 

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

306 """ 

307 data = [] 

308 for line in lines: 

309 data.extend(line.split()) 

310 return data 

311 

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

313 self.header = lines[:1] 

314 

315 data = lines_to_list(lines[1:]) 

316 

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

318 

319 self.Nelements = 1 

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

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

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

323 self.lattice = [data[3]] 

324 

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

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

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

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

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

330 

331 n = 9 + self.nrho 

332 self.embedded_data = np.array([np.float_(data[9:n])]) 

333 

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

335 self.nr]) 

336 

337 effective_charge = np.float_(data[n:n + self.nr]) 

338 # convert effective charges to rphi according to 

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

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

341 

342 self.density_data = np.array( 

343 [np.float_(data[n + self.nr:n + 2 * self.nr])]) 

344 

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

346 self.header = lines[:3] 

347 i = 3 

348 

349 data = lines_to_list(lines[i:]) 

350 

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

352 d = 1 

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

354 d += self.Nelements 

355 

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

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

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

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

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

361 

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

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

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

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

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

367 self.lattice = [] 

368 d += 5 

369 

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

371 for elem in range(self.Nelements): 

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

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

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

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

376 d += 4 

377 

378 self.embedded_data[elem] = np.float_( 

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

380 d += self.nrho 

381 self.density_data[elem] = np.float_(data[d:(d + self.nr)]) 

382 d += self.nr 

383 

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

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

386 self.nr]) 

387 

388 for i in range(self.Nelements): 

389 for j in range(i + 1): 

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

391 d += self.nr 

392 

393 elif self.form == 'fs': 

394 self.header = lines[:3] 

395 i = 3 

396 

397 data = lines_to_list(lines[i:]) 

398 

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

400 d = 1 

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

402 d += self.Nelements 

403 

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

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

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

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

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

409 

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

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

412 self.nr]) 

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

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

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

416 self.lattice = [] 

417 d += 5 

418 

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

420 for elem in range(self.Nelements): 

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

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

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

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

425 d += 4 

426 

427 self.embedded_data[elem] = np.float_( 

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

429 d += self.nrho 

430 self.density_data[elem, :, :] = np.float_( 

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

432 self.Nelements, self.nr]) 

433 d += self.nr * self.Nelements 

434 

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

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

437 self.nr]) 

438 

439 for i in range(self.Nelements): 

440 for j in range(i + 1): 

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

442 d += self.nr 

443 

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

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

446 

447 # choose the set_splines method according to the type 

448 if self.form == 'fs': 

449 self.set_fs_splines() 

450 else: 

451 self.set_splines() 

452 

453 if self.form == 'adp': 

454 self.read_adp_data(data, d) 

455 self.set_adp_splines() 

456 

457 def set_splines(self): 

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

459 # derivative functions) that define the potential 

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

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

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

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

464 

465 for i in range(self.Nelements): 

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

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

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

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

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

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

472 

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

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

475 

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

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

478 for i in range(self.Nelements): 

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

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

481 self.r[1:], 

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

483 

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

485 

486 if j != i: 

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

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

489 

490 def set_fs_splines(self): 

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

492 self.electron_density = np.empty( 

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

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

495 self.d_electron_density = np.empty( 

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

497 

498 for i in range(self.Nelements): 

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

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

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

502 for j in range(self.Nelements): 

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

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

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

506 self.electron_density[i, j]) 

507 

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

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

510 

511 for i in range(self.Nelements): 

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

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

514 self.r[1:], 

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

516 

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

518 

519 if j != i: 

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

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

522 

523 def set_adp_splines(self): 

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

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

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

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

528 

529 for i in range(self.Nelements): 

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

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

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

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

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

535 

536 # make symmetrical 

537 if j != i: 

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

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

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

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

542 

543 def read_adp_data(self, data, d): 

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

545 

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

547 # should be non symmetrical combinations of 2 

548 for i in range(self.Nelements): 

549 for j in range(i + 1): 

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

551 d += self.nr 

552 

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

554 # should be non symmetrical combinations of 2 

555 for i in range(self.Nelements): 

556 for j in range(i + 1): 

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

558 d += self.nr 

559 

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

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

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

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

564 """ 

565 

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

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

568 

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

570 assert self.nr % nc == 0 

571 assert self.nrho % nc == 0 

572 

573 for line in self.header: 

574 fd.write(line) 

575 

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

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

578 

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

580 (self.nrho, self.drho, self.nr, 

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

582 

583 # start of each section for each element 

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

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

586 

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

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

589 

590 for i in range(self.Nelements): 

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

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

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

594 np.savetxt(fd, 

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

596 nc), 

597 fmt=nc * [numformat]) 

598 if self.form == 'fs': 

599 for j in range(self.Nelements): 

600 np.savetxt(fd, 

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

602 self.nr // nc, nc), 

603 fmt=nc * [numformat]) 

604 else: 

605 np.savetxt(fd, 

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

607 self.nr // nc, nc), 

608 fmt=nc * [numformat]) 

609 

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

611 # as r*phi for alloy format 

612 for i in range(self.Nelements): 

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

614 np.savetxt(fd, 

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

616 nc), 

617 fmt=nc * [numformat]) 

618 

619 if self.form == 'adp': 

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

621 for i in range(self.Nelements): 

622 for j in range(i + 1): 

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

624 

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

626 for i in range(self.Nelements): 

627 for j in range(i + 1): 

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

629 

630 def update(self, atoms): 

631 # check all the elements are available in the potential 

632 self.Nelements = len(self.elements) 

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

634 unavailable = np.logical_not( 

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

636 

637 if np.any(unavailable): 

638 raise RuntimeError( 

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

640 

641 # cutoffs need to be a vector for NeighborList 

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

643 

644 # convert the elements to an index of the position 

645 # in the eam format 

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

647 for el in atoms.get_chemical_symbols()]) 

648 self.pbc = atoms.get_pbc() 

649 

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

651 # local electron density we cannot just calculate and use 

652 # one way neighbors 

653 self.neighbors = NeighborList(cutoffs, 

654 skin=self.parameters.skin, 

655 self_interaction=False, 

656 bothways=True) 

657 self.neighbors.update(atoms) 

658 

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

660 system_changes=all_changes): 

661 """EAM Calculator 

662 

663 atoms: Atoms object 

664 Contains positions, unit-cell, ... 

665 properties: list of str 

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

667 of 'energy', 'forces' 

668 system_changes: list of str 

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

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

671 'pbc', 'initial_charges' and 'initial_magmoms'. 

672 """ 

673 

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

675 

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

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

678 self.update(self.atoms) 

679 self.calculate_energy(self.atoms) 

680 

681 if 'forces' in properties: 

682 self.calculate_forces(self.atoms) 

683 

684 # check we have all the properties requested 

685 for property in properties: 

686 if property not in self.results: 

687 if property == 'energy': 

688 self.calculate_energy(self.atoms) 

689 

690 if property == 'forces': 

691 self.calculate_forces(self.atoms) 

692 

693 # we need to remember the previous state of parameters 

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

695# self.read_potential(potential) 

696 

697 def calculate_energy(self, atoms): 

698 """Calculate the energy 

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

700 the embedding energy of each atom into the electron cloud 

701 generated by its neighbors 

702 """ 

703 

704 pair_energy = 0.0 

705 embedding_energy = 0.0 

706 mu_energy = 0.0 

707 lam_energy = 0.0 

708 trace_energy = 0.0 

709 

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

711 if self.form == 'adp': 

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

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

714 

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

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

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

718 

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

720 atoms.positions[i]) 

721 

722 # calculate the distance to the nearest neighbors 

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

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

725 

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

727 for j_index in range(self.Nelements): 

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

729 if not use.any(): 

730 continue 

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

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

733 

734 if self.form == 'fs': 

735 density = np.sum( 

736 self.electron_density[j_index, 

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

738 else: 

739 density = np.sum( 

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

741 self.total_density[i] += density 

742 

743 if self.form == 'adp': 

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

745 r[nearest][use], 

746 rvec[nearest][use], 

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

748 

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

750 r[nearest][use], 

751 rvec[nearest][use], 

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

753 

754 # add in the electron embedding energy 

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

756 self.total_density[i]) 

757 

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

759 

760 if self.form == 'adp': 

761 mu_energy += np.sum(self.mu ** 2) / 2. 

762 lam_energy += np.sum(self.lam ** 2) / 2. 

763 

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

765 trace_energy -= np.sum(self.lam[i].trace() ** 2) / 6.