Many body models

Models that use the Embedded Atom Method (EAM) can be tabulated using potable. Embedded atom models take the general form

(11)\[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})\]

• Where:
  • \(\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.

In order to support the description of EAM when compared to pair-potential models, additional sections in the input file are required. These are:

• [EAM-Density]: defines EAM density functions.
• [EAM-Embed]: describe the model’s embedding functions.

As before, the pairwise section of the forcefield is specified in the [Pair] section of the input file.

The EAM sections of the input are defined in much the same way as the [Pair] section. Both [EAM-Density], [EAM-Embed] allow the use of multi-range potential definitions and potential modifiers.

As for pair only models, many-bodied force fields can specify [Potential-Form] and [Table-Form] sections can be provided if custom functional forms are required in [EAM-Density], [EAM-Embed] or [Pair].

[EAM-Density]

The density functions for embedded atom models are specified in this section. For the EAM form shown in equation (11) entries in this section take the form:

SPECIES : POTENTIAL_FORM PARAM_1 PARAM_2 ... PARAM_N

• Where:
  • SPECIES species for which density should be calculated
  • POTENTIAL_FORM PARAM_1 PARAM_2 ... PARAM_N potential form definition.

Density functions are tabulate in \(r_{ij}\) space, therefore their extent and resolution are controlled by the the dr, nr and cutoff fields in the [Tabulation] section in the same way as for [Pair] potentials.

Note

The SPECIES label can take the form ALPHA->BETA for Finnis-Sinclair tabulations. See Finnis Sinclair Models, below and [EAM-Density] in the reference section .

Example

A density function for silver may look something like this:

[EAM-Density]
Ag : as.exponential 4681.013008649 -6

This is taken from Sutton Ag EAM Example. As you can determine from the parameters given there, and just for fun, this could have been defined using potential-modifiers as this, which makes the original parameters self evident:

[EAM-Density]
Ag : pow(
         product(as.constant 4.09,
                 pow(
                     as.polynomial 0 1,
                     as.constant -1)),
         as.constant 6)

[EAM-Embed]

Embedding functions are defined in this section.

Entries have the following form:

SPECIES : POTENTIAL_FORM PARAM_1 PARAM_2 ... PARAM_N

Where:

• SPECIES is atomic type at which the surrounding electron density will be embedded using the specified potential form.
• POTENTIAL_FORM PARAM_1 ... : embedding functions instantiate potential forms in the same way as in the [Pair] section.

Note

Embedding functions are tabulated using rho values. The resolution and extent of functions in rho are defined by drho, nrho and cutoff_rho in the [Tabulation] section.

Example

This shows the embedding function used in the Sutton Ag EAM Example for Ag:

[EAM-Embed]
Ag : product(as.constant 2.5415e-3, as.sqrt -144.41)

Note the use of the product() modifier to apply a constant multiplication factor to the square root embedding function.

Standard EAM Examples

More complete examples of EAM tabulation are listed in the following table:

Example	Description
Sutton Ag EAM Example	An example of how to tabulate a single component EAM potential for Ag, to use in LAMMPS

Finnis-Sinclair Style EAM Models

A variation of the standard EAM is supported allowing different density functions to be specified for each pair of species. Before looking at this let’s have another look at original definition of the EAM given in (11):

\[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})\]

Here the many body term is:

\[F_\alpha \left( \sum_{j \neq i} \rho_\beta(r_{ij}) \right)\]

From this it can be seen that for any atom type \(\beta\) surrounding atom \(i\), the same density function is used for \(\beta\), no matter the species (\(\alpha\)) of the central atom. So in the standard EAM, the density due to a \(B\) atom neighbouring an \(A\) atom would be calculatd by \(\rho_B (r_{ij})\). Similarly a \(B\) atom next to another \(B\) atom would also have its density calculated using the same \(\rho_B (r_{ij})\) function.

By comparison, the EAM variant (referred to as Finnis-Sinclair by LAMMPS) has the following form:

(13)\[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})\]

The difference is subtle, but has important implications for the expressiveness of the potential model. The density function now becomes \(\rho_{\alpha\beta}(r_{ij})\) meaning it is now specific the types of the interacting species. So the density due to a \(B\) atom around an \(A\) atom would be given by \(\rho_{AB} (r_{ij})\) whilst a different function \(\rho_{BB} (r_{ij})\) would be used for \(B\) atoms around \(A\) atoms.

