index | modules | gitlab | page source

# EAM

## Introduction

The Embedded Atom Method (EAM) [1] is a classical potential which is good for modelling metals, particularly fcc materials. Because it is an equiaxial potential the EAM does not model directional bonds well. However, the Angular Dependent Potential (ADP) [2] which is an extended version of EAM is able to model directional bonds and is also included in the EAM calculator.

Generally all that is required to use this calculator is to supply a potential file or as a set of functions that describe the potential. The files containing the potentials for this calculator are not included but many suitable potentials can be downloaded from The Interatomic Potentials Repository Project at https://www.ctcms.nist.gov/potentials/

## Theory

A single element EAM potential is defined by three functions: the embedded energy, electron density and the pair potential. A two element alloy contains the individual three functions for each element plus cross pair interactions. The ADP potential has two additional sets of data to define the dipole and quadrupole directional terms for each alloy and their cross interactions.

The total energy \(E_{\rm tot}\) of an arbitrary arrangement of atoms is given by the EAM potential as

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

and

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

where \(F\) is an embedding function, namely the energy to embed an atom \(i\) in the combined electron density \(\bar\rho_i\) which is contributed from each of its neighbouring atoms \(j\) by an amount \(\rho(r_{ij})\), \(\phi(r_{ij})\) is the pair potential function representing the energy in bond \(ij\) which is due to the short-range electro-static interaction between atoms, and \(r_{ij}\) is the distance between an atom and its neighbour for that bond.

The ADP potential is defined as

\[E_\text{tot} = \sum_i F(\bar\rho_i) + \frac{1}{2}\sum_{i\ne j} \phi(r_{ij}) + \frac{1}{2} \sum_{i,\alpha} (\mu_i^\alpha)^2 + \frac{1}{2} \sum_{i,\alpha,\beta} (\lambda_i^{\alpha\beta})^2 - \frac{1}{6} \sum_i \nu_i^2\]

where \(\mu_i^\alpha\) is the dipole vector, \(\lambda_i^{\alpha\beta}\) is the quadrupole tensor and \(\nu_i\) is the trace of \(\lambda_i^{\alpha\beta}\).

The fs potential is defined as

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

where \(\alpha\) and \(\beta\) are element types of atoms. This form is similar to original EAM formula above, except that \(\rho\) and \(\phi\) are determined by element types.

## Running the Calculator

EAM calculates the cohesive atom energy and forces. Internally the potential functions are defined by splines which may be directly supplied or created by reading the spline points from a data file from which a spline function is created. The LAMMPS compatible `.alloy`, `.fs` and `.adp` formats are supported. The LAMMPS `.eam` format is slightly different from the `.alloy` format and is currently not supported.

For example:

    from ase.calculators.eam import EAM

    mishin = EAM(potential='Al99.eam.alloy')
    mishin.write_potential('new.eam.alloy')
    mishin.plot()

    slab.calc = mishin
    slab.get_potential_energy()
    slab.get_forces()

The breakdown of energy contribution from the indvidual components are stored in the calculator instance `.results['energy_components']`

## Arguments

Keyword | Description
--- | ---
`potential` | file of potential in `.eam`, `.alloy`, `.adp` or `.fs` format or file object (This is generally all you need to supply). For file object the `form` argument is required
`elements[N]` | array of N element abbreviations
`embedded_energy[N]` | arrays of embedded energy functions
`electron_density[N]` | arrays of electron density functions
`phi[N,N]` | arrays of pair potential functions
`d_embedded_energy[N]` | arrays of derivative embedded energy functions
`d_electron_density[N]` | arrays of derivative electron density functions
`d_phi[N,N]` | arrays of derivative pair potentials functions
`d[N,N], q[N,N]` | ADP dipole and quadrupole function
`d_d[N,N], d_q[N,N]` | ADP dipole and quadrupole derivative functions
`skin` | skin distance passed to NeighborList(). If no atom has moved more than the skin-distance since the last call to the `update()` method then the neighbor list can be reused. Defaults to 1.0.
`form` | the form of the potential `eam`, `alloy`, `adp` or `fs`. This will be determined from the file suffix or must be set if using equations or file object

## Additional parameters for writing potential files

The following parameters are only required for writing a potential in `.alloy`, `.adp` or `fs` format file.

Keyword | Description
--- | ---
`header` | Three line text header. Default is standard message.
`Z[N]` | Array of atomic number of each element
`mass[N]` | Atomic mass of each element
`a[N]` | Array of lattice parameters for each element
`lattice[N]` | Lattice type
`nrho` | No. of rho samples along embedded energy curve
`drho` | Increment for sampling density
`nr` | No. of radial points along density and pair potential curves
`dr` | Increment for sampling radius

## Special features

