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

2 

3# eam.py 

4# Embedded Atom Method Potential 

5# These routines integrate with the ASE simulation environment 

6# Paul White (Oct 2012) 

7# UNCLASSIFIED 

8# License: See accompanying license files for details 

9 

10import os 

11import numpy as np 

12 

13from ase.neighborlist import NeighborList 

14from ase.calculators.calculator import Calculator, all_changes 

15from scipy.interpolate import InterpolatedUnivariateSpline as spline 

16from ase.units import Bohr, Hartree 

17 

18 

19class EAM(Calculator): 

20 r""" 

21 

22 EAM Interface Documentation 

23 

24Introduction 

25============ 

26 

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

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

29an equiaxial potential the EAM does not model directional bonds 

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

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

32included in the EAM calculator. 

33 

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

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

36The files containing the potentials for this calculator are not 

37included but many suitable potentials can be downloaded from The 

38Interatomic Potentials Repository Project at 

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

40 

41Theory 

42====== 

43 

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

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

46element alloy contains the individual three functions for each element 

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

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

49each alloy and their cross interactions. 

50 

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

52given by the EAM potential as 

53 

54.. math:: 

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

56 

57and 

58 

59.. math:: 

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

61 

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

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

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

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

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

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

68atom and its neighbour for that bond. 

69 

70The ADP potential is defined as 

71 

72.. math:: 

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

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

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

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

77 

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

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

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

81 

82The fs potential is defined as 

83 

84.. math:: 

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

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

87 

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

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

90by element types. 

91 

92Running the Calculator 

93====================== 

94 

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

96potential functions are defined by splines which may be directly 

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

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

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

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

101supported. 

102 

103For example:: 

104 

105 from ase.calculators.eam import EAM 

106 

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

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

109 mishin.plot() 

110 

111 slab.calc = mishin 

112 slab.get_potential_energy() 

113 slab.get_forces() 

114 

115The breakdown of energy contribution from the indvidual components are 

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

117 

118Arguments 

119========= 

120 

121========================= ==================================================== 

122Keyword Description 

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

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

125 format or file object (This is generally all you need to supply). 

126 In case of file object the ``form`` argument is required 

127 

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

129 

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

131 

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

133 

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

135 

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

137 

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

139 

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

141 

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

143 

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

145 

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

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

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

149 list can be reused. Defaults to 1.0. 

150 

151``form`` the form of the potential ``eam``, ``alloy``, ``adp`` or 

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

153 or must be set if using equations or file object 

154 

155========================= ==================================================== 

156 

157 

158Additional parameters for writing potential files 

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

160 

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

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

163 

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

165Keyword Description 

166========================= ==================================================== 

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

168 

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

170 

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

172 

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

174 

175``lattice[N]`` Lattice type 

176 

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

178 

179``drho`` Increment for sampling density 

180 

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

182 potential curves 

183 

184``dr`` Increment for sampling radius 

185 

186========================= ==================================================== 

187 

188Special features 

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

190 

191``.plot()`` 

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

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

194 function requires the installation of the Matplotlib libraries. 

195 

196Notes/Issues 

197============= 

198 

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

200 small calculations or for creating new potentials by matching baseline 

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

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

203 with the ASE LAMMPS calculator interface. 

204 

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

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

207 potential will be determined from the file suffix. 

208 

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

210 

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

212 forces. 

213 

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

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

216 neighbor such as caused by evaluating a dimer. 

217 

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

219 

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

221 1285. 

222 

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

224 Acta Materialia 53 2005 4029--4041. 

225 

226 

227End EAM Interface Documentation 

228 """ 

229 

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

231 

232 default_parameters = dict( 

233 skin=1.0, 

234 potential=None, 

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

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

237 b'blank\n']) 

238 

239 def __init__(self, restart=None, 

240 ignore_bad_restart_file=Calculator._deprecated, 

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

242 

243 self.form = form 

244 

245 if 'potential' in kwargs: 

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

247 

248 Calculator.__init__(self, restart, ignore_bad_restart_file, 

249 label, atoms, **kwargs) 

250 

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

252 'cutoff', 'atomic_number', 'mass', 'a', 'lattice', 

253 'embedded_energy', 'electron_density', 'phi', 

254 # derivatives 

255 'd_embedded_energy', 'd_electron_density', 'd_phi', 

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

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

258 

259 # set any additional keyword arguments 

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

261 if arg in valid_args: 

262 setattr(self, arg, val) 

263 else: 

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

265 % (arg, valid_args)) 

