Energy

Overview

In RuNNer the total energy of the system consists of a short range contribution and of an electrostatic energy arising from the interaction of charges centered at the atoms,

\[ \begin{align} E_{\rm tot} = E_{\rm short} + E_{\rm elec} \label{eq:TotalEnergy} \end{align} \]

Short Range Energy Contribution

The short range contribution to the total energy is a sum over environment-dependent atomic energies. The atomic energies \(E_{\mu}\) are outputs of individual atomic NNs,

\[ \begin{align} E_{\rm short} = \sum_{\mu} E_{\mu} \label{eq:AtomicEnergy} \end{align} \]

The atomic energy expression is given here for an example 2-3-4-5-1 NN.

\[ \begin{align} E_{\mu} &= 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:EnergyExampleNN} \\ &= 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(x_j^1\right) \right) \right) \right) \notag \\ &= 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 y_j^1\right) \right) \right) \notag \\ &= 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(x_k^2\right) \right) \right) \notag \\ &= 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 y_k^2\right) \right) \notag \\ &= f^4\left(b_1^4 +\sum_{l=1}^5 a_{l1}^{34}\cdot f^3 \left(x_l^3\right) \right) \notag \\ &= f^4\left(b_1^4 +\sum_{l=1}^5 a_{l1}^{34} \cdot y_l^3\right) \notag \\ &= f^4\left(x_1^4\right) \notag \\ &= y_1^4 \notag \end{align} \]

Here we have introduced the following abbreviations:

• \(x_i^j\) is the value of node \(i\) in layer \(j\) before applying the activation function.

• \(y_i^j=f^k(x_i^j)\) is the value of node \(i\) in layer \(j\).

• \(G_i=y_i^0\) is the \(i\)th input node, i.e., symmetry function.

• \(a_{ij}^{kl}\) is the weight connecting node \(j\) in layer \(l\) with node \(i\) in layer \(k\).

• \(b_i^j\) is the bias weight acting on node \(i\) in layer \(j\).

• \(f^i\) is the activation function acting on the nodes in layer \(i\).

The output of the atomic NN for the short-range interaction is calculated in the subroutine calconeatom.f90.

Electrostatic Energy Contribution

The electrostatic energy contribution is calculated using Coulomb's law for non-periodic systems:

\[ \begin{align} E_{\rm elec}= \sum_i\sum_{j>i}\frac{Q_i Q_j}{r_{ij}} \label{eq:ElectrostaticCoulombsLaw} \end{align} \]

The atomic charges are fitted in the same way as the short range atomic energy contributions as a function of the local atomic environment.

In case of a periodic system a standard Ewald summation is used:

\[ \begin{align} E_{\rm elec}=E_{\rm real}+E_{\rm recip}+E_{\rm self} \label{eq:ElectrostaticEwaldSummation} \end{align} \]

The Ewald summation contains the real space term \(E_{\rm real}\), the reciprocal space term \(E_{\rm recip}\) and the self-energy correction term \(E_{\rm self}\).

\[ \begin{align} E_{\rm real} &=&&\frac{1}{2}\sum_{i=1}^N \sum_{j\ne i}^N Q_i Q_j \frac{{\rm erfc}\left(\alpha R_{ij}\right)}{R_{ij}} \label{eq:ElectrostaticEwaldReal} \\ E_{\rm recip} &=&&\frac{2\pi}{V}\sum_{k\ne0}^\infty \frac{4\pi^2}{k^2}{\rm exp}\left(-\frac{k^2}{4\alpha^2}\right) \notag \\ & &&\cdot \left[\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 \right] \label{eq:ElectrostaticEwaldRecip} \\ E_{\rm self} &=&&\frac{\alpha}{\sqrt{\pi}}\sum_{i=1}^NQ_i^2 \label{eq:ElectrostaticEwalSelf} \end{align} \]

In Eq. \(\eqref{eq:ElectrostaticEwaldReal}\) it is in principle assumed that the summation over the charges within the simulation cell is sufficient due to an appropriate choice of \(\alpha\). In RuNNer for the real space contribution a neighbor list is constructed based on the ewaldcutoff specified in input.nn. Then all real space interactions within this cutoff are included. Thus the definition of the unit cell is not important. \({\rm erfc}\) is the complementary error function.

When the electrostatic charges are fitted according to the local atomic environments, the constraint of overall charge neutrality has to be fulfilled.

Some general rules for testing the convergence of the parameters of the Ewald summation:

