Embedded Atom Method (EAM) Tabulation

An EAM model is defined by constructing instances of atsim.potentials.EAMPotential describing each species within the model. EAMPotential encapsulates the density and embedding functions specific to each species’ many bodied interactions. In addition the purely pairwise interactions within the EAM are defined using a list of atsim_potentials.Potential objects.

Once the EAM model has been described in terms of EAMPotential and Potential objects it can be tabulated for specific simulation codes. In addition to the differences in table files expected by different simulation codes, there are several variations on the embedded atom method, in order to support this variety, the atsim_potentials module contains several tabulation functions:

Function File-Format Simulation Code Example
writeFuncFL() funcfl LAMMPS Example 1: Ag in LAMMPS
writeSetFL() setfl LAMMPS Example 2a: Al-Cu in LAMMPS
writeTABEAM() TABEAM DL_POLY Example 2b: Al-Cu in LAMMPS
writeSetFLFinnisSinclair() setfl LAMMPS Example 3a: Al-Fe Finnis-Sinclair in LAMMPS
writeTABEAMFinnisSinclair() TABEAM DL_POLY Example 3b: Al-Fe Finnis-Sinclair in DL_POLY

Examples

Example 1: Using writeFuncFL() to Tabulate Ag Potential for LAMMPS

This example shows how to use writeFuncFL() function to tabulate an EAM model for the simulation of Ag metal. How to use this tabulation within LAMMPS will then be demonstrated. The final tabulation script can be found in eam_tabulate_example1.py.

Model Description

Within this example the Ag potential of Sutton will be tabulated [1]. Within the EAM the energy (\(E_i\)) of an atom \(i\) whose species is \(\alpha\) is given by:

\[E_i = F_\alpha \left( \sum_{j \neq i} \rho_\beta(r_{ij}) \right) + \frac{1}{2} \sum_{j \neq i} \phi_{\alpha \beta} (r_{ij})\]