266 

267 def set_form(self, name): 

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

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

270 

271 if extension == '.eam': 

272 self.form = 'eam' 

273 elif extension == '.alloy': 

274 self.form = 'alloy' 

275 elif extension == '.adp': 

276 self.form = 'adp' 

277 elif extension == '.fs': 

278 self.form = 'fs' 

279 else: 

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

281 

282 def read_potential(self, filename): 

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

284 and creates the interpolation functions from the data 

285 """ 

286 

287 if isinstance(filename, str): 

288 with open(filename) as fd: 

289 self._read_potential(fd) 

290 else: 

291 fd = filename 

292 self._read_potential(fd) 

293 

294 def _read_potential(self, fd): 

295 if self.form is None: 

296 self.set_form(fd.name) 

297 

298 lines = fd.readlines() 

299 

300 def lines_to_list(lines): 

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

302 """ 

303 data = [] 

304 for line in lines: 

305 data.extend(line.split()) 

306 return data 

307 

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

309 self.header = lines[:1] 

310 

311 data = lines_to_list(lines[1:]) 

312 

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

314 

315 self.Nelements = 1 

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

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

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

319 self.lattice = [data[3]] 

320 

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

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

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

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

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

326 

327 n = 9 + self.nrho 

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

329 

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

331 self.nr]) 

332 

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

334 # convert effective charges to rphi according to 

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

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

337 

338 self.density_data = np.array( 

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

340 

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

342 self.header = lines[:3] 

343 i = 3 

344 

345 data = lines_to_list(lines[i:]) 

346 

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

348 d = 1 

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

350 d += self.Nelements 

351 

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

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

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

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

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

357 

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

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

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

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

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

363 self.lattice = [] 

364 d += 5 

365 

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

367 for elem in range(self.Nelements): 

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

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

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

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

372 d += 4 

373 

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

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

376 d += self.nrho 

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

378 d += self.nr 

379 

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

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

382 self.nr]) 

383 

384 for i in range(self.Nelements): 

385 for j in range(i + 1): 

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

387 d += self.nr 

388 

389 elif self.form == 'fs': 

390 self.header = lines[:3] 

391 i = 3 

392 

393 data = lines_to_list(lines[i:]) 

394 

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

396 d = 1 

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

398 d += self.Nelements 

399 

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

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

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

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

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

405 

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

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

408 self.nr]) 

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

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

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

412 self.lattice = [] 

413 d += 5 

414 

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

416 for elem in range(self.Nelements): 

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

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

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

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

421 d += 4 

422 

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

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

425 d += self.nrho 

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

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

428 d += self.nr*self.Nelements 

429 

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

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

432 self.nr]) 

433 

434 for i in range(self.Nelements): 

435 for j in range(i + 1): 

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

437 d += self.nr 

438 

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

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

441 

442 # choose the set_splines method according to the type 

443 if self.form == 'fs': 

444 self.set_fs_splines() 

445 else: 

446 self.set_splines() 

447 

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

449 self.read_adp_data(data, d) 

450 self.set_adp_splines() 

451 

452 def set_splines(self): 

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

454 # derivative functions) that define the potential 

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

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

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

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

459 

460 for i in range(self.Nelements): 

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

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

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

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

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

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

467 

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

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

470 

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

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

473 for i in range(self.Nelements): 

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

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

476 self.r[1:], 

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

478 

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

480 

481 if j != i: 

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

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

484 

485 def set_fs_splines(self): 

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

