Science Foundation (CNS-1136075, DMS-1555072, DMS-1736364, CMMI-1634832, and
CMMI-1560834), Brookhaven National Laboratory subcontract 382247, ARO/MURI
grant W911NF-15-1-0562, and U.S. Department of Energy (DOE) Office of Science
Advanced Scientific Computing Research program DE-SC0021142.

## Credit author statement

Hugo Esquivel: Conceptualization, Data curation, Formal analysis,
Investigation, Methodology, Resources, Software, Validation, Visualization,
Writing–‘original draft’, Writing–‘review and editing’. Arun Prakash: Funding
acquisition, Project administration, Resources, Supervision, Writing–‘review
and editing’. Guang Lin: Funding acquisition, Project administration,
Resources, Supervision, Writing–‘review and editing’.

## References

  * [1] D. P. Kroese, T. Taimre, Z. I. Botev, Handbook of Monte Carlo methods, Vol. 706, John Wiley & Sons, 2013.
  * [2] G. Fishman, Monte Carlo: Concepts, algorithms, and applications, Springer Science & Business Media, 2013.
  * [3] R. Y. Rubinstein, D. P. Kroese, Simulation and the Monte Carlo method, Vol. 10, John Wiley & Sons, 2016.
  * [4] I. Babuška, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM Journal on Numerical Analysis 45 (3) (2007) 1005–1034.
  * [5] F. Nobile, R. Tempone, C. G. Webster, A sparse grid stochastic collocation method for partial differential equations with random input data, SIAM Journal on Numerical Analysis 46 (5) (2008) 2309–2345.
  * [6] J. Foo, X. Wan, G. E. Karniadakis, The multi-element probabilistic collocation method (ME-PCM): Error analysis and applications, Journal of Computational Physics 227 (22) (2008) 9572–9595.
  * [7] J. Foo, G. E. Karniadakis, Multi-element probabilistic collocation method in high dimensions, Journal of Computational Physics 229 (5) (2010) 1536–1557.
  * [8] A. L. Teckentrup, P. Jantsch, C. G. Webster, M. Gunzburger, A multilevel stochastic collocation method for partial differential equations with random input data, SIAM/ASA Journal on Uncertainty Quantification 3 (1) (2015) 1046–1074.
  * [9] J. D. Collins, W. T. Thomson, The eigenvalue problem for structural systems with statistical properties, AIAA journal 7 (4) (1969) 642–648.
  * [10] W. K. Liu, T. Belytschko, A. Mani, Probabilistic finite elements for nonlinear structural dynamics, Computer Methods in Applied Mechanics and Engineering 56 (1) (1986) 61–81.
  * [11] W. K. Liu, T. Belytschko, A. Mani, Random field finite elements, International journal for numerical methods in engineering 23 (10) (1986) 1831–1845.
  * [12] M. Kleiber, T. D. Hien, The stochastic finite element method: Basic perturbation technique and computer implementation, Wiley, 1992.
  * [13] M. Shinozuka, G. Deodatis, Response variability of stochastic finite element systems, Journal of Engineering Mechanics 114 (3) (1988) 499–519.
  * [14] F. Yamazaki, M. Shinozuka, G. Dasgupta, Neumann expansion for stochastic finite element analysis, Journal of Engineering Mechanics 114 (8) (1988) 1335–1354.
  * [15] G. Deodatis, Weighted integral method. I: Stochastic stiffness matrix, Journal of Engineering Mechanics 117 (8) (1991) 1851–1864.
  * [16] G. Deodatis, M. Shinozuka, Weighted integral method. II: Response variability and reliability, Journal of Engineering Mechanics 117 (8) (1991) 1865–1877.
  * [17] D. Xiu, Numerical Methods for Stochastic Computations: A Spectral Method Approach, Princeton University Press, 2010.
  * [18] O. P. Le Maître, O. M. Knio, Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics, Springer Science & Business Media, 2010.
  * [19] T. J. Sullivan, Introduction to Uncertainty Quantification, Vol. 63, Springer, 2015\.
  * [20] N. Wiener, The homogeneous chaos, American Journal of Mathematics 60 (4) (1938) 897–936.
  * [21] N. Wiener, A. Wintner, The discrete chaos, American Journal of Mathematics 65 (2) (1943) 279–298.
  * [22] N. Wiener, Nonlinear problems in random theory, MIT press, 1958.
  * [23] V. Volterra, M. Tomassetti, F. Zarlatti, Leçons sur les équations intégrales et les équations intégro-différentielles, Gauthier-Villars Paris, 1913.
  * [24] R. H. Cameron, W. T. Martin, The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals, Annals of Mathematics (1947) 385–392.
  * [25] D. Lucor, D. Xiu, G. Karniadakis, Spectral representations of uncertainty in simulations: Algorithms and applications, in: Proceedings of the International Conference on Spectral and High Order Methods (ICOSAHOM-01), Uppsala, Sweden, 2001.
  * [26] D. Xiu, G. E. Karniadakis, The Wiener-Askey polynomial chaos for stochastic differential equations, SIAM journal on scientific computing 24 (2) (2002) 619–644.
  * [27] R. Field, M. Grigoriu, A new perspective on polynomial chaos, Computational Stochastic Mechanics, Spanos & Deodatis eds, Millpress (2003) 199–205.
  * [28] K. Itô, Multiple Wiener integral, Journal of the Mathematical Society of Japan 3 (1) (1951) 157–169.
  * [29] H. Ogura, Orthogonal functionals of the Poisson process, IEEE Transactions on Information Theory 18 (4) (1972) 473–481.
  * [30] A. Segall, T. Kailath, Orthogonal functionals of independent-increment processes, IEEE Transactions on Information Theory 22 (3) (1976) 287–298.
  * [31] M. Schetzen, Nonlinear system modeling based on the Wiener theory, Proceedings of the IEEE 69 (12) (1981) 1557–1573.
  * [32] D. D. Engel, The multiple stochastic integral, Vol. 265, American Mathematical Soc., 1982.
  * [33] R. Ghanem, P. D. Spanos, Polynomial chaos in stochastic finite elements, Journal of Applied Mechanics 57 (1) (1990) 197–202.
  * [34] R. G. Ghanem, P. D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag New York, 1991.
  * [35] R. Ghanem, P. D. Spanos, A stochastic Galerkin expansion for nonlinear random vibration analysis, Probabilistic Engineering Mechanics 8 (3-4) (1993) 255–264.
  * [36] R. Ghanem, Stochastic finite elements with multiple random non-Gaussian properties, Journal of Engineering Mechanics 125 (1) (1999) 26–40.
  * [37] O. M. Knio, H. N. Najm, R. G. Ghanem, et al., A stochastic projection method for fluid flow. I: Basic formulation, Journal of computational Physics 173 (2) (2001) 481–511.
  * [38] O. P. Le Maıtre, M. T. Reagan, H. N. Najm, R. G. Ghanem, O. M. Knio, A stochastic projection method for fluid flow. II: Random process, Journal of computational Physics 181 (1) (2002) 9–44.
  * [39] D. Xiu, G. E. Karniadakis, Modeling uncertainty in flow simulations via generalized polynomial chaos, Journal of computational physics 187 (1) (2003) 137–167.
  * [40] B. V. Asokan, N. Zabaras, Using stochastic analysis to capture unstable equilibrium in natural convection, Journal of Computational Physics 208 (1) (2005) 134–153.
  * [41] O. Knio, O. Le Maitre, Uncertainty propagation in CFD using polynomial chaos decomposition, Fluid Dynamics Research 38 (9) (2006) 616.
  * [42] H. N. Najm, Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics, Annual review of fluid mechanics 41 (2009) 35–52.
  * [43] X. Wan, G. E. Karniadakis, An adaptive multi-element generalized polynomial chaos method for stochastic differential equations, Journal of Computational Physics 209 (2) (2005) 617–642.
  * [44] X. Wan, G. E. Karniadakis, Beyond Wiener-Askey expansions: Handling arbitrary PDFs, Journal of Scientific Computing 27 (1-3) (2006) 455–464.
  * [45] X. Wan, G. E. Karniadakis, Multi-element generalized polynomial chaos for arbitrary probability measures, SIAM Journal on Scientific Computing 28 (3) (2006) 901–928.
  * [46] X. Ma, N. Zabaras, An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations, Journal of Computational Physics 228 (8) (2009) 3084–3113.
  * [47] J. D. Jakeman, A. Narayan, D. Xiu, Minimal multi-element stochastic collocation for uncertainty quantification of discontinuous functions, Journal of Computational Physics 242 (2013) 790–808.
  * [48] X. Wan, G. E. Karniadakis, Long-term behavior of polynomial chaos in stochastic flow simulations, Computer methods in applied mechanics and engineering 195 (41-43) (2006) 5582–5596.
  * [49] N. Agarwal, N. R. Aluru, A domain adaptive stochastic collocation approach for analysis of MEMS under uncertainties, Journal of Computational Physics 228 (20) (2009) 7662–7688.
  * [50] G. Kewlani, K. Iagnemma, A multi-element generalized polynomial chaos approach to analysis of mobile robot dynamics under uncertainty, in: 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems, IEEE, 2009, pp. 1177–1182.
  * [51] P. Prempraneerach, F. S. Hover, M. S. Triantafyllou, G. E. Karniadakis, Uncertainty quantification in simulations of power systems: Multi-element polynomial chaos methods, Reliability Engineering & System Safety 95 (6) (2010) 632–646.
  * [52] S. Oladyshkin, H. Class, R. Helmig, W. Nowak, A concept for data-driven uncertainty quantification and its application to carbon dioxide storage in geological formations, Advances in Water Resources 34 (11) (2011) 1508–1518.
  * [53] G. Kewlani, J. Crawford, K. Iagnemma, A polynomial chaos approach to the analysis of vehicle dynamics under uncertainty, Vehicle system dynamics 50 (5) (2012) 749–774.
  * [54] T. P. Sapsis, P. F. Lermusiaux, Dynamically orthogonal field equations for continuous stochastic dynamical systems, Physica D: Nonlinear Phenomena 238 (23) (2009) 2347–2360.
  * [55] A. Papoulis, Probability, Random Variables and Stochastic Processes, McGraw-Hill, 1965.
  * [56] P. Holmes, J. L. Lumley, G. Berkooz, C. W. Rowley, Turbulence, coherent structures, dynamical systems and symmetry, Cambridge university press, 2012.
  * [57] M. Choi, T. P. Sapsis, G. E. Karniadakis, A convergence study for SPDEs using combined polynomial chaos and dynamically-orthogonal schemes, Journal of Computational Physics 245 (2013) 281–301.
  * [58] M. Ueckermann, P. F. Lermusiaux, T. P. Sapsis, Numerical schemes for dynamically orthogonal equations of stochastic fluid and ocean flows, Journal of Computational Physics 233 (2013) 272–294.
  * [59] M. Cheng, T. Y. Hou, Z. Zhang, A dynamically bi-orthogonal method for time-dependent stochastic partial differential equations. I: Derivation and algorithms, Journal of Computational Physics 242 (2013) 843–868.
  * [60] M. Cheng, T. Y. Hou, Z. Zhang, A dynamically bi-orthogonal method for time-dependent stochastic partial differential equations. II: Adaptivity and generalizations, Journal of Computational Physics 242 (2013) 753–776.
  * [61] E. Musharbash, F. Nobile, T. Zhou, Error analysis of the dynamically orthogonal approximation of time dependent random PDEs, SIAM Journal on Scientific Computing 37 (2) (2015) A776–A810.
  * [62] M. Gerritsma, J.-B. Van der Steen, P. Vos, G. Karniadakis, Time-dependent generalized polynomial chaos, Journal of Computational Physics 229 (22) (2010) 8333–8363.
  * [63] V. Heuveline, M. Schick, A hybrid generalized polynomial chaos method for stochastic dynamical systems, International Journal for Uncertainty Quantification 4 (1) (2014).
  * [64] D. M. Luchtenburg, S. L. Brunton, C. W. Rowley, Long-time uncertainty propagation using generalized polynomial chaos and flow map composition, Journal of Computational Physics 274 (2014) 783–802.
  * [65] H. C. Ozen, G. Bal, Dynamical polynomial chaos expansions and long time evolution of differential equations with random forcing, SIAM/ASA Journal on Uncertainty Quantification 4 (1) (2016) 609–635.
  * [66] P. J. Davis, P. Rabinowitz, Methods of numerical integration, Courier Corporation, 2007.
  * [67] E. Novak, K. Ritter, High dimensional integration of smooth functions over cubes, Numerische Mathematik 75 (1) (1996) 79–97.
  * [68] H.-J. Bungartz, M. Griebel, Sparse grids, Acta numerica 13 (2004) 147–269.
  * [69] G. Leobacher, F. Pillichshammer, Introduction to quasi-Monte Carlo integration and applications, Springer, 2014.
  * [70] W. Rudin, Real and complex analysis, McGraw-Hill, 1987.
  * [71] W. Cheney, D. Kincaid, Linear algebra: Theory and applications, 2nd Edition, Jones & Bartlett Learning, 2010.
  * [72] S. A. Smolyak, Quadrature and interpolation formulas for tensor products of certain classes of functions, in: Doklady Akademii Nauk, Vol. 148, Russian Academy of Sciences, 1963, pp. 1042–1045.
  * [73] J. C. Butcher, J. Butcher, The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods, Vol. 512, Wiley New York, 1987.
  * [74] The MathWorks Inc., MATLAB (2016).
