Coverage for /builds/debichem-team/python-ase/ase/calculators/eam.py: 82.84%

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.data import chemical_symbols 

18from ase.neighborlist import NeighborList 

19from ase.units import Bohr, Hartree 

20 

21 

22class EAM(Calculator): 

23 r""" 

24 

25 EAM Interface Documentation 

26 

27Introduction 

28============ 

29 

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

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

32an equiaxial potential the EAM does not model directional bonds 

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

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

35included in the EAM calculator. 

36 

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

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

39The files containing the potentials for this calculator are not 

40included but many suitable potentials can be downloaded from The 

41Interatomic Potentials Repository Project at 

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

43 

44Theory 

45====== 

46 

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

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

49element alloy contains the individual three functions for each element 

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

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

52each alloy and their cross interactions. 

53 

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

55given by the EAM potential as 

56 

57.. math:: 

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

59 

60and 

61 

62.. math:: 

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

64 

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

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

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

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

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

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

71atom and its neighbour for that bond. 

72 

73The ADP potential is defined as 

74 

75.. math:: 

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

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

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

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

80 

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

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

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

84 

85The fs potential is defined as 

86 

87.. math:: 

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

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

90 

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

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

93by element types. 

94 

95Running the Calculator 

96====================== 

97 

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

99potential functions are defined by splines which may be directly 

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

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

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

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

104supported. 

105 

106For example:: 

107 

108 from ase.calculators.eam import EAM 

109 

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

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

112 mishin.plot() 

113 

114 slab.calc = mishin 

115 slab.get_potential_energy() 

116 slab.get_forces() 

117 

118The breakdown of energy contribution from the indvidual components are 

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

120 

121Arguments 

122========= 

123 

124========================= ==================================================== 

125Keyword Description 

126========================= ==================================================== 

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

128 format or file object 

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

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

131 

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

133 

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

135 

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

137 

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

139 

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

141 

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

143 

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

145 

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

147 

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

149 

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

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

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

153 list can be reused. Defaults to 1.0. 

154 

155``form`` the form of the potential 

156 ``eam``, ``alloy``, ``adp`` or 

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

158 or must be set if using equations or file object 

159 

160========================= ==================================================== 

161 

162 

163Additional parameters for writing potential files 

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

165 

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

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

168 

169========================= ==================================================== 

170Keyword Description 

171========================= ==================================================== 

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

173 

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

175 

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

177 

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

179 

180``lattice[N]`` Lattice type 

181 

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

183 

184``drho`` Increment for sampling density 

185 

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

187 potential curves 

188 

189``dr`` Increment for sampling radius 

190 

191========================= ==================================================== 

192 

193Special features 

194================ 

195 

196``.plot()`` 

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

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

199 function requires the installation of the Matplotlib libraries. 

200 

201Notes/Issues 

202============= 

203 

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

205 small calculations or for creating new potentials by matching baseline 

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

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

208 with the ASE LAMMPS calculator interface. 

209 

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

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

212 potential will be determined from the file suffix. 

213 

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

215 

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

217 forces. 

218 

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

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

221 neighbor such as caused by evaluating a dimer. 

222 

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

224 

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

226 1285. 

227 

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

229 Acta Materialia 53 2005 4029--4041. 

230 

231 

232End EAM Interface Documentation 

233 """ 

234 

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

236 

237 default_parameters = dict( 

238 skin=1.0, 

239 potential=None, 

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

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

242 b'blank\n']) 

243 

