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 

12import numpy as np 

13 

14from ase.neighborlist import NeighborList 

15from ase.calculators.calculator import Calculator, all_changes 

16from scipy.interpolate import InterpolatedUnivariateSpline as spline 

17from ase.units import Bohr, Hartree 

18 

19 

20class EAM(Calculator): 

21 r""" 

22 

23 EAM Interface Documentation 

24 

25Introduction 

26============ 

27 

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

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

30an equiaxial potential the EAM does not model directional bonds 

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

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

33included in the EAM calculator. 

34 

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

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

37The files containing the potentials for this calculator are not 

38included but many suitable potentials can be downloaded from The 

39Interatomic Potentials Repository Project at 

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

41 

42Theory 

43====== 

44 

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

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

47element alloy contains the individual three functions for each element 

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

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

50each alloy and their cross interactions. 

51 

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

53given by the EAM potential as 

54 

55.. math:: 

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

57 

58and 

59 

60.. math:: 

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

62 

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

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

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

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

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

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

69atom and its neighbour for that bond. 

70 

71The ADP potential is defined as 

72 

73.. math:: 

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

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

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

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

78 

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

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

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

82 

83The fs potential is defined as 

84 

85.. math:: 

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

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

88 

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

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

91by element types. 

92 

93Running the Calculator 

94====================== 

95 

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

97potential functions are defined by splines which may be directly 

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

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

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

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

102supported. 

103 

104For example:: 

105 

106 from ase.calculators.eam import EAM 

107 

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

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

110 mishin.plot() 

111 

112 slab.calc = mishin 

113 slab.get_potential_energy() 

114 slab.get_forces() 

115 

116The breakdown of energy contribution from the indvidual components are 

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

118 

119Arguments 

120========= 

121 

122========================= ==================================================== 

123Keyword Description 

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

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

126 format or file object 

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

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

129 

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

131 

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

133 

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

135 

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

137 

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

139 

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

141 

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

143 

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

145 

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

147 

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

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

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

151 list can be reused. Defaults to 1.0. 

152 

153``form`` the form of the potential 

154 ``eam``, ``alloy``, ``adp`` or 

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

156 or must be set if using equations or file object 

157 

158========================= ==================================================== 

159 

160 

161Additional parameters for writing potential files 

162================================================= 

163 

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

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

166 

167========================= ==================================================== 

168Keyword Description 

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

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

171 

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

173 

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

175 

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

177 

178``lattice[N]`` Lattice type 

179 

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

181 

182``drho`` Increment for sampling density 

183 

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

185 potential curves 

186 

187``dr`` Increment for sampling radius 

188 

189========================= ==================================================== 

190 

191Special features 

192================ 

193 

194``.plot()`` 

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

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

197 function requires the installation of the Matplotlib libraries. 

198 

199Notes/Issues 

200============= 

201 

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

203 small calculations or for creating new potentials by matching baseline 

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

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

206 with the ASE LAMMPS calculator interface. 

207 

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

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

210 potential will be determined from the file suffix. 

211 

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

213 

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

215 forces. 

216 

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

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

219 neighbor such as caused by evaluating a dimer. 

220 

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

222 

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

224 1285. 

225 

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

227 Acta Materialia 53 2005 4029--4041. 

228 

229 

230End EAM Interface Documentation 

231 """ 

232 

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

234 

235 default_parameters = dict( 

236 skin=1.0, 

237 potential=None, 

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

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

240 b'blank\n']) 

241 