URL www.mathworks.com

## Errata

This updated version of the article fixes a typo found on Page 5888Also on
Page 5 of the journal article..

> Throughout this paper, we assume that the components of the $d$-tuple random
> variable $\xi=(\xi^{1},\ldots,\xi^{d})$ are pairwise mutually independent
> and that the random domain $\Xi$ is a hypercube of $d$ dimensions obtained
> by performing a $d$-fold Cartesian product of intervals
> $\bar{\Xi}_{i}:=\xi^{i}(\Omega)$.


# Accelerating GW calculations through machine learned dielectric matrices

Mario G. Zauchner  Department of Materials, Thomas Young Centre, Imperial
College London, South Kensington Campus, London, SW7 2AZ, United Kingdom
Andrew Horsfield  Department of Materials, Thomas Young Centre, Imperial
College London, South Kensington Campus, London, SW7 2AZ, United Kingdom
Johannes Lischner<EMAIL_ADDRESS>Department of Materials, Thomas
Young Centre, Imperial College London, South Kensington Campus, London, SW7
2AZ, United Kingdom

###### Abstract

The GW approach produces highly accurate quasiparticle energies, but its
application to large systems is computationally challenging, which can be
largely attributed to the difficulty in computing the inverse dielectric
matrix. To address this challenge, we develop a machine learning approach to
efficiently predict density-density response functions (DDRF) in materials.
For this, an atomic decomposition of the DDRF is introduced as well as the
neighbourhood density-matrix descriptor both of which transform in the same
way under rotations. The resulting DDRFs are then used to evaluate
quasiparticle energies via the GW approach. This technique is called the ML-GW
approach. To assess the accuracy of this method, we apply it to hydrogenated
silicon clusters and find that it reliably reproduces HOMO-LUMO gaps and
quasiparticle energy levels. The accuracy of the predictions deteriorates when
the approach is applied to larger clusters than those included in the training
set. These advances pave the way towards GW calculations of complex systems,
such as disordered materials, liquids, interfaces and nanoparticles.

##  I Introduction

Density functional theory (DFT)[1, 2] has shown tremendous success in the
calculation of electronic ground-state properties. However, it is well known
that band gaps of solids and HOMO-LUMO gaps of molecules are often
significantly underestimated when computed using Kohn-Sham (KS) eigenvalues
[3, 4]. In order to remedy this issue, the GW method [5, 6, 7] is often
employed in which a self-energy correction to the DFT KS energies is computed.
The resulting quasiparticle energies are in excellent agreement with
experimental measurements for a wide range of materials. However, the large
numerical effort required for GW calculations and its unfavorable scaling with
system size restrict applications to relatively small systems [8, 9]. The most
expensive step is the computation of the interacting density-density response
function (DDRF) which is closely related to the inverse dielectric matrix. In
particular, the non-interacting DDRF is typically computed by carrying out a
slowly-converging summation over all unoccupied states [10, 11, 8].
Afterwards, the non-interacting DDRF must be inverted to calculate the
interacting DDRF.

These difficulties have led to the development of model DDRFs (or model
dielectric matrices). For example, Hybertsen and Louie constructed a model
dielectric matrix based on the assumption that the local screening response of
the material is similar to that of a homogeneous medium with the same local
density [12]. A similar model was also proposed by Cappellini et al. [13, 14].
However, it has proven difficult to generalize these model dielectric
functions to highly non-uniform systems, such as isolated molecules or nano-
clusters whose screening properties are substantially different from uniform
systems. To overcome this limitation, Rohlfing [9] proposed to express the
dielectric matrix as a sum of atomic contributions attributing a density
response resulting from a Gaussian-shaped charge density to each atom. This
model dielectric matrix contains a number of parameters which need to be
determined, for example by comparison to calculated RPA dielectric functions.

In recent years, machine learning (ML) techniques have been widely adopted to
predict scalar properties of materials, such as the total energy. A key
ingredient in ML approaches is the descriptor which parametrizes the atomic
and chemical structure of the material. Many descriptors used in computational
chemistry are explicitly constructed to be invariant under rotations and
translations: for example, ACE [15], SOAP [16], the Coulomb matrix [17, 18],
bag-of-bonds [19] or fingerprint-based descriptors have been shown to be
reliable descriptors for the prediction of scalar quantities. When predicting
tensors or functions, however, it is no longer sufficient to employ a
rotationally invariant descriptor. To alleviate this problem, Grisafi et al.
[20] developed a symmetry-adapted version of the SOAP kernel which is
equivariant under rotations and was successfully used in the prediction of
polarizability tensors and first hyperpolarizabilities [21, 20], dipole
moments [22] and electronic densities [23]. Several other groups also explored
ML approaches for the electronic density including Brockherde et al.[24],
Alred et al. [25] and Chandrasekaran and co-workers [26]. Moreover, the
construction of group-equivariant neural networks, such as Clebsch-Gordan
networks [27, 28, 29], tensor-field networks [30] and spherical convolutional
neural networks (CNNs) [31, 32] have seen significant developments in recent
years and the implementation of these methods has been significantly
simplified by frameworks such as e3NN [33] developed by Geiger et al. [34],
thus providing promising alternatives to the symmetry-adapted SOAP for the
learning of functions.

To the best of our knowledge, however, there has been no attempt to develop ML
models for the prediction of non-local response functions, such as the DDRF.
Predicting such quantities is a formidable challenge: for example, the DDRF of
a small silicon cluster can be tens of gigabytes in size when represented in a
plane-wave basis even when a modest plane-wave cutoff is used. To address this
problem, we introduce a decomposition of the DDRF into atomic contributions
which can be predicted using ML techniques. To ensure that the ML model
appropriately incorporates the transformation properties of the DDRF, we also
develop a new descriptor called neighbourhood density-matrix (NDM) which
transforms in the same way as the DDRF under rotations and is used in
conjunction with a dense neural network to predict the atomic contributions to
the DDRF. We then use the ML DDRFs to carry out GW calculations of
hydrogenated silicon clusters. This approach which we refer to as the ML-GW
method produces accurate GW quasiparticle energies at a significantly reduced
computational cost compared to standard implementations.

##  II Results

###  II.1 Theoretical results

The GW method yields accurate quasiparticle energies by applying a self-energy
correction to the mean-field KS energy levels. The GW self-energy
$\Sigma(\mathbf{r,r^{\prime}},\omega)$ is calculated from the one-electron
Green’s function $G(\mathbf{r,r^{\prime}},\omega)$ and the screened Coulomb
interaction $W(\mathbf{r,r^{\prime}},\omega)$ according to [8, 7, 35]

$\Sigma(\mathbf{r,r^{\prime}},\omega)=\frac{i}{2\pi}\int
e^{-i\delta\omega^{\prime}}G(\mathbf{r,r^{\prime}},\omega+\omega^{\prime})W(\mathbf{r,r^{\prime}},\omega^{\prime})d\omega^{\prime}$
(1)

with $\delta$ denoting a positive infinitesimal. The screened Coulomb
interaction is in turn computed from the bare Coulomb interaction
$v(\mathbf{r,r^{\prime}})$ and the inverse dielectric matrix
$\epsilon^{-1}(\mathbf{r,r^{\prime}},\omega)$ via

$W(\mathbf{r,r^{\prime}},\omega)=\int\epsilon^{-1}(\mathbf{r},\mathbf{r}_{2},\omega)v(\mathbf{r}_{2},\mathbf{r^{\prime}})d\mathbf{r}_{2},$
(2)

which demonstrates that the dielectric matrix constitutes a key ingredient in
GW calculations. It can be obtained from the interacting DDRF
$\chi(\mathbf{r},\mathbf{r}^{\prime},\omega)$ according to

$\epsilon^{-1}(\mathbf{r,r^{\prime}},\omega)=\delta(\mathbf{r,r^{\prime}})+\int
v(\mathbf{r},\mathbf{r}_{2})\chi(\mathbf{r}_{2},\mathbf{r^{\prime}},\omega)d\mathbf{r}_{2}.$
(3)

In the remainder of this paper, we will assume that the frequency dependence
of the dielectric matrix can be approximated by the generalized plamon-pole
approximation (GPP) [7, 36, 37]. As a consequence, only the static DDRF
$\chi(\mathbf{r},\mathbf{r}^{\prime})\equiv\chi(\mathbf{r},\mathbf{r}^{\prime},\omega=0)$
needs to be determined.

Within the random-phase approximation (RPA), the interacting static DDRF is
given by

$\displaystyle\chi(\mathbf{r},\mathbf{r}^{\prime})=$
$\displaystyle\chi_{0}(\mathbf{r},\mathbf{r}^{\prime})+$ (4)
$\displaystyle\int
d\mathbf{r}_{1}d\mathbf{r}_{2}\chi_{0}(\mathbf{r},\mathbf{r}_{1})v(\mathbf{r}_{1},\mathbf{r}_{2})\chi(\mathbf{r}_{2},\mathbf{r}^{\prime})$
(5)

with $\chi_{0}(\mathbf{r},\mathbf{r}^{\prime})$ denoting the static non-
interacting DDRF which is typically computed as a sum over empty and occupied
states [10, 11] according to

$\displaystyle\chi_{0}(\mathbf{r},\mathbf{r}^{\prime})=$
$\displaystyle\sum_{ij}\frac{f_{i}(1-f_{j})}{\epsilon_{i}-\epsilon_{j}}\times$
(6)
$\displaystyle\left[\phi_{i}^{*}(\mathbf{r})\phi_{j}(\mathbf{r})\phi^{*}_{j}(\mathbf{r}^{\prime})\phi_{i}(\mathbf{r}^{\prime})+\text{c.c.}\right].$
(7)

