Stress

For periodic condensed systems also the stress can be calculated analytically.

In general, the stress tensor $\sigma$ contains a kinetic contribution and a static contribution.

\[\begin{aligned} \sigma = \sigma^{\rm kin} + \sigma^{\rm static} \end{aligned} \]

The kinetic stress $\sigma^{\rm kin}$ is calculated from the atomic velocities $v$ and masses $m$ in the molecular dynamics simulation according to

\[\begin{aligned} \sigma^{\rm kin}_{ij} = \frac{1}{V}\sum_{k=1}^{N}m_k v_{k,i}v_{k,j} \end{aligned} \]

This has to be done by the MD code and is not part of RuNNer.

The static stress $\sigma^{\rm static}$ can be calculated analytically from the NN. We define

\[ \begin{aligned} R_{ij}^{\alpha} = \alpha_i - \alpha_j \quad , \end{aligned} \]

$\alpha$ being the $x$, $y$ or $z$ coordinate. We first calculate atomic contributions to the stress tensor separately for the radial and the angular symmetry functions. The radial stress contribution of atom $i$ to the static stress matrix element $\sigma^{\rm static}_{i,\alpha\beta}$ is

\[ \begin{aligned} \sigma^{\rm static,rad}_{i,\alpha\beta} =\sum_{j=1}^N R_{ij}^{\alpha}\cdot F_{j}^{\beta} \label{eq:radstress} \end{aligned} \]

Here $F_j^{\beta}$ is the force acting on atom $j$ in direction $\beta$. The angular stress contribution of atom $i$ is

\[ \begin{aligned} \sigma^{\rm static,ang}_{i,\alpha\beta} = \sum_{j=1}^N R_{ij}^{\alpha}\cdot F_{j}^{\beta} + \sum_{k=1}^N R_{ik}^{\alpha}\cdot F_{k}^{\beta} \quad . \label{angstress} \end{aligned} \]

The final matrix element of the static stress tensor is then obtained by adding all atomic contributions. $$ \begin{aligned} \sigma^{\rm static}{\alpha\beta} =\sum^N \left( \sigma^{\rm static,rad}{i,\alpha\beta} +\sigma^{\rm static,ang} \right) \end{aligned} $$

In Eqs. $\eqref{eq:radstress} and $\eqref{eq:angstress} we have to take into account the mapping of the Cartesian Coordinates on the symmetry functions.

Short range stress

For the radial stress of the short range component we obtain

\[ \begin{alignat}{2} \sigma^{\rm short,rad}_{\alpha\beta} &=&&\sum_i \sum_j R_{ij}^{\alpha}\cdot F_{j}^{\beta} \notag \\ &=&&-\sum_i\sum_j R_{ij}^{\alpha} \cdot \frac{\partial E}{\partial \beta_{j}} \notag \\ &=&&-\sum_i\sum_j R_{ij}^{\alpha} \sum_k \frac{\partial E_k}{\partial \beta_j} \notag \\ &=&&-\sum_i\sum_j R_{ij}^{\alpha} \sum_k \sum_{\mu}\frac{\partial E_k}{\partial G_{k}^{\mu}} \cdot \frac{\partial G_{k}^{\mu}}{\partial \beta_j} \notag \\ &=&&-\sum_i\sum_j\sum_k \sum_{\mu} R_{ij}^{\alpha} \frac{\partial E_k}{\partial G_{k}^{\mu}} \cdot \frac{\partial G_{k}^{\mu}}{\partial\beta_j} \notag \\ &=&&-\sum_k \sum_{\mu} \frac{\partial E_k}{\partial G_{k}^{\mu}} \cdot\sum_i\sum_j R_{ij}^{\alpha} \frac{\partial G_{k}^{\mu}}{\partial\beta_j} \end{alignat} \]

and for the angular short range stress