242 def __init__(self, restart=None, 

243 ignore_bad_restart_file=Calculator._deprecated, 

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

245 

246 self.form = form 

247 

248 if 'potential' in kwargs: 

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

250 

251 Calculator.__init__(self, restart, ignore_bad_restart_file, 

252 label, atoms, **kwargs) 

253 

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

255 'cutoff', 'atomic_number', 'mass', 'a', 'lattice', 

256 'embedded_energy', 'electron_density', 'phi', 

257 # derivatives 

258 'd_embedded_energy', 'd_electron_density', 'd_phi', 

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

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

261 

262 # set any additional keyword arguments 

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

264 if arg in valid_args: 

265 setattr(self, arg, val) 

266 else: 

267 raise RuntimeError('unknown keyword arg "%s" : not in %s' 

268 % (arg, valid_args)) 

269 

270 def set_form(self, name): 

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

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

273 

274 if extension == '.eam': 

275 self.form = 'eam' 

276 elif extension == '.alloy': 

277 self.form = 'alloy' 

278 elif extension == '.adp': 

279 self.form = 'adp' 

280 elif extension == '.fs': 

281 self.form = 'fs' 

282 else: 

283 raise RuntimeError('unknown file extension type: %s' % extension) 

284 

285 def read_potential(self, filename): 

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

287 and creates the interpolation functions from the data 

288 """ 

289 

290 if isinstance(filename, str): 

291 with open(filename) as fd: 

292 self._read_potential(fd) 

293 else: 

294 fd = filename 

295 self._read_potential(fd) 

296 

297 def _read_potential(self, fd): 

298 if self.form is None: 

299 self.set_form(fd.name) 

300 

301 lines = fd.readlines() 

302 

303 def lines_to_list(lines): 

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

305 """ 

306 data = [] 

307 for line in lines: 

308 data.extend(line.split()) 

309 return data 

310 

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

312 self.header = lines[:1] 

313 

314 data = lines_to_list(lines[1:]) 

315 

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

317 

318 self.Nelements = 1 

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

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

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

322 self.lattice = [data[3]] 

323 

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

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

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

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

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

329 

330 n = 9 + self.nrho 

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

332 

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

334 self.nr]) 

335 

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

337 # convert effective charges to rphi according to 

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

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

340 

341 self.density_data = np.array( 

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

343 

344 elif self.form in ['alloy', 'adq']: 

345 self.header = lines[:3] 

346 i = 3 

347 

348 data = lines_to_list(lines[i:]) 

349 

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

351 d = 1 

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

353 d += self.Nelements 

354 

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

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

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

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

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

360 

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

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

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

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

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

366 self.lattice = [] 

367 d += 5 

368 

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

370 for elem in range(self.Nelements): 

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

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

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

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

375 d += 4 

376 

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

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

379 d += self.nrho 

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

381 d += self.nr 

382 

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

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

385 self.nr]) 

386 

387 for i in range(self.Nelements): 

388 for j in range(i + 1): 

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

390 d += self.nr 

391 

392 elif self.form == 'fs': 

393 self.header = lines[:3] 

394 i = 3 

395 

396 data = lines_to_list(lines[i:]) 

397 

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

399 d = 1 

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

401 d += self.Nelements 

402 

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

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

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

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

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

408 

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

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

411 self.nr]) 

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

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

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

415 self.lattice = [] 

416 d += 5 

417 

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

419 for elem in range(self.Nelements): 

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

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

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

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

424 d += 4 

425 

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

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

428 d += self.nrho 

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

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

431 self.Nelements, self.nr]) 

432 d += self.nr*self.Nelements 

433 

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

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

436 self.nr]) 

437 

438 for i in range(self.Nelements): 

439 for j in range(i + 1): 

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

441 d += self.nr 

442 

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

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

445 

446 # choose the set_splines method according to the type 

447 if self.form == 'fs': 

448 self.set_fs_splines() 

449 else: 

450 self.set_splines() 

451 

452 if (self.form == 'adp'): 

453 self.read_adp_data(data, d) 

454 self.set_adp_splines() 

455 

456 def set_splines(self): 

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

458 # derivative functions) that define the potential 

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

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

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

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

463 

464 for i in range(self.Nelements): 

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

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

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

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

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

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

471 

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

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

474 

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

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

477 for i in range(self.Nelements): 

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

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

480 self.r[1:], 

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

482 

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

484 

485 if j != i: 

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

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

488 

489 def set_fs_splines(self): 

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