Here, $\epsilon_{i}$, $f_{i}$ and $\phi_{i}(\mathbf{r})$ denote the orbital
energy, occupancy and wavefunctions of the KS state $i$.

Equations (5) and (7) highlight the two main challenges in computing the DDRF:
(1) the calculation of the non-interacting DDRF requires a summation of all
empty states which is slowly converging and (2) the calculation of the
interacting DDRF requires a matrix inversion which scales unfavorably with
system size.

####  II.1.1 Atomic decomposition of the density-density response function

In order to bypass the expensive computation of the DDRF and pave the way
towards a machine learning approach, we propose to express
$\chi(\mathbf{r,r^{\prime}})$ as a sum of atomic contributions
$\chi_{i}(\mathbf{r,r^{\prime}})$ according to

$\chi(\mathbf{r,r^{\prime}})=\sum_{i=1}^{N}\chi_{i}(\mathbf{r,r^{\prime}}),$
(8)

where $i$ labels atoms and $N$ is the total number of atoms.

How this partitioning is achieved is not immediately obvious. However, the
atomic contributions to the DDRF should have the following properties: (1) the
atomic contributions should be localized in the vicinity of the corresponding
atom, (2) they should retain the global symmetry of $\chi$, i.e.
$\chi(\mathbf{r,r^{\prime}})=\chi(\mathbf{r^{\prime},r})$, and (3) they should
integrate to zero, i.e.
$\int\chi(\mathbf{r},\mathbf{r}^{\prime})d\mathbf{r}=\int\chi(\mathbf{r},\mathbf{r}^{\prime})d\mathbf{r^{\prime}}=0$,
to ensure that the change in the charge density induced by a perturbing
potential is overall charge neutral [8].

We start by expressing the DDRF in a localized basis set of real orbitals
$\\{\phi_{\alpha_{a}}^{a}(\mathbf{r})\\}$, where $a$ labels the atom on which
the basis function is centered and $\alpha_{a}$ indexes the orbital on site
$a$. We should note that exploiting locality in GW calculations through local
orbital representations has been done before, even as early as the first
practical applications of the GW method by Strinati et al. [6]. In this basis
the DDRF is given by

$\chi(\mathbf{r},\mathbf{r}^{\prime})=\sum_{a,\alpha_{a}}\sum_{b,\alpha_{b}}\chi_{\alpha_{a}\alpha_{b}}^{ab}\phi_{\alpha_{a}}^{a}(\mathbf{r})\phi_{\alpha_{b}}^{b}(\mathbf{r}^{\prime}),$
(9)

where $\chi_{\alpha_{a}\alpha_{b}}^{ab}$ is a symmetric matrix. This
expression suggests the following decomposition of the DDRF into atomic
contributions

$\chi_{i}(\mathbf{r,r^{\prime}})=\frac{1}{2}\sum_{\alpha_{i}}\sum_{b,\alpha_{b}}\bigg{(}\chi^{ib}_{\alpha_{i}\alpha_{b}}\phi_{\alpha_{i}}^{i}(\mathbf{r})\phi_{\alpha_{b}}^{b}(\mathbf{r^{\prime}})\\\
+\chi^{bi}_{\alpha_{b}\alpha_{i}}\phi_{\alpha_{b}}^{b}(\mathbf{r})\phi_{\alpha_{i}}^{i}(\mathbf{r^{\prime}})\bigg{)}.$
(10)

We refer to the representation of the DDRF in the basis
$\\{\phi^{a}_{\alpha_{a}}(\mathbf{r})\\}$ as 2-center DDRF (2C-DDRF) because
it contains pairs of basis functions which are centered on different atoms.

Using the symmetry of $\chi^{iw}_{\alpha_{i}\alpha_{w}}$ and the fact that the
basis functions are real, it can be easily verified that
$\chi_{i}(\mathbf{r,r^{\prime}})=\chi_{i}(\mathbf{r^{\prime},r})$. We can also
ensure that $\int\chi_{i}(\mathbf{r},\mathbf{r}^{\prime})d\mathbf{r}=0$ by
removing all s-orbitals from the basis: see computational methods section for
details. The locality of $\chi_{i}(\mathbf{r,r^{\prime}})$ is directly
inherited from the corresponding properties of the full DDRF. In particular,
we have found that the expansion coefficients
$\chi^{iw}_{\alpha_{i}\alpha_{w}}$ decay rapidly as the distance between atom
$i$ and atom $w$ increases [38].

We stress that this atomic representation of the DDRF is exact, i.e.
$\sum_{i}\chi_{i}(\mathbf{r},\mathbf{r}^{\prime})$ reproduces the full
interacting DDRF when the local basis sets is complete. However, the atomic
contributions to the DDRF contain contributions from pairs of basis functions
which are centered on different atoms, see Eq. (10). These contributions are
difficult to learn using atom-centered descriptors.

To make progress, we exploit the localization of
$\chi_{i}(\mathbf{r},\mathbf{r}^{\prime})$ and expand it in terms of a set of
basis functions
$\psi^{i}_{nlm}(\mathbf{r})=Y_{lm}(\hat{\mathbf{r}})R_{n}(|\mathbf{r}|)$ (with
$Y_{lm}$ denoting the spherical harmonics and $R_{n}$ a set of radial
functions) which are all centered on atom $i$ according to

$\chi_{i}(\mathbf{r,r^{\prime}})=\sum_{nlm}\sum_{n^{\prime}l^{\prime}m^{\prime}}\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}\\\
Y_{lm}(\hat{\mathbf{r}})Y_{l^{\prime}m^{\prime}}^{*}(\hat{\mathbf{r}}^{\prime})R_{n}(|\mathbf{r}|)R_{n^{\prime}}^{*}(|\mathbf{r^{\prime}}|)$
(11)

with $\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}$ denoting the expansion
coefficients given by

$\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}=\int\int
d\mathbf{r}d\mathbf{r^{\prime}}\chi_{i}(\mathbf{r,r^{\prime}})R^{*}_{n}(|\mathbf{r}|)R_{n^{\prime}}(|\mathbf{r^{\prime}}|)\\\
Y^{*}_{lm}(\mathbf{\hat{r}})Y_{l^{\prime}m^{\prime}}(\mathbf{\hat{r}^{\prime}}).$
(12)

These coefficients can be learned using a neural network based on atom-
centered descriptors. We refer to the representation of the DDRF in the basis
$\\{\psi^{i}_{nlm}(\mathbf{r})\\}$ as 1-center DDRF (1C-DDRF) because it only
contains pairs of basis functions centered on the same atom.

####  II.1.2 Neighbourhood density-matrix descriptor

As discussed in the introduction, it is not appropriate to use a scalar
descriptor (such as the standard SOAP descriptor [39]) that is invariant under
rotations to develop a ML model for the DDRF: the behaviour of the atomic
DDRFs under rotations is determined by their analytical form: see Eq. (11). In
particular, we show in the Appendix that the coefficients of the atomic DDRF
transform according to

$\tilde{\chi}^{(i)}_{nlm_{1}n^{\prime}l^{\prime}m_{2}}=\sum_{m,m^{\prime}}D^{l}_{m_{1}m}(\hat{R})D^{l^{\prime}*}_{m_{2}m^{\prime}}(\hat{R})\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}},$
(13)

where $\tilde{\chi}^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}$ denote the
coefficients of the tranformed DDRF, $\hat{R}$ is a rotation and
$D^{l}_{mm^{\prime}}(\hat{R})$ is a Wigner D-matrix [40].

Next, we construct the NDM descriptor, which transforms under rotations in the
same way as the atomic DDRF. The starting point for such a descriptor is a
non-local extension of the smooth neighbourhood density of atom $i$ of species
$\eta$ employed in the SOAP descriptor [16], defined as

$\rho^{\eta}_{i}(\mathbf{r,r^{\prime}})=\sum_{k\in\eta}\sum_{l\in\eta}e^{-\alpha(\mathbf{r}-\mathbf{r}_{k})^{2}}e^{-\alpha(\mathbf{r}^{\prime}-\mathbf{r}_{l})^{2}},$
(14)

where $k$ and $l$ run over atoms in the neighbourhood of atom $i$ within a
cut-off radius $R_{cut}$ and $\alpha$ is a hyperparameter which describes the
size of an atom. The NDM is then expanded in a basis of spherical harmonics
and radial basis functions $R_{n}(|\mathbf{r}|)$ according to

$\rho^{\eta}_{i}(\mathbf{r,r^{\prime}})=\sum_{nlm}\sum_{n^{\prime}l^{\prime}m^{\prime}}\rho^{(i,\eta)}_{nlmn^{\prime}l^{\prime}m^{\prime}}\\\
Y_{lm}(\hat{\mathbf{r}})Y_{l^{\prime}m^{\prime}}^{*}(\hat{\mathbf{r}}^{\prime})R_{n}(|\mathbf{r}|)R_{n^{\prime}}^{*}(|\mathbf{r^{\prime}}|),$
(15)

with $\rho^{(i,\eta)}_{nlmn^{\prime}l^{\prime}m^{\prime}}$ being expansion
coefficients. The above equation shows that the NDM transforms in the same way
as the atomic DDRF: see Appendix for additional details. Therefore, we use the
expansion coefficients as a descriptor for learning the DDRF.

We note that the NDM can be written as the product of two neighbourhood
densities
$\rho^{\eta}_{i}(\mathbf{r})=\sum_{k\in\eta}\exp\\{-\alpha{(\mathbf{r}-\mathbf{r}_{k})^{2}}\\}$
according to

$\rho^{\eta}_{i}(\mathbf{r,r^{\prime}})=\rho^{\eta}_{i}(\mathbf{r})\rho^{\eta}_{i}(\mathbf{r}^{\prime}).$
(16)

Similar to the NDM, $\rho^{\eta}_{i}(\mathbf{r})$ can be expanded in a basis
of spherical harmonics and radial basis functions $R_{n}(|\mathbf{r}|)$ with
coefficients $\rho^{(i,\eta)}_{nlm}$. It follows that

$\rho^{(i,\eta)}_{nlmn^{\prime}l^{\prime}m^{\prime}}=\rho^{(i,\eta)}_{nlm}\rho^{(i,\eta)}_{n^{\prime}l^{\prime}m^{\prime}},$
(17)

which demonstrates that the coefficients of the neighbourhood density contain
the same information as the coefficients of the neighbourhood density matrix.
Indeed, we have found in our calculations that both types of coefficients
perform equally when used as descriptors to predict the atomic DDRFs. We
further note that the coefficients of the 3-body version of the SOAP
descriptor $d^{(\eta)}_{nn^{\prime}l}$ can be obtained from the NDM using

$d^{(\eta)}_{nn^{\prime}l}=\sum_{l^{\prime}mm^{\prime}}\sqrt{\frac{8\pi^{2}}{2l+1}}\rho^{(i,\eta)}_{nlm}\rho^{(i,\eta)}_{n^{\prime}l^{\prime}m^{\prime}}\delta_{ll^{\prime}}\delta_{mm^{\prime}},$
(18)

in the case where there is no coupling between different atomic species
$\eta$.

###  II.2 Machine learning

