Forces for the Embedded Atom Model

It is not easy to find the forces written down for embedded atom model.

The potential is always written as

U_i = F_i(\rho) + \frac{1}{2}\sum_{j \neq i} \phi(r_{ij})


\rho = \sum_{j\neq i} f_j(r_{ij})

Probably the reason for forces written down is that different fitting functions for F, \rho, \phi are used depending on the researcher.  For this reason, LAMMPS uses splines for their implementation.

As usual, the force is calculated by

\boldsymbol F_i = \nabla U_i

\varphi_{ij} =   \frac{\partial U}{\partial r_{ij}} = -\frac{\partial F(\rho_i)}{\partial \rho}\frac{\partial f_i}{\partial r_{ij}} \nabla_i r_{ij} -\frac{\partial F(\rho_j)}{\partial \rho}\frac{\partial f_j}{\partial r_{ij}}\nabla_i r_{ij}- \frac{1}{2}\frac{\partial \phi}{\partial r_{ij}}\nabla_i r_{ij}

Cai and Ye suggest the following forms will in the following forms (with derivatives added by me). Note that we have replaced r_{ij} with a more compact notation given by r.

f(r) = f_e \mathrm{exp}\left(-\chi (r- r_e) \right)

\frac{\partial f}{\partial r} = -f_e \chi \mathrm{exp}\left(-\chi (r - r_e) \right)

F(\rho) = -F_0\left[1- \mathrm{ln}\left(\frac{\rho}{\rho_e} \right)^n \right]\left(\frac{\rho}{\rho_e} \right)^n + F_1 \left(\frac{\rho}{\rho_e} \right)

\frac{\partial F}{\partial \rho} = -F_0n\left( \frac{\rho}{\rho_e} \right)^{n-1}\left(\frac{1}{\rho_e} - 1 \right) + F_0n^2 \frac{\rho^{n-1}}{\rho_e^n}\mathrm{ln}\left(\frac{\rho}{\rho_e} \right) + \frac{F_1}{\rho_e}


\phi(r) = -\alpha\left[1 + \beta\left(\frac{r}{r_a}-1 \right) \right]\mathrm{exp}\left(\beta\left(\frac{r}{r_a}-1 \right) \right)

\frac{\partial \phi}{\partial r} = - \alpha \beta \left[1 + \beta\left(\frac{r}{r_a} -1 \right)\right] \mathrm{exp}\left(\beta\left(\frac{r}{r_a} -1 \right)  \right) + -\alpha \frac{\beta}{r_a}\mathrm{exp}\left(\beta\left(\frac{r}{r_a} -1 \right)  \right)



Forces for the Embedded Atom Model

A Brief Introduction to Phase Field Modeling of Solidification