The potable tool allows this concept to be expressed by using a slightly different style of [EAM-Density] section. To define the two density functions described in the previous paragraph this would look like this:

[EAM-Density]
A->B : DENSITY_AB
B->B : DENSITY_BB

Where:

• DENSITY_AB and DENSITY_BB would be the two potential-forms for the ‘B surrounding A’ and ‘B surrounding B’ density functions.

Summarising, the [EAM-Density] section of Finnis-Sinclair style potable files has contains SPECIES labels of the form:

ALPHA->BETA

Where:

• ALPHA is the central atom species for the density function.
• BETA is the surrrounding atom species for the density function.

Finnis-Sinclair style models can be used with the following tabulation targets:

• excel_eam_fs
• DL_POLY_EAM_fs
• setfl_fs

Example

The following example aims to demonstrate the difference between the standard and Finnis-Sinclair potential models and how to tabulate them. For simplicity a ‘toy’ atomic configuration will be used, this is illustrated in Fig. 7 and can be downloaded as a LAMMPS file: toy_structure.lmpstruct (atom_style charge).

Fig. 7 The atomic configuration used in this example (viewed down z-axis).

The coordination environments of the A and B atoms in this structure are as follows:

Table 4 Coordination environment of the different atoms.

Central Atom	Surrounding Atoms	\(r_{ij}\)
A	4 × B	2.0
B	1 × A	2.0
	1 × B	4.0
	2 × B	\(2 \sqrt{2}\)

Standard EAM

For our toy example, let’s define the density due to an A atom as \(\rho_A(r_{ij}) = 2r_{ij}\) and that due to B as \(\rho_B(r_{ij}) = 3r_{ij}\). Using these we can extend our table to calculate the density expected around each type of atom:

Table 5 Density calculation

Central Atom	Surrounding Atoms	\(r_{ij}\)	\(\rho_A(r_{ij}) = 2r_{ij}\)	\(\rho_B(r_{ij}) = 3r_{ij}\)	Total \(\rho\)
A	4 × B	2.0		4 × 3 × 2.0 = 24.0	24.0
B	1 × A	2.0	1 × 2 × 2.0 = 4.0		4.0 + 12.0 + 16.971 = 32.971
	1 × B	4.0		1 × 3 × 4.0 = 12.0
	2 × B	\(2 \sqrt{2}\)		2 × 3 × \(2 \sqrt{2}\) = 16.971

Looking at the calculation in Table 5 we can see that the density around an A atom is 24.0 and that around B is 32.971.

Let’s confirm that this is the case by running a LAMMPS calculation that reproduces the hand calculation given in the table using the following potable input (standard_eam.aspot)

[Tabulation]
target : setfl
cutoff = 5.0
dr = 0.1
cutoff_rho = 50.0
drho = 0.1

[Species]
A.atomic_mass = 1
A.atomic_number = 1
B.atomic_mass = 2
B.atomic_number = 2

[EAM-Embed]
A = as.polynomial 0 1
B = as.zero

[EAM-Density]
A = as.polynomial 0 2
B = as.polynomial 0 3

[Pair]

Notes:

• The [Species] section has been included so we can use A and B as species labels rather than proper element names.
• The [Pair] section has been left empty as we haven’t defined any pair potentials to use with this model.

The [EAM-Density] section uses the as.polynomial potential form to give our \(\rho_A(r_{ij}) = 2r_{ij}\) and \(\rho_B(r_{ij}) = 3r_{ij}\) density functions:

[EAM-Density]
A = as.polynomial 0 2
B = as.polynomial 0 3

More unusual is the use of the as.polynomial and as.zero potential forms in the [EAM-Embed] section. Our example is in no way trying to reproduce the physics of atomic bonding but is instead showing how the EAM and its tabulations works. As a result the [EAM-Embed] section looks like this:

[EAM-Embed]
A = as.polynomial 0 1
B = as.zero

Here the embedding function for A has been set to as.polynomial 0 1. This is useful for debugging as this is an identity function meaning that the ‘energy’ calculated by LAMMPS for each A atom will instead be its summed density.