We apply our ML approach for predicting DDRFs to hydrogenated silicon clusters
and then use the DDRFs to calculate GW quasiparticle energies for these
systems. We refer to this technique as the ML-GW approach. The atomic
positions of the clusters were constructed as described in the methods section
and then relaxed using DFT.

To establish the accuracy of this approach, we first investigate the error in
the GW quasiparticle energies resulting from the expansion of the DDRF in
terms of the intermediate local basis
$\\{\phi^{a}_{\alpha_{a}}(\mathbf{r})\\}$: see Eq. (9). Fig. 1 compares the
HOMO-LUMO gaps obtained from mean-field DFT-PBE calculations, a standard
plane-wave G0W0 calculation using a generalized plasmon-pole approximation [7,
37] and a G0W0 calculation using the 2C-DDRF, where the DDRF is expanded in
terms of a modified version of the admm-2 basis set [41]: see methods section.
The DFT-PBE results show that the HOMO-LUMO gap decreases with increasing
cluster size from $E_{g}\approx 4.8$ eV for the smallest cluster containing 10
Si atoms to $E_{g}\approx 3$ eV for the biggest cluster with almost 60 Si
atoms. This decrease is a consequence of quantum confinement effects which are
less pronounced for bigger clusters. The plane-wave GW HOMO-LUMO gaps show a
similar trend as function of cluster size, but the gaps are larger than the
DFT-PBE gaps by several electron volts. Interestingly, the GW corrections are
larger for smaller clusters than for larger clusters. As a consequence, the
reduction in the GW HOMO-LUMO gaps as a function of cluster size is larger
compared to the DFT-PBE result: in particular, the gap is as large as 8.6 eV
for the smallest clusters and shrinks to 5.5 eV for the largest clusters
corresponding to a decrease of 3.1 eV (compared to a decrease of 1.8 eV in the
DFT-PBE HOMO-LUMO gap energies). Similar results were obtained by Chelikowsky
et al. [42] who also carried out GW calculations hydrogenated Si clusters. In
particular, they found that the HOMO-LUMO gap shrinks from $\sim 9$ eV for a
10 Si atom cluster to $\sim 6.5$ eV for a 47 Si atom cluster. The GW results
obtained with the 2C-DDRF are qualitatively similar to the plane-wave GW
results. However, the HOMO-LUMO gaps that are obtained with this approach are
consistently $\sim 0.4$ eV smaller than the plane-wave results. This is a
consequence of the incompleteness of the local basis set

Figure 1: HOMO-LUMO gaps of hydrogenated silicon clusters from DFT-PBE Kohn-
Sham eigenvalues, plane-wave G0W0 and G0W0 calculations using the 2C-DDRF, see
section "Atomic decomposition of the density-density response function".

Next, we determine the 1C-DDRF. For the basis set we use solid harmonic
Gaussians with optimized decay coefficients: see methods section. Fig. 2 (a)
compares the HOMO-LUMO gaps from G0W0 calculations with the 1C-DDRF to those
obtained with the 2C-DDRF and also to plane-wave G0W0 results. For small
clusters, the HOMO-LUMO gaps obtained with the 1C-DDRF are smaller than those
obtained with the 2C-DDRF, while the opposite behaviour is observed for larger
cluster. The largest difference between the two methods is obtained for
clusters containing $\sim 40$ Si atoms. The root-mean-square error (RMSE) of
the 1C-basis results relative to the 2C-basis results is 0.22 eV and the RMSE
relative to the plane-wave results is 0.45 eV for all clusters. Fig. 2 (b)
shows the HOMO and LUMO quasiparticle energies. It can be seen that better
agreement with the plane-wave result is obtained for the LUMO than for the
HOMO.

(a)

(b)

Figure 2: (a) HOMO-LUMO gaps of hydrogenated silicon clusters from plane-wave
G0W0 and G0W0 calculations using the 2C-DDRF and G0W0 calculations using the
1C-DDRF, see section "Atomic decomposition of the density-density response
function". (b) HOMO and LUMO energies of hydrogenated Si clusters.

Fig. 3 (a) shows the quasiparticle energy corrections of the ten lowest
conduction orbitals and the ten highest valence orbitals from plane-wave G0W0
and G0W0 with the 1C-DDRF. The corrections obtained with the 1C-DDRF follow a
similar trend as those obtained from the plane-wave calculation. For the
unoccupied states, the quantitative agreement is better than for the occupied
states, but the 1C-DDRF results for the unoccupied states are scattered over a
larger energy range than the plane-wave results. To analyze the errors that
arise from the use of the 1C-DDRF in more detail, Fig. 3 (b) shows a two-
dimensional histogram of the difference in QP corrections between plane-wave
G0W0 and G0W0 with the 1C-DDRF. For the occupied states the differences are
mostly smaller than 0.4 eV, while they are somewhat smaller for the unoccupied
states. The RMSE over all energy levels is 0.32 eV.

(a)

(b)

Figure 3: (a) Quasiparticle corrections from plane-wave G0W0 and G0W0 with the
1C-DDRF for the 10 highest valence orbitals and the 10 lowest conduction
orbitals of hydrogenated silicon clusters. (b) Histogram of difference in
quasiparticle corrections from plane-wave G0W0 and G0W0 calculations with the
1C-DDRF for the 10 highest valence orbitals and the 10 lowest conduction
orbitals of hydrogenated silicon clusters. The mean-field energies are
referenced to the middle of the mean-field HOMO-LUMO gap.

Now that we have established the accuracy of the method used to generate the
training set, we use a dense neural network (NN) in conjunction with the NDM
descriptor to generate the coefficients of the 1C-DDRF according to

$\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}=f(\rho^{(i,Si)}_{nlm},\rho^{(i,H)}_{nlm}),$
(19)

where $f$ is the neural network function. The hydrogen and silicon environment
descriptors are concatenated into a single vector before being fed into the
neural network. A separate network is trained for Si and H contributions to
the DDRF. The exact architecture of the network as well as the practical
computation of the atomic decomposition and the descriptors are described in
the Methods section. To generate the training data for the neural network, we
start from the set of relaxed hydrogenated Si clusters that were studied
above. From each relaxed cluster, we generate six new configurations by
randomly displacing the atoms with the magnitude of the displacements being
drawn from a uniform distribution with a maximum of 0.1 Å. For these clusters,
we then calculate the 1C-DDRF.

Once the neural network is trained on the 1C-DDRF of the randomly displaced
clusters, we use it to calculate the 1C-DDRFs of the relaxed clusters and then
determine quasiparticle energies via the ML-GW approach. Fig. 4 compares the
HOMO-LUMO gaps from ML-GW and GW with explicitly calculated 1C-DDRFs. Except
for the smallest cluster, the ML-GW method accurately reproduces the HOMO-LUMO
gaps of the explicit GW calculations. The worse performance for the smallest
cluster is a consequence of the training set which contains a large number of
bigger clusters containing atomic environments that differ from those found in
the smallest clusters. The overall RMSE of the ML-GW method relative to the
explicit GW with the 1C-basis is only 0.15 eV, but reduces to 0.06 eV when the
smallest cluster is excluded.

Fig. 5 shows the difference in QP corrections between ML-GW and GW with the
1C-DDRF for the 10 highest valence states and 10 lowest conduction states. ML-
GW produces QP shifts for both valence and conductions states within 0.1 eV
from the explicit G0W0 with the 1C-DDRF. The majority of valence states
exhibit a positive error, while for conduction states, the error is largely
negative.

Figure 4: HOMO-LUMO gaps of hydrogenated silicon clusters from plane-wave G0W0
and G0W0 calculations using the 1C-DDRF and ML-G0W0.  Figure 5: Histogram of
difference in quasiparticle corrections from G0W0 using the 1C-DDRF and
ML-G0W0 for the 10 highest valence orbitals and the 10 lowest conduction
orbitals of hydrogenated silicon clusters. The mean-field energies are
referenced to the middle of the mean-field HOMO-LUMO gap. The energies of the
smallest cluster were excluded.

Fig. 6 compares the ML-G0W0 QP corrections to plane-wave G0W0 results. As
expected, the differences are very similar to those between plane-wave G0W0
and the explicit G0W0 with the 1C-basis. In particular, the RMSE is 0.35 eV
for all clusters and reduces to 0.30 eV when the smallest cluster is excluded.
This results demonstrates that the key obstacle to improving the ML-GW
approach is the development of a better basis set.

Figure 6: Histogram of difference in quasiparticle corrections from plane-wave
G0W0 and ML-G0W0 DDRF for the 10 highest valence orbitals and the 10 lowest
conduction orbitals of hydrogenated silicon clusters. The mean-field energies
are referenced to the middle of the mean-field HOMO-LUMO gap. The energies of
the smallest cluster were excluded.

Finally, we test the ability of the ML-GW approach to predict the
quasiparticle energies of clusters which are larger than those included in the
training data. For this, we only include clusters with up to $N_{max}$ Si
atoms in the training set with $N_{max}$ being 60, 50 and 40. Again, the
training set only include clusters with randomly displaced atoms and the test
set consists of the relaxed clusters. The predicted ML-GW for the whole set of
relaxed clusters is shown in Fig. 7. From this graph, it is clear that the
accuracy of the prediction for the largest clusters deteriorates as $N_{max}$
is reduced: while for $N_{max}=60$, the gaps and QP corrections for clusters
with more than 60 Si atoms are still highly accurate, larger differences are
observed for $N_{max}=50$. For $N_{max}=40$, errors as larger as 1 eV are
obtained for the gaps of clusters with around 50 Si atoms. Fig. 7(f) shows
that the large error in the gaps are a consequence of having a negative error
in the QP shifts for occupied states and a positive error in the shift for
unoccupied states. In other words: instead of a cancellation, we get an
accumulation of errors when computing HOMO-LUMO gaps.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: HOMO-LUMO gaps (left panels) and errors in quasiparticles shifts
(right panels) from explicit G0W0 calculations with the 1C-DDRF and from
ML-G0W0 trained on clusters containing up to $N_{m}ax=60$ Si atoms (upper
panels), $N_{m}ax=50$ Si atoms (middle panels) and $N_{m}ax=40$ Si atoms
(lower panels). The red vertical line indicates $N_{m}ax$. The panels on the
right hand side only contain results for clusters with more Si atoms than
$N_{max}$. The mean-field energies are referenced to the middle of the mean-
field HOMO-LUMO gap.

##  III Discussion

We have developed a machine learning approach to predict the interacting
density-density response function (DDRF) of materials. To achieve this, we
introduce a decomposition of the DDRF into atomic contributions which form the
output of a neural network. We also introduce the NDM descriptor which is a
generalization of the widely used SOAP descriptor [16]: instead of
symmetrizing the descriptor using a Haar integral over a symmetry group [43],
we construct the tensor product of the expansion coefficients of the
neighbourhood density which transforms under rotation in the same way as the
atomic contributions to the DDRF. Thus, while not fully covariant, our
approach is able to distinguish between different orientations of a chemical
environment, which is a key requirement for predicting functions, such as the
DDRF.

The machine learning technique for DDRFs is then combined with the GW
approach. The resulting approach is called the ML-GW approach. We apply this
method to hydrogenated silicon clusters. The ML-GW approach reproduces HOMO-
LUMO gaps and quasiparticle energies of GW calculations using the explicitly
calculated 1C-DDRF, i.e. the DDRF in a pair basis where the basis functions of
each pair are centered on the same atom, with an accuracy of about 0.1 eV. The
accuracy of the results deteriorates when it is applied to clusters which are
larger than those included in the training set.