• Convergence of the Ewald cutoff: The Ewald cutoff just depends on the value of \(\alpha\). It has to be chosen sufficiently large so that all neighbors (the neighbor list is defined by the Ewald cutoff) which are interacting with the reference atom are in the neighbor list. Examples: For \(\alpha=0.5\) a cutoff of 7 Bohr seems to be converged. For \(\alpha=0.2\) a cutoff of 15 Bohr should be used. The larger the Ewald cutoff, the better (but also the more expensive). A smaller \(\alpha\) results in a broader Gaussian, thus a smaller \(\alpha\) requires a larger cutoff. The required size of the cutoff is independent of the actual system size. A given pair (\(\alpha\)/Ewald cutoff) will give the same real space energy independent of the system size (is size-consistent), but changes of course if the structure is changed.

• Convergence of \(\alpha\): A large \(\alpha\) (more localized Gaussian) reduces the costs of the real space part because the Ewald cutoff for the neighbor list can be reduced. For a given \(\alpha\) the energy is size-consistent. There is no optimum value for \(\alpha\), because the convergence depends on the pair (\(\alpha/k_{max}\)). A large \(\alpha\) means that a larger \(k_{max}\) is needed for the same level of convergence.

• Convergence of \(k_{max}\): The value of \(k_{max}\) defines the costs of the reciprocal energy part. The larger \(k_{max}\) the better (but also the more expensive). The energy with respect to \(k_{max}\) is not size consistent (!) but depends on the lattice vectors. A larger cell requires a larger value of k. If the cell is doubled in one direction, also the k-mesh should be doubled (note that this relation is inverse to k-points in DFT). A smaller value of \(\alpha\) requires a smaller k, because a smaller \(\alpha\) corresponds to a wider Gaussian, which covers more of the electrostatics in the real space part. For \(\alpha=0.5\) and a cell size with \(a_{lat}=40\) Bohr \(k_{max}=\) 15 seems to be converged. A cell with \(a_{lat}=80\) Bohr then needs \(k_{max}=\) 30.

The following combinations seem to work for bulk cells up to \(a_{lat}=40\) Bohr:

• \(\alpha=0.5\), Ewald cutoff = 7.0 Bohr, \(k_{max}=\) 30.

• \(\alpha=0.2\), Ewald cutoff = 15.0 Bohr, \(k_{max}=\) 10.

The balance between \(\alpha\) and \(k_{max}\) should be found to minimize the CPU time. CAUTION: Be careful when setting \(k_{max}\) for slab structures because of the large c-dimension!!! The scaling of the reciprocal energy part is cubically with respect to \(k_{max}\). Calculating \(E_{recip}\) for 512 atoms in a 40 Bohr box costs 25 s with \(\alpha=0.5\)/\(k_{max}=\) 30 and 1 s with \(\alpha=\) 0.2/\(k_{max}=\) 10.

van-der-Waals Dispersion Energy Contribution

In RuNNer, the symmetry function cutoff radius governs which energetic contributions are accounted for by the neural network potential. However, non-covalent interactions may extend well beyond typical choices for the cutoff. Especially van-der-Waals (vdW) interactions, comprised of electrostatic Keesom interaction, Debye polarization interaction, and London dispersion stand out, due to their \(R_{ij}^{-6}\) decline with distance.

Many efficient methods have been proposed to include dispersion interactions in DFT-based calculations, e.g. by introducing fixed pairwise interatomic coefficients \(C_{6ij}r^{-6}\) as in the Grimme DFT-D2 and DFT-D3 corrections. Tkatchenko and Scheffler (TS) expanded this scheme with on-the-fly coefficient calculations based on the atomic charge density, i.e. they introduced a dependency on the chemical environment around an atom.

RuNNer sets out to offer a framework to calculate the dispersion energy in its prediction mode. The current implementation offers Grimme-D2 and a simple, environment-dependent correction inspired by the TS method. The main difference is that the current implementation of the TS correction does not yet offer on-the-fly polarizability-based \(C_{6ij}\) (environment-dependent) coefficients. Instead, the factors must be fixed before a run for all possible combinations of element pairs in input.nn. Consequently, it is very similar to the Grimme DFT-D2 method with user-defined \(C_{6ij}\) coefficients.

In both schemes, the total dispersion energy correction term is

\[ \begin{align} E_{\mathrm{disp}}= -\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N} \sum_{\textbf{L}}\frac{C_{6ij}}{r_{ij,L}^6}\cdot f(r_{ij,L})\,. \label{eq:vdwEnergyExpression} \end{align} \]

Here, \(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\).

As dispersion interactions can only occur between non-overlapping charge densities, a Fermi-type screening function \(f\) is employed in \(\eqref{eq:vdwEnergyExpression}\) to constrain pairwise interactions to distances larger than the combined vdW radii \(R_{ij}\) of atoms \(i\) and \(j\),

