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. This is done by using the EAM_Tabulation objects from the atsim.potentials.eam_tabulation module:

Class	Format	Simulation Code	Example
`SetFL_EAMTabulation`	`setfl`, pair_style eam/alloy	LAMMPS	Example 1: Using SetFL_EAMTabulation to Tabulate Ag Potential for LAMMPS Example 2a: Tabulate Al-Cu Alloy Potentials Using SetFL_EAMTabulation for LAMMPS
`SetFL_FS_EAMTabulation`	pair_style eam/fs	LAMMPS	Example 3a: Tabulate Al-Fe Finnis-Sinclair Potentials Using SetFL_FS_EAMTabulation for LAMMPS
`TABEAM_EAMTabulation`	`TABEAM`	DL_POLY	Example 2b: Tabulate Al-Cu Alloy Potentials Using TABEAM_EAMTabulation for DL_POLY
`TABEAM_FinnisSinclair_EAMTabulation`	`EEAM TABEAM`	DL_POLY	Example 3b: Tabulate Al-Fe Finnis-Sinclair Potentials Using TABEAM_FinnisSinclair_EAMTabulation for DL_POLY
`Excel_EAMTabulation`	`.xlsx`
`Excel_FinnisSinclair_EAMTabulation`	`.xlsx`

Even though the use of EAM_Tabulation objects is preferred a legacy procedural interface is also provided. This is described here: Procedural EAM Tabulation.

Examples¶

Example 1: Using `SetFL_EAMTabulation` to Tabulate Ag Potential for LAMMPS¶

This example shows how to use the SetFL_EAMTabulation class 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_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, Potential
from atsim.potentials.eam_tabulation import SetFL_EAMTabulation


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:

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

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

Now all the components of the model have been defined a table file can be created in the setfl format. Before doing this, it is necessary to choose appropriate density and separation cut-offs together with the number of rows in the density/pair functions (nr) and embedding function (nrho) respectively:

Here 50000 density values will be tabulated to a cutoff of 50.0.

The pair-potential cut-off and the maximum \(r_{ij}\) value for the density function is 12 Å and both will have 12000 rows.

An instance of atsim.potentials.eam_tabulation.SetFL_EAMTabulation is created with the EAMPotential and Potential objects. This object is then used to tabulate the Ag potential by calling the write() method with the Ag.eam.alloy file object:

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

    cutoff_rho = 50.0
    nrho = 50000

    cutoff = 12.0
    nr = 12000

    tabulation = SetFL_EAMTabulation(
        pair_potentials,
        eam_potentials,
        cutoff, nr,
        cutoff_rho, nrho)

    with open("Ag.eam.alloy", 'w') as outfile:
        tabulation.write(outfile)

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

#! /usr/bin/env python
import math

from atsim.potentials import EAMPotential, Potential
from atsim.potentials.eam_tabulation import SetFL_EAMTabulation


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
    eam_potentials = [EAMPotential("Ag", 47, 107.8682, embed, density)]
    pair_potentials = [Potential('Ag', 'Ag', pair_AgAg)]

    cutoff_rho = 50.0
    nrho = 50000

    cutoff = 12.0
    nr = 12000

    tabulation = SetFL_EAMTabulation(
        pair_potentials,
        eam_potentials,
        cutoff, nr,
        cutoff_rho, nrho)

    with open("Ag.eam.alloy", 'w') as outfile:
        tabulation.write(outfile)


if __name__ == "__main__":
    main()

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

python eam_example1.py

Using the `Ag.eam.alloy` 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.alloy file you created previously. This describes a single FCC cell with a wildly inaccurate lattice parameter:

Title


atoms
atom types
0 5.000000 xlo xhi
0 5.000000 ylo yhi
0 5.000000 zlo zhi
000000 0.000000 0.000000 xy xz yz



Masses

107.86820000000000163709 #Ag

Atoms

0 1 0.000000 0.000000 0.000000 0.000000
0 1 0.000000 2.500000 2.500000 0.000000
0 1 0.000000 0.000000 2.500000 2.500000
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 example_eam_alloy_minimize.lmpin:

units metal
boundary p p p

atom_style full
read_data fcc.lmpstruct

pair_style eam/alloy
pair_coeff * * Ag.eam.alloy Ag

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/alloy command tells LAMMPS to use the EAM and expect pair_coeff commands mapping atom types to particular table files:

pair_style eam/alloy

The following pair_coeff directive indicates that the interaction between atom-type 1 (Ag) with itself should use the setfl formatted file contained within Ag.eam.alloy. The Ag label at the end of the line indicates that atom-type 1 should be associated with this label in the table file:

pair_coeff * * Ag.eam.alloy Ag

The example can then be run by invoking LAMMPS:

lammps -in example_eam_alloy_minimize.lmpin

