Forces

Overview

The force acting on coordinate $$\alpha$$ is given by

\begin{align} F_{\alpha} &= - \frac{\partial E}{\partial \alpha} \label{eq:ForceExpression}\\ &= - \frac{\partial E_{\rm short} + E_{\rm el}}{\partial \alpha} \notag \\ &= - \frac{\partial E_{\rm short}}{\partial \alpha} - \frac{\partial E_{\rm el}}{\partial \alpha} \notag \\ &= - \sum_i \frac{\partial E_{\rm i,short}}{\partial \alpha} - \sum_i \sum_{j>i} \frac{\partial \frac{Q_i Q_j}{r_{ij}}}{\partial \alpha} \notag \\ &= - \sum_i \sum_j\frac{\partial E_{i,\mathrm{short}}}{\partial G_j^i} \frac{\partial G_j^i}{\partial \alpha} - \sum_i \sum_{j>i} \frac{ \frac{\partial Q_i Q_j}{\partial \alpha} \cdot r_{ij} - Q_i Q_j \frac{\partial r_{ij}}{\partial \alpha} }{r_{ij}^2} \notag \\ &= - \sum_i \sum_j \frac{\partial E_{i,\mathrm{short}}}{\partial G_j^i} \frac{\partial G_j^i}{\partial \alpha} - \sum_i \sum_{j>i} \frac{1}{r_{ij}^2} \cdot \left( \left( \frac{\partial Q_i}{\partial \alpha}Q_j + Q_i\frac{\partial Q_j}{\partial \alpha} \right) \cdot r_{ij}- Q_iQ_j \cdot \frac{\partial r_{ij}}{\partial \alpha} \right) \notag \\ &= - \sum_i \sum_j \frac{\partial E_{i,\mathrm{short}}}{\partial G_j^i} \frac{\partial G_j^i}{\partial \alpha} - \sum_i \sum_{j>i} \frac{1}{r_{ij}^2} \cdot \left( \left( \sum_k \frac{\partial Q_i}{\partial G_k^i} \frac{\partial G_k^i}{\partial \alpha}Q_j + Q_i \sum_k \frac{\partial Q_j}{\partial G_k^j} \frac{\partial G_k^j}{\partial \alpha} \right) \cdot r_{ij}- Q_i Q_j \cdot \frac{\partial r_{ij}}{\partial \alpha} \right) \notag \end{align}

The calculation of the derivative $$\frac{\partial G_j^i}{\partial \alpha}$$ (dsfuncdxyz in the code) is explained in another chapter.

The calculation of $$\frac{\partial E_i}{\partial G_j^i}$$ (deshortdsfunc in the code) is illustrated on the example of a 2-3-4-5-1 NN. Its output is given by

\begin{align} E &= f^4 \left(b_1^4 + \sum_{l=1}^5 a_{l1}^{34} \cdot f^3 \left(b_l^3 + \sum_{k=1}^4 a_{kl}^{23} \cdot f^2 \left(b_k^2 + \sum_{j=1}^3 a_{jk}^{12} \cdot f^1 \left(b_j^1 + \sum_{i=1}^2 a_{ij}^{01} \cdot G_i\right)\right)\right)\right)\,. \label{eq:ForceExampleNN} \end{align}

We can calculate $$\frac{\partial E_i}{\partial G_j^i}$$ individually for each atom. We will therefore leave out the index $$i$$ and consider $$E$$ to be an atomic short range energy contribution here.