244 def __init__(self, restart=None, 

245 ignore_bad_restart_file=Calculator._deprecated, 

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

247 

248 self.form = form 

249 

250 if 'potential' in kwargs: 

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

252 

253 Calculator.__init__(self, restart, ignore_bad_restart_file, 

254 label, atoms, **kwargs) 

255 

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

257 'cutoff', 'atomic_number', 'mass', 'a', 'lattice', 

258 'embedded_energy', 'electron_density', 'phi', 

259 # derivatives 

260 'd_embedded_energy', 'd_electron_density', 'd_phi', 

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

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

263 

264 # set any additional keyword arguments 

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

266 if arg in valid_args: 

267 setattr(self, arg, val) 

268 else: 

269 raise RuntimeError( 

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

271 

272 def set_form(self, name): 

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

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

275 

276 if extension == '.eam': 

277 self.form = 'eam' 

278 elif extension == '.alloy': 

279 self.form = 'alloy' 

280 elif extension == '.adp': 

281 self.form = 'adp' 

282 elif extension == '.fs': 

283 self.form = 'fs' 

284 else: 

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

286 

287 def read_potential(self, filename): 

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

289 and creates the interpolation functions from the data 

290 """ 

291 

292 if isinstance(filename, str): 

293 with open(filename) as fd: 

294 self._read_potential(fd) 

295 else: 

296 fd = filename 

297 self._read_potential(fd) 

298 

299 def _read_potential(self, fd): 

300 if self.form is None: 

301 self.set_form(fd.name) 

302 

303 lines = fd.readlines() 

304 

305 def lines_to_list(lines): 

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

307 """ 

308 data = [] 

309 for line in lines: 

310 data.extend(line.split()) 

311 return data 

312 

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

314 self.header = lines[:1] 

315 

316 data = lines_to_list(lines[1:]) 

317 

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

319 

320 self.Nelements = 1 

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

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

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

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

325 self.lattice = [data[3]] 

326 

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

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

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

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

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

332 

333 n = 9 + self.nrho 

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

335 

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

337 self.nr]) 

338 

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

340 # convert effective charges to rphi according to 

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

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

343 

344 self.density_data = np.array( 

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

346 

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

348 self.header = lines[:3] 

349 i = 3 

350 

351 data = lines_to_list(lines[i:]) 

352 

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

354 d = 1 

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

356 d += self.Nelements 

357 

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

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

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

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

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

363 

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

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

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

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

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

369 self.lattice = [] 

370 d += 5 

371 

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

373 for elem in range(self.Nelements): 

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

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

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

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

378 d += 4 

379 

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

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

382 d += self.nrho 

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

384 d += self.nr 

385 

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

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

388 self.nr]) 

389 

390 for i in range(self.Nelements): 

391 for j in range(i + 1): 

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

393 d += self.nr 

394 

395 elif self.form == 'fs': 

396 self.header = lines[:3] 

397 i = 3 

398 

399 data = lines_to_list(lines[i:]) 

400 

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

402 d = 1 

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

404 d += self.Nelements 

405 

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

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

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

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

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

411 

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

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

414 self.nr]) 

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

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

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

418 self.lattice = [] 

419 d += 5 

420 

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

422 for elem in range(self.Nelements): 

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

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

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

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

427 d += 4 

428 

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

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

431 d += self.nrho 

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

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

434 self.Nelements, self.nr]) 

435 d += self.nr * self.Nelements 

436 

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

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

439 self.nr]) 

440 

441 for i in range(self.Nelements): 

442 for j in range(i + 1): 

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

444 d += self.nr 

445 

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

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

448 

449 # choose the set_splines method according to the type 

450 if self.form == 'fs': 

451 self.set_fs_splines() 

452 else: 

453 self.set_splines() 

454 

455 if self.form == 'adp': 

456 self.read_adp_data(data, d) 

457 self.set_adp_splines() 

458 

459 def set_splines(self): 

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

461 # derivative functions) that define the potential 

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

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

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

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

466 

467 for i in range(self.Nelements): 

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

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

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

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

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

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

474 

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

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

477 

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

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

480 for i in range(self.Nelements): 

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

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

483 self.r[1:], 

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

485 

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

487 

488 if j != i: 

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

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

491 

492 def set_fs_splines(self): 

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