However, the error of ML-GW is significantly larger when compared to standard
plane-wave GW results: HOMO-LUMO gaps are reproduced to within 0.5 eV, but the
error reduces to 0.4 eV when the smallest cluster is excluded from the test
set. These errors are comparable to those obtained by Rohlfing in his GW
calculations for silane using a model dielectric function [9].

These findings demonstrate that the main challenge towards improving the ML-GW
method is the construction of better local basis sets for the DDRF. The basis
used for the 2C-DDRF can be improved straightforwardly by using larger basis
sets, such as aug-admm-2, admm-3 or aug-admm-3 [41]. However, it is more
difficult to increase the basis used for the 1C-DDRF as this leads to linear
dependencies which deteriorate the predictive accuracy of the neural network.
This was also observed by Grisafi et. al. [23] when predicting the expansion
coefficients of the electronic density using the symmetry adapted SOAP kernel
[20]. In the future, we plan to explore the use of orthogonal radial basis
sets, such as Laguerre polynomials instead of solid harmonic Gaussians.

We expect that the ML-GW method can be applied to calculate quasiparticle
energies in systems that have so far been out of reach for standard
implementations. Examples include disordered materials, liquids, interfaces or
nanoparticles. It could also be combined with on-the-fly machine learning
methods [44] to perform GW calculations on molecular-dynamics snapshots to
determine finite-temperature quasiparticle energies.

##  IV Methods

###  IV.1 Data generation

The atomic structures of the hydrogenated silicon clusters were obtained in
the same way as described by Zauchner et. al. [45]: starting from the
Si123H100 cluster of the silicon Quantum Dot data set [46], we remove the
silicon atom furthest from the centre of the cluster, terminate the dangling
bonds with hydrogen atoms and relax the resulting structure using DFT. The
process is repeated until only 10 silicon atoms remain. In addition, we also
include the clusters with less than 123 Si atoms from the silicon Quantum Dot
data set. From this set of silicon clusters, only clusters with less than 60
silicon atoms were used in the training set for DDRF prediction. From each
cluster with less than 60 silicon atoms, we created six additional clusters in
which random displacements were added to the atomic positions. The magnitudes
of the displacements were drawn from a normal distribution with mean of 0 Å
and standard deviation of 0.1 Å. Finally, calculations were also carried out
for clusters with between 60 and 70 silicon atoms. These clusters are not part
of the training set, but are used to test the extrapolation capacity of the ML
approach.

###  IV.2 DFT and GW calculations

The DDRF and QP corrections were calculated using the BerkeleyGW software
package [47, 48]. Mean-field DFT calculations were performed using the Quantum
Espresso code [49, 50]. Norm-conserving pseudopotentials from the Quantum
Espresso Pseudopotential Library were used. The parameters of the DFT
calculations were the same as those used by Zauchner et. al. [45]: a plane-
wave cut-off of 65 Ry, and a supercell with sufficient vacuum to avoid
interactions between periodic images. For the calculation of the DDRF a total
of 1000 Kohn-Sham states were used in the summation. Also, a plane-wave cut-
off of 6 Ry and a truncated Coulomb interaction were used. The QP corrections
were calculated using the generalized plasmon-pole approximation (GPP) [7], an
explicit sum over 1000 Kohn-Sham states and also a static remainder correction
[51]. To calculate the HOMO and LUMO energies, the vacuum level was determined
by averaging the electrostatic potential over the faces of the supercell.

###  IV.3 Projection onto intermediate basis

We first use BerkeleyGW to calculate the inverse dielectric matrix
$\epsilon^{-1}_{\mathbf{GG^{\prime}}}$ in a plane-wave basis [48]. From this,
we determine the interacting DDRF via

$\chi_{\mathbf{GG^{\prime}}}=(\epsilon_{\mathbf{GG^{\prime}}}-\delta_{\mathbf{GG^{\prime}}})/v_{\mathbf{G}}$
(20)

with $v_{\mathbf{G}}$ being the Fourier transform of the truncated Coulomb
interaction.

Next, the DDRF in real space is obtained as

$\chi(\mathbf{r,r^{\prime}})=\frac{1}{V}\sum_{\mathbf{G,G^{\prime}}}e^{i\mathbf{G\cdot
r}}\chi_{\mathbf{GG^{\prime}}}e^{-i\mathbf{G^{\prime}\cdot r^{\prime}}},$ (21)

where $V$ is the volume of the supercell.

Starting from a set of real atom-centred basis functions
${\phi}^{i}_{\alpha_{i}}(\mathbf{r})$, where $\alpha_{i}$ labels the basis
function on atom $i$, we construct an orthogonal basis set
$\tilde{\phi}^{i}_{\alpha_{i}}(\mathbf{r})$

$\tilde{\phi}^{i}_{\alpha_{i}}(\mathbf{r})=\sum_{k}\sum_{\alpha_{k}}A_{ik}^{\alpha_{i}\alpha_{k}}{\phi}^{k}_{\alpha_{k}}(\mathbf{r}),$
(22)

where $A_{ik}^{\alpha_{i}\alpha_{k}}$ is the matrix of eigenvectors of the
overlap matrix. The coefficients of the DDRF when expanded in the
orthogonalized basis are

$\tilde{\chi}_{\alpha_{i}\alpha_{j}}^{ij}=\frac{1}{V}\sum_{\mathbf{G,G^{\prime}}}\chi_{\mathbf{G,G^{\prime}}}\\\
\times\int_{-\infty}^{{\infty}}\tilde{\phi}^{i}_{\alpha_{i}}(\mathbf{r})e^{i\mathbf{G\cdot
r}}d\mathbf{r}\int_{-\infty}^{{\infty}}e^{-i\mathbf{G^{\prime}\cdot
r^{\prime}}}\tilde{\phi}^{j}_{\alpha_{j}}(\mathbf{r^{\prime}})d\mathbf{r^{\prime}},$
(23)

where, due to the localised nature of the basis functions, we extended the
integral from an integral over the supercell to an integral over all space.
These integrals are proportional to the Fourier transforms of the basis
functions (or their complex conjugates).

We then transform back to the non-orthogonal localised basis set using Eq.
(22) to find

$\chi(\mathbf{r,r^{\prime}})=\sum_{\alpha_{i}\alpha_{j}}\sum_{ij}\tilde{\chi}_{\alpha_{i}\alpha_{j}}^{ij}\tilde{\phi}^{i}_{\alpha_{i}}(\mathbf{r})\tilde{\phi}^{j}_{\alpha_{j}}(\mathbf{r^{\prime}})=\\\
\sum_{\alpha_{k}\alpha_{l}}\sum_{kl}\sum_{\alpha_{i}\alpha_{j}}\sum_{ij}A_{ik}^{\alpha_{i}\alpha_{k}}A_{jl}^{\alpha_{i}\alpha_{k}}\tilde{\chi}_{\alpha_{i}\alpha_{j}}^{ij}{\phi}^{k}_{\alpha_{k}}(\mathbf{r}){\phi}^{l}_{\alpha_{l}}(\mathbf{r^{\prime}})\\\
=\sum_{\alpha_{k}\alpha_{l}}\sum_{kl}{\chi}_{\alpha_{k}\alpha_{l}}^{kl}{\phi}^{k}_{\alpha_{k}}(\mathbf{r}){\phi}^{l}_{\alpha_{l}}(\mathbf{r}^{\prime}),$
(24)

where we defined

${\chi}_{\alpha_{k}\alpha_{l}}^{kl}=\sum_{\alpha_{i}\alpha_{j}}\sum_{ij}\tilde{\chi}_{\alpha_{i}\alpha_{j}}^{ij}A_{ik}^{\alpha_{i}\alpha_{k}}A_{jl}^{\alpha_{i}\alpha_{k}}.$
(25)

The basis functions we employed are the real solid harmonic Gaussians as
defined in LibInt2 [52]

$\phi_{lm}(r,\theta,\phi)=N_{l}(\beta)r^{l}e^{-\beta
r^{2}}R_{lm}(\theta,\phi),$ (26)

where $\beta$ is a decay parameter, $N_{l}(\beta)$ is a normalization factor
and $R_{lm}$ are the real spherical harmonics given by [53]

$\displaystyle R_{lm}(\theta,\phi)=$ (27)
$\displaystyle\begin{cases}\frac{i}{\sqrt{2}}\left(Y_{l-|m|}(\theta,\phi)-(-1)^{m}Y_{l|m|}(\theta,\phi)\right)\text{
if }m<0\\\ Y_{lm}(\theta,\phi)\text{ if }m=0\\\
\frac{1}{\sqrt{2}}\left(Y_{l-|m|}(\theta,\phi)+(-1)^{m}Y_{l|m|}(\theta,\phi)\right)\text{
if }m>0,\end{cases}$ (28)

where $Y_{lm}(\theta,\phi)$ are the complex spherical harmonics with the
Condon-Shortley phase convention. Kuang and Lin showed that the Fourier
transform of the complex solid harmonic Gaussians is again a solid harmonic
Gaussian [54]

$\frac{1}{(2\pi)^{3/2}}\int d\mathbf{r}e^{-i\mathbf{G\cdot
r}}N_{l}(\beta)r^{l}e^{-\beta r^{2}}Y_{lm}(\hat{\mathbf{r}})\\\
=(-i)^{l}\tilde{N}_{l}(\beta)G^{l}e^{-G^{2}/(4\beta)}Y_{lm}(\hat{\mathbf{G}}),$
(29)

with $\tilde{N}_{l}(\beta)=N_{l}(\beta)/(2\beta)^{3/2}$. The Fourier transform
of the real solid harmonic Gaussians can then be easily computed using Eq.
(28).

The basis set used in this work is a modified version of the admm-2 basis set
[41] (see Appendix for details), in which the s-orbitals were removed and
contracted Gaussians were uncontracted into individual basis functions.
Removing the s-orbitals ensures that $\int
d\mathbf{r}\chi(\mathbf{r},\mathbf{r}^{\prime})=0$ since only the Fourier
transform of s-orbitals has a $\mathbf{G}=0$ contribution.

###  IV.4 Projection onto atomic basis

The fully atom-centred basis set also consists of solid harmonic Gaussians.
The basis set was constructed following the same procedure as in the DScribe
library [55], where individual basis functions are given by

$\psi_{nlm}(r,\theta,\phi)=N_{l}(\beta_{nl})r^{l}e^{-\beta_{nl}r^{2}}R_{lm}(\theta,\phi),$
(30)

where the basis set is truncated at a maximum angular momentum $l_{max}$ and a
maximum principal quantum number $n_{max}$. For silicon atoms we use
$l_{max}=n_{max}=4$. For hydrogen atoms we use $l_{max}=n_{max}=3$.

The exponents $\beta_{nl}$ are constructed such that the corresponding basis
functions decay to zero at a cutoff radius $R_{n}$, i.e.
$\beta_{nl}=-\ln(\frac{T}{R_{n}^{l}})/R_{n}^{2}$ with $T=10^{-3}$ Ål being a
threshold parameter. The cutoff radius $R_{n}=R_{i}+(R_{o}-R_{i})/n$ lies
between an inner radius $R_{i}$ and an outer radius $R_{o}$. For hydrogen
atoms, we used $R_{i}=0.1$ Å and $R_{o}=3.0$ Å and for silicon atoms, we used
$R_{i}=1.0$ Å and $R_{o}=8.0$ Å. Both $R_{i}$ and $R_{o}$ were optimized to
minimize linear dependencies in the basis set as such dependencies
significantly deteriorate the accuracy of the neural network predictions. A
similar observation was made by Grisafi et. al. [23] when learning electron
densities, although a different approach was taken to remedy this issue in
their work.