Note

  • \(\rho_\beta(r_{ij})\) is the density function which gives the electron density for atom \(j\) with species \(\beta\) as a function of its separation from atom \(i\), \(r_{ij}\).
  • The electron density for atom \(i\) is obtained by summing over the density (\(\rho_\beta (r_{ij}\)) contributions due to its neighbours.
  • The embedding function \(F_\alpha(\rho)\) is used to calculate the many-bodied energy contribution from this summed electron density.
  • The sum \(\frac{1}{2} \sum_{j \neq i} \phi_{\alpha \beta} (r_{ij})\) gives the pair-potential contribution to atom \(i\)‘s energy.
  • \(\phi_{\alpha \beta} (r_{ij})\) are simply pair potentials that describe the energy between two atoms as a function of their separation.

The embedding function used by Sutton is:

\[F_\alpha (\rho) = - \sqrt{\rho}\]

and the density function is:

\[\rho_\beta (r_{ij}) = \left( \frac{a}{r_{ij}} \right)^m\]

whilst pair interactions are given by:

\[\phi_{\alpha \beta} (r_{ij}) = \left( \frac{b}{r_{ij}} \right)^n\]

The model parameters are given as:

Parameter Value
\(m\) 6
\(n\) 12
\(a\) \(2.928323832 \text{Å} \text{eV}^{\frac{1}{3}}\)
\(b\) \(2.485883762 \text{eV}^\frac{1}{12} \text{Å}\)

Define the Model

It is now necessary to describe the model in python code. Hard-coding the model parameters from the previous table, embed() and density() functions can be defined for \(F_{\text{Ag}} (\rho)\) and \(\rho_{\text{Ag}}\) respectively:

import math
from atsim.potentials import EAMPotential
from atsim.potentials import Potential

def embed(rho):
  return -math.sqrt(rho)


def density(rij):
  if rij == 0:
    return 0.0
  return (2.928323832 / rij) ** 6.0

The embedding and density functions should then be wrapped in an EAMPotential object to create a single item list:

  # Create EAMPotential
  eamPotentials = [ EAMPotential("Ag", 47, 107.8682, embed, density) ]

Similarly the pair potential component, \(\phi_{\text{Ag}-\text{Ag}} (r_{ij}\), of the model can easily be defined as:

def pair_AgAg(rij):
    if rij == 0:
      return 0.0
    return (2.485883762/rij) ** 12

This can then be wrapped in a atsim.potentials.Potential object to create a list of pair potentials.

  pairPotentials = [ Potential('Ag', 'Ag', pair_AgAg) ]

Note

writeFuncFL() only accepts a single Potential object and this should be the X-X interaction (where X is the species for which the funcfl tabulation is being created). Within the ‘pair’-potential is tabulated as

\(\sqrt{\frac{\phi(r_{ij}) r_{ij}}{27.2 \times 0.529}}\)

The numerical constants (27.2 and 0.529) convert from eV and Å into the units of Hartree and Bohr radius used by the funcfl format. The square rooting of the potential function is important: the simulation code effectively reconstitutes a pair potential by multiplying two of these tabulated square-rooted functions (one for each species in each interacting pair) together. If atoms \(i\) and \(j\) in an interacting pair, have the same species then effectively the original pair-potential is obtained (albeit multiplied by \(r_{ij}\)).

By comparison, if multiple funcfl files are used to define multiple species within a simulation (e.g. for alloy systems), then the pair potential functions of each species are effectively ‘mixed’ when they are multiplied together for heterogeneous atom pairs. If more control is required, with pair-potential functions specific to distinct pairs of species being necessary, then the setfl format produced by the writeSetFL() and writeSetFLFinnisSinclair() functions should be used instead.

Now all the components of the model have been defined a table file can be created in the funcfl format. Before doing this, it is necessary to choose appropriate density and separation cut-offs together with \(d r_{ij}\) and \(d \rho\) increments for the density/pair functions and embedding function respectively:

  • Here a \(d \rho\) value of 0.001 will be used and 50000 density values tabulated.
  • This means the maximum density that can be accepted by the embedding function is \(49999 \times 0.001 = 49.999\)
  • \(dr = 0.001\) Å using 12000 rows.
  • The pair-potential cut-off and the maximum \(r_{ij}\) value for the density function is therefore 11.999 Å.

Invoking the writeFuncFL() function with these values and the EAMPotential and Potential objects, can be used to tabulate the Ag potential into the Ag.eam file:

  nrho = 50000
  drho = 0.001

  nr = 12000
  dr = 0.001

  from atsim.potentials import writeFuncFL

  with open("Ag.eam", 'wb') as outfile:
    writeFuncFL(
            nrho, drho,
            nr, dr,
            eamPotentials,
            pairPotentials,
            out= outfile,
            title='Sutton Chen Ag')

Putting this together the following script is obtained (this script can also be downloaded eam_tabulate_example1.py:

#! /usr/bin/env python
import math
from atsim.potentials import EAMPotential
from atsim.potentials import Potential

def embed(rho):
  return -math.sqrt(rho)


def density(rij):
  if rij == 0:
    return 0.0
  return (2.928323832 / rij) ** 6.0


def pair_AgAg(rij):
    if rij == 0:
      return 0.0
    return (2.485883762/rij) ** 12


def main():
  # Create EAMPotential
  eamPotentials = [ EAMPotential("Ag", 47, 107.8682, embed, density) ]
  pairPotentials = [ Potential('Ag', 'Ag', pair_AgAg) ]

  nrho = 50000
  drho = 0.001

  nr = 12000
  dr = 0.001

  from atsim.potentials import writeFuncFL

  with open("Ag.eam", 'wb') as outfile:
    writeFuncFL(
            nrho, drho,
            nr, dr,
            eamPotentials,
            pairPotentials,
            out= outfile,
            title='Sutton Chen Ag')

if __name__ == "__main__":
  main()

Running this script will produce a table file named Ag.eam in the same directory as the script:

python eam_tabulate_example1.py

Using the Ag.eam file within LAMMPS

This section of the example will now demonstrate how the table file can be used used to perform a static energy minimisation of an FCC Ag structure in LAMMPS.

Place the following in a file called fcc.lmpstruct in the same directory as the Ag.eam file you created previously. This describes a single FCC cell with a wildly inaccurate lattice parameter:

Title


4 atoms
1 atom types
0.0 5.000000 xlo xhi
0.0 5.000000 ylo yhi
0.0 5.000000 zlo zhi
0.000000 0.000000 0.000000 xy xz yz



Masses

1 107.86820000000000163709 #Ag

Atoms

1 0 1 0.000000 0.000000 0.000000 0.000000
2 0 1 0.000000 2.500000 2.500000 0.000000
3 0 1 0.000000 0.000000 2.500000 2.500000
4 0 1 0.000000 2.500000 0.000000 2.500000


The following LAMMPS input file describes a minimisation run. The lines describing potentials are highlighted. Put its contents in a file called example1_minimize.lmpin:

units metal
boundary p p p

atom_style full
read_data fcc.lmpstruct

pair_style eam
pair_coeff 1 1 Ag.eam

fix 1 all box/relax x 0.0 y 0.0 z 0.0

minimize 0.0 1.0e-8 1000 100000 

The pair_style eam command tells LAMMPS to use the EAM and expect pair_coeff commands mapping atom types to particular table files:

pair_style eam

The following pair_coeff directive indicates that the interaction between atom-type 1 (Ag) with itself should use the funcfl formatted file contained within Ag.eam:

pair_coeff 1 1 Ag.eam

The example can then be run by invoking LAMMPS:

lammps -in example1_minimize.lmpin

Example 2a: Tabulate Al-Cu Alloy Potentials Using writeSetFL() for LAMMPS

Within the following example the process required to generate and use a setfl file that tabulates the Al-Cu alloy model of Zhou et al [2]. By comparison to the funcfl format, setfl allows multiple elements to be given in the same file and additionally pair-potentials for particular pairs of interacting species can be specified (funcfl relies on the simulation code to ‘mix’ pair-potentials within alloy systems). The eam_tabulate_example2a.py gives a complete example of how the Zhou model can be tabulated.

Model Description

The model makes use of the EAM as described above (see Example 1 Model Description) . The density function, \(\rho_\beta (r_{ij})\) for an atom \(j\) of species \(\beta\) separated from atom \(i\) by \(r_{ij}\) is:

\[\rho_\beta (r_{ij}) = \frac{f_e \exp \left[ -\omega (r_{ij} / r_e - 1) \right] }{1+(r_{ij}/r_e - \lambda)^{20}}\]

where \(f_e\), \(r_e\), \(\omega\) and \(\lambda\) are parameters specific to species \(\beta\). The pair-potential function acting between species \(\alpha\)\(\beta\) is obtained by combining the density functions of the interacting species:

\[\phi_{\alpha\beta} (r_{ij}) = \frac{1}{2} \left[ \frac{\rho_\beta(r_{ij})}{\rho_\alpha(r_{ij})} \phi_{\alpha\alpha}(r_{ij}) + \frac{\rho_\alpha(r_{ij})}{\rho_\beta(r_{ij})} \phi_{\beta\beta}(r_{ij}) \right]\]

The homogeneous pair-interactions, \(\phi_{\alpha\alpha}(r_{ij})\) and \(\phi_{\beta\beta}(r_{ij})\) have the form:

\[\phi_{\alpha \alpha}(r_{ij}) = \frac{A \exp \left[ -\gamma (r_{ij} / r_e - 1) \right] }{1+(r_{ij}/r_e - \kappa)^{20}} - \frac{B \exp \left[ -\omega (r_{ij} / r_e - 1) \right] }{1+(r_{ij}/r_e - \lambda)^{20}}\]

again, \(A\), \(B\), \(\gamma\), \(\omega\), \(\kappa\) and \(\omega\) are parameters specific to the species \(\alpha\).

The embedding function for each species, \(F_\alpha (\rho)\), is defined over three density ranges using the following:

\[\begin{split}F_\alpha (\rho) = \left\{ \begin{array}{ll} \sum_{i=0}^3 F_{ni}\left( \frac{\rho}{\rho_n} - 1\right)^i & \rho < \rho_n, & \rho_n = 0.85 \rho_e \\ \sum_{i=0}^3 F_{i}\left( \frac{\rho}{\rho_e} - 1\right)^i & \rho_n \leq \rho < \rho_0, & \rho_0 = 1.15 \rho_e \\ F_e \left[1-\eta\ln\left(\frac{\rho}{\rho_s}\right)\right]\left(\frac{\rho}{\rho_s}\right)^\eta &\rho_0 \leq \rho & \\ \end{array} \right.\end{split}\]

The model parameters for Cu and Al are given in the following table:

Parameter Cu Al
\(r_e\) 2.556162 2.863924
\(f_e\) 1.554485 1.403115
\(\rho_e\) 21.175871 20.418205
\(\rho_s\) 21.175395 23.195740
\(\gamma\) 8.127620 6.613165
\(\omega\) 4.334731 3.527021
\(A\) 0.396620 0.314873
\(B\) 0.548085 0.365551
\(\kappa\) 0.308782 0.379846
\(\lambda\) 0.756515 0.759692
\(F_{n0}\) -2.170269 -2.807602
\(F_{n1}\) -0.263788 -0.301435
\(F_{n2}\) 1.088878 1.258562
\(F_{n3}\) -0.817603 -1.247604
\(F_0\) -2.19 -2.83
\(F_1\) 0 0
\(F_2\) 0.561830 0.622245
\(F_3\) -2.100595 -2.488244
\(\eta\) 0.310490 0.785902
\(F_e\) -2.186568 -2.824528

Note

The Al \(A\) value is given as 0.134873 in Zhou’s original Phys. Rev. B paper. However parameter file provided by Zhou for this model, at http://www.ctcms.nist.gov/potentials/Zhou04.html gives the parameter as 0.314873. It is this latter value that is used here.

In addition the final term of the embedding function has been modified to match that used in fortran tabulation code also provided at http://www.ctcms.nist.gov/potentials/Zhou04.html

Define the Model

A series of python functions are defined to describe the embedding, density and pair interaction functions. To encourage code re-use a number of function factories are defined. Using the parameters passed to them they return specialised functions appropriate for the parameters. The given factory functions make use of python’s support for closures in their implementation.

The makeFunc() factory function is used to define density functions. As this functional form is also used as a component of the pair-potentials makeFunc() is re-used within the makePairPotAA() factory function.

def makeFunc(a, b, r_e, c):
  # Creates functions of the form used for density function.
  # Functional form also forms components of pair potential.
  def func(r):
    return (a * math.exp(-b*(r/r_e -1)))/(1+(r/r_e - c)**20.0)
  return func

The following factory returns the functions used to describe the homogeneous Al-Al and Cu-Cu pair-potential interactions:

def makePairPotAA(A, gamma, r_e, kappa,
                  B, omega, lamda):
  # Function factory that returns functions parameterised for homogeneous pair interactions
  f1 = makeFunc(A, gamma, r_e, kappa)
  f2 = makeFunc(B, omega, r_e, lamda)
  def func(r):
    return f1(r) - f2(r)
  return func

Whilst makePairPotAB() describes the Al-Cu pair-potential:

def makePairPotAB(dens_a, phi_aa, dens_b, phi_bb):
  # Function factory that returns functions parameterised for heterogeneous pair interactions
  def func(r):
    return 0.5 * ( (dens_b(r)/dens_a(r) * phi_aa(r)) + (dens_a(r)/dens_b(r) * phi_bb(r)) )
  return func

The makeEmbed() function describes the embedding function:

def makeEmbed(rho_e, rho_s, F_ni, F_i, F_e, eta):
  # Function factory returning parameterised embedding function.
  rho_n = 0.85*rho_e
  rho_0 = 1.15*rho_e

  def e1(rho):
    return sum([F_ni[i] * (rho/rho_n - 1)**float(i) for i in range(4)])

  def e2(rho):
    return sum([F_i[i] * (rho/rho_e - 1)**float(i) for i in range(4)])

  def e3(rho):
    return F_e * (1.0 - eta*math.log(rho/rho_s)) * (rho/rho_s)**eta

  def func(rho):
    if rho < rho_n:
      return e1(rho)
    elif rho_n <= rho < rho_0:
      return e2(rho)
    return e3(rho)
  return func

Lists of EAMPotential and Potential objects are created and returned as a tuple by the makePotentialObjects() function within eam_tabulate_example2a.py. Before invoking the factory functions we just defined, the model parameters are assigned to easily identifiable variables within this function:

def makePotentialObjects():
  # Potential parameters
  r_eCu     = 2.556162
  f_eCu     = 1.554485
  gamma_Cu  = 8.127620
  omega_Cu  = 4.334731
  A_Cu      = 0.396620
  B_Cu      = 0.548085
  kappa_Cu  = 0.308782
  lambda_Cu = 0.756515

  rho_e_Cu  = 21.175871
  rho_s_Cu  = 21.175395
  F_ni_Cu   = [-2.170269, -0.263788, 1.088878, -0.817603]
  F_i_Cu    = [-2.19, 0.0, 0.561830, -2.100595]
  eta_Cu     = 0.310490
  F_e_Cu    = -2.186568

  r_eAl     = 2.863924
  f_eAl     = 1.403115
  gamma_Al  = 6.613165
  omega_Al  = 3.527021
  # A_Al      = 0.134873
  A_Al      = 0.314873
  B_Al      = 0.365551
  kappa_Al  = 0.379846
  lambda_Al = 0.759692

  rho_e_Al  =  20.418205
  rho_s_Al  =  23.195740
  F_ni_Al   = [-2.807602, -0.301435, 1.258562, -1.247604]
  F_i_Al    = [-2.83, 0.0, 0.622245, -2.488244]
  eta_Al     = 0.785902
  F_e_Al    = -2.824528

Now the functions required by the EAMPotential instances for Al and Cu can be created:

  # Define the density functions
  dens_Cu   = makeFunc(f_eCu, omega_Cu, r_eCu, lambda_Cu )
  dens_Al   = makeFunc(f_eAl, omega_Al, r_eAl,  lambda_Al )

  # Finally, define embedding functions for each species
  embed_Cu = makeEmbed(rho_e_Cu, rho_s_Cu, F_ni_Cu, F_i_Cu, F_e_Cu, eta_Cu)
  embed_Al = makeEmbed(rho_e_Al, rho_s_Al, F_ni_Al, F_i_Al, F_e_Al, eta_Al)

Now these are wrapped up in EAMPotential objects to give the eamPotentials list:

  # Wrap them in EAMPotential objects
  eamPotentials = [
    EAMPotential("Al", 13, 26.98, embed_Al, dens_Al),
    EAMPotential("Cu", 29, 63.55, embed_Cu, dens_Cu)]

Similarly, using the makePairPotAA() and makePairPotAB() function factories the Potential objects required for the tabulation are defined:

  # Define pair functions
  pair_CuCu = makePairPotAA(A_Cu, gamma_Cu, r_eCu, kappa_Cu,
                           B_Cu, omega_Cu, lambda_Cu)

  pair_AlAl = makePairPotAA(A_Al, gamma_Al, r_eAl, kappa_Al,
                           B_Al, omega_Al, lambda_Al)

  pair_AlCu = makePairPotAB(dens_Cu, pair_CuCu, dens_Al, pair_AlAl)

  # Wrap them in Potential objects
  pairPotentials = [
      Potential('Al', 'Al', pair_AlAl),
      Potential('Cu', 'Cu', pair_CuCu),
      Potential('Al', 'Cu', pair_AlCu)]

Now we have all the objects required for writeSetFL(). The next excerpt call makeObjects() to get the EAM and pair-potential objects before invoking the tabulation function, writing the data into a file called Zhou_AlCu.setfl:

def main():
  eamPotentials, pairPotentials = makePotentialObjects()

  # Perform tabulation
  # Make tabulation
  nrho = 2000
  drho = 0.05

  nr = 2000
  dr = 0.003

  with open("Zhou_AlCu.setfl", 'wb') as outfile:
    writeSetFL(
      nrho, drho,
      nr, dr,
      eamPotentials,
      pairPotentials,
      out= outfile,
      comments = ['Zhou Al Cu', "", ""]) # <-- Note: title lines given as list of three strings

Putting this all together gives the following script (which can also be downloaded using the following link eam_tabulate_example2a.py:). Running this (python eam_tabulate_example2a.py) produces the Zhou_AlCu.setfl file in current working directory.

#! /usr/bin/env python

from atsim.potentials import writeSetFL
from atsim.potentials import Potential
from atsim.potentials import EAMPotential

import math

def makeFunc(a, b, r_e, c):
  # Creates functions of the form used for density function.
  # Functional form also forms components of pair potential.
  def func(r):
    return (a * math.exp(-b*(r/r_e -1)))/(1+(r/r_e - c)**20.0)
  return func


def makePairPotAA(A, gamma, r_e, kappa,
                  B, omega, lamda):
  # Function factory that returns functions parameterised for homogeneous pair interactions
  f1 = makeFunc(A, gamma, r_e, kappa)
  f2 = makeFunc(B, omega, r_e, lamda)
  def func(r):
    return f1(r) - f2(r)
  return func


def makePairPotAB(dens_a, phi_aa, dens_b, phi_bb):
  # Function factory that returns functions parameterised for heterogeneous pair interactions
  def func(r):
    return 0.5 * ( (dens_b(r)/dens_a(r) * phi_aa(r)) + (dens_a(r)/dens_b(r) * phi_bb(r)) )
  return func

def makeEmbed(rho_e, rho_s, F_ni, F_i, F_e, eta):
  # Function factory returning parameterised embedding function.
  rho_n = 0.85*rho_e
  rho_0 = 1.15*rho_e

  def e1(rho):
    return sum([F_ni[i] * (rho/rho_n - 1)**float(i) for i in range(4)])

  def e2(rho):
    return sum([F_i[i] * (rho/rho_e - 1)**float(i) for i in range(4)])

  def e3(rho):
    return F_e * (1.0 - eta*math.log(rho/rho_s)) * (rho/rho_s)**eta

  def func(rho):
    if rho < rho_n:
      return e1(rho)
    elif rho_n <= rho < rho_0:
      return e2(rho)
    return e3(rho)
  return func

def makePotentialObjects():
  # Potential parameters
  r_eCu     = 2.556162
  f_eCu     = 1.554485
  gamma_Cu  = 8.127620
  omega_Cu  = 4.334731
  A_Cu      = 0.396620
  B_Cu      = 0.548085
  kappa_Cu  = 0.308782
  lambda_Cu = 0.756515

  rho_e_Cu  = 21.175871
  rho_s_Cu  = 21.175395
  F_ni_Cu   = [-2.170269, -0.263788, 1.088878, -0.817603]
  F_i_Cu    = [-2.19, 0.0, 0.561830, -2.100595]
  eta_Cu     = 0.310490
  F_e_Cu    = -2.186568

  r_eAl     = 2.863924
  f_eAl     = 1.403115
  gamma_Al  = 6.613165
  omega_Al  = 3.527021
  # A_Al      = 0.134873
  A_Al      = 0.314873
  B_Al      = 0.365551
  kappa_Al  = 0.379846
  lambda_Al = 0.759692

  rho_e_Al  =  20.418205
  rho_s_Al  =  23.195740
  F_ni_Al   = [-2.807602, -0.301435, 1.258562, -1.247604]
  F_i_Al    = [-2.83, 0.0, 0.622245, -2.488244]
  eta_Al     = 0.785902
  F_e_Al    = -2.824528

  # Define the density functions
  dens_Cu   = makeFunc(f_eCu, omega_Cu, r_eCu, lambda_Cu )
  dens_Al   = makeFunc(f_eAl, omega_Al, r_eAl,  lambda_Al )

  # Finally, define embedding functions for each species
  embed_Cu = makeEmbed(rho_e_Cu, rho_s_Cu, F_ni_Cu, F_i_Cu, F_e_Cu, eta_Cu)
  embed_Al = makeEmbed(rho_e_Al, rho_s_Al, F_ni_Al, F_i_Al, F_e_Al, eta_Al)

  # Wrap them in EAMPotential objects
  eamPotentials = [
    EAMPotential("Al", 13, 26.98, embed_Al, dens_Al),
    EAMPotential("Cu", 29, 63.55, embed_Cu, dens_Cu)]

  # Define pair functions
  pair_CuCu = makePairPotAA(A_Cu, gamma_Cu, r_eCu, kappa_Cu,
                           B_Cu, omega_Cu, lambda_Cu)

  pair_AlAl = makePairPotAA(A_Al, gamma_Al, r_eAl, kappa_Al,
                           B_Al, omega_Al, lambda_Al)

  pair_AlCu = makePairPotAB(dens_Cu, pair_CuCu, dens_Al, pair_AlAl)

  # Wrap them in Potential objects
  pairPotentials = [
      Potential('Al', 'Al', pair_AlAl),
      Potential('Cu', 'Cu', pair_CuCu),
      Potential('Al', 'Cu', pair_AlCu)]

  return eamPotentials, pairPotentials

def main():
  eamPotentials, pairPotentials = makePotentialObjects()

  # Perform tabulation
  # Make tabulation
  nrho = 2000
  drho = 0.05

  nr = 2000
  dr = 0.003

  with open("Zhou_AlCu.setfl", 'wb') as outfile:
    writeSetFL(
      nrho, drho,
      nr, dr,
      eamPotentials,
      pairPotentials,
      out= outfile,
      comments = ['Zhou Al Cu', "", ""]) # <-- Note: title lines given as list of three strings

if __name__ == '__main__':
  main()

Using the Zhou_AlCu.setfl file within LAMMPS

Within LAMMPS the setfl files generated by writeSetFL() are used with the eam/alloy pair_style. The pair_coeff directive used with this pair_style effectively maps LAMMPS species numbers to the element names within the table file.

Single Element Systems

Assuming a LAMMPS system containing only Al (i.e. Al is species 1) then the pair_style and pair_coeff directives would be given as:

pair_style eam/alloy
pair_coeff * * Zhou_AlCu.setfl Al

Likewise if a copper system was being simulated:

pair_style eam/alloy
pair_coeff * * Zhou_AlCu.setfl Cu
Mixed Al-Cu System

For an Al-Cu system where Al is species 1 and Cu species 2 then the directives would be:

pair_style eam/alloy
pair_coeff * * Zhou_AlCu.setfl Al Cu

Or if Cu was 1 and Al 2:

pair_style eam/alloy
pair_coeff * * Zhou_AlCu.setfl Cu Al

Example 2b: Tabulate Al-Cu Alloy Potentials Using writeTABEAM() for DL_POLY

The tabulation script used with Example 2a can be easily modified to produce the TABEAM format expected by the DL_POLY simulation code. See the tabulation script for this example: eam_tabulate_example2b.py.

The EAMPotential and Potential lists are created in exactly the same way as Example 2a, however rather than calling writeSetFL() the main() function is modified to use the DL_POLY specific writeTABEAM() function instead and to write into a file named TABEAM. The main() function of eam_tabulate_example2b.py is now given:

def main():
  eamPotentials, pairPotentials = makePotentialObjects()

  # Perform tabulation
  # Make tabulation
  nrho = 2000
  drho = 0.05

  nr = 2000
  dr = 0.003

  with open("TABEAM", 'wb') as outfile:
    writeTABEAM(
      nrho, drho,
      nr, dr,
      eamPotentials,
      pairPotentials,
      out= outfile)

Using the TABEAM file with DL_POLY

Running eam_tabulate_example2b.py will create a file named TABEAM in the working directory. This should be copied into the simulation directory containing the DL_POLY input files (CONTROL, CONFIG and FIELD).

The following should be added at the bottom of the FIELD file:

metal 3
Al Al eam
Cu Cu eam
Al Cu eam

Example 3a: Tabulate Al-Fe Finnis-Sinclair Potentials Using writeSetFLFinnisSinclair() for LAMMPS

This example will show how to reproduce the EAM model described by Mendelev et al. for Fe segregation at grain boundaries within Al [3]. As a result this example effectively shows how to reproduce the AlFe_mm.eam.fs file provided with the LAMMPS source distribution using the writeSetFLFinnisSinclair() function.

The file format created by writeSetFLFinnisSinclair() is supported by the LAMMPS pair_style eam/fs command. This adds an additional level of flexibility in comparison to the eam/alloy style; when calculating the density surrounding an atom with species \(\alpha\), each neighbouring atom’s contribution to the density is calculated as a function of its separation from the central atom using \(\rho_{\alpha\beta}(r_{ij})\). This means that the density function is now specific to both the central atom species, \(\alpha\) and that of the surrounding atom, \(\beta\). By comparison when using eam/alloy tabulations the same \(\rho_\beta(r_{ij})\) function is used, no matter the type of the central atom. This means that the equation describing eam/fs style models becomes:

\[E_i = F_\alpha \left( \sum_{j \neq i} \rho_{\alpha\beta}(r_{ij}) \right) + \frac{1}{2} \sum_{j \neq i} \phi_{\alpha \beta} (r_{ij})\]

Here a binary Al, Fe, model is being described and the resultant eam/fs file should contain definitions for the following:

  • Pair-Potentials: \(\phi_{\text{Al}\text{Al}}(r_{ij})\), \(\phi_{\text{Fe}\text{Fe}}(r_{ij})\) and \(\phi_{\text{Al}\text{Fe}}(r_{ij})\).
  • Embedding-Functions: \(F_{\text{Al}}(\rho)\) and \(F_{\text{Fe}}(\rho)\).
  • Density-Functions: \(\rho_{\text{Al}\text{Al}}(r_{ij})\), \(\rho_{\text{Fe}\text{Fe}}(r_{ij})\), \(\rho_{\text{Al}\text{Fe}}(r_{ij})\) and \(\rho_{\text{Fe}\text{Al}}(r_{ij})\).

From this it can be seen that, when using eam/fs style potentials, the density functions must have both the \(\alpha\beta\) and \(\beta\alpha\) interactions specified to writeSetFLFinnisSinclair().

Although both the \(\alpha\beta\) and \(\beta\alpha\) can be described using eam/fs files, the Mendelev model used in this example uses the same density function for both Al-Fe and Fe-Al cross density functions [3].

Using writeSetFLFinnisSinclair to Tabulate the Model

As in previous examples it is necessary to define pair, density and embedding functions in python code that are then wrapped in EAMPotential and Potential objects to be passed to the tabulation function. For brevity only the names of the functions, as defined in the attached example file (eam_tabulate_example3a.py) are now given:

  • Pair-Potentials:
    • def ppfuncAlAl(r): - Al-Al pair-potential \(\phi_{\text{Al}\text{Al}}(r_{ij})\).
    • def ppfuncAlFe(r): - Al-Fe pair-potential \(\phi_{\text{Al}\text{Fe}}(r_{ij})\).
    • def ppfuncFeFe(r): - Fe-Fe pair-potential \(\phi_{\text{Fe}\text{Fe}}(r_{ij})\).
  • Embedding-Functions:
    • def AlEmbedFunction(rho): - Al embedding function \(F_{\text{Al}}(\rho)\).
    • def FeEmbedFunction(rho): - Fe embedding function \(F_{\text{Fe}}(\rho)\).
  • Density-Functions:
    • def AlAlDensityFunction(r): - Al density function \(\rho_{\text{Al}\text{Al}}(r_{ij})\).
    • def FeFeDensityFunction(r): - Fe density function \(\rho_{\text{Al}\text{Al}}(r_{ij})\).
    • def FeAlDensityFunction(r): - Al-Fe density function \(\rho_{\text{Al}\text{Fe}}(r_{ij})\).

Note

The functional forms used within the Mendelev paper [3] are somewhat long, and including their implementations here would detract from the readability of this example. However, they are included in the attached python file: eam_tabulate_example3a.py.

These functions are used within the main() function of the eam_tabulate_example3a.py file which is now shown:

def main():
  # Define list of pair potentials
  pairPotentials = [
    Potential('Al', 'Al', ppfuncAlAl),
    Potential('Al', 'Fe', ppfuncAlFe),
    Potential('Fe', 'Fe', ppfuncFeFe)]

  # Assemble the EAMPotential objects
  eamPotentials = [
    #Al
    EAMPotential('Al', 13, 26.98154, AlEmbedFunction,
      { 'Al' : AlAlDensityFunction,
        'Fe' : FeAlDensityFunction },
      latticeConstant = 4.04527,
      latticeType = 'fcc'),
    #Fe
    EAMPotential('Fe', 26, 55.845, FeEmbedFunction,
      { 'Al': FeAlDensityFunction,
        'Fe' : FeFeDensityFunction},
      latticeConstant = 2.855312,
      latticeType = 'bcc') ]

  # Number of grid points and cut-offs for tabulation.
  nrho = 10000
  drho = 3.00000000000000E-2
  nr   = 10000
  dr   = 6.50000000000000E-4

  with open("Mendelev_Al_Fe.eam.fs", "wb") as outfile:
    writeSetFLFinnisSinclair(
      nrho, drho,
      nr, dr,
      eamPotentials,
      pairPotentials,
      outfile)
  1. Breaking main() into its components, first a list of Potential objects is created, this is common with the other tabulation methods already discussed:
  # Define list of pair potentials
  pairPotentials = [
    Potential('Al', 'Al', ppfuncAlAl),
    Potential('Al', 'Fe', ppfuncAlFe),
    Potential('Fe', 'Fe', ppfuncFeFe)]
  1. Next, the EAMPotential objects for Al and Fe are instantiated. This is where the use of writeSetFLFinnisSinclair() differs from writeSetFL(), as a dictionary of density functions is passed to the constructor instead of the single function used previously (see highlighted lines):
  # Assemble the EAMPotential objects
  eamPotentials = [
    #Al
    EAMPotential('Al', 13, 26.98154, AlEmbedFunction,
      { 'Al' : AlAlDensityFunction,
        'Fe' : FeAlDensityFunction },
      latticeConstant = 4.04527,
      latticeType = 'fcc'),
    #Fe
    EAMPotential('Fe', 26, 55.845, FeEmbedFunction,
      { 'Al': FeAlDensityFunction,
        'Fe' : FeFeDensityFunction},
      latticeConstant = 2.855312,
      latticeType = 'bcc') ]
  1. The density function dictionary keys refer to the \(\beta\) species in each \(\alpha\beta\) pair. This means that:
  • for the Al EAMPotential instance:
    • \(\rho_{\text{Al}\text{Al}}\) = AlAlDensityFunction(),
    • \(\rho_{\text{Al}\text{Fe}}\) = FeAlDensityFunction().
  • for the Fe EAMPotential instance:
    • \(\rho_{\text{Fe}\text{Al}}\) = FeAlDensityFunction(),
    • \(\rho_{\text{Fe}\text{Fe}}\) = FeFeDensityFunction().
  1. Finally, having defined the list of EAMPotential instances the writeSetFLFinnisSinclair() function is called, in this case writing the data to Mendelev_Al_Fe.eam.fs in the current directory:
  # Number of grid points and cut-offs for tabulation.
  nrho = 10000
  drho = 3.00000000000000E-2
  nr   = 10000
  dr   = 6.50000000000000E-4

  with open("Mendelev_Al_Fe.eam.fs", "wb") as outfile:
    writeSetFLFinnisSinclair(
      nrho, drho,
      nr, dr,
      eamPotentials,
      pairPotentials,
      outfile)

Using the Mendelev_Al_Fe.eam.fs file within LAMMPS

For a binary system where Al and Fe have IDs of 1 and 2 the Mendelev_Al_Fe.eam.fs file is specified to LAMMPS as follows:

pair_style eam/fs
pair_coeff * * Mendelev_Al_Fe.eam.fs Al Fe

Example 3b: Tabulate Al-Fe Finnis-Sinclair Potentials Using writeTABEAMFinnisSinclair() for DL_POLY

Using exactly the same model definition as for Example 3a, the Al-Fe model can be re-tabulated for DL_POLY with minimal modification to the main() function. The modified version of the tabulation script can be found in eam_tabulate_example3b.py.

The main() function is given below:

def main():
  # Define list of pair potentials
  pairPotentials = [
    Potential('Al', 'Al', ppfuncAlAl),
    Potential('Al', 'Fe', ppfuncAlFe),
    Potential('Fe', 'Fe', ppfuncFeFe)]

  # Assemble the EAMPotential objects
  eamPotentials = [
    #Al
    EAMPotential('Al', 13, 26.98154, AlEmbedFunction,
      { 'Al' : AlAlDensityFunction,
        'Fe' : FeAlDensityFunction },
      latticeConstant = 4.04527,
      latticeType = 'fcc'),
    #Fe
    EAMPotential('Fe', 26, 55.845, FeEmbedFunction,
      { 'Al': FeAlDensityFunction,
        'Fe' : FeFeDensityFunction},
      latticeConstant = 2.855312,
      latticeType = 'bcc') ]

  # Number of grid points and cut-offs for tabulation.
  nrho = 10000
  drho = 3.00000000000000E-2
  nr   = 10000
  dr   = 6.50000000000000E-4
  cutoff = 6.5

  with open("TABEAM", "wb") as outfile:
    writeTABEAMFinnisSinclair(
      nrho, drho,
      nr, dr,
      eamPotentials,
      pairPotentials,
      outfile)

Excluding the import statement at the top of the file, only two lines have been changed (highlighted). The first changes the filename to TABEAM whilst the second tells python to call writeTABEAMFinnisSinclair() instead of writeSetFLFinnisSinclair():

  with open("TABEAM", "wb") as outfile:
    writeTABEAMFinnisSinclair(
      nrho, drho,
      nr, dr,
      eamPotentials,
      pairPotentials,
      outfile)

That’s it, nothing else has changed.

Using the TABEAM file with DL_POLY

Running eam_tabulate_example3b.py produces a file names TABEAM within the working directory. This should be placed in the same directory as the other DL_POLY input files (CONTROL, CONFIG and FIELD). Then the following should be added to the end of the FIELD file:

metal 3
Al Al eeam
Fe Fe eeam
Al Fe eeam

Note

The Extended EAM (eeam) variant of the TABEAM file generated here is only supported in DL_POLY versions >= 4.05.

Footnotes

[1]A.P. Sutton, and J. Chen, “Long-range Finnis-Sinclair potentials”, Philos. Mag. Lett. 61 (1990) 139.
[2]
  1. Zhou, R. Johnson and H. Wadley, “Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers”, Phys. Rev. B. 69 (2004) 144113.
[3](1, 2, 3) M.I. Mendelev, D.J. Srolovitz, G.J. Ackland, and S. Han, “Effect of Fe Segregation on the Migration of a Non-Symmetric Σ5 Tilt Grain Boundary in Al”, J. Mater. Res. 20 (2011) 208.