\[ \begin{align} f(r_{ij,L})= s_6 \cdot \left(1 + \mathrm{e}^{-d\left( \frac{r_{ij}}{s_{\mathrm{R}} R_{ij}} - 1 \right)}\right)^{-1}\,, \label{eq:vdwScreeningFunction} \end{align} \]

where three scaling parameters \(s_6\), \(d\), and \(s_{\mathrm{R}}\) define the shape of the function.

For the calculation of vdW dispersion force components or the vdW stress tensor, see vdW Dispersion Forces and vdW Stress, respectively.

Parameter choice

In the current implementation, RuNNer requires the user to define the \(C_{6ij}\) coefficients for each possible pair of elements in the structure, as well as the atomic vdW radii \(R_{i}\). Moreover, the parameters of the screening function \(s_6\), \(d\), and \(s_{\mathrm{R}}\) need to be chosen. Furthermore, the user needs to define an additional cutoff radius for the dispersion interactions.

• Choice of \(C_{6ij}\) coefficients: For same element interactions, \(C_{6ii}\) coefficients may be found in the literature. The TS implementation in VASP utilises the values determined by Grimme et al. ¹.

  The definition of reasonable heterogeneous coefficients \(C_{6ij}\), required for the current TS implementation, may be more challenging. This is because one needs to find a value which captures all possible atomic environments in the system in a meaningful way. The Grimme DFT-D2 combination rule,

  \[ \begin{align} C_{6ij}=\sqrt{C_{6ii}C_{6jj}}\,, \label{eq:vdwCombinationRule} \end{align} \]

  may give a valuable first estimate of the coefficients for the system at hand. For fairly homogeneous and ordered systems with little deviations in the atomic environments (of same-element species), one may also carry out a DFT Single Point calculation for a comparable system which applies the TS dispersion correction and use the mean value of all coefficients as an approximation of the real coefficient.

  Be aware of the correct units

  RuNNer requires atomic units for all input parameters. For example, the coefficients by Grimme et al ¹ are given in J nm\(^{-6}\) mol\(^{-1}\) which can be converted by

  \[ \begin{align*} 1\,\mathrm{J}\,\mathrm{nm}^{-6}\,\mathrm{mol}^{-1} = 17.34527759\,\mathrm{E}_{\mathrm{h}}\,\mathrm{a}_0^6\,\mathrm{atom}^{-1}\,. \end{align*} \]

• Choice of atomic vdW radii \(R_{i}\): Values for \(R_{i}\) may be found in the literature¹. Be aware of the correct units: RuNNer requires atomic units for all input parameters.

• Choice of screening function parameters \(s_6\), \(d\), and \(s_{\mathrm{R}}\): The shape of the screening function depends strongly on the functional used in the reference DFT calculations. In their original paper, TS determined the parameters for the PBE functional.²

  Parameter	Symbol	Value
  Global scaling	\(s_6\)	0.94 (0.94)
  Exponential scaling	\(R\)	0.75 (0.75)
  Damping	\(d\)	12 < \(d\) < 45 (20.0)

  The value in brackets shows the VASP default for the PBE functional.

• Choice of vdW cutoff radius: As the implementation is very fast, the cutoff radius may be set to a very large value. For most systems, the VASP default of 50.0Å (\(=94.5\;\mathrm{a}_{0}\)) should be a reasonable choice.

Code implementation

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

Name	Data Type	Description
vdwenergy	REAL*8	The total vdw dispersion correction energy of the structure.
atomvdw(:)	REAL*8	Dispersion energy contribution of each atom.

predict calls predictionvdw, which initiates a timer and goes on to invoke distribute_vdw_calculation.

This sub-module splits the structure into separate blocks of atoms which can be either handled consecutively (serial mode) or handed off to different MPI processes (parallel mode).

For each block of atoms, the neighbour list is determined in neighbor_vdw. Next, the main vdw calculation routine calc_tkatchenko_scheffler starts which iterates over all atoms (loop_block_of_atoms) and all their neighbours (loop_neighbors), summing up the pairwise dispersion energy according to \(\eqref{eq:vdwEnergyExpression}\) and storing it in atomvdw(:). Be aware that the neighbor list already accounts for periodic images in the vdW cutoff radius and prevents double counting. Furthermore, the contributions of the individual atoms in the block are summed up and stored in vdwenergy_block, which in turn is added to vdwenergy in the distribution routine.

---

1. S. Grimme, J. Comput. Chem. 27 (15), 1787-1799 (2006). ↩↩↩

2. A. Tkatchenko, M. Scheffler, Phys. Rev. Lett. 102 (7), 073005 (2009). ↩