Notes on MEAM Potential

The modified embedded-atom-method (MEAM) potential extends the conventional embedded-atom-method by including directional bonding [1]. This allows the application to both metallic and covalent systems, as well as mixed systems.

The energy function for the Modified Embedded Atom Method as presented e.g. in [2] is given as the sum of some pair term and some embedding term as follows:

\[E= \sum_i \left\{ F_i(\bar{\rho}_i) + \frac{1}{2} \sum_{j\ne i} \Phi_{ij}(r_{ij})\right\}\]

The embedding function is of the form \(F_i(\bar{\rho}_i) = A_i E_i^0\bar{\rho}_i \log \bar{\rho}_i\) for \(\bar{\rho}_i > 0\) and for \(\bar{\rho}_i \le 0\) is \(F_i(\bar{\rho}_i) = 0\) if embedding_negative is 0 or \(F_i(\bar{\rho}_i) = -A_i E_i^0 \bar{\rho}_i\) if it is 1, where \(A_i\) and \(E_i^0\) are element dependent parameters. The background electron density \(\bar{\rho}_i\) is given by:

\[\bar{\rho}_i = \frac{\bar{\rho}_i^{(0)}}{\rho_i^0} G_i (\Gamma_i)\]

where

\[\Gamma_i = \sum_{k=1}^3 t_i^{(k)}\left(\frac{\bar{\rho}_i^{(k)}}{\bar{\rho}_i^{(0)}}\right)^2\]

and \(G_i\) is one of

\[\begin{split}G_i(\Gamma) = \begin{cases} \sqrt{1+\Gamma} & \texttt{gamma}=0\\ \exp(\Gamma / 2) & \texttt{gamma}=1\\ \text{sign}(1+\Gamma) \sqrt{|1+\Gamma|} & \texttt{gamma}=2\\ \frac{2}{1+\exp(-\Gamma)} & \texttt{gamma}=3\\ \sqrt{1+\Gamma} & \texttt{gamma}=4 \end{cases}\end{split}\]

With the element dependent constant \(\rho_i^0\), the number of nearest neighbors in the reference structure \(Z_{i0}\) and \(\Gamma_i^\text{ref}\) \(\Gamma\) evaluated in the reference structure, the composition dependent density scaling can be chosen to be either

\[\begin{equation} \label{meamNormalScaling} \rho_i^0 = \rho_{i0}Z_{i0}G_i(\Gamma_i^\text{ref}) \end{equation}\]

or \(\rho_i^0 = \rho_{i0} Z_{i0}\). Note that for \gamma=0 and \gamma=2 \(G_i(\Gamma_i^\text{ref})\) is always set to 1. The partial electron densities are defined as

\[\begin{split}\begin{eqnarray*} \bar{\rho}_i^{(0)} &=& \sum_{j\ne i} \rho_j^{a(0)}(r_{ij})S_{ij} \\ (\bar{\rho}_i^{(1)})^2 &=& \sum_{\alpha = 1}^3 \left[\sum_{j \ne i} \rho_j^{a(1)} \frac{r_{ij\alpha}}{r_{ij}} S_{ij} \right]^2 \\ (\bar{\rho}_i^{(2)})^2 &=& \sum_{\alpha =1}^3 \sum_{\beta =1}^3 \left[\sum_{j\ne i} \rho_j^{a(2)} \frac{r_{ij\alpha}r_{ij\beta}}{r_{ij}^2} S_{ij} \right]^2 - \frac{1}{3} \left[ \sum_{j\ne i} \rho_j^{a(2)}S_{ij}\right]^2\\ (\bar{\rho}_i^{(3)})^2 &=& \sum_{\alpha =1}^3 \sum_{\beta =1}^3 \sum_{\gamma =1}^3 \left[\sum_{j\ne i} \rho_j^{a(3)} \frac{r_{ij\alpha}r_{ij\beta}r_{ij\gamma}}{r_{ij}^3} S_{ij}\right]^2 - \frac{3}{5}\sum_{\alpha =1}^3\left[\sum_{j \ne i} \rho_j^{a(3)} \frac{r_{ij\alpha}}{r_{ij}} S_{ij}\right]^2 \end{eqnarray*}\end{split}\]

with the atomic electron densities given by

\[\rho_i^{a(k)}(r_{ij}) = \rho_{i0} \exp \left[-\beta_i^{(k)} \left(\frac{r_{ij}}{r_i^0} - 1 \right)\right]\]

Finally, the average weighting factors are given by

\[\begin{split}t_i^{(k)} = \begin{cases} \frac{1}{\bar{\rho}_i^{(0)}} \sum_{j\ne i} t_{0,j}^{(k)} \rho_j^{a(0)} S_{ij} & \text{wf-mixing=0} \\ \frac{ \sum_{j\ne i} t_{0,j}^{(k)} \rho_j^{a(0)} S_{ij}} { \sum_{j\ne i} \left(t_{0,j}^{(k)}\right)^2 \rho_j^{a(0)} S_{ij}}& \text{wf-mixing=1} \\ t_{0,i}^{(k)} & \text{wf-mixing=2} \end{cases}\end{split}\]

with element dependent parameters \(t_{0,i}^{(k)}\).

The pair term \(\Phi_{ij}\) is defined as