\begin{align} \frac{\partial E}{\partial G_i} &= \frac{\partial f^4(x_1^4)}{\partial G_i} \cdot \sum_{l=1}^5 a_{l1}^{34} \cdot \frac{\partial}{\partial G_i} f^3 \left( b_l^3 +\sum_{k=1}^4 a_{kl}^{23} \cdot f^2 \left( b_k^2 + \sum_{j=1}^3 a_{jk}^{12} \cdot f^1 \left( b_j^1 + \sum_{i=1}^2 a_{ij}^{01} \cdot G_i \right) \right) \right) \label{eq:ForceExampleNNDerivative} \\ &= \frac{\partial f^4(x_1^4)}{\partial G_i} \cdot \sum_{l=1}^5 a_{l1}^{34} \cdot \frac{\partial f^3(x_l^3)}{\partial G_i} \cdot\sum_{k=1}^4 a_{kl}^{23} \cdot\frac{\partial}{\partial G_i} f^2 \left( b_k^2 + \sum_{j=1}^3 a_{jk}^{12} \cdot f^1 \left( b_j^1 + \sum_{i=1}^2 a_{ij}^{01} \cdot G_i \right) \right) \notag \\ &= \frac{\partial f^4(x_1^4)}{\partial G_i} \cdot \sum_{l=1}^5 a_{l1}^{34} \cdot \frac{\partial f^3(x_l^3)}{\partial G_i} \cdot \sum_{k=1}^4 a_{kl}^{23} \cdot \frac{\partial f^2(x_k^2)}{\partial G_i} \cdot \sum_{j=1}^3 a_{jk}^{12} \cdot\frac{\partial }{\partial G_i} f^1 \left( b_j^1 + \sum_{i=1}^2 a_{ij}^{01}\cdot G_i \right) \notag \\ &= \frac{\partial f^4(x_1^4)}{\partial G_i} \cdot \sum_{l=1}^5 a_{l1}^{34} \cdot \frac{\partial f^3(x_l^3)}{\partial G_i} \cdot \sum_{k=1}^4 a_{kl}^{23} \cdot \frac{\partial f^2(x_k^2)}{\partial G_i} \cdot \sum_{j=1}^3 a_{jk}^{12} \cdot \frac{\partial f^1(x_j^1)}{\partial G_i} a_{ij}^{01} \notag \end{align}

For the short range contribution to the forces we use an auxiliary array tempderivative in the code. It has the dimensions tempderivative(num_layersshort, maxnodes_short, num_funcvalues). This array is used to calculate $$\frac{\partial E_j}{\partial G_i}$$ step by step. For tempderivative we will use the short form $$T_{L,j,i}$$, where $$L$$ is the index of the target layer, $$i$$ the index of the symmetry function (input node) and $$j$$ the index for the node in the target layer. We then have

\begin{align} T_{1,j,i} &= \frac{\partial f^1(x_j^1)}{\partial G_i}\cdot a_{ij}^{01} \label{eq:ForceExampleNNCodeDerivative} \\ T_{2,k,i} &= \frac{\partial f^2(x_k^2)}{\partial G_i} \cdot \sum_{j=1}^3 a_{jk}^{12} T_{1,j,i} \notag \\ T_{3,l,i} &= \frac{\partial f^3(x_l^3)}{\partial G_i} \cdot \sum_{k=1}^4 a_{kl}^{23} T_{2,k,i} \notag \\ T_{4,1,i} &= \frac{\partial f^4(x_1^4)}{\partial G_i} \cdot \sum_{l=1}^5 a_{l1}^{34} T_{3,l,i} \notag \\ \frac{\partial E_n}{\partial G_i} &= T_{4,1,i} \notag \end{align}

Electrostatic Forces

The calculation of the atomic forces requires the calculation of the derivative of the short range energy $$F_{\rm short}$$ and the derivative of the Ewald energy $$F_{\rm Ewald}$$ with respect to the Cartesian coordinates $$\alpha$$ of each atom $$i$$.

\begin{align} F_{\alpha_i,} &= F_{\alpha_i,{\rm short}} + F_{\alpha_i,{\rm Ewald}} \label{eq:ForceElectrostaticExpression} \end{align}

The derivative of the short range energy is

