Force and stress calculation

Theory

Consider a material containing \(N\) atoms has a total potential energy \(\cal E\). The atomistic environment of each atom can be described by \(M\) fingerprints \(G_{im}\) (\(i=1,2,\dots,N\) and \(m=1,2,\dots,M\)), which are used to determine the potential energy \(E_{i}\) of \(i\)-th atom. The atomic forces acting on each atom can be calculated as:

(1)\[f_{j\alpha}=-\frac{\partial {\cal E}}{\partial r_{j\alpha}}= - \sum_{i=1}^{N} \frac{\partial E_i (G_{i1}, G_{i2},\dots,G_{iM})}{\partial r_{j\alpha}} = - \sum_{i=1}^{N} \sum_{m=1}^{M} {\frac{\partial E_i}{\partial G_{im}}} {\frac{\partial G_{im}}{\partial r_{j\alpha}}}\]

where \(r_{j\alpha}\) (\(\alpha=1,2,3\)) are the Cartesian coordinates of \(j\)-th atom (\(j=1,2,\dots,N\)).

The first Piola-Kirchhoff (PK) stress tensor can be calculated as the work conjugate of deformation gradient tensor

(2)\[\begin{equation} P_{\alpha\beta} = \frac{1}{V_0} \frac{\partial {\cal E}}{\partial F_{\alpha\beta}}, \end{equation}\]

where \(V_0\) is the volume of the reference configuration and \(F_{\alpha\beta}\) is the deformation gradient tensor. Similar to the atomic force calculation, the potential energy can be written as the sum of atomic potential energies, so the stress can be written as

(3)\[\begin{equation} P_{\alpha\beta} = \frac{1}{V_0} \sum_{i=1}^{N} \sum_{m=1}^{M} \frac{\partial E_i}{\partial G_{im}} \frac{\partial G_{im} } {\partial F_{\alpha \beta}}. \end{equation}\]

The fingerprint \(G_{im}\) is determined by the coordinates of the atoms inside the cutoff of atom \(i\). Therefore,

(4)\[\begin{equation} \frac{\partial G_{im} } {\partial F_{\alpha \beta}} = \sum_{j \in \text{NB}_{i}} \sum_{\gamma=1}^{3} \frac{\partial G_{im}}{\partial r_{j\gamma}} \frac{\partial r_{j\gamma}}{\partial F_{\alpha \beta}}, \end{equation}\]

where atom \(j\) is inside the neighbor list of atom \(i\). By definition, the deformation gradient maps the atom positionfrom the reference configuration to the current configuration.

(5)\[\begin{equation} r_{j\gamma} = \sum_{\beta=1}^{3}F_{\gamma\beta}R_{j\beta} , \end{equation}\]

where \(R_{j\beta}\) is the coordinates of atom \(j\) in the reference configuration. Then, the stress in Eq. (3) can be written as

(6)\[\begin{equation} P_{\alpha\beta} = \frac{1}{V_0} \sum_{i=1}^{N} \sum_{m=1}^{M} \sum_{j \in \text{NB}_i} \frac{\partial E_i}{\partial G_{im}} \frac{\partial G_{im} } {\partial r_{j \alpha}} R_{j\beta}. \end{equation}\]

Cauchy stress can be further calculated by

(7)\[\begin{equation} \sigma_{\alpha\gamma} = \det ({\bf F})^{-1} \sum_{\beta=1}^{3}P_{\alpha\beta} F_{\gamma \beta}. \end{equation}\]

Substitute Eq. (6) into Eq. (7), and then apply \(V = V_0 \det ({\bf F})\), which is the volume in current configuration, as well as \(r_{j\gamma} = \sum_{\beta =1}^{3}F_{\gamma\beta}R_{j\beta}\), we get

(8)\[\begin{equation} \sigma_{\alpha\beta} = \frac{1}{V} \sum_{i=1}^{N} \sum_{m=1}^{M} \sum_{j \in \text{NB}_i} \frac{\partial E_i}{\partial G_{im}} \frac{\partial G_{im} } {\partial r_{j \alpha}} r_{j\beta} \end{equation}\]

after replacing \(\gamma\) with \(\beta\).

Implementation

