Procedural EAM Tabulation

As an alternative to the EAM_Tabulation objects described here Embedded Atom Method (EAM) Tabulation a procedural interface is also provided for EAM tabulation using the following 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.

The same model as used for the SetFL_EAMTabulation example (Example 1: Using SetFL_EAMTabulation to Tabulate Ag Potential for LAMMPS) and is described the same way in python. In terms of the code, the only significant difference between the object based example and this one, is the use of writeFuncFL() to tabulate the model into a file. The output format used in this example is also different, it uses the simpler funcfl format. Each funcfl file contains a single species, making alloy systems less convenient. Further more, alloy models are simulated by combining the funcfl files using pre-determined mixing rules meaning there is much less control over the specific interactions between the various elements in the alloy. To prodce the same setfl files as produced by the SetFL_EAMTabulation class, the writeSetFL() function can be used (an example of which is given below).

The embed() and density() functions are 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 potentialsPotential 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", 'w') 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", 'w') 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.

This example is almost entirely the same as that given for the object based interface (Example 2a: Tabulate Al-Cu Alloy Potentials Using SetFL_EAMTabulation for LAMMPS) with the only difference being the use of the writeSetFL() function for the final tabulation. For a description of the Zhou model and how it is coded in python please see here.

Putting everything 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.eam.alloy 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.eam.alloy", 'w') 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()

See also

• See Using the Zhou_AlCu.eam.alloy file within LAMMPS for details of how to use the tabulation file with LAMMPS.

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 by using the writeTABEAM(). See the tabulation script for this example: eam_tabulate_example2b.py.

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

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

  nr = 2000
  dr = 0.003

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

See also

• See the object oriented version of this example Example 2b: Tabulate Al-Cu Alloy Potentials Using TABEAM_EAMTabulation for DL_POLY.

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 example uses the writeSetFLFinnisSinclair() function to produce files supported by the LAMMPS pair_style eam/fs command.

The potential model and definition of potential objects is detailed in Example 3b: Tabulate Al-Fe Finnis-Sinclair Potentials Using TABEAM_FinnisSinclair_EAMTabulation for DL_POLY which uses a tabulation class but is otherwise very similar to this example. 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:

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", "w") as outfile:
    writeSetFLFinnisSinclair(
      nrho, drho,
      nr, dr,
      eamPotentials,
      pairPotentials,
      outfile)

The full tabulation script can be downloaded as eam_tabulate_example3a.py.

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", "w") 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", "w") as outfile:
    writeTABEAMFinnisSinclair(
      nrho, drho,
      nr, dr,
      eamPotentials,
      pairPotentials,
      outfile)

Footnotes

[1]	A.P. Sutton, and J. Chen, “Long-range Finnis-Sinclair potentials”, Philos. Mag. Lett. 61 (1990) 139 doi:10.1080/09500839008206493.

[2]	Zhou, R. Johnson and H. Wadley, “Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers”, Phys. Rev. B. 69 (2004) 144113.

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