\[ \begin{alignat}{2} \sigma^{\rm short,ang}_{\alpha\beta} &=&&\sum_i \sum_j R_{ij}^{\alpha}\cdot F_{j}^{\beta} + \sum_i\sum_m R_{im}^{\alpha}\cdot F_m^{\beta} \notag \\ &=&&-\sum_k \sum_{\mu} \frac{\partial E_k}{\partial G_{k}^{\mu}}\cdot \left( \sum_i\sum_j R_{ij}^{\alpha} \frac{\partial G_{k}^{\mu}}{\partial\beta_j} + \sum_i\sum_m R_{im}^{\alpha} \frac{\partial G_{k}^{\mu}}{\partial\beta_m} \right) \end{alignat} \]

In RuNNer the short range stress is calculated using the quantities strs$(3,3,\mu,k)$=$\sum_i\sum_j R_{ij}^{\alpha}\frac{\partial G_k^{\mu}}{\partial \beta_j}$ and deshortdsfunc$(k,\mu)$$=\frac{\partial E_k}{\partial G_k^{\mu}}$.

Electrostatic stress

The electrostatic stress is not yet implemented.

van-der-Waals Stress

For a general introduction to dispersion interactions in RuNNer, see here.

In the original derivation of the Tkatchenko-Scheffler (TS) dispersion correction scheme, as presented by Bucko et al.,[^1] the expression for the gradient of the dispersion energy $E_{\mathrm{disp}}$ with respect to the lattice vector components $h^{\alpha,\beta}$ takes the form

\[\begin{aligned} \frac{\partial E_{\mathrm{disp}}}{\partial h^{\alpha,\beta}} = &-\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{\textbf{L}} \left[ -6\frac{C_{6ij}}{r_{ij,L}^7} f(r_{ij,L}) \frac{\partial r_{ij,L}}{h^{\alpha,\beta}} + \frac{C_{6ij}}{r_{ij,L}^6} \frac{\partial f(r_{ij,L})}{\partial r_{ij,L}} \frac{\partial r_{ij,L}}{h^{\alpha,\beta}} \right] \notag \\ &-\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{\textbf{L}} \left[ \frac{C_{6ij}}{r_{ij,L}^6} \frac{\partial f(r_{ij,L})}{\partial R_{ij}} \frac{\partial R_{ij}}{h^{\alpha,\beta}} + \frac{f(r_{ij,L})}{r_{ij,L}^6} \frac{\partial C_{6ij}}{\partial h^{\alpha,\beta}} \right]\,. \label{eq:vdwStressExpression} \end{aligned} \]

Please note the formal similarity of this equation to $\eqref{eq:vdwForceExpression}$.

The derivative of the screening function $f$ with respect to the interatomic distance $r_{ij,L}$ is given in $\eqref{eq:vdWScreeningDerivative}$. The derivative of the interatomic distance with respect to $h^{\alpha,\beta}$ is

\[ \begin{aligned} \frac{\partial r_{ij,L}}{\partial h^{\alpha,\beta}} = \frac{r_{ij,L}^{\beta}}{r_{ij,L}}\sum_{\gamma} \left(h^{\alpha, \beta}\right)^{-1}r_{ij,L}^{\gamma}\,. \label{eq:vdWDistDerivativeWRTLattice} \end{aligned} \]

Code implementation

In RuNNer prediction mode, the predict routine collects all stress contributions. Relevant variables are:

Name	Data Type	Description
`vdwstress(3, 3)`	REAL*8	The vdw contribution to the stress tensor components.

The stress is calculated alongside the energies. The main vdw calculation routine calc_tkatchenko_scheffler iterates over all atoms (loop_block_of_atoms) in one block of atoms and all their neighbours (loop_neighbors). For each pair, the contribution to the stress tensor is computed according to the equations given above and stored in vdwstress(3, 3).

[^1] T. Bucko, S. Lebègue, J. Hafner, J. G. Angtn, Phys. Rev. B 87 (6) 064110 (2013)