\begin{align} F_{\alpha_i,{\rm short}} &= -\frac{\partial E}{\partial \alpha_i} \notag \\ &= -\frac{\partial \sum_j E_j}{\partial \alpha_i} \notag \\ &= -\sum_j \frac{\partial E_j}{\partial \alpha_i} \notag \\ &= -\sum_j \sum_{\mu} \frac{\partial E_j}{\partial G_j^{\mu}} \cdot\frac{\partial G_j^{\mu}}{\partial \alpha_i} \end{align}

The derivative of the electrostatic energy depends on the presence of periodic boundary conditions, i.e., if the standard Coulomb law or the Ewald summation is used.

The Coulomb Case

The derivative of Coulomb's law is given by

\begin{align} F_{\alpha_{\mu},{\rm Coulomb}} &= -\frac{\partial E_{\rm Coulomb}}{\partial \alpha_{\mu}} \notag \\ &= -\frac{\partial}{\partial \alpha_{\mu}} \frac{1}{2}\sum_i\sum_{j\ne i}\frac{Q_i Q_j}{r_{ij}} \notag \\ &= -\sum_i \sum_{j\ne i}\frac{1}{2r_{ij}^2} \left[ \frac{\partial }{\partial \alpha_{\mu}}\left(Q_i Q_j\right)\cdot r_{ij} -Q_iQ_j\frac{\partial}{\partial \alpha_{\mu}} r_{ij} \right] \notag \\ &= -\sum_i \sum_{j\ne i}\frac{1}{2r_{ij}^2} \left[ \frac{\partial Q_i}{\partial \alpha_{\mu}}Q_j r_{ij} + Q_i\frac{\partial Q_j}{\partial \alpha_{\mu}}r_{ij} - Q_iQ_j\frac{\partial r_{ij}}{\partial \alpha_{\mu}} \right] \label{eq:ForceElectrostaticCoulomb} \end{align}

The first two terms can be combined by splitting the three sums and renaming the indices in the second sum:

\begin{align} F_{\alpha_{\mu},{\rm Coulomb}} &= -\sum_i \sum_{j\ne i} \left( \frac{Q_i}{r_{ij}}\frac{\partial Q_j}{\partial \alpha_{\mu}} \right) + \sum_i \sum_{j\ne i} \left( \frac{Q_iQ_j}{2r_{ij}^2}\frac{\partial r_{ij}}{\partial \alpha_{\mu}} \right) \notag \\ &= \sum_{i}\frac{Q_i}{r_{ij}}\sum_{j\ne i} \left[ \left(\frac{-\partial Q_j}{\partial \alpha_{\mu}}\right) + \frac{1}{2}\frac{Q_j}{r_{ij}} \frac{\partial r_{ij}}{\partial \alpha_{\mu}} \right] \end{align}

With

\begin{align} \frac{\partial Q_i}{\partial \alpha_{\mu}} = \sum_k\frac{\partial Q_i}{\partial G_{ik}} \frac{\partial G_{ik}}{\partial \alpha_{\mu}} \end{align}

we obtain

\begin{align} F_{\alpha_{\mu},{\rm Coulomb}} &= \sum_i \sum_{j\ne i} \frac{Q_i}{r_{ij}} \cdot\left[ \frac{1}{2}\frac{Q_j}{r_{ij}} \frac{\partial r_{ij}}{\partial \alpha_{\mu}} - \sum_k \frac{\partial Q_j}{\partial G_{jk}} \frac{\partial G_{jk}}{\partial \alpha_{\mu}} \right] \notag \\ &= \sum_j \sum_{i\ne j} \frac{Q_i}{r_{ij}} \cdot\left[ \frac{1}{2}\frac{Q_j}{r_{ij}} \frac{\partial r_{ij}}{\partial \alpha_{\mu}} - \sum_k \frac{\partial Q_j}{\partial G_{jk}} \frac{\partial G_{jk}}{\partial \alpha_{\mu}} \right] \end{align}

