4GHDNNPs
Overview
The total energy expression of 4GHDNNP consists of two parts, the shortrange energy and longrange electrostatic energy,
The electrostatic energy is calculated from the atomic charges, which are obtained from a charge equilibration scheme (Qeq) based on environmentdependent electronegativities that are predicted by atomic neural networks. The shortrange energy is a sum of atomic energies computed by atomic neural networks, which compared to 2DHDNNPs have an additional input node providing the atomic charge. With this additional information the shortrange atomic neural networks are able to accurately predict energetic changes due to modifications in the local electronic structure resulting from longrange charge transfer. Atomic forces are calculated by taking the negative analytic derivatives of the total energy in Eq. \(\eqref{eq:4GTotalEnergy}\) with respect to the atomic positions.
Charge equilibration
The charge equilibration method is based on the idea that the electrons in a system are distributed in a way that minimizes the total energy. This energy does not only include the electrostatic Coulomb terms but also a term that describes the local energy caused by some amount of charge on an atom. This local energy is usually described using a Taylor series up to second order. The expansion factor for the linear term is called the electronegativity (\(\chi_i\)), while the second order factor is called the atomic hardness (\(J_i\)). The charge on each atom is Gaussian distributed with width \(\sigma_i\). This results in the energy expression for the charge equilibration scheme
\(N_{\mathrm{at}}\) is the number of atoms in the system. \(E_{\mathrm{elec}}\) is the electrostatic energy resulting from the Gaussian charge distributions,
and
As only linear and quadratic terms in \(Q_i\) appear in \(E_{\mathrm{elec}}\) and \(E_{\mathrm{Qeq}}\), they can be expressed using matrix notation
and
with \(\mathbf{Q}\) being a column vector containing the atomic charges \(Q_i\), being a column vector of the electronegativities and \(\mathbf{A}\) being the matrix
It has to be noted that because the total energy of a continuous charge distribution is always positive, the matrix \(\mathbf{A}\) is positive definite, if the \(J_i>0\). The \(Q_i\) are now chosen, so that they minimize the energy \(E_{\mathrm{Qeq}}\) under the additional constraint of total charge conservation.
To solve this minimization problem, we set the derivatives with respect to the charges to zero,
Including the constraint of total charge conservation using a Lagrange multiplier \(\lambda\) we end up with the system of linear equations
which we can rewrite for simplicity as
where \(\mathbf{A'}\) and \(\mathbf{Q'}\) represent the (\(N_{\mathrm{at}}+1\times N_{\mathrm{at}}+1\)) matrix and column vector in the left hand side respectively, while \(\mathbf{b}\) is the column vector on the right hand side. The electronegativities \(\chi_i\) are predicted by neural networks, for each atom individually, depending on the local chemical environments. The hardness values of \(J_i\) are constant for a given element, and they are also optimized during the training of the neural networks.
Derivatives used for the calculation of the forces and for the neural network training
In this section we provide further details about some of the derivatives which are required for the calculation of the atomic forces as well as for the gradientbased optimization of the neural network parameters. The atomic force component \(F_{r_\alpha}\) is given as the negative derivative of the energy with respect to the atomic coordinate \(r_{\mathrm{\alpha}}\),
In this equation the partial derivatives of the atomic charges with respect to the atomic positions appear. As we will show in section , the calculation of these terms can be avoided for the determination of the forces. During the training phase, however, these derivatives as well as the derivatives of the charges w.r.t. the electronegativity and hardness are needed.
Calculation of the charge derivative with respect to the Cartesian coordinates
To calculate the \(\frac{\partial Q_i}{\partial r_{\mathrm{\alpha}}}\) we take the derivative with respect to the spatial coordinate \(r_{\mathrm{\alpha}}\) of the Qeq in Eq.$\eqref{eq:qeq_prime}. Reordering the terms yields
as well as a Lagrange multiplier that ensures \(\sum_i \frac{\partial Q_i}{\partial r_\alpha} = 0\). To obtain all the 3\(N_{\mathrm{at}}\) required derivatives we will have to solve 3 \(N_{\mathrm{at}}\) linear equation systems of size \(N_{\mathrm{at}}\)+1. However, this can be avoided as explained in the following section, which allows the calculation of the total force by only solving one linear equation system.
## Calculation of the charge derivative with respect to electronegativity and hardness
With a similar procedure we can calculate the equation systems for the derivatives w.r.t. the electronegativity and hardness
and
\(\mathbf{\delta}_i\) is a vector filled with zeros except entry \(i\), which is one. As before, a Lagrange multiplier will be necessary, to ensure that the sum of the derivatives adds up to zero.
An efficient method for the force computation
In the last section we showed how the forces can be calculated using the partial derivatives of the charges w.r.t. the atomic coordinates. This calculation is computationally expensive, since 3 \(N_{\mathrm{at}}\) linear equation systems need to be solved. This can be avoided by exploiting a method that allows the calculation of the forces by solving only one linear equation system instead.
The total energy of the 4GHDNNP is a function of the atomic coordinates (\(\V{R}\)) and the charges (\(\mathbf{Q}\)), which also depend on the atomic coordinates,
We now define an auxiliary function \(L\) with
Here \(\sum_{j=1}^{N_{\mathrm{at}}+1} A{'}_{ij} Q{'}_j  b_i\) are the differences of the left hand side minus the right hand sides of Eq. \eqref{eq:qeq_prime}, which were solved to determine the charges \(Q_i\). These terms are therefore always zero, making \(L\) equal to \(E_{\mathrm{total}}\). We now choose \(\mathbf{\lambda}\) such that the partial derivatives \(\frac{\partial L}{\partial Q{'}_i}\) are zero,
For this we solve the linear equation system
Note that \(\mathbf{A{'}}\) is a symmetric matrix. We now turn to the derivative \(\frac{dL}{dr_{\alpha}}\), which is equal to \(\frac{dE_{\mathrm{total}}}{dr_{\alpha}}\).
Rearranging the equation yields
The term \(\frac{\partial E_{\mathrm{total}}}{\partial Q{'}_i} + \sum_j A{'}_{ij} \lambda_j\) is zero by definition of \(\mathbf{\lambda}\) and can therefore be omitted, which leads to the expression
Charge equilibration for periodic systems
The Qeq equations for periodic boundary conditions are essentially identical to the corresponding equations for free boundary conditions, and the main difference is the calculation of the matrix \(\mathbf{A}\). Because of the periodic boundary conditions we have to resort to an Ewald summation to calculate the electrostatic interaction energy.
The basic idea of Ewald summation is, that by placing Gaussian charges of the opposite sign on each of the point charges, the remaining electrostatic interaction becomes shortranged. This shortranged energy can then be calculated in real space (\(E_{\mathrm{real}}\)). We then subtract the interaction energy of the auxiliary Gaussian charges again to obtain the desired total energy of the point charges. This interaction energy of the Gaussians can be efficiently calculated in reciprocal space, resulting in the energies \(E_{\mathrm{recip}}\) and \(E_{\mathrm{self}}\). The electrostatic energy of \(N_{\mathrm{at}}\) point charges can be hence calculated as
The real space part is given by
Here, \(N_{\mathrm{neig}}\) indicates, that the sum goes over all neighbouring atoms withing the real space cutoff radius \(r_{\mathrm{cut}}\). \(r_{ij}\) is the distance between atoms \(i\) and \(j\). The reciprocal space part is
with
\(V\) being the volume of the unit cell and the sum going over all reciprocal lattice points inside reciprocal space cutoff radius \(r_{\mathrm{cut}}^{\mathrm{recip}}\). Finally, the selfinteraction correction is
In these equations \(\eta\) is the standard deviation of the Gaussian charges, which are placed on the point charges to remove the longrange interactions.
Since we use Gaussian charge distributions for the charge equilibration process, the following terms have to be added that account for the different interaction in the shortrange part as well as for the self interaction of the Gaussian charges.
Here \(E_{\mathrm{elec}}^{\mathrm{pc}}\) is the electrostatic energy of the point charges as given above.
The important observation is that the total energy expression of the Ewald summation contains only terms of the form \(\frac{1}{2} e_{ij} Q_i Q_j\). By calculating the individual coefficients \(e_{ij}\) we can therefore construct the matrix \(\mathbf{E}\), so that
Including the terms for the hardness and adding the electronegativity results in a formalism equivalent to that of the Qeq method for free boundary condition.
Calculation of of the derivatives of the coefficient matrix
The differentiation of the above equation allows us to calculate the derivatives \(\frac{dA_{ij}}{dr_{\alpha}}\). Explicit calculation of these derivatives however can be computationally expensive, as they are quite numerous (\(3 N_{\mathrm{at}}^3\)). Most of these coefficients however are zero, and the matrix is very sparse. As only the product \(\sum_{j=1} \frac{dA_{ij}}{dr_{\alpha}} Q_j\) is ever needed in our computations, explicit calculation can be avoided and only the non zero terms have to be considered.
As \(\sum_{j=1} \frac{dA_{ij}}{dr_{\alpha}} Q_j\) has to be calculated anyways for the efficient computations of the forces, we can also use it to calculate the electrostatic forces.
Calculation of the electrostatic energies and forces boils down to the following.
 Calculate symmetry functions.
 Calculate environment depend electronegativities using a first set of neural networks.
 Construct the matrix \(\mathbf{A}\)
 Calculate partial charges \(Q_i\) by solving the system of \(N_{\mathrm{at}}+1\) linear equations (Eq. \(\eqref{eq:qeq}\))

Feed the atomic charges into the atomic NNs to calculate the short range energy and forces

Use the efficient method (see above) to calculate the total force by solving one more \(N_{\mathrm{at}} + 1\) dimensional linear equation system.
The only difference for the periodic case is in the step 3 and 6, where Ewald summation has to be used.
Shortrange part
In the short range neural network, we also include the atomic charge via an additional input neuron, such that the atomic energy contribution also depends on global charge distributions. The expression of atomic energies is very similar to the 2GHDNNP and it can be expressed as a function of symmetry functions and atomic charges. The atomic forces can be calculated by taking the derivatives of the energy with respect to the atomic positions
where \(F_{\mathrm{\alpha}}\) and \(G_{j,k}\) represent the force component \(\mathrm{\alpha}\) and the \(k^{\mathrm{th}}\) symmetry function of atom \(j\) respectively. In addition, \(N_{\mathrm{SF},i}\), \(N_{\mathrm{neig},i}\) equal to number of symmetry functions and neighbors of atom \(i\) and \(N_{\mathrm{at}}\) the total number of atom. Note that \(N_{\mathrm{neig},i}\) includes the atom \(i\) itself. If the method described in section is used, the last term, which includes the partial derivatives of the atomic charges w.r.t. the atomic coordinates, can be avoided.