`.plot()`

  Plots the individual functions. This may be called from multiple EAM potentials to compare the shape of the individual curves. This function requires the installation of the Matplotlib libraries.

## Notes/Issues

- Although currently not fast, this calculator can be good for trying small calculations or for creating new potentials by matching baseline data such as from DFT results. The format for these potentials is compatible with LAMMPS and so can be used either directly by LAMMPS or with the ASE LAMMPS calculator interface.

- Supported formats are the LAMMPS `.alloy` and `.adp`. The `.eam` format is currently not supported. The form of the potential will be determined from the file suffix.

- Any supplied values will override values read from the file.

- The derivative functions, if supplied, are only used to calculate forces.

- There is a bug in early versions of scipy that will cause eam.py to crash when trying to evaluate splines of a potential with one neighbor such as caused by evaluating a dimer.

[1]

M.S. Daw and M.I. Baskes, Phys. Rev. Letters 50 (1983) 1285.

[2]

Y. Mishin, M.J. Mehl, and D.A. Papaconstantopoulos, Acta Materialia 53 2005 4029–4041.

## Example:

    import numpy as np
    from scipy.interpolate import InterpolatedUnivariateSpline as spline

    from ase.build import bulk
    from ase.calculators.eam import EAM

    def test_eam(testdir):
        # test to generate an EAM potential file using a simplified
        # approximation to the Mishin potential Al99.eam.alloy data

        cutoff = 6.28721

        n = 21
        rs = np.arange(0, n) * (cutoff / n)
        rhos = np.arange(0, 2, 2. / n)

        # generated from
        # mishin = EAM(potential='../potentials/Al99.eam.alloy')
        # m_density = mishin.electron_density[0](rs)
        # m_embedded = mishin.embedded_energy[0](rhos)
        # m_phi = mishin.phi[0,0](rs)

        m_density = np.array([2.78589606e-01, 2.02694937e-01, 1.45334053e-01,
                              1.06069912e-01, 8.42517168e-02, 7.65140344e-02,
                              7.76263116e-02, 8.23214224e-02, 8.53322309e-02,
                              8.13915861e-02, 6.59095390e-02, 4.28915711e-02,
                              2.27910928e-02, 1.13713167e-02, 6.05020311e-03,
                              3.65836583e-03, 2.60587564e-03, 2.06750708e-03,
                              1.48749693e-03, 7.40019174e-04, 6.21225205e-05])

        m_embedded = np.array([1.04222211e-10, -1.04142633e+00, -1.60359806e+00,
                               -1.89287637e+00, -2.09490167e+00, -2.26456628e+00,
                               -2.40590322e+00, -2.52245359e+00, -2.61385603e+00,
                               -2.67744693e+00, -2.71053295e+00, -2.71110418e+00,
                               -2.69287013e+00, -2.68464527e+00, -2.69204083e+00,
                               -2.68976209e+00, -2.66001244e+00, -2.60122024e+00,
                               -2.51338548e+00, -2.39650817e+00, -2.25058831e+00])

        m_phi = np.array([6.27032242e+01, 3.49638589e+01, 1.79007014e+01,
                          8.69001383e+00, 4.51545250e+00, 2.83260884e+00,
                          1.93216616e+00, 1.06795515e+00, 3.37740836e-01,
                          1.61087890e-02, -6.20816372e-02, -6.51314297e-02,
                          -5.35210341e-02, -5.20950200e-02, -5.51709524e-02,
                          -4.89093894e-02, -3.28051688e-02, -1.13738785e-02,
                          2.33833655e-03, 4.19132033e-03, 1.68600692e-04])

        m_densityf = spline(rs, m_density)
        m_embeddedf = spline(rhos, m_embedded)
        m_phif = spline(rs, m_phi)

        a = 4.05  # Angstrom lattice spacing
        al = bulk('Al', 'fcc', a=a)

        mishin_approx = EAM(
            elements=['Al'], embedded_energy=np.array([m_embeddedf]),
            electron_density=np.array([m_densityf]),
            phi=np.array([[m_phif]]), cutoff=cutoff, form='alloy',
            # the following terms are only required to write out a file
            Z=[13], nr=n, nrho=n, dr=cutoff / n, drho=2. / n,
            lattice=['fcc'], mass=[26.982], a=[a])

        al.calc = mishin_approx
        mishin_approx_energy = al.get_potential_energy()

        mishin_approx.write_potential('Al99-test.eam.alloy')

        mishin_check = EAM(potential='Al99-test.eam.alloy')
        al.calc = mishin_check
        mishin_check_energy = al.get_potential_energy()

        print('Cohesive Energy for Al = ', mishin_approx_energy, ' eV')

        error = (mishin_approx_energy - mishin_check_energy) / mishin_approx_energy
        print('read/write check error = ', error)

        assert abs(error) < 1e-4