This can be parallelized over blocks of atoms $$j$$ as long as all charges $$Q_i$$ on all atoms are known.

The derivatives of the Coulomb energy with respect to the atomic positions are calculated in subroutine getcoulombforces.f90. The force $$-\frac{\partial E_{\rm Coulomb}}{\partial \alpha_i}$$ is nnewaldforce(max_num_atoms,3), the partial derivatives $$\frac{\partial Q_i}{\partial \alpha}$$ dchargedxyz(max_num_atoms, max_num_atoms,3), $$\frac{\partial Q_i}{\partial G_j}$$ is dchargedsfunc(max_num_atoms, num_funcvaluese), $$\frac{\partial G_j}{\partial \alpha}$$ is dsfuncdxyze(num_funcvaluese, max_num_atoms,max_num_atoms,3).

The Ewald Case

The derivative of the Ewald energy is given by

\begin{align} F_{\alpha_i,{\rm Ewald}} &= -\frac{\partial E_{\rm Ewald}}{\partial \alpha_i} \nonumber\\ &= -\frac{\partial E_{\rm real}}{\partial \alpha_i}-\frac{\partial E_{\rm recip}}{\partial \alpha_i}-\frac{\partial E_{\rm self}}{\partial \alpha_i} \end{align}

We have to take into account that the Ewald energy contributions $$E_{\rm real}$$ and $$E_{\rm recip}$$ directly depend on the coordinates $$\alpha_i$$ via the bond lengths $$R_{ij}$$ and that all three Ewald energy contributions depend on the $$\alpha_i$$ indirectly via the charges $$Q_i$$, since the charges are determined by the local atomic environments. Consequently, we have the derivatives

\begin{alignat}{2} -\frac{\partial E_{\rm real}}{\partial \alpha_{\mu}} &=&&-\frac{\partial}{\partial \alpha_{\mu}} \frac{1}{2}\sum_{i=1}^N \sum_{j\ne i}^N \left[ Q_i Q_j \frac{{\rm erfc}\left(\alpha R_{ij}\right)}{R_{ij}} \right] \notag \\ &=&&-\frac{1}{2}\sum_i^N \sum_{j\ne i}^N \left[ \frac{\partial}{\partial \alpha_{\mu}} \left( Q_i Q_j \frac{{\rm erfc}\left(\alpha R_{ij}\right)}{R_{ij}} \right) \right] \notag \\ &=&&-\frac{1}{2}\sum_i^N \sum_{j\ne i}^N \left[ \frac{\partial Q_i Q_j}{\partial \alpha_{\mu}} \frac{{\rm erfc}\left(\alpha R_{ij}\right)}{R_{ij}} + Q_i Q_j \frac{\partial}{\partial \alpha_{\mu}} \frac{{\rm erfc} \left(\alpha R_{ij}\right)}{R_{ij}} \right] \notag \\ &=&&-\frac{1}{2}\sum_i^N \sum_{j\ne i}^N \Big[ \left( Q_j\frac{\partial Q_i}{\partial \alpha_{\mu} } + Q_i\frac{\partial Q_j}{\partial \alpha_{\mu} } \right) \cdot \frac{{\rm erfc} \left(\alpha R_{ij}\right)}{R_{ij}} + Q_i Q_j \bigg( \frac{1}{R_{ij}} \frac{\partial}{\partial\alpha_{\mu}}{\rm erfc} \left(\alpha R_{ij}\right) \notag \\ & &&-{\rm erfc} \left(\alpha R_{ij}\right) \frac{1}{R_{ij}^2}\frac{\partial R_{ij}}{\partial\alpha_{\mu}} \bigg) \Big] \notag \\ &=&&-\frac{1}{2}\sum_i^N \sum_{j\ne i}^N \Big[ \left( Q_j\sum_k\frac{\partial Q_i}{\partial G_i^k} \frac{\partial G_i^k}{\partial \alpha_{\mu}} + Q_i\sum_k\frac{\partial Q_j}{\partial G_j^k} \frac{\partial G_j^k}{\partial\alpha_{\mu}} \right) \notag \\ & &&\cdot \frac{{\rm erfc} \left(\alpha R_{ij}\right)}{R_{ij}} + Q_i Q_j \bigg( \frac{1}{R_{ij}} \frac{\partial}{\partial\alpha_{\mu}}{\rm erfc} \left(\alpha R_{ij}\right) -{\rm erfc}\left(\alpha R_{ij}\right) \frac{1}{R_{ij}^2}\frac{\partial R_{ij}}{\partial\alpha_{\mu}} \bigg) \Big] \notag \\ &=&&-\frac{1}{2} \left[ \sum_i^N \sum_{j\ne i}^N Q_j \sum_k \frac{\partial Q_i}{\partial G_i^k} \frac{\partial G_i^k}{\partial\alpha_{\mu}} + \sum_i^N \sum_{j\ne i}^N Q_i \sum_k \frac{\partial Q_j}{\partial G_j^k} \frac{\partial G_j^k}{\partial \alpha_{\mu}} \right] \cdot \frac{{\rm erfc}\left(\alpha R_{ij}\right)}{R_{ij}} \notag \\ & &&-\frac{1}{2}\sum_i^N \sum_{j\ne i}^N Q_iQ_j \cdot \left[ \frac{1}{R_{ij}} \frac{\partial}{\partial\alpha_{\mu}}{\rm erfc} \left(\alpha R_{ij}\right) -{\rm erfc}\left(\alpha R_{ij}\right) \cdot \frac{1}{R_{ij}^2} \frac{\partial R_{ij}}{\partial \alpha_{\mu}} \right] \end{alignat}