In order to compute the coefficients of the atomic contributions to the DDRF
in the fully atom-centered basis the same procedure as in the intermediate
basis was used: the basis was first orthogonalized by computing the
eigenvectors of the overlap matrix. Then the atomic DDRFs in the intermediate
basis were projected onto the orthogonalized fully-atom centred basis with
overlaps between the different basis functions being computed using LibInt-2
[52]. Then the atomic DDRFs were transformed back to the non-orthogonal basis
produced the desired coefficients
$\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}$.

###  IV.5 Descriptors

The basis set for neighbourhood densities was generated using the same
procedure as for the fully atom-centered basis for the DDRF. However,
s-orbitals were not removed and the basis functions of the admm-2 basis set
were not included. We used $R_{i}=1.0$ Å for both hydrogen and silicon atoms
and $R_{o}=4.0$ Å for hydrogen atoms and $R_{o}=9.0$ Å for silicon atoms. The
exponents of the Gaussians in Eq. (14) were set such that the standard
deviation of the Gaussians is 0.5 Å. LibInt-2 [52] was again used to compute
the required integrals for the projection.

###  IV.6 Neural network

A dense neural network with four hidden layers with 2000, 1500, 1000 and 2000
nodes respectively was constructed for both silicon and hydrogen atoms. Each
layer uses a Leaky-ReLu activation function with a leak parameter of 0.1. The
output layer was further symmetrized by adding its transpose. The loss used
was the mean-squared error between the predicted and true expansion
coefficients $\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}$. The neural
network was trained on the perturbed clusters for 20,000 epochs. We found that
adding dropout to the layers does not significantly improve the quasiparticle
energies resulting from the predictions which is likely due to the similarity
between the atomic environments in the training and test set.

##  V Data availability statement

The data that support the findings of this study are available upon reasonable
request from the authors.

##  VI Acknowledgements

This work was supported through a studentship in the Centre for Doctoral
Training on Theory and Simulation of Materials at Imperial College London
funded by the EPSRC (EP/L015579/1). We acknowledge the Thomas Young Centre
under grant number TYC-101. This work used the ARCHER2 UK National
Supercomputing Service via J.L.’s membership of the HEC Materials Chemistry
Consortium of UK, which is funded by EPSRC (EP/L000202).

## References

## References

  * [1] Hohenberg, P. & Kohn, W.  Inhomogeneous electron gas.  _Phys. Rev._ 136, B864–B871 (1964).  URL https://link.aps.org/doi/10.1103/PhysRev.136.B864.
  * [2] Kohn, W. & Sham, L. J.  Self-consistent equations including exchange and correlation effects.  _Phys. Rev._ 140, A1133–A1138 (1965).  URL https://link.aps.org/doi/10.1103/PhysRev.140.A1133.
  * [3] Sham, L. J. & Schlüter, M.  Density-functional theory of the energy gap.  _Phys. Rev. Lett._ 51, 1888–1891 (1983).  URL https://link.aps.org/doi/10.1103/PhysRevLett.51.1888.
  * [4] Schultz, P. A.  Theory of defect levels and the “band gap problem” in silicon.  _Phys. Rev. Lett._ 96, 246401 (2006).  URL https://link.aps.org/doi/10.1103/PhysRevLett.96.246401.
  * [5] Hedin, L.  Effect of electronic correlation on band structure of solids 30, 231–231 (1965).
  * [6] Strinati, G., Mattausch, H. J. & Hanke, W.  Dynamical aspects of correlation corrections in a covalent crystal.  _Phys. Rev. B_ 25, 2867–2888 (1982).  URL https://link.aps.org/doi/10.1103/PhysRevB.25.2867.
  * [7] Hybertsen, M. S. & Louie, S. G.  Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies.  _Phys. Rev. B_ 34, 5390–5413 (1986).  URL https://link.aps.org/doi/10.1103/PhysRevB.34.5390.
  * [8] Onida, G., Reining, L. & Rubio, A.  Electronic excitations: density-functional versus many-body green’s-function approaches.  _Rev. Mod. Phys._ 74, 601–659 (2002).  URL https://link.aps.org/doi/10.1103/RevModPhys.74.601.
  * [9] Rohlfing, M.  Electronic excitations from a perturbative $\text{LDA}+gdw$ approach.  _Phys. Rev. B_ 82, 205127 (2010).  URL https://link.aps.org/doi/10.1103/PhysRevB.82.205127.
  * [10] Adler, S. L.  Quantum theory of the dielectric constant in real solids.  _Phys. Rev._ 126, 413–420 (1962).  URL https://link.aps.org/doi/10.1103/PhysRev.126.413.
  * [11] Wiser, N.  Dielectric constant with local field effects included.  _Phys. Rev._ 129, 62–69 (1963).  URL https://link.aps.org/doi/10.1103/PhysRev.129.62.
  * [12] Hybertsen, M. S. & Louie, S. G.  Model dielectric matrices for quasiparticle self-energy calculations.  _Phys. Rev. B_ 37, 2733–2736 (1988).  URL https://link.aps.org/doi/10.1103/PhysRevB.37.2733.
  * [13] Cappellini, G., Sole R, D., Reining, L. & Bechstedt, F.  Model dielectric function for semiconductors.  _Physical review. B, Condensed matter_ 47, 9892–9895 (1993).
  * [14] Bechstedt, F., Sole], R. D., Cappellini, G. & Reining, L.  An efficient method for calculating quasiparticle energies in semiconductors.  _Solid State Communications_ 84, 765 – 770 (1992).  URL http://www.sciencedirect.com/science/article/pii/003810989290476P.
  * [15] Drautz, R.  Atomic cluster expansion for accurate and transferable interatomic potentials.  _Phys. Rev. B_ 99, 014104 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.99.014104.
  * [16] Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  _Phys. Rev. Lett._ 104, 136403 (2010).  URL https://link.aps.org/doi/10.1103/PhysRevLett.104.136403.
  * [17] Hansen, K. _et al._ Assessment and validation of machine learning methods for predicting molecular atomization energies.  _Journal of Chemical Theory and Computation_ 9, 3404–3419 (2013).  URL https://doi.org/10.1021/ct400195d.  PMID: 26584096, eprint https://doi.org/10.1021/ct400195d.
  * [18] Rupp, M., Tkatchenko, A., Müller, K.-R. & von Lilienfeld, O. A.  Fast and accurate modeling of molecular atomization energies with machine learning.  _Phys. Rev. Lett._ 108, 058301 (2012).  URL https://link.aps.org/doi/10.1103/PhysRevLett.108.058301.
  * [19] Hansen, K. _et al._ Machine Learning Predictions of Molecular Properties: Accurate Many-Body Potentials and Nonlocality in Chemical Space.  _The Journal of Physical Chemistry Letters_ 6, 2326–2331 (2015).  URL https://dx.doi.org/10.1021/acs.jpclett.5b00831.
  * [20] Grisafi, A., Wilkins, D. M., Csányi, G. & Ceriotti, M.  Symmetry-adapted machine learning for tensorial properties of atomistic systems.  _Phys. Rev. Lett._ 120, 036002 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.036002.
  * [21] Wilkins, D. M. _et al._ Accurate molecular polarizabilities with coupled cluster theory and machine learning.  _Proceedings of the National Academy of Sciences_ 116, 3401–3406 (2019).  URL https://www.pnas.org/content/116/9/3401.  eprint https://www.pnas.org/content/116/9/3401.full.pdf.
  * [22] Veit, M., Wilkins, D. M., Yang, Y., DiStasio, R. A. & Ceriotti, M.  Predicting molecular dipole moments by combining atomic partial charges and atomic dipoles.  _The Journal of Chemical Physics_ 153, 024113 (2020).  URL https://doi.org/10.1063/5.0009106.  eprint https://doi.org/10.1063/5.0009106.
  * [23] Grisafi, A. _et al._ Transferable machine-learning model of the electron density.  _ACS Central Science_ 5, 57–64 (2019).  URL https://doi.org/10.1021/acscentsci.8b00551.  eprint https://doi.org/10.1021/acscentsci.8b00551.
  * [24] Brockherde, F. _et al._ Bypassing the Kohn-Sham equations with machine learning.  _Nature Communications_ 8 (2017).  URL https://dx.doi.org/10.1038/s41467-017-00839-3.
  * [25] Alred, J. M., Bets, K. V., Xie, Y. & Yakobson, B. I.  Machine learning electron density in sulfur crosslinked carbon nanotubes.  _Composites Science and Technology_ 166, 3–9 (2018).  URL https://dx.doi.org/10.1016/j.compscitech.2018.03.035.
  * [26] Chandrasekaran, A. _et al._ Solving the electronic structure problem with machine learning.  _npj Computational Materials_ 5 (2019).  URL https://dx.doi.org/10.1038/s41524-019-0162-7.
  * [27] Kondor, R., Lin, Z. & Trivedi, S.  Clebsch–gordan nets: a fully fourier space spherical convolutional neural network.  In Bengio, S. _et al._ (eds.) _Advances in Neural Information Processing Systems_ , vol. 31 (Curran Associates, Inc., 2018).  URL https://proceedings.neurips.cc/paper/2018/file/a3fc981af450752046be179185ebc8b5-Paper.pdf.
  * [28] Kondor, R. & Trivedi, S.  On the generalization of equivariance and convolution in neural networks to the action of compact groups.  In Dy, J. & Krause, A. (eds.) _Proceedings of the 35th International Conference on Machine Learning_ , vol. 80 of _Proceedings of Machine Learning Research_ , 2747–2755 (PMLR, 2018).  URL https://proceedings.mlr.press/v80/kondor18a.html.
  * [29] Anderson, B., Hy, T. S. & Kondor, R.  Cormorant: Covariant molecular neural networks.  In Wallach, H. _et al._ (eds.) _Advances in Neural Information Processing Systems_ , vol. 32 (Curran Associates, Inc., 2019).  URL https://proceedings.neurips.cc/paper/2019/file/03573b32b2746e6e8ca98b9123f2249b-Paper.pdf.
  * [30] Thomas, N. _et al._ Tensor field networks: Rotation- and translation-equivariant neural networks for 3d point clouds.  _CoRR_ abs/1802.08219 (2018).  URL http://arxiv.org/abs/1802.08219.  eprint 1802.08219.
  * [31] Cohen, T. S., Geiger, M., Köhler, J. & Welling, M.  Spherical cnns.  _CoRR_ abs/1801.10130 (2018).  URL http://arxiv.org/abs/1801.10130.  eprint 1801.10130.
  * [32] Cohen, T. & Welling, M.  Group equivariant convolutional networks.  In Balcan, M. F. & Weinberger, K. Q. (eds.) _Proceedings of The 33rd International Conference on Machine Learning_ , vol. 48 of _Proceedings of Machine Learning Research_ , 2990–2999 (PMLR, New York, New York, USA, 2016).  URL https://proceedings.mlr.press/v48/cohenc16.html.
  * [33] Lapchevskyi, K. _et al._ Euclidean neural networks (e3nn) v1.0, version v1.0 (2020).  URL https://www.osti.gov//servlets/purl/1770279.
  * [34] Geiger, M. & Smidt, T.  e3nn: Euclidean neural networks (2022).  URL https://arxiv.org/abs/2207.09453.
  * [35] Hybertsen, M. S. & Louie, S. G.  First-principles theory of quasiparticles: Calculation of band gaps in semiconductors and insulators.  _Phys. Rev. Lett._ 55, 1418–1421 (1985).  URL https://link.aps.org/doi/10.1103/PhysRevLett.55.1418.
  * [36] Lischner, J., Sharifzadeh, S., Deslippe, J., Neaton, J. B. & Louie, S. G.  Effects of self-consistency and plasmon-pole models on $gw$ calculations for closed-shell molecules.  _Phys. Rev. B_ 90, 115130 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevB.90.115130.
  * [37] Sharifzadeh, S., Tamblyn, I., Doak, P., Darancet, P. T. & Neaton, J. B.  Quantitative molecular orbital energies within a g0w0 approximation.  _The European Physical Journal B_ 85, 323 (2012).  URL https://doi.org/10.1140/epjb/e2012-30206-0.
  * [38] Mussard, B. & Ángyán, J. G.  Relationships between charge density response functions, exchange holes and localized orbitals.  _Computational and Theoretical Chemistry_ 1053, 44 – 52 (2015).  URL http://www.sciencedirect.com/science/article/pii/S2210271X14004952.  Special Issue: Understanding structure and reactivity from topology and beyond.
  * [39] Bartók, A. P., Kondor, R. & Csányi, G.  On representing chemical environments.  _Phys. Rev. B_ 87, 184115 (2013).  URL https://link.aps.org/doi/10.1103/PhysRevB.87.184115.
  * [40] Rose, M.  _Elementary Theory of Angular Momentum_.  Dover Books on Physics Series (Dover Publications, Incorporated, 2013).  URL https://books.google.at/books?id=1c7lngEACAAJ.
  * [41] Kumar, C. _et al._ Accelerating kohn-sham response theory using density fitting and the auxiliary-density-matrix method.  _International Journal of Quantum Chemistry_ 118, e25639 (2018).  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/qua.25639.  eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.25639.
  * [42] Tiago, M. L. & Chelikowsky, J. R.  Optical excitations in organic molecules, clusters, and defects studied by first-principles green’s function methods.  _Phys. Rev. B_ 73, 205334 (2006).  URL https://link.aps.org/doi/10.1103/PhysRevB.73.205334.
  * [43] Langer, M. F., Goeßmann, A. & Rupp, M.  Representations of molecules and materials for interpolation of quantum-mechanical simulations via machine learning.  _npj Computational Materials_ 8, 41 (2022).  URL https://doi.org/10.1038/s41524-022-00721-x.
  * [44] Li, Z., Kermode, J. R. & De Vita, A.  Molecular dynamics with on-the-fly machine learning of quantum-mechanical forces.  _Phys. Rev. Lett._ 114, 096405 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevLett.114.096405.
  * [45] Zauchner, M. G., Forno, S. D., Cśanyi, G., Horsfield, A. & Lischner, J.  Predicting polarizabilities of silicon clusters using local chemical environments.  _Machine Learning: Science and Technology_ 2, 045029 (2021).  URL https://dx.doi.org/10.1088/2632-2153/ac2cfe.
  * [46] Barnard, A. W. & Hugh.  Silicon quantum dot data set.  _CSIRO_ v2. (2015).  URL https://doi.org/10.4225/08/5721BB609EDB0.
  * [47] Hybertsen, M. S. & Louie, S. G.  Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies.  _Phys. Rev. B_ 34, 5390–5413 (1986).  URL https://link.aps.org/doi/10.1103/PhysRevB.34.5390.
  * [48] Deslippe, J. _et al._ Berkeleygw: A massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures.  _Computer Physics Communications_ 183, 1269 – 1289 (2012).  URL http://www.sciencedirect.com/science/article/pii/S0010465511003912.
  * [49] Giannozzi, P. _et al._ Advanced capabilities for materials modelling with quantum espresso.  _Journal of Physics: Condensed Matter_ 29, 465901 (2017).  URL http://dx.doi.org/10.1088/1361-648X/aa8f79.
  * [50] Giannozzi, P. _et al._ QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials.  _Journal of Physics: Condensed Matter_ 21, 395502 (2009).  URL https://doi.org/10.1088%2F0953-8984%2F21%2F39%2F395502.
  * [51] Deslippe, J., Samsonidze, G., Jain, M., Cohen, M. L. & Louie, S. G.  Coulomb-hole summations and energies for g w calculations with limited number of empty orbitals: A modified static remainder approach.  _Physical Review B_ 87, 165124 (2013).
  * [52] Valeev, E. F.  Libint: A library for the evaluation of molecular integrals of many-body operators over gaussian functions.  http://libint.valeyev.net/ (2022).  Version 2.8.0.
  * [53] Schlegel, H. B. & Frisch, M. J.  Transformation between cartesian and pure spherical harmonic gaussians.  _International Journal of Quantum Chemistry_ 54, 83–87 (1995).
  * [54] Kuang, J. & Lin, C. D.  Molecular integrals over spherical gaussian-type orbitals: I.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ 30, 2529–2548 (1997).  URL https://doi.org/10.1088/0953-4075/30/11/007.
  * [55] Himanen, L. _et al._ DScribe: Library of descriptors for machine learning in materials science.  _Computer Physics Communications_ 247, 106949 (2020).  URL https://doi.org/10.1016/j.cpc.2019.106949.
  * [56] Jeevanjee, N.  _An Introduction to Tensors and Group Theory for Physicists_.  An Introduction to Tensors and Group Theory for Physicists (Birkhäuser Boston, 2011).  URL https://books.google.at/books?id=O2zuQNwRIHsC.

