Skip to content



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)