494 self.electron_density = np.empty( 

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

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

497 self.d_electron_density = np.empty( 

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

499 

500 for i in range(self.Nelements): 

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

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

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

504 for j in range(self.Nelements): 

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

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

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

508 self.electron_density[i, j]) 

509 

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

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

512 

513 for i in range(self.Nelements): 

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

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

516 self.r[1:], 

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

518 

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

520 

521 if j != i: 

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

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

524 

525 def set_adp_splines(self): 

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

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

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

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

530 

531 for i in range(self.Nelements): 

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

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

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

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

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

537 

538 # make symmetrical 

539 if j != i: 

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

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

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

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

544 

545 def read_adp_data(self, data, d): 

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

547 

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

549 # should be non symmetrical combinations of 2 

550 for i in range(self.Nelements): 

551 for j in range(i + 1): 

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

553 d += self.nr 

554 

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

556 # should be non symmetrical combinations of 2 

557 for i in range(self.Nelements): 

558 for j in range(i + 1): 

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

560 d += self.nr 

561 

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

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

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

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

566 """ 

567 

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

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

570 

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

572 assert self.nr % nc == 0 

573 assert self.nrho % nc == 0 

574 

575 for line in self.header: 

576 fd.write(line) 

577 

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

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

580 

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

582 (self.nrho, self.drho, self.nr, 

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

584 

585 # start of each section for each element 

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

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

588 

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

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

591 

592 for i in range(self.Nelements): 

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

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

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

596 np.savetxt(fd, 

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

598 nc), 

599 fmt=nc * [numformat]) 

600 if self.form == 'fs': 

601 for j in range(self.Nelements): 

602 np.savetxt(fd, 

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

604 self.nr // nc, nc), 

605 fmt=nc * [numformat]) 

606 else: 

607 np.savetxt(fd, 

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

609 self.nr // nc, nc), 

610 fmt=nc * [numformat]) 

611 

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

613 # as r*phi for alloy format 

614 for i in range(self.Nelements): 

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

616 np.savetxt(fd, 

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

618 nc), 

619 fmt=nc * [numformat]) 

620 

621 if self.form == 'adp': 

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

623 for i in range(self.Nelements): 

624 for j in range(i + 1): 

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

626 

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

628 for i in range(self.Nelements): 

629 for j in range(i + 1): 

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

631 

632 def update(self, atoms): 

633 # check all the elements are available in the potential 

634 self.Nelements = len(self.elements) 

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

636 unavailable = np.logical_not( 

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

638 

639 if np.any(unavailable): 

640 raise RuntimeError( 

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

642 

643 # cutoffs need to be a vector for NeighborList 

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

645 

646 # convert the elements to an index of the position 

647 # in the eam format 

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

649 for el in atoms.get_chemical_symbols()]) 

650 self.pbc = atoms.get_pbc() 

651 

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

653 # local electron density we cannot just calculate and use 

654 # one way neighbors 

655 self.neighbors = NeighborList(cutoffs, 

656 skin=self.parameters.skin, 

657 self_interaction=False, 

658 bothways=True) 

659 self.neighbors.update(atoms) 

660 

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

662 system_changes=all_changes): 

663 """EAM Calculator 

664 

665 atoms: Atoms object 

666 Contains positions, unit-cell, ... 

667 properties: list of str 

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

669 of 'energy', 'forces' 

670 system_changes: list of str 

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

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

673 'pbc', 'initial_charges' and 'initial_magmoms'. 

674 """ 

675 

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

677 

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

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

680 self.update(self.atoms) 

681 self.calculate_energy(self.atoms) 

682 

683 if 'forces' in properties: 

684 self.calculate_forces(self.atoms) 

685 

686 # check we have all the properties requested 

687 for property in properties: 

688 if property not in self.results: 

689 if property == 'energy': 

690 self.calculate_energy(self.atoms) 

691 

692 if property == 'forces': 

693 self.calculate_forces(self.atoms) 

694 

695 # we need to remember the previous state of parameters 

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

697# self.read_potential(potential) 

698 

699 def calculate_energy(self, atoms): 

700 """Calculate the energy 

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

702 the embedding energy of each atom into the electron cloud 

703 generated by its neighbors 

704 """ 

705 

706 pair_energy = 0.0 

707 embedding_energy = 0.0 

708 mu_energy = 0.0 

709 lam_energy = 0.0 

710 trace_energy = 0.0 

711 

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

713 if self.form == 'adp': 

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

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

716 

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

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

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

720 

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

722 atoms.positions[i]) 

723 

724 # calculate the distance to the nearest neighbors 

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

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

727 

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

729 for j_index in range(self.Nelements): 

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

731 if not use.any(): 

732 continue 

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

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

735 

736 if self.form == 'fs': 

737 density = np.sum( 

738 self.electron_density[j_index, 

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

740 else: 

741 density = np.sum( 

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

743 self.total_density[i] += density 

744 

745 if self.form == 'adp': 

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

747 r[nearest][use], 

748 rvec[nearest][use], 

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

750 

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

752 r[nearest][use], 

753 rvec[nearest][use], 

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

755 

756 # add in the electron embedding energy 

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

758 self.total_density[i]) 

759 

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

761 

762 if self.form == 'adp': 

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

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