##  VII Appendix A

###  VII.1 Transformation of atomic density-density response function and NDM
under rotation

We study the transformation properties of the atomic contributions to the DDRF
under the action of a representation of the rotation group $D(\hat{R})$ [56,
40]:

$D({\hat{R}})\otimes
D^{\dagger}({\hat{R}})\chi_{i}(\mathbf{r,r^{\prime}})=\chi_{i}(\hat{R}^{-1}\mathbf{r},\hat{R}^{-1}\mathbf{r^{\prime}}).$
(31)

The transformation is defined by the tensor product representation
$D({\hat{R}})\otimes D^{\dagger}({\hat{R}})$ of two SO(3) representations.
Using the expansion of $\chi_{i}$ in the $\psi$-basis, we find

$\chi_{i}(\hat{R}^{-1}\mathbf{r},\hat{R}^{-1}\mathbf{r^{\prime}})=\sum_{nlm}\sum_{n^{\prime}l^{\prime}m^{\prime}}\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}\\\
R_{n}(r)R_{n^{\prime}}^{*}(r^{\prime})Y_{lm}(\hat{R}^{-1}\mathbf{\hat{r}})Y_{l^{\prime}m^{\prime}}^{*}(\hat{R}^{-1}\mathbf{\hat{r}^{\prime}})\\\
\\\
=\sum_{nlm}\sum_{n^{\prime}l^{\prime}m^{\prime}}\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}R_{n}(r)R_{n^{\prime}}^{*}(r^{\prime})\\\
\sum_{m_{1},m_{2}}D^{l*}_{mm_{1}}(\hat{R}^{-1})D^{l^{\prime}}_{m^{\prime}m_{2}}(\hat{R}^{-1})Y_{lm_{1}}(\mathbf{\hat{r}})Y_{l^{\prime}m_{2}}^{*}(\mathbf{\hat{r}^{\prime}}).$
(32)

Exchanging the summations of $m,m^{\prime}$ and $m_{1},m_{2}$ and using the
unitarity property of the Wigner-D matrices $D^{l}(\hat{R})$ [40] yields

$\chi_{i}(\hat{R}^{-1}\mathbf{r},\hat{R}^{-1}\mathbf{r^{\prime}})=\\\
\sum_{nlm_{1}}\sum_{n^{\prime}l^{\prime}m_{2}}R_{n}(r)R_{n^{\prime}}^{*}(r^{\prime})Y_{lm_{1}}(\mathbf{\hat{r}})Y_{l^{\prime}m_{2}}^{*}(\mathbf{\hat{r}^{\prime}})\\\
\sum_{m,m^{\prime}}D^{l}_{m_{1}m}(\hat{R})D^{l^{\prime}*}_{m_{2}m^{\prime}}(\hat{R})\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}.$
(33)

Thus, the coefficients of the atomic DDRF transform under rotation as (with
the tilde denoting the rotated coefficients)

$\tilde{\chi}^{(i)}_{nlm_{1}n^{\prime}l^{\prime}m_{2}}=\sum_{m,m^{\prime}}D^{l}_{m_{1}m}(\hat{R})D^{l^{\prime}*}_{m_{2}m^{\prime}}(\hat{R})\chi^{(i)}_{nlmn^{\prime}l^{\prime}m^{\prime}}$
(34)

or in matrix notation

$\tilde{\chi}^{(i)}_{nn^{\prime}ll^{\prime}}=\mathbf{D}^{l}(\hat{R})\mathbf{\chi}^{(i)}_{nn^{\prime}ll^{\prime}}\mathbf{D}^{l^{\prime}\dagger}(\hat{R}).$
(35)

By going through the same steps, it can be shown that the expansion
coefficients of neighbourhood density matrix obey

$\tilde{\rho}^{(i,\eta)}_{nlm_{1}n^{\prime}l^{\prime}m_{2}}=\sum_{m,m^{\prime}}D^{l}_{m_{1}m}(\hat{R})D^{l^{\prime}*}_{m_{2}m^{\prime}}(\hat{R})\rho^{(i,\eta)}_{nlmn^{\prime}l^{\prime}m^{\prime}}$
(36)

under rotation.

##  VIII Appendix B

###  VIII.1 Exponents of modified admm-2 basis

The exponents of the modified admm-2 basis sets for hydrogen and silicon atoms
are shown in Tables 1 and 2. The basis functions were normalized to 1.

Table 1: Hydrogen basis $l$ |  $\beta$ [$1/a_{0}^{2}$]
---|---
1 | 1.0
1 | 0.457639
Table 2: Silicon basis $l$ |  $\beta$ [$1/a_{0}^{2}$]
---|---
1 | 13.8028
1 | 59.9261
1 | 4.34446
1 | 0.267360
1 | 0.0765250
2 | 0.45


