r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# flake8: noqa

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

4# eam.py

5# Embedded Atom Method Potential

6# These routines integrate with the ASE simulation environment

7# Paul White (Oct 2012)

8# UNCLASSIFIED

11import os

12import numpy as np

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

20class EAM(Calculator):

21 r"""

23 EAM Interface Documentation

25Introduction

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

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.

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

39Interatomic Potentials Repository Project at

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

42Theory

43======

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

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

50each alloy and their cross interactions.

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

53given by the EAM potential as

55.. math::

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

58and

60.. math::

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

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.

71The ADP potential is defined as

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

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}.

83The fs potential is defined as

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})

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.

93Running the Calculator

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

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.

104For example::

106 from ase.calculators.eam import EAM

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

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

110 mishin.plot()

112 slab.calc = mishin

113 slab.get_potential_energy()

114 slab.get_forces()

116The breakdown of energy contribution from the indvidual components are

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

119Arguments

120=========

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

123Keyword Description

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

125potential 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

130elements[N] array of N element abbreviations

132embedded_energy[N] arrays of embedded energy functions

134electron_density[N] arrays of electron density functions

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

138d_embedded_energy[N] arrays of derivative embedded energy functions

140d_electron_density[N] arrays of derivative electron density functions

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

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

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

148skin 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.

153form 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

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

161Additional parameters for writing potential files

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

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

165.alloy, .adp or fs format file.

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

168Keyword Description

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

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

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

174mass[N] Atomic mass of each element

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

178lattice[N] Lattice type

180nrho No. of rho samples along embedded energy curve

182drho Increment for sampling density

184nr No. of radial points along density and pair

185 potential curves

187dr Increment for sampling radius

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

191Special features

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

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.

199Notes/Issues

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

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.

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.

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

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

215 forces.

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.

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

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

224 1285.

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

227 Acta Materialia 53 2005 4029--4041.

230End EAM Interface Documentation

231 """

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

235 default_parameters = dict(

236 skin=1.0,

237 potential=None,

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

240 b'blank\n'])

242 def __init__(self, restart=None,

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

246 self.form = form

248 if 'potential' in kwargs:

252 label, atoms, **kwargs)

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')

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))

270 def set_form(self, name):

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

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

274 if extension == '.eam':

275 self.form = 'eam'

276 elif extension == '.alloy':

277 self.form = 'alloy'

280 elif extension == '.fs':

281 self.form = 'fs'

282 else:

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

287 and creates the interpolation functions from the data

288 """

290 if isinstance(filename, str):

291 with open(filename) as fd:

293 else:

294 fd = filename

298 if self.form is None:

299 self.set_form(fd.name)

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

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

314 data = lines_to_list(lines[1:])

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

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]]

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])

330 n = 9 + self.nrho

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

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

334 self.nr])

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)

341 self.density_data = np.array(

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

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

346 i = 3

348 data = lines_to_list(lines[i:])

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

351 d = 1

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

353 d += self.Nelements

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])

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

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

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

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])

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

392 elif self.form == 'fs':

394 i = 3

396 data = lines_to_list(lines[i:])

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

399 d = 1

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

401 d += self.Nelements

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])

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

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

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

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])

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

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

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

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()

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)

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])

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

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

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)

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

485 if j != i:

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

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

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)

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])

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

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

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)

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

518 if j != i:

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

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

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)

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])

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]

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

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

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 """

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

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

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

569 assert self.nr % nc == 0

570 assert self.nrho % nc == 0

573 fd.write(line)

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

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

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

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

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

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)

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

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

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])

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])

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])

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])

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]))

636 if np.any(unavailable):

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

638 elements[unavailable])

640 # cutoffs need to be a vector for NeighborList

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

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()

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)

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

659 system_changes=all_changes):

660 """EAM Calculator

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 """

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

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)

680 if 'forces' in properties:

681 self.calculate_forces(self.atoms)

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)

689 if property == 'forces':

690 self.calculate_forces(self.atoms)

692 # we need to remember the previous state of parameters

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

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 """

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

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

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

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

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())

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

719 atoms.positions[i])

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

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.

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

744 r[nearest][use],

745 rvec[nearest][use],

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

749 r[nearest][use],

750 rvec[nearest][use],

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

753 # add in the electron embedding energy

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

755 self.total_density[i])

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

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

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

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.

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

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

774 energy = 0.0

775 for i in components.keys():

776 energy += components[i]

778 self.energy_free = energy

779 self.energy_zero = energy

781 self.results['energy_components'] = components

782 self.results['energy'] = energy

784 def calculate_forces(self, atoms):

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

787 self.update(atoms)

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

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]

798 d_embedded_energy_i = self.d_embedded_energy[

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

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

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

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

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)))

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

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)

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)

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)

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

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.

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

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

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

875 # calculate the dipole contribution

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

878 return mu # sign to agree with lammps

881 # slow way of calculating the quadrupole contribution

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

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]

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

892 def deriv(self, spline):

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

894 def d_spline(aspline):

895 return spline(aspline, 1)

897 return d_spline

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

900 """Plot the individual curves"""

902 import matplotlib.pyplot as plt

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

905 nrow = 2

907 nrow = 3

908 else:

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

911 if hasattr(self, 'r'):

912 r = self.r

913 else:

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

916 if hasattr(self, 'rho'):

917 rho = self.rho

918 else:

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

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)

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)

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

942 plt.subplot(nrow, 2, 5)

943 self.multielem_subplot(r, self.d,

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

946 plt.subplot(nrow, 2, 6)

947 self.multielem_subplot(r, self.q,

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

950 plt.plot()

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()

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()