491 self.electron_density = np.empty( 

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

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

494 self.d_electron_density = np.empty( 

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

496 

497 for i in range(self.Nelements): 

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

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

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

501 for j in range(self.Nelements): 

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

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

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

505 self.electron_density[i, j]) 

506 

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

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

509 

510 for i in range(self.Nelements): 

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

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

513 self.r[1:], 

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

515 

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

517 

518 if j != i: 

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

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

521 

522 def set_adp_splines(self): 

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

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

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

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

527 

528 for i in range(self.Nelements): 

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

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

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

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

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

534 

535 # make symmetrical 

536 if j != i: 

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

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

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

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

541 

542 def read_adp_data(self, data, d): 

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

544 

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

546 # should be non symmetrical combinations of 2 

547 for i in range(self.Nelements): 

548 for j in range(i + 1): 

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

550 d += self.nr 

551 

552 self.q_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.q_data[j, i] = data[d:d + self.nr] 

557 d += self.nr 

558 

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

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

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

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

563 """ 

564 

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

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

567 

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

569 assert self.nr % nc == 0 

570 assert self.nrho % nc == 0 

571 

572 for line in self.header: 

573 fd.write(line) 

574 

575 fd.write('{0} '.format(self.Nelements).encode()) 

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

577 

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

579 (self.nrho, self.drho, self.nr, 

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

581 

582 # start of each section for each element 

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

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

585 

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

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

588 

589 for i in range(self.Nelements): 

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

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

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

593 np.savetxt(fd, 

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

595 nc), 

596 fmt=nc * [numformat]) 

597 if self.form == 'fs': 

598 for j in range(self.Nelements): 

599 np.savetxt(fd, 

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

601 self.nr // nc, nc), 

602 fmt=nc * [numformat]) 

603 else: 

604 np.savetxt(fd, 

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

606 self.nr // nc, nc), 

607 fmt=nc * [numformat]) 

608 

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

610 # as r*phi for alloy format 

611 for i in range(self.Nelements): 

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

613 np.savetxt(fd, 

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

615 nc), 

616 fmt=nc * [numformat]) 

617 

618 if self.form == 'adp': 

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

620 for i in range(self.Nelements): 

621 for j in range(i + 1): 

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

623 

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

625 for i in range(self.Nelements): 

626 for j in range(i + 1): 

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

628 

629 def update(self, atoms): 

630 # check all the elements are available in the potential 

631 self.Nelements = len(self.elements) 

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

633 unavailable = np.logical_not( 

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

635 

636 if np.any(unavailable): 

637 raise RuntimeError('These elements are not in the potential: %s' % 

638 elements[unavailable]) 

639 

640 # cutoffs need to be a vector for NeighborList 

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

642 

643 # convert the elements to an index of the position 

644 # in the eam format 

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

646 for el in atoms.get_chemical_symbols()]) 

647 self.pbc = atoms.get_pbc() 

648 

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

650 # local electron density we cannot just calculate and use 

651 # one way neighbors 

652 self.neighbors = NeighborList(cutoffs, 

653 skin=self.parameters.skin, 

654 self_interaction=False, 

655 bothways=True) 

656 self.neighbors.update(atoms) 

657 

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

659 system_changes=all_changes): 

660 """EAM Calculator 

661 

662 atoms: Atoms object 

663 Contains positions, unit-cell, ... 

664 properties: list of str 

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

666 of 'energy', 'forces' 

667 system_changes: list of str 

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

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

670 'pbc', 'initial_charges' and 'initial_magmoms'. 

671 """ 

672 

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

674 

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

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

677 self.update(self.atoms) 

678 self.calculate_energy(self.atoms) 

679 

680 if 'forces' in properties: 

681 self.calculate_forces(self.atoms) 

682 

683 # check we have all the properties requested 

684 for property in properties: 

685 if property not in self.results: 

686 if property == 'energy': 

687 self.calculate_energy(self.atoms) 

688 

689 if property == 'forces': 

690 self.calculate_forces(self.atoms) 

691 

692 # we need to remember the previous state of parameters 

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

694# self.read_potential(potential) 

695 

696 def calculate_energy(self, atoms): 

697 """Calculate the energy 

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

699 the embedding energy of each atom into the electron cloud 

700 generated by its neighbors 

701 """ 

702 

703 pair_energy = 0.0 

704 embedding_energy = 0.0 

705 mu_energy = 0.0 

706 lam_energy = 0.0 

707 trace_energy = 0.0 

708 

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

710 if (self.form == 'adp'): 

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

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

713 

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

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

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

717 

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

719 atoms.positions[i]) 

720 

721 # calculate the distance to the nearest neighbors 

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

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

724 

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

726 for j_index in range(self.Nelements): 

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

728 if not use.any(): 

729 continue 

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

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

732 

733 if self.form == 'fs': 

734 density = np.sum( 

735 self.electron_density[j_index, 

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

737 else: 

738 density = np.sum( 

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

740 self.total_density[i] += density 

741 

742 if self.form == 'adp': 

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

744 r[nearest][use], 

745 rvec[nearest][use], 

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

747 

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

749 r[nearest][use], 

750 rvec[nearest][use], 

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

752 

753 # add in the electron embedding energy 

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

755 self.total_density[i]) 

756 

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

758 

759 if self.form == 'adp': 

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

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

762 

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

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

765 

766 adp_result = dict(adp_mu=mu_energy, 

767 adp_lam=lam_energy, 

768 adp_trace=trace_energy) 

769 components.update(adp_result) 

770 

771 self.positions = atoms.positions.copy() 

772 self.cell = atoms.get_cell().copy() 

773 

774 energy = 0.0 

775 for i in components.keys(): 

776 energy += components[i] 

777 

778 self.energy_free = energy 

779 self.energy_zero = energy 

780 

781 self.results['energy_components'] = components 

782 self.results['energy'] = energy 

783 

784 def calculate_forces(self, atoms): 

785 # calculate the forces based on derivatives of the three EAM functions 

786 

787 self.update(atoms) 

788 self.results['forces'] = np.zeros((len(atoms), 3)) 

789 

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

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

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

793 # create a vector of relative positions of neighbors 

794 rvec = atoms.positions[neighbors] + offset - atoms.positions[i] 

795 r = np.sqrt(np.sum(np.square(rvec), axis=1)) 

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

797 

798 d_embedded_energy_i = self.d_embedded_energy[ 

799 self.index[i]](self.total_density[i]) 

800 urvec = rvec.copy() # unit directional vector 

801 

802 for j in np.arange(len(neighbors)): 

803 urvec[j] = urvec[j] / r[j] 

804 

805 for j_index in range(self.Nelements): 

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

807 if not use.any(): 

808 continue 

809 rnuse = r[nearest][use] 

810 density_j = self.total_density[neighbors[nearest][use]] 

811 if self.form == 'fs': 

812 scale = (self.d_phi[self.index[i], j_index](rnuse) + 

813 (d_embedded_energy_i * 

814 self.d_electron_density[j_index, 

815 self.index[i]](rnuse)) + 

816 (self.d_embedded_energy[j_index](density_j) * 

817 self.d_electron_density[self.index[i], 

818 j_index](rnuse))) 

819 else: 

820 scale = (self.d_phi[self.index[i], j_index](rnuse) + 

821 (d_embedded_energy_i * 

822 self.d_electron_density[j_index](rnuse)) + 

823 (self.d_embedded_energy[j_index](density_j) * 

824 self.d_electron_density[self.index[i]](rnuse))) 

825 

826 self.results['forces'][i] += np.dot(scale, urvec[nearest][use]) 

827 

828 if (self.form == 'adp'): 

829 adp_forces = self.angular_forces( 

830 self.mu[i], 

831 self.mu[neighbors[nearest][use]], 

832 self.lam[i], 

833 self.lam[neighbors[nearest][use]], 

834 rnuse, 

835 rvec[nearest][use], 

836 self.index[i], 

837 j_index) 

838 

839 self.results['forces'][i] += adp_forces 

840 

841 def angular_forces(self, mu_i, mu, lam_i, lam, r, rvec, form1, form2): 

842 # calculate the extra components for the adp forces 

843 # rvec are the relative positions to atom i 

844 psi = np.zeros(mu.shape) 

845 for gamma in range(3): 

846 term1 = (mu_i[gamma] - mu[:, gamma]) * self.d[form1][form2](r) 

847 

848 term2 = np.sum((mu_i - mu) * 

849 self.d_d[form1][form2](r)[:, np.newaxis] * 

850 (rvec * rvec[:, gamma][:, np.newaxis] / 

851 r[:, np.newaxis]), axis=1) 

852 

853 term3 = 2 * np.sum((lam_i[:, gamma] + lam[:, :, gamma]) * 

854 rvec * self.q[form1][form2](r)[:, np.newaxis], 

855 axis=1) 

856 term4 = 0.0 

857 for alpha in range(3): 

858 for beta in range(3): 

859 rs = rvec[:, alpha] * rvec[:, beta] * rvec[:, gamma] 

860 term4 += ((lam_i[alpha, beta] + lam[:, alpha, beta]) * 

861 self.d_q[form1][form2](r) * rs) / r 

862 

863 term5 = ((lam_i.trace() + lam.trace(axis1=1, axis2=2)) * 

864 (self.d_q[form1][form2](r) * r + 

865 2 * self.q[form1][form2](r)) * rvec[:, gamma]) / 3. 

866 

867 # the minus for term5 is a correction on the adp 

868 # formulation given in the 2005 Mishin Paper and is posted 

869 # on the NIST website with the AlH potential 

870 psi[:, gamma] = term1 + term2 + term3 + term4 - term5 

871 

872 return np.sum(psi, axis=0) 

873 

874 def adp_dipole(self, r, rvec, d): 

875 # calculate the dipole contribution 

876 mu = np.sum((rvec * d(r)[:, np.newaxis]), axis=0) 

877 

878 return mu # sign to agree with lammps 

879 

880 def adp_quadrupole(self, r, rvec, q): 

881 # slow way of calculating the quadrupole contribution 

882 r = np.sqrt(np.sum(rvec ** 2, axis=1)) 

883 

884 lam = np.zeros([rvec.shape[0], 3, 3]) 

885 qr = q(r) 

886 for alpha in range(3): 

887 for beta in range(3): 

888 lam[:, alpha, beta] += qr * rvec[:, alpha] * rvec[:, beta] 

889 

890 return np.sum(lam, axis=0) 

891 

892 def deriv(self, spline): 

893 """Wrapper for extracting the derivative from a spline""" 

894 def d_spline(aspline): 

895 return spline(aspline, 1) 

896 

897 return d_spline 

898 

899 def plot(self, name=''): 

900 """Plot the individual curves""" 

901 

902 import matplotlib.pyplot as plt 

903 

904 if self.form == 'eam' or self.form == 'alloy' or self.form == 'fs': 

905 nrow = 2 

906 elif self.form == 'adp': 

907 nrow = 3 

908 else: 

909 raise RuntimeError('Unknown form of potential: %s' % self.form) 

910 

911 if hasattr(self, 'r'): 

912 r = self.r 

913 else: 

914 r = np.linspace(0, self.cutoff, 50) 

915 

916 if hasattr(self, 'rho'): 

917 rho = self.rho 

918 else: 

919 rho = np.linspace(0, 10.0, 50) 

920 

921 plt.subplot(nrow, 2, 1) 

922 self.elem_subplot(rho, self.embedded_energy, 

923 r'$\rho$', r'Embedding Energy $F(\bar\rho)$', 

924 name, plt) 

925 

926 plt.subplot(nrow, 2, 2) 

927 if self.form == 'fs': 

928 self.multielem_subplot( 

929 r, self.electron_density, 

930 r'$r$', r'Electron Density $\rho(r)$', name, plt, half=False) 

931 else: 

932 self.elem_subplot( 

933 r, self.electron_density, 

934 r'$r$', r'Electron Density $\rho(r)$', name, plt) 

935 

936 plt.subplot(nrow, 2, 3) 

937 self.multielem_subplot(r, self.phi, 

938 r'$r$', r'Pair Potential $\phi(r)$', name, plt) 

939 plt.ylim(-1.0, 1.0) # need reasonable values 

940 

941 if self.form == 'adp': 

942 plt.subplot(nrow, 2, 5) 

943 self.multielem_subplot(r, self.d, 

944 r'$r$', r'Dipole Energy', name, plt) 

945 

946 plt.subplot(nrow, 2, 6) 

947 self.multielem_subplot(r, self.q, 

948 r'$r$', r'Quadrupole Energy', name, plt) 

949 

950 plt.plot() 

951 

952 def elem_subplot(self, curvex, curvey, xlabel, ylabel, name, plt): 

953 plt.xlabel(xlabel) 

954 plt.ylabel(ylabel) 

955 for i in np.arange(self.Nelements): 

956 label = name + ' ' + self.elements[i] 

957 plt.plot(curvex, curvey[i](curvex), label=label) 

958 plt.legend() 

959 

960 def multielem_subplot(self, curvex, curvey, xlabel, 

961 ylabel, name, plt, half=True): 

962 plt.xlabel(xlabel) 

963 plt.ylabel(ylabel) 

964 for i in np.arange(self.Nelements): 

965 for j in np.arange((i + 1) if half else self.Nelements): 

966 label = name + ' ' + self.elements[i] + '-' + self.elements[j] 

967 plt.plot(curvex, curvey[i, j](curvex), label=label) 

968 plt.legend()