In Tensorflow, it is more convenient and efficient to implement the above summation in terms of matrix multiplication. Therefore, we define matrix

\[[dEdG]_i = \begin{pmatrix} \dfrac{\partial E_i}{\partial G_{i1}} & \dfrac{\partial E_i}{\partial G_{i2}} & \dots & \dfrac{\partial E_i}{\partial G_{iM}} \end{pmatrix}\]

and

\[\begin{split}[dGdr]_{ij} = \begin{bmatrix} \dfrac{\partial G_{i1}}{\partial r_{j1}} & \dfrac{\partial G_{i1}}{\partial r_{j2}} & \dfrac{\partial G_{i1}}{\partial r_{j3}} \\ \rule{0pt}{22pt} \dfrac{\partial G_{i2}}{\partial r_{j1}} & \dfrac{\partial G_{i2}}{\partial r_{j2}} & \dfrac{\partial G_{i2}}{\partial r_{j3}} \\ \rule{0pt}{15pt} \vdots & \vdots & \vdots \\\rule{0pt}{15pt} \dfrac{\partial G_{iM}}{\partial r_{j1}} & \dfrac{\partial G_{iM}}{\partial r_{j2}} & \dfrac{\partial G_{iM}}{\partial r_{j3}} \\ \end{bmatrix}\end{split}\]

Next, we use an example to demonstrate the matrix multiplicaiton process. Assume there are 8 pairs of derivative (\(dGdr\)) data in one structure that contains 6 atoms. The center atoms IDs for thes 8 pairs are

\[\begin{bmatrix} 1&2&2&4&5&0&1&2 \end{bmatrix}\]

The neighbor atoms IDs are

\[\begin{equation} \begin{bmatrix} 2&3&4&2&1&4&0&4 \end{bmatrix} \end{equation}\]

The atomic forces can be computed as:

\[\begin{split}\begin{equation} \begin{bmatrix} \left[dEdG\right]_1 \\ \rule{0pt}{10pt} \left[dEdG\right]_2 \\ \rule{0pt}{10pt} \left[dEdG\right]_2 \\ \rule{0pt}{10pt} \left[dEdG\right]_4 \\ \rule{0pt}{10pt} \left[dEdG\right]_5 \\ \rule{0pt}{10pt} \left[dEdG\right]_0 \\ \rule{0pt}{10pt} \left[dEdG\right]_1 \\ \rule{0pt}{10pt} \left[dEdG\right]_2 \\ \end{bmatrix} \begin{bmatrix} \left[dGdr\right]_{12} \\ \rule{0pt}{10pt} \left[dGdr\right]_{23} \\ \rule{0pt}{10pt} \left[dGdr\right]_{24} \\ \rule{0pt}{10pt} \left[dGdr\right]_{42} \\ \rule{0pt}{10pt} \left[dGdr\right]_{51} \\ \rule{0pt}{10pt} \left[dGdr\right]_{04} \\ \rule{0pt}{10pt} \left[dGdr\right]_{10} \\ \rule{0pt}{10pt} \left[dGdr\right]_{24} \\ \end{bmatrix} = \begin{bmatrix} (f_2^x)_1 & (f_2^y)_1 & (f_2^z)_1 \\ \rule{0pt}{10pt} (f_3^x)_2 & (f_3^y)_2 & (f_3^z)_2 \\ \rule{0pt}{10pt} (f_4^x)_2 & (f_4^y)_2 & (f_4^z)_2 \\ \rule{0pt}{10pt} (f_2^x)_4 & (f_2^y)_4 & (f_2^z)_4 \\ \rule{0pt}{10pt} (f_1^x)_5 & (f_1^y)_5 & (f_1^z)_5 \\ \rule{0pt}{10pt} (f_4^x)_0 & (f_4^y)_0 & (f_4^z)_0 \\ \rule{0pt}{10pt} (f_0^x)_1 & (f_0^y)_1 & (f_0^z)_1 \\ \rule{0pt}{10pt} (f_4^x)_2 & (f_4^y)_2 & (f_4^z)_2 \end{bmatrix} \end{equation}\end{split}\]