The first two terms can be combined by renaming the indices of summation. We then get

\begin{alignat}{2} -\frac{\partial E_{\rm real}}{\partial\alpha_{\mu}} &=&&-\sum_i^N Q_i \sum_{j\ne i}^N \sum_k \frac{\partial Q_j}{\partial G_j^k} \frac{\partial G_j^k}{\partial \alpha_{\mu}} \cdot \frac{{\rm erfc}\left(\alpha R_{ij}\right)}{R_{ij}} \notag \\ & &&-\sum_i^N Q_i \sum_{j\ne i}^N \frac{Q_j}{2}\cdot\left[ \frac{1}{R_{ij}}\frac{\partial}{\partial\alpha_{\mu}} {\rm erfc}\left(\alpha R_{ij}\right) - {\rm erfc}\left(\alpha R_{ij}\right) \cdot \frac{1}{R_{ij}^2}\frac{\partial R_{ij}}{\partial \alpha_{\mu}} \right] \notag \\ &=&&-\sum_i^N Q_i \cdot \sum_{j\ne i}^N \Big[ \frac{{\rm erfc}\left(\alpha R_{ij}\right)}{R_{ij}} \sum_k \left[ \frac{\partial Q_j}{\partial G_j^k} \frac{\partial G_j^k}{\partial \alpha_{\mu}} \right] \notag \\ & &&+ \frac{Q_j}{2R_{ij}} \left( \frac{\partial}{\partial\alpha_{\mu}} {\rm erfc}\left(\alpha R_{ij}\right) - {\rm erfc}\left(\alpha R_{ij}\right) \cdot \frac{1}{R_{ij}} \frac{\partial R_{ij}}{\partial \alpha_{\mu}} \right) \Big] \end{alignat}
\begin{alignat}{2} -\frac{\partial E_{\rm recip}}{\partial \alpha_i} &=&&-\frac{2\pi}{V}\sum_{k\ne 0}^\infty \frac{{\rm exp}\left(-\frac{k^2}{4\alpha^2}\right)}{k^2} \cdot\frac{\partial}{\partial \alpha_i} \Bigg( \left( \sum_{j=1}^N Q_j\:{\rm cos}\left(\mathbf{k} \mathbf{R}_j\right) \right)^2 +\left( \sum_{j=1}^N Q_j\:{\rm sin} \left(\mathbf{k}\mathbf{R}_j\right) \right)^2 \Bigg) \notag \\ &=&&-\frac{2\pi}{V}\sum_{k\ne 0}^\infty \frac{{\rm exp}\left(-\frac{k^2}{4\alpha^2}\right)}{k^2}\cdot \Bigg( \frac{\partial}{\partial \alpha_i} \left( \sum_{j=1}^N Q_j\:{\rm cos}\left(\mathbf{k} \mathbf{R}_j\right) \right)^2 +\frac{\partial}{\partial \alpha_i} \left( \sum_{j=1}^N Q_j\:{\rm sin}\left(\mathbf{k} \mathbf{R}_j\right) \right)^2 \Bigg) \notag \\ &=&&-\frac{2\pi}{V}\sum_{k\ne 0}^\infty \frac{{\rm exp}\left(-\frac{k^2}{4\alpha^2}\right)}{k^2}\cdot \Bigg( 2\left( \sum_{j=1}^N Q_j\:{\rm cos}\left(\mathbf{k}\mathbf{R}_j\right) \right) \left( \sum_{j=1}^N \frac{\partial}{\partial \alpha_i} \left( Q_j\:{\rm cos}\left(\mathbf{k} \mathbf{R}_j\right) \right) \right) \notag \\ & &&+2\left( \sum_{j=1}^N Q_j\:{\rm sin}\left(\mathbf{k}\mathbf{R}_j\right) \right) \left( \sum_{j=1}^N\frac{\partial}{\partial \alpha_i} \left( Q_j\:{\rm sin}\left(\mathbf{k} \mathbf{R}_j\right) \right) \right) \Bigg) \notag \\ &=&&-\frac{2\pi}{V}\sum_{k\ne 0}^\infty \frac{{\rm exp}\left(-\frac{k^2}{4\alpha^2}\right)}{k^2}\cdot \Bigg( 2\left( \sum_{j=1}^N Q_j\:{\rm cos}\left(\mathbf{k}\mathbf{R}_j\right) \right) \notag \\ & &&\left( \sum_{j=1}^N \left[ {\rm cos} \left(\mathbf{k} \mathbf{R}_j\right) \frac{\partial Q_j}{\partial\alpha_i}-Q_j {\rm sin} \left(\mathbf{k} \mathbf{R}_j\right) \frac{\partial}{\partial \alpha_i} \left(\mathbf{k} \mathbf{R}_j\right) \right] \right) \notag \\ & &&+2\left( \sum_{j=1}^N Q_j\:{\rm sin}\left(\mathbf{k}\mathbf{R}_j\right) \right) \left( \sum_{j=1}^N \left[ {\rm sin}\left(\mathbf{k} \mathbf{R}_j\right) \frac{\partial Q_j}{\partial\alpha_i} + Q_j \:{\rm cos}\left(\mathbf{k}\mathbf{R}_j\right) \frac{\partial}{\partial\alpha_i} \left(\mathbf{k} \mathbf{R}_j\right) \right] \right) \Bigg) \end{alignat}