The embedding function for B has been set to as.zero meaning the density around these atoms will not contribute to the value output by LAMMPS. We will turn the B embedding function back on later in the example, however for our first run the potable tabulation should result in the density surrounding the single A atom in Fig. 7 being produced. According to the calculation in Table 5 this should be 24.0.

Run potable on standard_eam.aspot to produce a table file named standard_eam.eam

potable standard_eam.aspot standard_eam.eam

Now download the structure file (toy_structure.lmpstruct) and the following LAMMPS input script (standard_evaluate.lmpin) into the same directory as your table:

units metal
atom_style charge

read_data toy_structure.lmpstruct

pair_style eam/alloy
pair_coeff * * standard_eam.eam A B

run 0

print ENERGY:$(pe)

Now let’s run this through LAMMPS to evaluate our table file for the example structure:

lammps -in standard_evaluate.lmpin  -log standard_evaluate.lmpout

This should produce output similar to this:

LAMMPS (3 Mar 2020)
units metal
atom_style charge

read_data toy_structure.lmpstruct
  orthogonal box = (0 0 0) to (20 20 20)
  1 by 1 by 1 MPI processor grid
  reading atoms ...
  5 atoms
  read_data CPU = 0.004225 secs

pair_style eam/alloy
pair_coeff * * standard_eam.eam A B

run 0
WARNING: No fixes defined, atoms won't move (../verlet.cpp:52)
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 7.1
  ghost atom cutoff = 7.1
  binsize = 3.55, bins = 6 6 6
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair eam/alloy, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d/newton
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.436 | 3.436 | 3.436 Mbytes
Step Temp E_pair E_mol TotEng Press 
       0            0           24            0           24   -1602.1765 
Loop time of 1e-06 on 1 procs for 0 steps with 5 atoms

100.0% CPU use with 1 MPI tasks x no OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0          | 0          | 0          |   0.0 |  0.00
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 0          | 0          | 0          |   0.0 |  0.00
Output  | 0          | 0          | 0          |   0.0 |  0.00
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 1e-06      |            |       |100.00

Nlocal:    5 ave 5 max 5 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:    0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:    10 ave 10 max 10 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 10
Ave neighs/atom = 2
Neighbor list builds = 0
Dangerous builds = 0

print ENERGY:$(pe)
print ENERGY:24.000000000000014211
ENERGY:24.000000000000014211
Total wall time: 0:00:00

As you can see from the penultimate line we obtain our expected value of 24.0 for the density of the A atom:

ENERGY:24.000000000000014211

Now change the [EAM-Embed] section as follows, to turn off the A density and turn on the B density:

[EAM-Embed]
A = as.zero
B = as.polynomial 0 1

Running the same process as above we get the following ENERGY line in our output:

ENERGY:131.88225099390862738

Remembering that we have four B atoms in our system, each with an expected density of 32.971 then it can be seen that this equates well to what we would expect: 4×32.971 = 131.884.

Finnis-Sinclair

We will now extend this example to use different density functions for the AB and BA interactions. The following functions will be used:

• \(\rho_{AB} (r_{ij}) = 3 r_{ij}\)
• \(\rho_{BB} (r_{ij}) = 5 r_{ij}\)
• \(\rho_{BA} (r_{ij}) = 2 r_{ij}\)

Using these functions will produce a density around the A atom that is the same as in the standard EAM example above. This is because the \(\rho_{AB} (r_{ij})\) function is the same as the \(\rho_{A} (r_{ij})\) used earlier. However, the function for calculating the B density surrounding a B atom is now different from the earlier example (\(\rho_{BB} (r_{ij}) = 5 r_{ij}\)). As the \(\rho_{BA} (r_{ij})\) function also matches the equivalent \(\rho_{A}\) function from earlier, it is only the change introduced by the new \(\rho_{BB} (r_{ij}) = 5 r_{ij}\) function for B-B pairs, that will change the LAMMPS output for this example. Due to the single A atom there are no density contribution from A-A pairs in this system.

The density calculation for the Finnis-Sinclair model can now be performed by hand as shown in Table 6.

Table 6 Density calculation

