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:
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, :).