In RuNNer each process has access to the charges of all atoms. Therefore, for a given $$\mathbf{k}$$ vector the sums over all atoms can be precalculated as

\begin{align} C_1 = 2\left( \sum_{j=1}^N Q_j\:{\rm cos}\left(\mathbf{k} \mathbf{R}_j\right) \right) \end{align}

and

\begin{align} C_2 = 2\left( \sum_{j=1}^N Q_j\:{\rm sin}\left(\mathbf{k} \mathbf{R}_j\right) \right) \end{align}

Further, we introduce the abbreviation

\begin{align} C_3 = -\frac{2\pi}{V}\frac{{\rm exp} \left(-\frac{k^2}{4\alpha^2}\right)}{k^2} \end{align}

We then have

\begin{alignat}{2} -\frac{\partial E_{\rm recip}}{\partial \alpha_i} &=&&\sum_{k\ne 0}^\infty C_3 \cdot \Bigg( C_1 \cdot \left(\sum_{j=1}^N \left[ {\rm cos}\left(\mathbf{k} \mathbf{R}_j\right) \frac{\partial Q_j}{\partial \alpha_i}-Q_j{\rm sin} \left(\mathbf{k} \mathbf{R}_j\right) \frac{\partial}{\partial \alpha_i} \left(\mathbf{k} \mathbf{R}_j\right) \right] \right) \notag \\ & &&+ C_2\cdot \left(\sum_{j=1}^N \left[ {\rm sin} \left(\mathbf{k} \mathbf{R}_j\right) \frac{\partial Q_j}{\partial \alpha_i} +Q_j \:{\rm cos}\left(\mathbf{k}\mathbf{R}_j\right) \frac{\partial}{\partial\alpha_i} \left(\mathbf{k} \mathbf{R}_j\right) \right] \right) \Bigg) \end{alignat}