Central Atom	Surrounding Atoms	\(r_{ij}\)	\(\rho_{AB}(r_{ij}) = 3r_{ij}\)	\(\rho_{BA}(r_{ij}) = 2r_{ij}\)	\(\rho_{BB}(r_{ij}) = 5r_{ij}\)	Total \(\rho\)
A	4 × B	2.0	4 × 3 × 2.0 = 24.0			24.0
B	1 × A	2.0		1 × 2 × 2.0 = 4.0		4.0 + 20.0 + 28.284 = 52.284
	1 × B	4.0			1 × 5 × 4.0 = 20.0
	2 × B	\(2 \sqrt{2}\)			2 × 5 × \(2 \sqrt{2}\) = 28.284

This can be described in the potable format as follows (finnis_sinclair_eam.aspot):

[Tabulation]
target : setfl_fs
cutoff = 5.0
dr = 0.1
cutoff_rho = 50.0
drho = 0.1

[Species]
A.atomic_mass = 1
A.atomic_number = 1
B.atomic_mass = 2
B.atomic_number = 2

[EAM-Embed]
A = as.zero
B = as.polynomial 0 1

[EAM-Density]
A->B = as.polynomial 0 3
B->A = as.polynomial 0 2
B->B = as.polynomial 0 5

[Pair]

Note that the [Tabulation] section now specifies a Finnis-Sinclair compatible tabulation target:

[Tabulation]
target : setfl_fs

The important differences between this file and the standard EAM example can be found in the [EAM-Density]

[EAM-Density]
A->B = as.polynomial 0 3
B->A = as.polynomial 0 2
B->B = as.polynomial 0 5

The first entry in this section defines the \(\rho_{AB}(r_{ij}) = 3r_{ij}\) function and it should be apparent how the remaining density functions are specified.

Tabulate the finnis_sinclair_eam.aspot file:

potable finnis_sinclair_eam.aspot finnis_sinclair.eam.fs

Copy the toy_structure.lmpstruct structure file and the following LAMMPs input file (finnis_sinclair_evaluate.lmpin) to the same directory as your tabulation:

units metal
atom_style charge

read_data toy_structure.lmpstruct

pair_style eam/fs
pair_coeff * * finnis_sinclair.eam.fs A B

run 0

print ENERGY:$(pe)

Running LAMMPS to obtain the total B density:

lammps -in finnis_sinclair_evaluate.lmpin -log finnis_sinclair_evaluate.lmpout

gives the following output:

ENERGY:209.13708498984715334

Referring back to Table 6 and remembering that there are four B atoms in the system the value expected from the hand calculation was 52.284 × 4 = 209.136. This very closely matches the value obtained from LAMMPS.

ADP Style EAM Models

Added in: 0.4.0

The potable tool allows EAM models using the angular dependent potential (ADP) extension to be tabulated.

(14)\[\begin{split}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})+ \frac{1}{2} \sum_s (\mu_i^s)^2 + \frac{1}{2} \sum_{s,t} (\lambda_i^{st})^2 - \frac{1}{6} \nu_i^2 \\ \mu_i^s & = \sum_{j\neq i}u_{\alpha\beta}(r_{ij})r_{ij}^s\\ \lambda_i^{st} & = \sum_{j\neq i}w_{\alpha\beta}(r_{ij})r_{ij}^sr_{ij}^t\\ \nu_i & = \sum_s\lambda_i^{ss}\end{split}\]

ADP extends the standard EAM (11) with additional dipole (\(\mu_i^s\)) and quadrupole terms (\(\lambda_i^{st}\)). The dipoles and quadrupoles are defined via \(u_{\alpha\beta}(r_{ij})\) and \(w_{\alpha\beta}(r_{ij})\) functions respectively.

As the \(u\) and \(w\) functions are specified for pairs of species they are defined in a potable file in the same way as pair potentials (see [Pair] section). The \(u_{\alpha\beta}(r_{ij})\) functions are given in the section named [EAM-ADP-Dipole] and the \(w_{\alpha\beta}(r_{ij})\) functions appear in [EAM-ADP-Quadrupole]. Other than these new sections ADP models are defined in the same way as the standard EAM (see Many body models). Consequently the potable file for an ADP model minimally contains the sections:

• [Tabulation]
  • For ADP the target value should be set as eam_adp.
• [EAM-Density]
• [EAM-Embed]
• [Pair]
• [EAM-ADP-Dipole]
• [EAM-ADP-Quadrupole]