487 self.electron_density = np.empty( 

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

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

490 self.d_electron_density = np.empty( 

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

492 

493 for i in range(self.Nelements): 

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

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

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

497 for j in range(self.Nelements): 

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

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

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

501 self.electron_density[i, j]) 

502 

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

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

505 

506 for i in range(self.Nelements): 

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

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

509 self.r[1:], 

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

511 

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

513 

514 if j != i: 

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

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

517 

518 def set_adp_splines(self): 

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

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

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

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

523 

524 for i in range(self.Nelements): 

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

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

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

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

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

530 

531 # make symmetrical 

532 if j != i: 

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

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

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

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

537 

538 def read_adp_data(self, data, d): 

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

540 

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

542 # should be non symmetrical combinations of 2 

543 for i in range(self.Nelements): 

544 for j in range(i + 1): 

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

546 d += self.nr 

547 

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

553 d += self.nr 

554 

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

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

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

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

559 """ 

560 

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

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

563 

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

565 assert self.nr % nc == 0 

566 assert self.nrho % nc == 0 

567 

568 for line in self.header: 

569 fd.write(line) 

570 

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

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

573 

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

575 (self.nrho, self.drho, self.nr, 

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

577 

578 # start of each section for each element 

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

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

581 

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

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

584 

585 for i in range(self.Nelements): 

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

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

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

589 np.savetxt(fd, 

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

591 nc), 

592 fmt=nc * [numformat]) 

593 if self.form == 'fs': 

594 for j in range(self.Nelements): 

595 np.savetxt(fd, 

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

597 self.nr // nc, nc), 

598 fmt=nc * [numformat]) 

599 else: 

600 np.savetxt(fd, 

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

602 self.nr // nc, nc), 

603 fmt=nc * [numformat]) 

604 

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

606 # as r*phi for alloy format 

607 for i in range(self.Nelements): 

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

609 np.savetxt(fd, 

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

611 nc), 

612 fmt=nc * [numformat]) 

613 

614 if self.form == 'adp': 

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

616 for i in range(self.Nelements): 

617 for j in range(i + 1): 

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

619 

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

621 for i in range(self.Nelements): 

622 for j in range(i + 1): 

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

624 

625 def update(self, atoms): 

626 # check all the elements are available in the potential 

627 self.Nelements = len(self.elements) 

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

629 unavailable = np.logical_not( 

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

631 

632 if np.any(unavailable): 

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

634 elements[unavailable]) 

635 

636 # cutoffs need to be a vector for NeighborList 

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

638 

639 # convert the elements to an index of the position 

640 # in the eam format 

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

642 for el in atoms.get_chemical_symbols()]) 

643 self.pbc = atoms.get_pbc() 

644 

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

646 # local electron density we cannot just calculate and use 

647 # one way neighbors 

648 self.neighbors = NeighborList(cutoffs, 

649 skin=self.parameters.skin, 

650 self_interaction=False, 

651 bothways=True) 

652 self.neighbors.update(atoms) 

653 

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

655 system_changes=all_changes): 

656 """EAM Calculator 

657 

658 atoms: Atoms object 

659 Contains positions, unit-cell, ... 

660 properties: list of str 

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

662 of 'energy', 'forces' 

663 system_changes: list of str 

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

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

666 'pbc', 'initial_charges' and 'initial_magmoms'. 

667 """ 

668 

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

670 

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

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

673 self.update(self.atoms) 

674 self.calculate_energy(self.atoms) 

675 

676 if 'forces' in properties: 

677 self.calculate_forces(self.atoms) 

678 

679 # check we have all the properties requested 

680 for property in properties: 

681 if property not in self.results: 

682 if property == 'energy': 

683 self.calculate_energy(self.atoms) 

684 

685 if property == 'forces': 

686 self.calculate_forces(self.atoms) 

687 

688 # we need to remember the previous state of parameters 

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

690# self.read_potential(potential) 

691 

692 def calculate_energy(self, atoms): 

693 """Calculate the energy 

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

695 the embedding energy of each atom into the electron cloud 

696 generated by its neighbors 

697 """ 

698 

699 pair_energy = 0.0 

700 embedding_energy = 0.0 

701 mu_energy = 0.0 

702 lam_energy = 0.0 

703 trace_energy = 0.0 

704 

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

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

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

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

709 

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

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

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

713 

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

715 atoms.positions[i]) 

716 

717 # calculate the distance to the nearest neighbors 

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

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

720 

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

722 for j_index in range(self.Nelements): 

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

724 if not use.any(): 

725 continue 

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

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

728 

729 if self.form == 'fs': 

730 density = np.sum( 

731 self.electron_density[j_index, self.index[i]](r[nearest][use])) 

732 else: 

733 density = np.sum( 

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

735 self.total_density[i] += density 

736 

737 if self.form == 'adp': 

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

739 r[nearest][use], 

740 rvec[nearest][use], 

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

742 

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

744 r[nearest][use], 

745 rvec[nearest][use], 

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

747 

748 # add in the electron embedding energy 

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

750 self.total_density[i]) 

751 

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

753 

754 if self.form == 'adp': 

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

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

757 

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

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

760 

761 adp_result = dict(adp_mu=mu_energy, 

762 adp_lam=lam_energy, 

763 adp_trace=trace_energy) 

764 components.update(adp_result) 

765 

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

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

768 

769 energy = 0.0 

770 for i in components.keys(): 

771 energy += components[i] 

772 

773 self.energy_free = energy 

774 self.energy_zero = energy 

775 

776 self.results['energy_components'] = components 

777 self.results['energy'] = energy 

778 

779 def calculate_forces(self, atoms): 

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

781 

782 self.update(atoms) 

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

784 

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

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

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

788 # create a vector of relative positions of neighbors 

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

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

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

792 

793 d_embedded_energy_i = self.d_embedded_energy[ 

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

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

796 

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

798 urvec[j] = urvec[j] / r[j] 

799 

800 for j_index in range(self.Nelements): 

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

802 if not use.any(): 

803 continue 

804 rnuse = r[nearest][use] 

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

806 if self.form == 'fs': 

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

808 (d_embedded_energy_i * 

809 self.d_electron_density[j_index, self.index[i]](rnuse)) + 

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

811 self.d_electron_density[self.index[i], j_index](rnuse))) 

812 else: 

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

814 (d_embedded_energy_i * 

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

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

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

818 

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

820 

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

822 adp_forces = self.angular_forces( 

823 self.mu[i], 

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

825 self.lam[i], 

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

827 rnuse, 

828 rvec[nearest][use], 

829 self.index[i], 

830 j_index) 

831 

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

833 

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

835 # calculate the extra components for the adp forces 

836 # rvec are the relative positions to atom i 

837 psi = np.zeros(mu.shape) 

838 for gamma in range(3): 

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

840 

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

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

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

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

845 

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

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

848 axis=1) 

849 term4 = 0.0 

850 for alpha in range(3): 

851 for beta in range(3): 

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

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

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

855 

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

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

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

859 

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

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

862 # on the NIST website with the AlH potential 

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

864 

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

866 

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

868 # calculate the dipole contribution 

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

870 

871 return mu # sign to agree with lammps 

872 

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

874 # slow way of calculating the quadrupole contribution 

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

876 

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

878 qr = q(r) 

879 for alpha in range(3): 

880 for beta in range(3): 

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

882 

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

884 

885 def deriv(self, spline): 

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

887 def d_spline(aspline): 

888 return spline(aspline, 1) 

889 

890 return d_spline 

891 

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

893 """Plot the individual curves""" 

894 

895 import matplotlib.pyplot as plt 

896 

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

898 nrow = 2 

899 elif self.form == 'adp': 

900 nrow = 3 

901 else: 

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

903 

904 if hasattr(self, 'r'): 

905 r = self.r 

906 else: 

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

908 

909 if hasattr(self, 'rho'): 

910 rho = self.rho 

911 else: 

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

913 

914 plt.subplot(nrow, 2, 1) 

915 self.elem_subplot(rho, self.embedded_energy, 

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

917 name, plt) 

918 

919 plt.subplot(nrow, 2, 2) 

920 if self.form == 'fs': 

921 self.multielem_subplot(r, self.electron_density, 

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

923 else: 

924 self.elem_subplot(r, self.electron_density, 

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

926 

927 plt.subplot(nrow, 2, 3) 

928 self.multielem_subplot(r, self.phi, 

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

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

931 

932 if self.form == 'adp': 

933 plt.subplot(nrow, 2, 5) 

934 self.multielem_subplot(r, self.d, 

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

936 

937 plt.subplot(nrow, 2, 6) 

938 self.multielem_subplot(r, self.q, 

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

940 

941 plt.plot() 

942 

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

944 plt.xlabel(xlabel) 

945 plt.ylabel(ylabel) 

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

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

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

949 plt.legend() 

950 

951 def multielem_subplot(self, curvex, curvey, xlabel, ylabel, name, plt, half=True): 

952 plt.xlabel(xlabel) 

953 plt.ylabel(ylabel) 

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

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

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

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

958 plt.legend()