For the derivative of the self-interaction part we have:

\begin{alignat}{2} -\frac{\partial E_{\rm self}}{\partial \alpha_i} &=&&-\frac{\alpha}{\sqrt{\pi}}\sum_{i=1}^N \frac{\partial Q_i^2}{\partial \alpha_i} \notag \\ &=&&-\frac{\alpha}{\sqrt{\pi}}\sum_{i=1}^N 2Q_i \frac{\partial Q_i}{\partial \alpha_i} \notag \\ &=&&-\frac{\alpha}{\sqrt{\pi}}\sum_{i=1}^N 2Q_i \sum_j \frac{\partial Q_i}{\partial {G_i^j}'} \frac{\partial{G_i^j}'}{\partial \alpha_i} \end{alignat}

For the derivatives in RuNNer the following variables are used in subroutine getewaldforces.f90:

• $$\frac{\partial Q_i}{\partial G_i^k}$$ = dchargedsfunc$$(i,k)$$
• $$\frac{\partial Q_i}{\partial \alpha_j}$$ = dchargedxyz$$(i,j,3)$$
• $$\frac{\partial G_j^k}{\partial \alpha_i}$$ = dsfuncdxyze$$(k,j,i,3)$$

The final electrostatic force $$F_{\alpha_i,Ewald}$$ is nnewaldforce$$(i,3)$$.

There are two different ways to derive the derivative.

Example for the self-energy:

\begin{alignat}{2} -\frac{\partial E_{\rm self}}{\partial \alpha_i} &=&&-\frac{\alpha}{\sqrt{\pi}} \sum_i \frac{\partial}{\partial \alpha_i}Q_i^2 \notag \\ &=&&-\frac{\alpha}{\sqrt{\pi}}\sum_i 2Q_i \frac{\partial Q_i}{\partial \alpha_i} \notag \end{alignat}
\begin{alignat}{2} -\frac{\partial E_{\rm self}}{\partial \alpha_i} &=&&-\sum_i \frac{\partial E_{\rm self}}{\partial Q_i} \frac{\partial Q_i}{\partial \alpha_i} \notag \\ &=&&\sum_i -\frac{\alpha}{\sqrt{\pi}}2Q_i \frac{\partial Q_i}{\partial \alpha} \notag \\ &=&& -\frac{\alpha}{\sqrt{\pi}} \sum_i 2Q_i \frac{\partial Q_i}{\partial \alpha} \notag \end{alignat}

because

\begin{align} \frac{\partial E_{\rm self}}{\partial Q_i} = \frac{\alpha}{\sqrt{\pi}}\sum_{j=1}^N \frac{\partial}{\partial Q_j}Q_i^2 = \frac{\alpha}{\sqrt{\pi}}2 Q_i \end{align}