where \((f_j^x)_i\) is the force on \(j\)-th atom contributed by atom \(i\). Then we need to compute the total force acting on atom \(j\), \(F_j\), by summing up all the contributions for other atoms. This is done by using one-hot matrix generated by the neighbor atom IDs:

\[\begin{split}\begin{equation} -\begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ \rule{0pt}{10pt} 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ \rule{0pt}{10pt} 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \rule{0pt}{10pt} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \rule{0pt}{10pt} 0 & 0 & 1 & 0 & 0 & 1 & 0 & 1 \\ \rule{0pt}{10pt} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} (f_2^x)_1 & (f_2^y)_1 & (f_2^z)_1 \\ \rule{0pt}{10pt} (f_3^x)_2 & (f_3^y)_2 & (f_3^z)_2 \\ \rule{0pt}{10pt} (f_4^x)_2 & (f_4^y)_2 & (f_4^z)_2 \\ \rule{0pt}{10pt} (f_2^x)_4 & (f_2^y)_4 & (f_2^z)_4 \\ \rule{0pt}{10pt} (f_1^x)_5 & (f_1^y)_5 & (f_1^z)_5 \\ \rule{0pt}{10pt} (f_4^x)_0 & (f_4^y)_0 & (f_4^z)_0 \\ \rule{0pt}{10pt} (f_0^x)_1 & (f_0^y)_1 & (f_0^z)_1 \\ \rule{0pt}{10pt} (f_4^x)_2 & (f_4^y)_2 & (f_4^z)_2 \end{bmatrix} = \begin{bmatrix} F_0^x & F_0^y & F_0^z \\ \rule{0pt}{10pt} F_1^x & F_1^y & F_1^z \\ \rule{0pt}{10pt} F_2^x & F_2^y & F_2^z \\ \rule{0pt}{10pt} F_3^x & F_3^y & F_3^z \\ \rule{0pt}{10pt} F_4^x & F_4^y & F_4^z \\ \rule{0pt}{10pt} F_5^x & F_5^y & F_5^z \end{bmatrix} \end{equation}\end{split}\]

Define the coordinates of neighboring atom

\[\begin{split}\begin{equation} [r]_j = \begin{bmatrix} r_{jx}\\ r_{jy}\\ r_{jz} \end{bmatrix} \end{equation}\end{split}\]

The stress can be calculated as

\[\begin{split}\begin{equation} \begin{bmatrix} [r]_2 \\ [r]_3 \\ [r]_4 \\ [r]_2 \\ [r]_1 \\ [r]_4 \\ [r]_0 \\ [r]_4 \end{bmatrix} \begin{bmatrix} [(f_2^x)_1 & (f_2^y)_1 & (f_2^z)_1] \\ \rule{0pt}{10pt} [(f_3^x)_2 & (f_3^y)_2 & (f_3^z)_2] \\ \rule{0pt}{10pt} [(f_4^x)_2 & (f_4^y)_2 & (f_4^z)_2]\\ \rule{0pt}{10pt} [(f_2^x)_4 & (f_2^y)_4 & (f_2^z)_4] \\ \rule{0pt}{10pt} [(f_1^x)_5 & (f_1^y)_5 & (f_1^z)_5] \\ \rule{0pt}{10pt} [(f_4^x)_0 & (f_4^y)_0 & (f_4^z)_0] \\ \rule{0pt}{10pt} [(f_0^x)_1 & (f_0^y)_1 & (f_0^z)_1] \\ \rule{0pt}{10pt} [(f_4^x)_2 & (f_4^y)_2 & (f_4^z)_2] \end{bmatrix} = \begin{bmatrix} [\sigma_{\alpha\beta}]_1 \\ [\sigma_{\alpha\beta}]_2 \\ [\sigma_{\alpha\beta}]_3 \\ [\sigma_{\alpha\beta}]_4 \\ [\sigma_{\alpha\beta}]_5 \\ [\sigma_{\alpha\beta}]_6 \\ [\sigma_{\alpha\beta}]_7 \\ [\sigma_{\alpha\beta}]_8 \\ \end{bmatrix} \end{equation}\end{split}\]

The final stress is

\[\begin{equation} \sigma_{\alpha\beta} = \frac{1}{V} \sum_{k=1}^{8} [\sigma_{\alpha\beta}]_k \end{equation}\]