Example 2a: Tabulate Al-Cu Alloy Potentials Using `SetFL_EAMTabulation` 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]. In comparison to the previous example this example contains density and embedding functions for multiple elements and includes pair-potentials specific to pairs of interacting species. The eam_example2a.py file 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_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

    # 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
    eam_potentials = [
        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
    pair_potentials = [
        Potential('Al', 'Al', pair_AlAl),
        Potential('Cu', 'Cu', pair_CuCu),
        Potential('Al', 'Cu', pair_AlCu)]

    return eam_potentials, pair_potentials

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:

    eam_potentials = [
        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
    pair_potentials = [
        Potential('Al', 'Al', pair_AlAl),
        Potential('Cu', 'Cu', pair_CuCu),
        Potential('Al', 'Cu', pair_AlCu)]

Now we have all the objects required for SetFL_EAMTabulation. The next excerpt calls makeObjects() to get the EAM and pair-potential objects before creating the tabulation object, and invoking its write() method to write the data into a file called Zhou_AlCu.eam.alloy:

def main():
    eam_potentials, pair_potentials = makePotentialObjects()

    # Perform tabulation
    # Make tabulation
    cutoff_rho = 100.0
    nrho = 2000

    cutoff = 6.0
    nr = 2000

    tabulation = SetFL_EAMTabulation(
        pair_potentials,
        eam_potentials,
        cutoff, nr,
        cutoff_rho, nrho
    )

    with open("Zhou_AlCu.eam.alloy", 'w') as outfile:
        tabulation.write(outfile)

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

#! /usr/bin/env python

import math

from atsim.potentials import EAMPotential, Potential
from atsim.potentials.eam_tabulation import SetFL_EAMTabulation


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
    eam_potentials = [
        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
    pair_potentials = [
        Potential('Al', 'Al', pair_AlAl),
        Potential('Cu', 'Cu', pair_CuCu),
        Potential('Al', 'Cu', pair_AlCu)]

    return eam_potentials, pair_potentials


def main():
    eam_potentials, pair_potentials = makePotentialObjects()

    # Perform tabulation
    # Make tabulation
    cutoff_rho = 100.0
    nrho = 2000

    cutoff = 6.0
    nr = 2000

    tabulation = SetFL_EAMTabulation(
        pair_potentials,
        eam_potentials,
        cutoff, nr,
        cutoff_rho, nrho
    )

    with open("Zhou_AlCu.eam.alloy", 'w') as outfile:
        tabulation.write(outfile)


if __name__ == '__main__':
    main()

Using the `Zhou_AlCu.eam.alloy` file within LAMMPS¶

Within LAMMPS the setfl files generated by SetFL_EAMTabulation 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.eam.alloy Al

Likewise if a copper system was being simulated:

pair_style eam/alloy
pair_coeff * * Zhou_AlCu.eam.alloy 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.eam.alloy Al Cu

Or if Cu was 1 and Al 2:

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

Example 2b: Tabulate Al-Cu Alloy Potentials Using `TABEAM_EAMTabulation` 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_example2b.py.

The EAMPotential and Potential lists are created in exactly the same way as Example 2a, however rather than creating an instance of SetFL_EAMTabulation in the main() function it is modified to use the DL_POLY specific TABEAM_EAMTabulation class instead and to write into a file named TABEAM. The main() function of eam_example2b.py is now given:

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

    # Perform tabulation
    # Make tabulation
    cutoff_rho = 100.0
    nrho = 2000

    cutoff = 6
    nr = 2000

    tabulation = TABEAM_EAMTabulation(
        pairPotentials, eamPotentials, cutoff, nr, cutoff_rho, nrho)

    with open("TABEAM", 'w') as outfile:
        tabulation.write(outfile)

Using the `TABEAM` file with DL_POLY ¶

Running eam_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 `SetFL_FS_EAMTabulation` 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 SetFL_FS_EAMTabulation class.

Using `SetFL_FS_EAMTabulation` 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_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_example3a.py.

These functions are used within the main() function of the eam_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.
    cutoff_rho = 300.0
    nrho = 10000

    cutoff = 6.5
    nr = 10000

    tabulation = SetFL_FS_EAMTabulation(
        pairPotentials, eamPotentials, cutoff, nr, cutoff_rho, nrho)

    with open("Mendelev_Al_Fe.eam.fs", "w") as outfile:
        tabulation.write(outfile)

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

Next, the EAMPotential objects for Al and Fe are instantiated. This is where a Finnis-Sinclair model differs from the standard EAM seen earlier. Instead of a single density callable being passed to the EAMPotential constructor a dictionary of density functions is passed instead (see highlighted lines):

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

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

Finally, having defined the list of EAMPotential instances these are passed to the constructor of SetFL_FS_EAMTabulation, in this case writing the data to Mendelev_Al_Fe.eam.fs in the current directory:

    tabulation = SetFL_FS_EAMTabulation(
        pairPotentials, eamPotentials, cutoff, nr, cutoff_rho, nrho)

    with open("Mendelev_Al_Fe.eam.fs", "w") as outfile:
        tabulation.write(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 `TABEAM_FinnisSinclair_EAMTabulation` 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_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.
    cutoff_rho = 300.0
    nrho = 10000

    cutoff = 6.5
    nr = 10000

    tabulation = TABEAM_FinnisSinclair_EAMTabulation(
        pairPotentials, eamPotentials, cutoff, nr, cutoff_rho, nrho)

    with open("TABEAM", "w") as outfile:
        tabulation.write(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 create an object of TABEAM_FinnisSinclair_EAMTabulation instead of atsim.potentials.eam_tabulation.SetFL_FS_EAMTabulation:

    tabulation = TABEAM_FinnisSinclair_EAMTabulation(
        pairPotentials, eamPotentials, cutoff, nr, cutoff_rho, nrho)

    with open("TABEAM", "w") as outfile:
        tabulation.write(outfile)

That’s it, nothing else has changed.

Using the `TABEAM` file with DL_POLY¶

Running eam_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 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]	(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.

Table of Contents

Previous topic

Next topic

This Page

Embedded Atom Method (EAM) Tabulation¶

Examples¶

Example 1: Using `SetFL_EAMTabulation` to Tabulate Ag Potential for LAMMPS¶

Model Description¶

Define the Model¶

Using the `Ag.eam.alloy` file within LAMMPS¶

Example 2a: Tabulate Al-Cu Alloy Potentials Using `SetFL_EAMTabulation` for LAMMPS¶

Model Description¶

Define the Model¶

Using the `Zhou_AlCu.eam.alloy` file within LAMMPS¶

Single Element Systems¶

Mixed Al-Cu System¶

Example 2b: Tabulate Al-Cu Alloy Potentials Using `TABEAM_EAMTabulation` for DL_POLY¶

Using the `TABEAM` file with DL_POLY ¶

Example 3a: Tabulate Al-Fe Finnis-Sinclair Potentials Using `SetFL_FS_EAMTabulation` for LAMMPS¶

Using `SetFL_FS_EAMTabulation` to Tabulate the Model¶

Using the `Mendelev_Al_Fe.eam.fs` file within LAMMPS¶

Example 3b: Tabulate Al-Fe Finnis-Sinclair Potentials Using `TABEAM_FinnisSinclair_EAMTabulation` for DL_POLY¶

Using the `TABEAM` file with DL_POLY¶

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

Embedded Atom Method (EAM) Tabulation¶

Examples¶

Example 1: Using SetFL_EAMTabulation to Tabulate Ag Potential for LAMMPS¶

Model Description¶

Define the Model¶

Using the Ag.eam.alloy file within LAMMPS¶

Example 2a: Tabulate Al-Cu Alloy Potentials Using SetFL_EAMTabulation for LAMMPS¶

Model Description¶

Define the Model¶

Using the Zhou_AlCu.eam.alloy file within LAMMPS¶

Single Element Systems¶

Mixed Al-Cu System¶

Example 2b: Tabulate Al-Cu Alloy Potentials Using TABEAM_EAMTabulation for DL_POLY¶

Using the TABEAM file with DL_POLY¶

Example 3a: Tabulate Al-Fe Finnis-Sinclair Potentials Using SetFL_FS_EAMTabulation for LAMMPS¶

Using SetFL_FS_EAMTabulation to Tabulate the Model¶

Using the Mendelev_Al_Fe.eam.fs file within LAMMPS¶

Example 3b: Tabulate Al-Fe Finnis-Sinclair Potentials Using TABEAM_FinnisSinclair_EAMTabulation for DL_POLY¶

Using the TABEAM file with DL_POLY¶

Example 1: Using `SetFL_EAMTabulation` to Tabulate Ag Potential for LAMMPS¶

Using the `Ag.eam.alloy` file within LAMMPS¶

Example 2a: Tabulate Al-Cu Alloy Potentials Using `SetFL_EAMTabulation` for LAMMPS¶

Using the `Zhou_AlCu.eam.alloy` file within LAMMPS¶

Example 2b: Tabulate Al-Cu Alloy Potentials Using `TABEAM_EAMTabulation` for DL_POLY¶

Using the `TABEAM` file with DL_POLY ¶

Example 3a: Tabulate Al-Fe Finnis-Sinclair Potentials Using `SetFL_FS_EAMTabulation` for LAMMPS¶

Using `SetFL_FS_EAMTabulation` to Tabulate the Model¶

Using the `Mendelev_Al_Fe.eam.fs` file within LAMMPS¶

Example 3b: Tabulate Al-Fe Finnis-Sinclair Potentials Using `TABEAM_FinnisSinclair_EAMTabulation` for DL_POLY¶

Using the `TABEAM` file with DL_POLY¶