Caution: The self-energy is a particularly easy case, because there is just $$Q_i$$ depending on $$\alpha_i$$, and no $$r_{ij}$$.

van-der-Waals Dispersion Forces

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

\begin{align} \frac{\partial E_{\mathrm{disp}}}{\partial r_{k}^{\alpha}} = &-\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}}{r_{k}^{\alpha}} + \frac{C_{6ij}}{r_{ij,L}^6}\frac{\partial f(r_{ij,L})}{\partial r_{ij,L}}\frac{\partial r_{ij,L}}{r_{k}^{\alpha}} \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}}{r_{k}^{\alpha}} + \frac{f(r_{ij,L})}{r_{ij,L}^6}\frac{\partial C_{6ij}}{\partial r_{k}^{\alpha}} \right]\,. \label{eq:vdwForceExpression} \end{align}

As before, $$N$$ denotes the number of atoms in the structure, $$\textbf{L}$$ is the unit cell translation vector, and $$r_{ij,L}$$ stands for the interatomic pairwise distance. Note that $$i\neq j$$ for $$\textbf{L}=0$$.

The derivative of the screening function $$f$$ with respect to the interatomic distance $$r_{ij,L}$$ is

\begin{align} \frac{\partial f(r_{ij,L})}{\partial r_{ij,L}} = \frac{f^2(r_{ij,L})}{s_6}\frac{d}{s_{\mathrm{R}}R_{ij}}\mathrm{e}^{-d\left(\frac{r_{ij}}{s_{\mathrm{R}} R_{ij}} - 1\right)}\,. \label{eq:vdWScreeningDerivative} \end{align}

Furthermore, the derivative of the interatomic distance with respect to the Cartesian components $$r_{k}^{\alpha}$$ of the center atom coordinates is given by

\begin{align} \frac{\partial r_{ij,L}}{r_{k}^{\alpha}} = \frac{| r_{ij,L}^{\alpha}|}{r_{ij,L}}(\delta_{ik}-\delta_{kj})\,. \label{eq:vdWDistDerivativeWRTCart} \end{align}

In this equation, $$\delta$$ denotes Kronecker's delta function.

Two realisations help simplify $$\eqref{eq:vdwForceExpression}$$: First, as the current implementation of the TS scheme in RuNNer does not support charge-density dependent $$C_{6ij}$$ coefficients or vdW radii, the third and fourth term of $$\eqref{eq:vdwForceExpression}$$ both equal zero.

Second, with regard to the computational implementation it is easier to change ones point of view from calculating the derivative of the total dispersion energy to calculating the derivative of the atomic dispersion energy. As the condition in $$\eqref{eq:vdWDistDerivativeWRTCart}$$ can only be non-zero in two cases, either $$k=i$$ or $$k=j$$ (but not both at the same time), this effectively reduces equation eq.$$\eqref{eq:vdwForceExpression}$$ to

\begin{align} \frac{\partial E_{\mathrm{k,disp}}}{\partial r_{k}^{\alpha}} = &\sum_{l}^{N_{\mathrm{nei}}} \left[ -6\frac{C_{6kl}}{r_{kl}^7} f(r_{kl})\frac{| r_{kl}^{\alpha}|}{r_{kl}} + \frac{C_{6kl}}{r_{kl}^6}\frac{\partial f(r_{kl})}{\partial r_{kl}}\frac{| r_{kl}^{\alpha}|}{r_{kl}} \right]\,,\label{eq:vdwForceExpressionShortened} \end{align}

where $$N_{\mathrm{nei}}$$ now denotes the number of neighbours of atom $$k$$.

Code implementation

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

Name Data Type Description
vdwforces(3, :) REAL*8 The vdw dispersion force components of each atom.

The forces are 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), summing up the pairwise dispersion force components according to $$\eqref{eq:vdwForceExpression}$$ and storing them in vdwforces(3, :).