\[\begin{split}\begin{eqnarray} \Phi_{ij}(r_{ij}) &=& \bar{\Phi}_{ij}(r_{ij}) S_{ij}\\ \bar{\Phi}_{ij}(r_{ij}) &=& \frac{1}{Z_{ij0}} [2E_i^u(r_{ij}) - F_i(\hat{\rho}_i (r_{ij})) - F_j(\hat{\rho}_j(r_{ij})] \end{eqnarray}\end{split}\]

where \(Z_{ij0}\) is the number of nearest neighbors in the reference structure, \(F_i(\hat{\rho}_{ij})\) the embedding function evaluated in the reference structure and

\[\begin{split}E_i^u(r_{ij}) = \begin{cases} -E_{ij}^c \left(1+a_{ij}^{*} + a_{ij}^{(3)} {a_{ij}^{*}}^3 \frac{r^0_{ij}}{r_{ij}}\right) \exp(-a_{ij}^*) & \text{erose=0} \\ -E_{ij}^c \left(1+a_{ij}^* + \left(-\text{attrac}+\frac{\text{repuls}}{r_{ij}}\right){a_{ij}^*}^3\right)\exp(-a_{ij}^*) & \text{erose=1}\\ -E_{ij}^c\left(1+a_{ij}^* + a_{ij}^{(3)}{a_{ij}^*}^3\right)\exp(-a_{ij}^*) & \text{erose=2} \end{cases}\end{split}\]

with

\[a_{ij}^*(r_{ij}) = \alpha_{ij} \left(\frac{r_{ij}}{r^0_{ij}} - 1\right)\]

and

\[\begin{split}a_{ij}^{(3)} = \begin{cases} \text{repuls} & a_{ij}^* < 0\\ \text{attrac} & a_{ij}^* \ge 0 \end{cases}\end{split}\]

Some reference structures use a more complex form of \(\Phi_{ij}(r_{ij})\). For small \(r_{ij}\), this pair term is blended with a ZBL potential by default.newline Last, the screening function \(S_{ij}\) is defined by

\[\begin{split}\begin{eqnarray*} S_{ij}&=&\bar{S}_{ij} f_c\left(\frac{r_c - r_{ij}}{\Delta r}\right)\\ \bar{S}_{ij} &=& \prod_{k\ne i,j} S_{ikj}\\ S_{ikj} &=& f_c\left(\frac{C_{ikj}-C_{\text{min},ikj}}{C_{\text{max},ikj} -C_{\text{min},ikj}}\right)\\ C_{ikj}&=& 1+2\frac{r_{ij}^2 r_{ik}^2 + r_{ij}^2r_{jk}^2-r_{ij}^4}{r_{ij}^4-(r_{ik}^2-r_{jk}^2)^2}\\ f_c(x)&=&\begin{cases} 1 & x \ge 1 \\ [1-(1-x)^4]^2 & 0<x<1\\ 0 & x \le 0 \end{cases} \end{eqnarray*}\end{split}\]

where \(r_c\) is the global cutoff radius, \(\Delta r\) is the length of a smoothing region for \(r\) near \(r_c\) and \(C_\text{min}\) and \(C_\text{max}\) are element-triple dependent parameters.

Note that the cutoff radius for this term needs to be larger than \(r_c\), depending on \(C_{\text{max}}\), because even triples with some particle distance \(> r_c\) may have a screening effect.

The parameter lattice defines the reference structure for the given particle type. The parameters alpha, beta_k, re, Ec, scaling_factor, weighting_factor_k, rho, \gamma, attrac and repuls correspond to the variables \(\alpha\), \(\beta^{(k)}\), \(r^0\), \(E^c\), \(A\), \(t_0^{(k)}\), \(\rho_{0}\), attrac and repuls from the previous section respectively. nn2 enables the second nearest neighbor formulation as described in cite{lee2000second}, zbl enables ZBL blending for small distances if set to 1.

re and Ec correspond to \(r^0\) and \(E^c\). Note that all element_data entries must be given before the first element_pair entry.newline If augment_1st is set to 1, we set

\[t_{0,i}^{(1)} = t_{0,i}^{(1)} + \frac{3}{5} t_{0,i}^{(3)}\]

for all \(i\). This is usually only needed for older parameter sets, modern ones already include this correction. The default density_scaling uses \(\rho_i^0=\rho_{i0}Z_{i0}G_i(\Gamma_i^{ref})\), setting \(G_i(\Gamma_i^{ref})=1\) uses \(\rho_i^0=\rho_{i0}Z_{i0}\) instead.

The latticeType parameter must be one of the following strings:

• dia = diamond

• fcc = face centered cubic

• bcc = body centered cubic

• dim = dimer

• b1 = rock salt

• hcp = hexagonal close-packed

• c11 = MoSi2 structure

• l12 = Cu3Au structure

• b2 = CsCl structure

[1]

M. I. Baskes. Modified embedded-atom potentials for cubic materials and impurities. Phys. Rev. B, 46:2727–2742, Aug 1992. URL: http://link.aps.org/doi/10.1103/PhysRevB.46.2727, doi:10.1103/PhysRevB.46.2727.

[2]

PM Gullet, G Wagner, and A Slepoy. Numerical tools for atomistic simulations. SANDIA Report, 8782:2003, 2003.