This is copy and pasted from my Master’s thesis, and is a very short introduction into phase field modeling of solidification. It is my notes on Boettingers 2002 review article Phase-Field Simulation of Solidification
A free energy functional governing solid-liquid phase transition is defined as follows
F = \int_V \left[ f(\phi, c, T) - \frac{\epsilon_c^2}{2}\vert \nabla c\vert^2 - \frac{\epsilon_\phi^2}{2}\vert \nabla \phi\vert^2 \right] \mathrm{dV}
where f, c and \phi are the free energy density, concentration and phase field, respectively with \epsilon_c and \epsilon_\phi are the associated gradient entropy coefficients. The free energy density f contains at least one minima for each phase (a two phase system is typically a double well). The time dependent changes can be derived from variational techniques in combination of mass continuity \cite{Boettinger2002} and are given by
\frac{\partial \phi}{\partial t} = -M_{\phi}\left(\frac{\partial f}{\partial \phi} - \epsilon_\phi^2 \nabla^2 \phi \right)
\frac{\partial c}{\partial t} = \nabla \cdot \left(M_{c} c(1-c) \nabla \left(\frac{\partial f}{\partial c} - \epsilon_c^2 \nabla^2 c \right)\right) \label{eq:phase_conserved}
where M_{\phi} and M_{c} are positive mobilities related to the material kinetics and the diffusion process. The difference between the two equations is due to the fact that c (in the bottom equation) is conserved while \phi (the top equation) is not. If the solid and liquid free energies are available the free energy density can be constructed as
f(\phi, c, T) = \left[(1-c)W_A + cW_B\right] g(\phi) + (1-p(\phi))f^S(c,t) + p(\phi) f^L(c,T)
where g(\phi) is a double well, W_A is the energy barrier between material A liquid phase and solid phase, W_B is the energy barrier between material B liquid and solid phase, p(\phi) is the interpolation function–typically \phi or C_1 \int g(\phi) \mathrm{d}\phi where C_1 is a constant. f^S(c,t) is the solid phase free energy and f^L(c,t) is the liquid phase free energy. For pure elements, the following simplification is frequently made. Take the solid state as the standard state f_S^A(T) = 0 and expand the difference between the liquid and solid free energies around the melting point
f_L^A(T) - f_S^A(T) = \frac{L(T_M^A - T)}{T_M^A}
where T_M^A is the melting temperature of the material. This gives the free energy
f_A(\phi,T) = W_A g(\phi) + L_A \frac{T-M^A -T}{T_M^A} p(\phi)
The free enthalpy which is used to derive the heat equation is described as
h = h_0 + C_p T + L \phi
This gives our heat and phase field equations for a pure element as follows
C_p \frac{\partial T}{\partial t} + L \frac{\partial \phi}{\partial t} = \nabla \cdot (k\nabla T)
\frac{\partial \phi}{\partial t} = M_\phi \epsilon_\phi\left(\nabla^2 \phi - \frac{2 W_A}{\epsilon_\phi} g'(\phi)\right) - \frac{M_\phi L}{T_M}(T_M - T)
The reader is directed to the review paper from Boettinger et al. for further detail and examples.

A Brief Introduction to Phase Field Modeling of Solidification

Laplacians, Gradients and Divergence

I can never remember the Laplacians, gradients and divergence formulas in spherical and cylindrical coordinates. Usually I use Griffiths Introduction to Electrodynamics, however I don’t always carry the book with me. As such I am leaving them here as a useful place for me to reference. While there are many places to reference they are usually not front in center.

I hope you find them useful too.


Cylindrical Polar Coordinates

\nabla U = \frac{\partial U}{\partial r}\boldsymbol{e}_r + \frac{1}{r}\frac{\partial U}{\partial \phi} \boldsymbol{e}_\phi + \frac{\partial U}{\partial z} \boldsymbol{e}_z 

\nabla \cdot \boldsymbol{F} = \frac{1}{r}\left[\frac{\partial(rF_r)}{\partial r} + \frac{\partial(rF_\phi)}{\partial \phi}  + \frac{\partial(rF_z)}{\partial z} \right]

\nabla^2  U = \frac{1}{r}\frac{\partial}{\partial r}r \frac{\partial U}{\partial r} +\frac{1}{r^2} \frac{\partial^2 U}{\partial \phi^2}  + \frac{\partial^2U}{\partial z^2}

Spherical Polar Coordinates

\nabla U = \frac{\partial U}{\partial r}\boldsymbol{e}_r + \frac{1}{r}\frac{\partial U}{\partial \theta} \boldsymbol{e}_\theta +  \frac{1}{r \mathrm{sin}\theta}\frac{\partial U}{\partial \phi} \boldsymbol{e}_\phi

\nabla \cdot \boldsymbol{F} = \frac{1}{r^2 \mathrm{sin}\theta}\left[\frac{\partial (r^2 \mathrm{sin} \theta F_r)}{\partial r} + \frac{\partial (\mathrm{sin}\theta F_\theta)}{\partial \theta} + \frac{\partial(rF_\phi)}{\partial \phi} \right]

\nabla^2  U = \frac{1}{r^2}\frac{\partial}{\partial r}r^2 \frac{\partial(U)}{\partial r} + \frac{1}{\mathrm{sin}\theta} \frac{\partial}{\partial \theta}\left(\mathrm{sin}\theta \frac{\partial}{\partial \theta}\right) +\frac{1}{\mathrm{sin}^2\theta} \frac{\partial^2 U}{\partial \phi^2}  

Laplacians, Gradients and Divergence