* [56] R. Chisnell, J. S. Helton, D. E. Freedman, D. K. Singh, R. I. Bewley, D. G. Nocera, and Y. S. Lee, “Topological magnon bands in a kagome lattice ferromagnet,” Phys. Rev. Lett., vol. 115, p. 147201, Sep 2015.
  * [57] Y. Su and X. R. Wang, “Chiral anomaly of Weyl magnons in stacked honeycomb ferromagnets,” Phys. Rev. B, vol. 96, p. 104437, Sep 2017.
  * [58] R. Seshadri and D. Sen, “Topological magnons in a kagome-lattice spin system with $xxz$ and Dzyaloshinskii-Moriya interactions,” Phys. Rev. B, vol. 97, p. 134411, Apr 2018.
  * [59] A. Mook, J. Henk, and I. Mertig, “Magnon waveguide with nanoscale confinement constructed from topological magnon insulators,” Phys. Rev. B, vol. 91, p. 174409, May 2015.
  * [60] M. Hirschberger, R. Chisnell, Y. S. Lee, and N. P. Ong, “Thermal Hall effect of spin excitations in a kagome magnet,” Phys. Rev. Lett., vol. 115, p. 106603, Sep 2015.
  * [61] A. Mook, J. Henk, and I. Mertig, “Spin dynamics simulations of topological magnon insulators: From transverse current correlation functions to the family of magnon Hall effects,” Phys. Rev. B, vol. 94, p. 174444, Nov 2016.
  * [62] Y. Su, X. S. Wang, and X. R. Wang, “Magnonic Weyl semimetal and chiral anomaly in pyrochlore ferromagnets,” Phys. Rev. B, vol. 95, p. 224403, Jun 2017.
  * [63] X. Cao, K. Chen, and D. He, “Magnon Hall effect on the lieb lattice,” J. Phys.: Condens. Matter, vol. 27, no. 16, p. 166003, 2015.
  * [64] G. Wirth, M. Ölschläger, and A. Hemmerich, “Evidence for orbital superfluidity in the p-band of a bipartite optical square lattice,” Nat. Phys., vol. 7, pp. 147–153, Dec. 2010.
  * [65] P. S. Kumar, I. F. Herbut, and R. Ganesh, “Dirac hamiltonians for bosonic spectra,” Phys. Rev. Res., vol. 2, p. 033035, Jul 2020.
  * [66] W. Li, D. N. Sheng, C. S. Ting, and Y. Chen, “Fractional quantum spin Hall effect in flat-band checkerboard lattice model,” Phys. Rev. B, vol. 90, p. 081102, Aug 2014.
  * [67] S. M. Rezende, R. L. Rodríguez-Suárez, and A. Azevedo, “Diffusive magnonic spin transport in antiferromagnetic insulators,” Phys. Rev. B, vol. 93, p. 054412, Feb 2016.
  * [68] R. Cheng, S. Okamoto, and D. Xiao, “Spin Nernst effect of magnons in collinear antiferromagnets,” Phys. Rev. Lett., vol. 117, no. 21, p. 217202, 2016.
  * [69] V. A. Zyuzin and A. A. Kovalev, “Magnon spin Nernst effect in antiferromagnets,” Phys. Rev. Lett., vol. 117, p. 217203, Nov 2016.
  * [70] H. Zhang and R. Cheng, “A perspective on magnon spin Nernst effect in antiferromagnets,” Appl. Phys. Lett., vol. 120, p. 090502, Feb. 2022.
  * [71] Y. Lu, X. Guo, V. Koval, and C. Jia, “Topological thermal Hall effect driven by spin-chirality fluctuations in frustrated antiferromagnets,” Phys. Rev. B, vol. 99, p. 054409, Feb 2019.
  * [72] P. Laurell and G. A. Fiete, “Magnon thermal Hall effect in kagome antiferromagnets with Dzyaloshinskii-Moriya interactions,” Phys. Rev. B, vol. 98, p. 094419, Sep 2018.
  * [73] H. Doki, M. Akazawa, H.-Y. Lee, J. H. Han, K. Sugii, M. Shimozawa, N. Kawashima, M. Oda, H. Yoshida, and M. Yamashita, “Spin thermal hall conductivity of a kagome antiferromagnet,” Phys. Rev. Let., vol. 121, Aug. 2018.
  * [74] K. Li, “Thermal Hall conductivity with sign change in the Heisenberg–Kitaev kagome magnet,” Commun. Theor. Phys., vol. 75, p. 015702, Dec. 2022.
  * [75] A. Mook, J. Henk, and I. Mertig, “Thermal Hall effect in noncollinear coplanar insulating antiferromagnets,” Phys. Rev. B, vol. 99, p. 014427, Jan 2019.
  * [76] M. Hirschberger, J. W. Krizan, R. J. Cava, and N. P. Ong, “Large thermal Hall conductivity of neutral spin excitations in a frustrated quantum magnet,” Science, vol. 348, pp. 106–109, Apr. 2015.
  * [77] V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, “Antiferromagnetic spintronics,” Rev. Mod. Phys., vol. 90, p. 015005, Feb 2018.
  * [78] T. Jungwirth, J. Sinova, A. Manchon, X. Marti, J. Wunderlich, and C. Felser, “The multiple directions of antiferromagnetic spintronics,” Nat. Phys., vol. 14, pp. 200–203, Mar. 2018.
  * [79] R. A. Duine, K.-J. Lee, S. S. P. Parkin, and M. D. Stiles, “Synthetic antiferromagnetic spintronics,” Nat. Phys., vol. 14, pp. 217–219, Mar. 2018.
  * [80] O. Gomonay, V. Baltz, A. Brataas, and Y. Tserkovnyak, “Antiferromagnetic spin textures and dynamics,” Nat. Phys., vol. 14, pp. 213–216, Mar. 2018.
  * [81] J. Zelezny, P. Wadley, K. Olejnik, A. Hoffmann, and H. Ohno, “Spin transport and spin torque in antiferromagnetic devices,” Nat. Phys., vol. 14, pp. 220–228, Mar. 2018.
  * [82] P. Nemec, M. Fiebig, T. Kampfrath, and A. V. Kimel, “Antiferromagnetic opto-spintronics,” Nat. Phys., vol. 14, pp. 229–241, Mar. 2018.
  * [83] L. Smejkal, Y. Mokrousov, B. Yan, and A. H. MacDonald, “Topological antiferromagnetic spintronics,” Nat. Phys., vol. 14, pp. 242–251, Mar. 2018.
  * [84] M. B. Jungfleisch, W. Zhang, and A. Hoffmann, “Perspectives of antiferromagnetic spintronics,” Phys. Lett. A, vol. 382, pp. 865–871, Apr. 2018.
  * [85] H. Chen, Q. Niu, and A. H. MacDonald, “Anomalous Hall effect arising from noncollinear antiferromagnetism,” Phys. Rev. Lett., vol. 112, p. 017205, Jan 2014.
  * [86] J. Kübler and C. Felser, “Non-collinear antiferromagnets and the anomalous Hall effect,” Europhys. Lett., vol. 108, p. 67001, Dec. 2014.
  * [87] A. S. T. Pires, “Topological magnons in the antiferromagnetic checkerboard lattice,” Physica E, vol. 118, p. 113899, Apr. 2020.
  * [88] A. S. T. Pires, “Dynamical spin conductivity of the Heisenberg antiferromagnetic checkerboard lattice,” Braz. J. Phys., vol. 51, pp. 429–433, Apr. 2021.
  * [89] P. Le and M. Yarmohammadi, “Magnonic heat transport in the Lieb lattice,” J. Magn. Magn. Mater., vol. 469, pp. 623–628, Jan. 2019.
  * [90] P. G. de Oliveira and A. S. T. Pires, “Transversal transport of magnons in a modified Lieb lattice,” Physica B, vol. 654, p. 414721, Apr. 2023.
  * [91] A. Bhattacharya and B. Pal, “Flat bands and nontrivial topological properties in an extended Lieb lattice,” Phys. Rev. B, vol. 100, Dec. 2019.
  * [92] Y. Aharonov and A. Casher, “Topological quantum effects for neutral particles,” Phys. Rev. Lett., vol. 53, pp. 319–321, Jul 1984.
  * [93] G. D. Mahan, Many-particle physics.  Physics of Solids and Liquids, New York, NY: Springer, 2 ed., Dec. 2012\.
  * [94] D. Malz, J. Knolle, and A. Nunnenkamp, “Topological magnon amplification,” Nat. Commun., vol. 10, Sept. 2019.
  * [95] D. Bhowmick, H. Sun, B. Yang, and P. Sengupta, “Tuning bulk topological magnon properties with light-induced magnons,” Phys. Rev. B, vol. 107, June 2023\.
  * [96] R. Kubo, M. Toda, and N. Hashitsume, Statistical physics II.  Springer Series in Solid-State Sciences, Berlin, Germany: Springer, 1985 ed., June 1985.
  * [97] A. S. T. Pires, Theoretical tools for spin models in magnetic systems.  2053-2563, IOP Publishing, 2021.
  * [98] J. Colpa, “Diagonalization of the quadratic boson hamiltonian,” Physica A, vol. 93, pp. 327–353, Sept. 1978.
  * [99] H. Katsura, N. Nagaosa, and A. V. Balatsky, “Spin current and magnetoelectric effect in noncollinear magnets,” Phys. Rev. Lett., vol. 95, p. 057205, Jul 2005.
  * [100] A. Collins, J. McEvoy, D. Robinson, C. J. Hamer, and Z. Weihong, “Quantum spin model with frustration on the Union Jack lattice,” Phys. Rev. B, vol. 73, p. 024407, Jan 2006.
  * [101] R. F. Bishop, P. H. Y. Li, D. J. J. Farnell, and C. E. Campbell, “Magnetic order on a frustrated spin-1/2 Heisenberg antiferromagnet on the Union Jack lattice,” Phys. Rev. B, vol. 82, p. 024416, Jul 2010.
  * [102] W. Zheng, J. Oitmaa, and C. J. Hamer, “Phase diagram of the frustrated Heisenberg antiferromagnet on the Union Jack lattice,” Phys. Rev. B, vol. 75, p. 184418, May 2007.
  * [103] R. F. Bishop and P. H. Y. Li, “A frustrated quantum spin-s model on the Union Jack lattice with spins s $>$ 1/2,” Eur. Phys. J. B, vol. 81, pp. 37–48, Apr. 2011.
  * [104] H. Kondo and Y. Akagi, “Nonlinear magnon spin Nernst effect in antiferromagnets and strain-tunable pure spin current,” Phys. Rev. Res., vol. 4, p. 013186, Mar 2022.
  * [105] J.-M. Hou and W. Chen, “Hidden symmetry and protection of Dirac points on the honeycomb lattice,” Sci. Rep., vol. 5, p. 17571, Dec. 2015.
  * [106] R. Hidalgo-Sacoto, R. I. Gonzalez, E. E. Vogel, S. Allende, J. D. Mella, C. Cardenas, R. E. Troncoso, and F. Munoz, “Magnon valley Hall effect in ${C}r{I}_{3}$-based van der Waals heterostructures,” Phys. Rev. B, vol. 101, p. 205425, May 2020.
  * [107] Z. Wang, S. Zhou, W. Chen, and F.-C. Zhang, “$t-j$ model on the effective brick-wall lattice for the recently discovered high-temperature superconductor ${B}a_{2}{C}u{O}_{3+\delta}$,” Phys. Rev. B, vol. 101, p. 180509, May 2020.
  * [108] H. Kondo, Y. Akagi, and H. Katsura, “$\mathbb{Z}_{2}$ topological invariant for magnon spin Hall systems,” Phys. Rev. B, vol. 99, p. 041110(R), Jan. 2019\.
  * [109] H. Wioland, F. G. Woodhouse, J. Dunkel, and R. E. Goldstein, “Ferromagnetic and antiferromagnetic order in bacterial vortex lattices,” Nat. Phys., vol. 12, pp. 341–345, Jan. 2016.
  * [110] M. Olschlager, G. Wirth, T. Kock, and A. Hemmerich, “Topologically induced avoided band crossing in an optical checkerboard lattice,” Phys. Rev. Lett., vol. 108, p. 075302, Feb 2012.
  * [111] F. Azizi and H. Rezania, “Spin transport properties of anisotropic Heisenberg antiferromagnet on honeycomb lattice in the presence of magnetic field,” Eur. Phys. J. B, vol. 93, Feb. 2020.
  * [112] F. Azizi and H. Rezania, “Dynamical and static spin structure factors of Heisenberg antiferromagnet on honeycomb lattice in the presence of Dzyaloshinskii-Moriya interaction (preprint),” 2021.