Forces
\[ \begin{align}\begin{aligned}\newcommand{\pdev}[2]{\frac{\partial #1}{\partial #2}}\\\begin{split}E(\{ \mathbf{r}_{i} \}_{i})
&= \sum_{i} E_{i}(\{ \mathbf{r}_{j} \}_{j}) \\
\mathbf{F}_{i}
&= \pdev{ E(\{ \mathbf{r}_{j} \}_{j}) }{ \mathbf{r}_{i} }\end{split}\end{aligned}\end{align} \]
We can obtain the following gradients (here we call them pair forces) via autograd:
\[ \begin{align}\begin{aligned}\newcommand{\pdev}[2]{\frac{\partial #1}{\partial #2}}\\\mathbf{F}_{ji} = -\pdev{ E_{i}(\{ \mathbf{r}_{k} \}_{k}) }{ \mathbf{r}_{j} }.\end{aligned}\end{align} \]
With the translation symmetry of \(E_{i}(\{ \mathbf{r}_{k} \}_{k})\), we can write it as a function of displacements \(\mathbf{r}_{ij} = \mathbf{r}_{j} - \mathbf{r}_{i}\),
\[E_{i}(\{ \mathbf{r}_{k} \}_{k})
= \Psi( \{ \mathbf{r}_{ij} \}_{j \in \mathcal{N}_{i} } ).\]
Then, we can prove a relation similar to the Newton’s thrid law,
\[ \begin{align}\begin{aligned}\newcommand{\pdev}[2]{\frac{\partial #1}{\partial #2}}\\\begin{split}\mathbf{F}_{ii}
&= -\pdev{ E_{i}(\{ \mathbf{r}_{j} \}_{j}) }{ \mathbf{r}_{i} } \\
&= \sum_{j \in \mathcal{N}_{i}} \pdev{\Psi}{\mathbf{r}_{ij}}(\{ \mathbf{r}_{ik} \}_{k \in \mathcal{N}_{i}}) \\
&= \sum_{j \in \mathcal{N}_{i}} \pdev{ \Psi(\{ \mathbf{r}_{ik} \}_{k \in \mathcal{N}_{i}}) }{\mathbf{r}_{j}} \\
&= -\sum_{j \in \mathcal{N}_{i}} \mathbf{F}_{ji}. \\\end{split}\end{aligned}\end{align} \]
The forces can be written by the pair forces as
\[\begin{split}\mathbf{F}_{i}
&= \mathbf{F}_{ii} + \sum_{ j \in \mathcal{N}_{i} } \mathbf{F}_{ij} \\
&= \sum_{ j \in \mathcal{N}_{i} } \left( \mathbf{F}_{ij} - \mathbf{F}_{ji} \right).\end{split}\]
Virial stress
Ref. [BuvckoHAngyan05, TPM09]
Consider fractional coordinates \(\mathbf{r}_{i} = \mathbf{A} \mathbf{s}_{i}^{0}\) with basis vectors \(\mathbf{A} = \left( \mathbf{a}_{1}, \mathbf{a}_{2}, \mathbf{a}_{3} \right)\)
\[\begin{split}A_{\beta p}( \epsilon )
&:= \sum_{\alpha} A_{\alpha p} (\delta_{\alpha\beta} + \epsilon_{\alpha\beta}) \\
\tilde{E}( \epsilon )
&:= E( \{ \mathbf{r}_{i} \leftarrow \mathbf{A}(\epsilon)\mathbf{s}_{i}^{0} \} )\end{split}\]
We define a virial stress tensor as
\[ \begin{align}\begin{aligned}\newcommand{\pdev}[2]{\frac{\partial #1}{\partial #2}}\\\begin{split}\sigma_{\alpha\beta}
&:= \lim_{ \epsilon \to \mathbf{O} } -\frac{1}{V} \pdev{ \tilde{E}(\epsilon) }{ \epsilon_{\alpha\beta} } \\
&= \lim_{ \epsilon \to \mathbf{O} } - \frac{1}{V} \sum_{i\gamma} \pdev{E}{r_{i\gamma}}( \{ \mathbf{A}(\epsilon)\mathbf{s}_{j}^{0} \} ) \pdev{}{\epsilon_{\alpha\beta}} \sum_{p} A_{\gamma p}(\epsilon) s_{ip}^{0} \\
&= \lim_{ \epsilon \to \mathbf{O} } - \frac{1}{V} \sum_{i p} \pdev{E}{r_{i\beta}}( \{ \mathbf{A}(\epsilon)\mathbf{s}_{j}^{0} \} ) A_{\alpha p} s_{ip}^{0} \\
&= \frac{1}{V} \sum_{i} r_{i\alpha} F_{i\beta} .\end{split}\end{aligned}\end{align} \]
Rewrite the virial stress tensor with the pair forces:
\[\begin{split}\mathbf{\sigma}V
&= \sum_{i} \mathbf{r}_{i} \otimes \mathbf{F}_{i} \\
&= \sum_{i} \sum_{j \in \mathcal{N}_{i}}
\mathbf{r}_{i} \otimes \left( \mathbf{F}_{ij} - \mathbf{F}_{ji} \right) \\
&= \frac{1}{2} \sum_{i} \sum_{j \in \mathcal{N}_{i}}
\left(
\mathbf{r}_{i} \otimes \left( \mathbf{F}_{ij} - \mathbf{F}_{ji} \right)
+ \mathbf{r}_{j} \otimes \left( \mathbf{F}_{ji} - \mathbf{F}_{ij} \right) \\
\right) \\
&= \frac{1}{2} \sum_{i} \sum_{j \in \mathcal{N}_{i}}
\mathbf{r}_{ij} \otimes \left( \mathbf{F}_{ji} - \mathbf{F}_{ij} \right) \\\end{split}\]