<figure><img src="image_1.png" /></figure>

HAL
open science

# GW density matrix in solids as an estimate of self-consistent GW total energies

Hassan Denawi, Fabien Bruneval, Marc Torrent, Mauricio Rodriguez-Mayorga

► To cite this version:

Hassan Denawi, Fabien Bruneval, Marc Torrent, Mauricio Rodriguez-Mayorga. GW density matrix in solids as an estimate of self-consistent GW total energies. Physical Review B, 2023, 108, pp.125107. 10.1103/PhysRevB.108.125107. cea-04603990

HAL Id: cea-04603990

https://cea.hal.science/cea-04603990

Submitted on 6 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# GW density matrix in solids as an estimate of self-consistent GW total energies

Hassan Denawi, Fabien Bruneval, Marc Torrent, Mauricio Rodriguez-Mayorga

► To cite this version:

Hassan Denawi, Fabien Bruneval, Marc Torrent, Mauricio Rodriguez-Mayorga. GW density matrix in solids as an estimate of self-consistent GW total energies. Physical Review B, 2023, 108, pp.125107. 10.1103/PhysRevB.108.125107. cea-04603990

HAL Id: cea-04603990

https://cea.hal.science/cea-04603990v1

Submitted on 6 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# GW density matrix for estimation of self-consistent GW total energies in solids

Adam Hassan Denawi and Fabien Bruneval

Université Paris-Saclay, CEA, Service de Recherche en Corrosion et Comportement des Matériaux, SRMP, 91191 Gif-sur-Yvette, France

Marc Torrent

CEA, DAM, DIF, F-91297 Arpajon, France

and Université Paris-Saclay, CEA, Laboratoire Matière en Conditions Extrêmes, 91680 Bruyères-le-Châtel, France

Mauricio Rodríguez-Mayorga

Université Paris-Saclay, CEA, Service de Recherche en Corrosion et Comportement des Matériaux, SRMP, 91191 Gif-sur-Yvette, France

and Department of Physical Chemistry, University of Alicante, E-03080 Alicante, Spain

<figure><img src="image_3.png" /></figure>

(Received 21 April 2023; revised 26 June 2023; accepted 24 August 2023; published 5 September 2023)

The GW approximation is a well-established method for calculating ionization potentials and electron affinities in solids and molecules. For numerous years, obtaining self-consistent GW total energies in solids has been a challenging objective that is not accomplished yet. However, it was shown recently that the linearized GW density matrix permits a reliable prediction of the self-consistent GW total energy for molecules [F. Bruneval, M. Rodríguez-Mayorga, P. Rinke, and M. Dvorak, J. Chem. Theory Comput. 17, 2126 (2021)] for which self-consistent GW energies are available. Here we implement, test, and benchmark the linearized GW density matrix for several solids. We focus on the total energy, lattice constant, and bulk modulus obtained from the GW density matrix and compare our findings to more traditional results obtained within the random-phase approximation (RPA). We conclude on the improved stability of the total energy obtained from the linearized GW density matrix with respect to the mean-field starting point. We bring compelling clues that the RPA and the GW density matrix total energies are certainly close to the self-consistent GW total energy in solids if we use hybrid functionals with enriched exchange as a starting point.

DOI: 10.1103/PhysRevB.108.125107

## I. INTRODUCTION

While a few self-consistent GW (scGW) calculations for the band gaps of real solids are available [1–3], scGW total energies are still not available today, certainly because of their high computational cost. However, there exist hints that scGW could be accurate: First, the results on the homogeneous electron gas are extremely good [4–6]; second, the random-phase approximation (RPA) which derives from the same family has been shown to yield total energies capable of describing the tenuous van der Waals interactions [7–13]. Unfortunately, scGW calculations are very involved in real solids. That is why it would be highly desirable to obtain scGW quality energies without actually performing the cumbersome self-consistency.

Pursuing the quest for a non-self-consistent approximation to scGW, a series of studies have been published in the early 2000s [14–19]. More recently, some of us proposed an alternative non-self-consistent total energy expression based on the GW density matrix, labeled $\gamma^{GW}$ [20,21]. Benchmarks on small molecules for which scGW calculations are possible [22,23] confirmed the remarkable properties of the $\gamma^{GW}$ total energies: Although it is evaluated non-self-consistently using a generalized Kohn-Sham (gKS) input, the resulting total energy remains quite insensitive to the gKS choice and approximates very well the reference scGW total energies.

In this context, this paper focuses on the evaluation of total energies in solids. We port to the solid systems the $\gamma^{GW}$ total energy with the sensible prospect that it will remain a good approximation to the scGW total energy. In doing so, we obtain the correlated density matrix $\gamma^{GW}$ as a physically meaningful intermediate object with unique properties due to correlation.

When considering solid systems several technical questions have to be addressed. First, the closed formulas obtained for finite systems [20,24] have to be adapted for numerical efficiency. Second, the pseudopotentials [25] that are customary in the plane-wave basis codes are typically designed to be used in conjunction with standard semilocal approximations to density-functional theory (DFT). Therefore, it is necessary to investigate which type of pseudopotential is suitable for obtaining consistently GW-type total energies.

With this, we will study the performance of the $\gamma^{GW}$ total energy for solids. We will compare it to the popular RPA total energy which may be derived as the GW approximation of the Klein functional [14,26,27]. As all these calculations are performed as a one-shot procedure, memory about the mean-field starting point is present. We will particularly investigate this issue by varying the content of exact exchange $\alpha$ in the hybrid functional PBEh($\alpha$), zero being the standard Perdew-Burke-Ernzerhof (PBE) functional [28] and 0.25 yielding regular PBE0 [29]. It is therefore always necessary to specify the
# GW density matrix for estimation of self-consistent GW total energies in solids

Adam Hassan Denawi and Fabien Bruneval

Université Paris-Saclay, CEA, Service de Recherche en Corrosion et Comportement des Matériaux, SRMP, 91191 Gif-sur-Yvette, France

Marc Torrent

CEA, DAM, DIF, F-91297 Arpajon, France

and Université Paris-Saclay, CEA, Laboratoire Matière en Conditions Extrêmes, 91680 Bruyères-le-Châtel, France

Mauricio Rodríguez-Mayorga

Université Paris-Saclay, CEA, Service de Recherche en Corrosion et Comportement des Matériaux, SRMP, 91191 Gif-sur-Yvette, France

and Department of Physical Chemistry, University of Alicante, E-03080 Alicante, Spain

<figure><img src="image_3.png" /></figure>

(Received 21 April 2023; revised 26 June 2023; accepted 24 August 2023; published 5 September 2023)

The GW approximation is a well-established method for calculating ionization potentials and electron affinities in solids and molecules. For numerous years, obtaining self-consistent GW total energies in solids has been a challenging objective that is not accomplished yet. However, it was shown recently that the linearized GW density matrix permits a reliable prediction of the self-consistent GW total energy for molecules [F. Bruneval, M. Rodríguez-Mayorga, P. Rinke, and M. Dvorak, J. Chem. Theory Comput. 17, 2126 (2021)] for which self-consistent GW energies are available. Here we implement, test, and benchmark the linearized GW density matrix for several solids. We focus on the total energy, lattice constant, and bulk modulus obtained from the GW density matrix and compare our findings to more traditional results obtained within the random-phase approximation (RPA). We conclude on the improved stability of the total energy obtained from the linearized GW density matrix with respect to the mean-field starting point. We bring compelling clues that the RPA and the GW density matrix total energies are certainly close to the self-consistent GW total energy in solids if we use hybrid functionals with enriched exchange as a starting point.

DOI: 10.1103/PhysRevB.108.125107

## I. INTRODUCTION

While a few self-consistent GW (scGW) calculations for the band gaps of real solids are available [1–3], scGW total energies are still not available today, certainly because of their high computational cost. However, there exist hints that scGW could be accurate: First, the results on the homogeneous electron gas are extremely good [4–6]; second, the random-phase approximation (RPA) which derives from the same family has been shown to yield total energies capable of describing the tenuous van der Waals interactions [7–13]. Unfortunately, scGW calculations are very involved in real solids. That is why it would be highly desirable to obtain scGW quality energies without actually performing the cumbersome self-consistency.

Pursuing the quest for a non-self-consistent approximation to scGW, a series of studies have been published in the early 2000s [14–19]. More recently, some of us proposed an alternative non-self-consistent total energy expression based on the GW density matrix, labeled $\gamma^{GW}$ [20,21]. Benchmarks on small molecules for which scGW calculations are possible [22,23] confirmed the remarkable properties of the $\gamma^{GW}$ total energies: Although it is evaluated non-self-consistently using a generalized Kohn-Sham (gKS) input, the resulting total energy remains quite insensitive to the gKS choice and approximates very well the reference scGW total energies.

In this context, this paper focuses on the evaluation of total energies in solids. We port to the solid systems the $\gamma^{GW}$ total energy with the sensible prospect that it will remain a good approximation to the scGW total energy. In doing so, we obtain the correlated density matrix $\gamma^{GW}$ as a physically meaningful intermediate object with unique properties due to correlation.

When considering solid systems several technical questions have to be addressed. First, the closed formulas obtained for finite systems [20,24] have to be adapted for numerical efficiency. Second, the pseudopotentials [25] that are customary in the plane-wave basis codes are typically designed to be used in conjunction with standard semilocal approximations to density-functional theory (DFT). Therefore, it is necessary to investigate which type of pseudopotential is suitable for obtaining consistently GW-type total energies.

With this, we will study the performance of the $\gamma^{GW}$ total energy for solids. We will compare it to the popular RPA total energy which may be derived as the GW approximation of the Klein functional [14,26,27]. As all these calculations are performed as a one-shot procedure, memory about the mean-field starting point is present. We will particularly investigate this issue by varying the content of exact exchange $\alpha$ in the hybrid functional PBEh($\alpha$), zero being the standard Perdew-Burke-Ernzerhof (PBE) functional [28] and 0.25 yielding regular PBE0 [29]. It is therefore always necessary to specify the
exact procedure used to obtain a one-shot energy. We do so here using the @ notation (e.g., RPA@PBE stands for the RPA total energy evaluated with self-consistent PBE inputs).

The paper is organized as follows: In Sec. II, we recapitulate the theoretical foundation for the GW density matrix and derive the working equations. In Sec. III, we detail the technical aspects of the implementation in a plane-wave code and assess the pseudopotential choice. Section IV shows some of the unique properties of the GW density matrix, exemplified with bulk silicon. In Sec. V, we benchmark the total energies obtained with the different approximations with a test set of seven standard covalent semiconductors and one layered material. Finally, Sec. VI concludes our paper.

## II. THEORY AND WORKING FORMULAS

### A. Green's-function-derived density matrix

In the many-body perturbation theory, the central quantity, namely the one-particle Green's function $G(\mathbf{r}, \mathbf{r}', \omega)$, contains a great deal of information. In particular, by virtue of the Galistkii-Midgal formula [30], it is sufficient to calculate the total energy of an electronic system.

Also it straightforwardly yields the density matrix $\gamma(\mathbf{r}, \mathbf{r}')$:

$$ \gamma(\mathbf{r}, \mathbf{r}') = -\frac{i}{2\pi} \int_{-\infty}^{+\infty} d\omega e^{i\eta\omega} G(\mathbf{r}, \mathbf{r}', \omega), \quad (1) $$

where $\eta$ is a vanishing positive real number that enforces the sensitivity to the occupied manifold of the time-ordered Green's function $G$. Hence, the electronic density can be obtained as the diagonal: $\rho(\mathbf{r}) = \gamma(\mathbf{r}, \mathbf{r}')$.

Therefore, the Green's-function methods, such as the GW approximation, can give access to an approximate density matrix.

### B. Linearized Dyson equation

In the many-body perturbation theory, the overall strategy is to connect the exact Green's function $G$ to a known Green's function $G_0$. The connection between the two is ensured by the complicated self-energy $\Sigma_{xc}$ that is in charge of all the correlation effects.

The expression of $G_0$ that is derived from a mean-field approach [Kohn-Sham (KS), Hartree-Fock (HF), etc.], is simple:

$$ G_0(\mathbf{r}, \mathbf{r}', \omega) = 2 \sum_{\mathbf{k}\mathbf{i}} \frac{\varphi_{\mathbf{k}\mathbf{i}}(\mathbf{r}) \varphi_{\mathbf{k}\mathbf{i}}^*(\mathbf{r}')}{\omega - \epsilon_{\mathbf{k}\mathbf{i}} \pm i\eta}, \quad (2) $$

where $\varphi_{\mathbf{k}\mathbf{i}}(\mathbf{r})$ and $\epsilon_{\mathbf{k}\mathbf{i}}$ are the mean-field wave functions and eigenvalue for state $i$ at the $\mathbf{k}$ point and the small positive $\eta$ ensures the correct location of the poles for a time-ordered function (above the real axis for occupied states $\epsilon_{\mathbf{k}\mathbf{i}} < \mu$ and below the real axis for empty states $\epsilon_{\mathbf{k}\mathbf{i}} > \mu$, $\mu$ being the Fermi level). Spin-restricted calculations are assumed here and the factor 2 accounts for it.

Then the connection from $G_0$ to $G$ is made with the so-called Dyson equation:

$$ G = G_0 + G_0(\Sigma_{xc} - V_{xc})G, \quad (3) $$

where $V_{xc}$ is the exchange-correlation operator (possibly including nonlocal exchange) and the space and frequency indices have been omitted for conciseness.

The self-energy $\Sigma_{xc}$ is itself a functional of the exact $G$. When approximating $\Sigma_{xc}$ and $G$, only a self-consistent solution ensures the conservation of the electron count [31]. In particular, for the GW approximation of $\Sigma_{xc}$ that is most often evaluated with a one-shot procedure, the violation of electron conservation is well documented [18,21,32,33].

In a previous study of ours [21], it was demonstrated analytically and verified numerically that linearizing the Dyson equation completely cures the problem of electron count conservation in the GW approximation. The linearized Dyson equation (LDE) reads

$$ G = G_0 + G_0(\Sigma_{xc} - V_{xc})G_0, \quad (4) $$

where the last $G$ in Eq. (3) had been simply replaced by $G_0$. The LDE is customary in the context of the Sham-Schüter equation [34].

This electron-conserving equation is then applied with the GW approximation to $\Sigma_{xc}$.

### C. GW self-energy-based density matrix

The GW approximation [35] is simply sketched here, since it has been the subject of numerous detailed reviews [36–38].

The screened Coulomb interaction $W$ is defined with a Dyson-like equation:

$$ W = v + v\chi_0 W, \quad (5) $$

where $\chi_0 = -2iGG$ is the non-interacting polarizability and $v$ is the usual bare Coulomb interaction.

The GW self-energy then reads

$$ \Sigma_{xc} = iGW. \quad (6) $$

It is convenient to decompose the self-energy into pure exchange and correlation. The exchange part $\Sigma_x$ is static, whereas the correlation part carries the frequency dependence $\Sigma_c(\omega)$. These quantities are routinely obtained with a one-shot procedure in standard periodic codes [39,40]

Some general properties of the density matrix are detailed in Appendix A. For instance, it is demonstrated that the density matrix can be fully characterized with a single $\mathbf{k}$-point index within the first Brillouin zone, even though it is a function of two spatial indices.

Now let us focus on the GW density matrix. It is handy to project into the mean-field orbitals $|\mathbf{k}\mathbf{i}\rangle$, which form a valid orthogonal basis:

$$ \begin{aligned} \gamma_{\mathbf{k}\mathbf{i}\mathbf{j}} &= \langle \mathbf{k}\mathbf{i} | \gamma | \mathbf{k}\mathbf{j} \rangle \\ &= \int d\mathbf{r} d\mathbf{r}' \varphi_{\mathbf{k}\mathbf{i}}^*(\mathbf{r}) \gamma(\mathbf{r}, \mathbf{r}') \varphi_{\mathbf{k}\mathbf{j}}(\mathbf{r}'). \end{aligned} \quad (7) $$

The first term on the right-hand side of Eq. (4) is $G_0$. Let us insert it in Eq. (1) and project on the orbital basis to obtain the spin-summed density matrix elements:

$$ \gamma_{\mathbf{k}\mathbf{i}\mathbf{j}}^{\text{gKS}} = 2\delta_{\mathbf{i}\mathbf{j}}\theta(\mu - \epsilon_{\mathbf{k}\mathbf{i}}). \quad (8) $$

This expression has been obtained by closing the contour in the upper part of the complex plane so that only the poles located above the real axis have survived.
exact procedure used to obtain a one-shot energy. We do so here using the @ notation (e.g., RPA@PBE stands for the RPA total energy evaluated with self-consistent PBE inputs).

The paper is organized as follows: In Sec. II, we recapitulate the theoretical foundation for the GW density matrix and derive the working equations. In Sec. III, we detail the technical aspects of the implementation in a plane-wave code and assess the pseudopotential choice. Section IV shows some of the unique properties of the GW density matrix, exemplified with bulk silicon. In Sec. V, we benchmark the total energies obtained with the different approximations with a test set of seven standard covalent semiconductors and one layered material. Finally, Sec. VI concludes our paper.

## II. THEORY AND WORKING FORMULAS

### A. Green's-function-derived density matrix

In the many-body perturbation theory, the central quantity, namely the one-particle Green's function $G(\mathbf{r}, \mathbf{r}', \omega)$, contains a great deal of information. In particular, by virtue of the Galistkii-Midgal formula [30], it is sufficient to calculate the total energy of an electronic system.

Also it straightforwardly yields the density matrix $\gamma(\mathbf{r}, \mathbf{r}')$:

$$ \gamma(\mathbf{r}, \mathbf{r}') = -\frac{i}{2\pi} \int_{-\infty}^{+\infty} d\omega e^{i\eta\omega} G(\mathbf{r}, \mathbf{r}', \omega), \quad (1) $$

where $\eta$ is a vanishing positive real number that enforces the sensitivity to the occupied manifold of the time-ordered Green's function $G$. Hence, the electronic density can be obtained as the diagonal: $\rho(\mathbf{r}) = \gamma(\mathbf{r}, \mathbf{r}')$.

Therefore, the Green's-function methods, such as the GW approximation, can give access to an approximate density matrix.

### B. Linearized Dyson equation

In the many-body perturbation theory, the overall strategy is to connect the exact Green's function $G$ to a known Green's function $G_0$. The connection between the two is ensured by the complicated self-energy $\Sigma_{xc}$ that is in charge of all the correlation effects.

The expression of $G_0$ that is derived from a mean-field approach [Kohn-Sham (KS), Hartree-Fock (HF), etc.], is simple:

$$ G_0(\mathbf{r}, \mathbf{r}', \omega) = 2 \sum_{\mathbf{k}\mathbf{i}} \frac{\varphi_{\mathbf{k}\mathbf{i}}(\mathbf{r}) \varphi_{\mathbf{k}\mathbf{i}}^*(\mathbf{r}')}{\omega - \epsilon_{\mathbf{k}\mathbf{i}} \pm i\eta}, \quad (2) $$

where $\varphi_{\mathbf{k}\mathbf{i}}(\mathbf{r})$ and $\epsilon_{\mathbf{k}\mathbf{i}}$ are the mean-field wave functions and eigenvalue for state $i$ at the $\mathbf{k}$ point and the small positive $\eta$ ensures the correct location of the poles for a time-ordered function (above the real axis for occupied states $\epsilon_{\mathbf{k}\mathbf{i}} < \mu$ and below the real axis for empty states $\epsilon_{\mathbf{k}\mathbf{i}} > \mu$, $\mu$ being the Fermi level). Spin-restricted calculations are assumed here and the factor 2 accounts for it.

Then the connection from $G_0$ to $G$ is made with the so-called Dyson equation:

$$ G = G_0 + G_0(\Sigma_{xc} - V_{xc})G, \quad (3) $$

where $V_{xc}$ is the exchange-correlation operator (possibly including nonlocal exchange) and the space and frequency indices have been omitted for conciseness.

The self-energy $\Sigma_{xc}$ is itself a functional of the exact $G$. When approximating $\Sigma_{xc}$ and $G$, only a self-consistent solution ensures the conservation of the electron count [31]. In particular, for the GW approximation of $\Sigma_{xc}$ that is most often evaluated with a one-shot procedure, the violation of electron conservation is well documented [18,21,32,33].

In a previous study of ours [21], it was demonstrated analytically and verified numerically that linearizing the Dyson equation completely cures the problem of electron count conservation in the GW approximation. The linearized Dyson equation (LDE) reads

$$ G = G_0 + G_0(\Sigma_{xc} - V_{xc})G_0, \quad (4) $$

where the last $G$ in Eq. (3) had been simply replaced by $G_0$. The LDE is customary in the context of the Sham-Schüter equation [34].

This electron-conserving equation is then applied with the GW approximation to $\Sigma_{xc}$.

### C. GW self-energy-based density matrix

The GW approximation [35] is simply sketched here, since it has been the subject of numerous detailed reviews [36–38].

The screened Coulomb interaction $W$ is defined with a Dyson-like equation:

$$ W = v + v\chi_0 W, \quad (5) $$

where $\chi_0 = -2iGG$ is the non-interacting polarizability and $v$ is the usual bare Coulomb interaction.

The GW self-energy then reads

$$ \Sigma_{xc} = iGW. \quad (6) $$

It is convenient to decompose the self-energy into pure exchange and correlation. The exchange part $\Sigma_x$ is static, whereas the correlation part carries the frequency dependence $\Sigma_c(\omega)$. These quantities are routinely obtained with a one-shot procedure in standard periodic codes [39,40]

Some general properties of the density matrix are detailed in Appendix A. For instance, it is demonstrated that the density matrix can be fully characterized with a single $\mathbf{k}$-point index within the first Brillouin zone, even though it is a function of two spatial indices.

Now let us focus on the GW density matrix. It is handy to project into the mean-field orbitals $|\mathbf{k}\mathbf{i}\rangle$, which form a valid orthogonal basis:

$$ \begin{aligned} \gamma_{\mathbf{k}\mathbf{i}\mathbf{j}} &= \langle \mathbf{k}\mathbf{i} | \gamma | \mathbf{k}\mathbf{j} \rangle \\ &= \int d\mathbf{r} d\mathbf{r}' \varphi_{\mathbf{k}\mathbf{i}}^*(\mathbf{r}) \gamma(\mathbf{r}, \mathbf{r}') \varphi_{\mathbf{k}\mathbf{j}}(\mathbf{r}'). \end{aligned} \quad (7) $$

The first term on the right-hand side of Eq. (4) is $G_0$. Let us insert it in Eq. (1) and project on the orbital basis to obtain the spin-summed density matrix elements:

$$ \gamma_{\mathbf{k}\mathbf{i}\mathbf{j}}^{\text{gKS}} = 2\delta_{\mathbf{i}\mathbf{j}}\theta(\mu - \epsilon_{\mathbf{k}\mathbf{i}}). \quad (8) $$

This expression has been obtained by closing the contour in the upper part of the complex plane so that only the poles located above the real axis have survived.
A similar approach can be used for the static terms in the right-hand part of Eq. (4), $G_0(\Sigma_x - V_{xc})G_0$:

$$ \Delta\gamma_{\mathbf{k}ij}^{\text{HF}} = 2\theta(\mu - \epsilon_{\mathbf{k}i})\theta(\epsilon_{\mathbf{k}j} - \mu)\frac{\langle\mathbf{k}i|\Sigma_x - V_{xc}|\mathbf{k}j\rangle}{\epsilon_{\mathbf{k}i} - \epsilon_{\mathbf{k}j}}. \quad (9) $$

We denote it with a “HF” superscript because this contribution to the (spin-summed) linearized density matrix is obtained from a pure exchange self-energy; thus, it vanishes when the HF approximation is employed to obtain the mean-field orbitals (i.e., $\Delta\gamma^{\text{HF}} = 0$ for $\gamma^{\text{gKS}} = \gamma^{\text{HF}}$).

For the last term in Eq. (4), $G_0\Sigma_c G_0$, the self-energy $\Sigma_c$ has a frequency dependence and therefore the calculations cannot be performed analytically in contrast with the two previous terms; it is convenient to perform the integration along the imaginary axis, so as to keep some distance with the poles of $G_0$ and of $\Sigma_c$. Closing the contour, we transform the real-axis integration of Eq. (1) into

$$ \Delta\gamma_{\mathbf{k}ij}^{\text{GW}} = \frac{1}{\pi} \int_{-\infty}^{+\infty} d\omega \frac{\langle\mathbf{k}i|\Sigma_c(\mu + i\omega)|\mathbf{k}j\rangle}{(\mu + i\omega - \epsilon_{\mathbf{k}i})(\mu + i\omega - \epsilon_{\mathbf{k}j})}. \quad (10) $$

The complete spin-summed linearized GW density matrix finally reads

$$ \gamma^{\text{GW}} = \gamma^{\text{gKS}} + \Delta\gamma^{\text{HF}} + \Delta\gamma^{\text{GW}}. \quad (11) $$

Lastly, the corresponding electronic density is $\rho^{\text{GW}}(\mathbf{r}) = \gamma^{\text{GW}}(\mathbf{r}, \mathbf{r})$.

### D. Total energies from the GW density matrix

In previous studies [20,21], we introduced a new total energy functional:

$$ E_{\text{total}}^{\text{GW}} = T[\gamma^{\text{GW}}] + V_{\text{ne}}[\rho^{\text{GW}}] + E_{\text{H}}[\rho^{\text{GW}}] \\ + E_{\text{x}}[\gamma^{\text{GW}}] + E_{\text{c}}[G_0] + V_{\text{nn}}, \quad (12) $$

where the kinetic energy $T$, the electron-nucleus interaction $V_{ne}$, the Hartree energy $E_H$, and the exchange energy $E_x$ are evaluated with the $\gamma^{\text{GW}}$ density matrix. Klimeš et al. [41] also used the GW (or RPA) density matrix to improve subparts of the energy. Just the correlation energy $E_c$ cannot be calculated with $\gamma^{\text{GW}}$ and is pragmatically obtained from the Galitskii-Migdal equation [21]:

$$ E_c[G_0] = \frac{1}{4\pi} \int_{-\infty}^{+\infty} d\omega \text{Tr}\{v\chi_0(i\omega) - v\chi(i\omega)\}, \quad (13) $$

where $\chi = \chi_0 + \chi_0 v \chi$ is the RPA polarizability.

This one-shot energy expression has the desirable property that all the input quantities conserve the number of electrons. Of course, being a one-shot total energy, it keeps a dependence with respect to the starting point. This will be studied in detail in Sec. V.

### E. RPA total energy

For completeness, we report here without derivation the RPA expression for the total energy as we will extensively compare $E_{\text{total}}^{\text{RPA}}$ and $E_{\text{total}}^{\gamma^{\text{GW}}}$ in the following.

The RPA correlation is defined as [42]

$$ \Phi_c[G_0] = \frac{1}{4\pi} \int_{-\infty}^{+\infty} d\omega \text{Tr}\{v\chi_0(i\omega) + \ln[1 - v\chi_0(i\omega)]\}. \quad (14) $$

By construction, $\Phi_c[G_0]$ contains the correlation part of the kinetic energy. The total one-shot energy expression reads

$$ E_{\text{total}}^{\text{RPA}} = T[\gamma^{\text{gKS}}] + V_{\text{ne}}[\rho^{\text{gKS}}] + E_{\text{H}}[\rho^{\text{gKS}}] \\ + E_{\text{x}}[\gamma^{\text{gKS}}] + \Phi_c[G_0] + V_{\text{nn}}. \quad (15) $$

The one-shot RPA total energy is known to have a noticeable starting point dependence [15,21,43,44]. We will quantify this in Sec. V in comparison with $E_{\text{total}}^{\gamma^{\text{GW}}}$.

## III. IMPLEMENTATION AND COMPUTATIONAL DETAILS

The linearized GW density matrix in periodic systems has not been studied before to the best of our knowledge. It should be noted though that the linearized GW density matrix appears as an intermediate quantity in the RPA forces derived by Ramberger et al. [45].

In this section, we provide a detailed description of our implementation of $\gamma^{\text{GW}}$ in the ABINIT code [40]. We also highlight the key technical aspects that are crucial for producing accurate results.

### A. Implementation in a periodic plane-wave approach

ABINIT is a standard plane-wave-based DFT code. The core electrons are frozen and hidden in a pseudopotential. While the Kohn-Sham part of ABINIT is able to use the more accurate and smoother projector augmented-wave atomic datasets [46–48], the extension to GW is very delicate [49,50]. As of today, the GW part of ABINIT is fully validated only for regular norm-conserving pseudopotentials [51,52].

The existing implementation in ABINIT provides us with $\langle\mathbf{k}i|\Sigma_c(\mu + i\omega)|\mathbf{k}j\rangle$ for any value of $\omega$. From this starting point, we have then implemented a Gauss-Legendre quadra-ture to perform the integral in Eq. (10). The symmetry relation

$$ \langle\mathbf{k}i|\Sigma_c(\mu - i\omega)|\mathbf{k}j\rangle = \langle\mathbf{k}j|\Sigma_c(\mu + i\omega)|\mathbf{k}i\rangle^* \quad (16) $$

is employed to limit the integration from zero to $+\infty$. A grid with typically 50 to 120 grid points is sufficient to ensure a very accurate convergence: We monitor the electron count deviation, which is always kept below $10^{-3}$. In the future, grid design could be optimized to minimize the computational burden [33,53].

The static term from Eq. (9) has been implemented as well, for any type of exchange-correlation potential $V_{xc}$, including those based on hybrid functionals. Note that for a Hartree-Fock mean-field starting point, the static term $\Sigma_x - V_{xc}$ vanishes. Furthermore, it is clear from Eq. (9) that the linearized density matrix is limited to systems with a finite band gap, or else diverging denominators would occur.

The matrix representation of $\gamma^{\text{GW}}$ is obtained on the gKS states $|\mathbf{k}i\rangle$ for $i \le N_b$. We then diagonalize it to obtain the
A similar approach can be used for the static terms in the right-hand part of Eq. (4), $G_0(\Sigma_x - V_{xc})G_0$:

$$ \Delta\gamma_{\mathbf{k}ij}^{\mathrm{HF}} = 2\theta(\mu - \epsilon_{\mathbf{k}i})\theta(\epsilon_{\mathbf{k}j} - \mu)\frac{\langle\mathbf{k}i|\Sigma_x - V_{xc}|\mathbf{k}j\rangle}{\epsilon_{\mathbf{k}i} - \epsilon_{\mathbf{k}j}}. \quad (9) $$

We denote it with a “HF” superscript because this contribution to the (spin-summed) linearized density matrix is obtained from a pure exchange self-energy; thus, it vanishes when the HF approximation is employed to obtain the mean-field orbitals (i.e., $\Delta\gamma^{\mathrm{HF}} = 0$ for $\gamma^{\mathrm{gKS}} = \gamma^{\mathrm{HF}}$).

For the last term in Eq. (4), $G_0\Sigma_c G_0$, the self-energy $\Sigma_c$ has a frequency dependence and therefore the calculations cannot be performed analytically in contrast with the two previous terms; it is convenient to perform the integration along the imaginary axis, so as to keep some distance with the poles of $G_0$ and of $\Sigma_c$. Closing the contour, we transform the real-axis integration of Eq. (1) into

$$ \Delta\gamma_{\mathbf{k}ij}^{\mathrm{GW}} = \frac{1}{\pi} \int_{-\infty}^{+\infty} d\omega \frac{\langle \mathbf{k}i | \Sigma_c (\mu + i\omega) | \mathbf{k}j \rangle}{(\mu + i\omega - \epsilon_{\mathbf{k}i})(\mu + i\omega - \epsilon_{\mathbf{k}j})}. \quad (10) $$

The complete spin-summed linearized GW density matrix finally reads

$$ \gamma^{\mathrm{GW}} = \gamma^{\mathrm{gKS}} + \Delta\gamma^{\mathrm{HF}} + \Delta\gamma^{\mathrm{GW}}. \quad (11) $$

Lastly, the corresponding electronic density is $\rho^{\mathrm{GW}}(\mathbf{r}) = \gamma^{\mathrm{GW}}(\mathbf{r}, \mathbf{r})$.

### D. Total energies from the GW density matrix

In previous studies [20,21], we introduced a new total energy functional:

$$ E_{\mathrm{total}}^{\gamma^{\mathrm{GW}}} = T[\gamma^{\mathrm{GW}}] + V_{\mathrm{ne}}[\rho^{\mathrm{GW}}] + E_{\mathrm{H}}[\rho^{\mathrm{GW}}] \\ + E_{\mathrm{x}}[\gamma^{\mathrm{GW}}] + E_{\mathrm{c}}[G_0] + V_{\mathrm{nn}}, \quad (12) $$

where the kinetic energy $T$, the electron-nucleus interaction $V_{\mathrm{ne}}$, the Hartree energy $E_{\mathrm{H}}$, and the exchange energy $E_{\mathrm{x}}$ are evaluated with the $\gamma^{\mathrm{GW}}$ density matrix. Klimeš et al. [41] also used the GW (or RPA) density matrix to improve subparts of the energy. Just the correlation energy $E_{\mathrm{c}}$ cannot be calculated with $\gamma^{\mathrm{GW}}$ and is pragmatically obtained from the Galitskii-Migdal equation [21]:

$$ E_c[G_0] = \frac{1}{4\pi} \int_{-\infty}^{+\infty} d\omega \mathrm{Tr}\{v\chi_0(i\omega) - v\chi(i\omega)\}, \quad (13) $$

where $\chi = \chi_0 + \chi_0 v \chi$ is the RPA polarizability.

This one-shot energy expression has the desirable property that all the input quantities conserve the number of electrons. Of course, being a one-shot total energy, it keeps a dependence with respect to the starting point. This will be studied in detail in Sec. V.

### E. RPA total energy

For completeness, we report here without derivation the RPA expression for the total energy as we will extensively compare $E_{\mathrm{total}}^{\mathrm{RPA}}$ and $E_{\mathrm{total}}^{\gamma^{\mathrm{GW}}}$ in the following.

The RPA correlation is defined as [42]

$$ \Phi_c[G_0] = \frac{1}{4\pi} \int_{-\infty}^{+\infty} d\omega \mathrm{Tr}\{v\chi_0(i\omega) + \ln[1 - v\chi_0(i\omega)]\}. \quad (14) $$

By construction, $\Phi_c[G_0]$ contains the correlation part of the kinetic energy. The total one-shot energy expression reads

$$ E_{\mathrm{total}}^{\mathrm{RPA}} = T[\gamma^{\mathrm{gKS}}] + V_{\mathrm{ne}}[\rho^{\mathrm{gKS}}] + E_{\mathrm{H}}[\rho^{\mathrm{gKS}}] \\ + E_{\mathrm{x}}[\gamma^{\mathrm{gKS}}] + \Phi_c[G_0] + V_{\mathrm{nn}}. \quad (15) $$

The one-shot RPA total energy is known to have a noticeable starting point dependence [15,21,43,44]. We will quantify this in Sec. V in comparison with $E_{\mathrm{total}}^{\gamma^{\mathrm{GW}}}$.

## III. IMPLEMENTATION AND COMPUTATIONAL DETAILS

The linearized GW density matrix in periodic systems has not been studied before to the best of our knowledge. It should be noted though that the linearized GW density matrix appears as an intermediate quantity in the RPA forces derived by Ramberger et al. [45].

In this section, we provide a detailed description of our implementation of $\gamma^{\mathrm{GW}}$ in the ABINIT code [40]. We also highlight the key technical aspects that are crucial for producing accurate results.

### A. Implementation in a periodic plane-wave approach

ABINIT is a standard plane-wave-based DFT code. The core electrons are frozen and hidden in a pseudopotential. While the Kohn-Sham part of ABINIT is able to use the more accurate and smoother projector augmented-wave atomic datasets [46–48], the extension to GW is very delicate [49,50]. As of today, the GW part of ABINIT is fully validated only for regular norm-conserving pseudopotentials [51,52].

The existing implementation in ABINIT provides us with $\langle \mathbf{k}i | \Sigma_c (\mu + i\omega) | \mathbf{k}j \rangle$ for any value of $\omega$. From this starting point, we have then implemented a Gauss-Legendre quadra-ture to perform the integral in Eq. (10). The symmetry relation

$$ \langle \mathbf{k}i | \Sigma_c (\mu - i\omega) | \mathbf{k}j \rangle = \langle \mathbf{k}j | \Sigma_c (\mu + i\omega) | \mathbf{k}i \rangle^* \quad (16) $$

is employed to limit the integration from zero to $+\infty$. A grid with typically 50 to 120 grid points is sufficient to ensure a very accurate convergence: We monitor the electron count deviation, which is always kept below $10^{-3}$. In the future, grid design could be optimized to minimize the computational burden [33,53].

The static term from Eq. (9) has been implemented as well, for any type of exchange-correlation potential $V_{xc}$, including those based on hybrid functionals. Note that for a Hartree-Fock mean-field starting point, the static term $\Sigma_x - V_{xc}$ vanishes. Furthermore, it is clear from Eq. (9) that the linearized density matrix is limited to systems with a finite band gap, or else diverging denominators would occur.

The matrix representation of $\gamma^{\mathrm{GW}}$ is obtained on the gKS states $|\mathbf{k}i\rangle$ for $i \le N_b$. We then diagonalize it to obtain the
<table><caption>TABLE I. Energy cutoff for the wave functions ($E_{\text{cut}}^{\varphi}$, in hartree), energy cutoff for the plane-wave expansion of the screened Coulomb interaction $W$ ($E_{\text{cut}}^{\text{W}}$, in hartree), number of bands ($N_b$) in $G_0$ and in Eqs. (9) and (10), number of pure imaginary frequencies ($N_{\omega_W}$) used for $W(i\omega)$, and number of frequency points ($N_{\omega_{\Sigma}}$), used for the quadrature along the imaginary axis in Eq. (10).</caption><thead><tr><th></th><th>Si</th><th>C</th><th>SiC</th><th>zb-BN</th><th>AlP</th><th>AlAs</th><th>Ge</th><th>h-BN</th></tr></thead><tbody><tr><td>$E_{\text{cut}}^{\varphi}$</td><td>12</td><td>25</td><td>50</td><td>40</td><td>35</td><td>35</td><td>30</td><td>55</td></tr><tr><td>$E_{\text{cut}}^{\text{W}}$</td><td>8</td><td>15</td><td>15</td><td>15</td><td>15</td><td>10</td><td>10</td><td>25</td></tr><tr><td>$N_b$</td><td>150</td><td>175</td><td>175</td><td>175</td><td>175</td><td>175</td><td>175</td><td>1400</td></tr><tr><td>$N_{\omega_W}$</td><td>120</td><td>120</td><td>120</td><td>120</td><td>120</td><td>120</td><td>120</td><td>40</td></tr><tr><td>$N_{\omega_{\Sigma}}$</td><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td><td></td></tr></tbody></table>

natural orbitals in the gKS basis:

$$ \sum_{j=1}^{N_b} \gamma_{\mathbf{k}i j} U_{\mathbf{k}j \lambda} = n_{\mathbf{k} \lambda} U_{\mathbf{k}i \lambda}, \qquad (17) $$

where $n_{\mathbf{k}\lambda}$, the eigenvalues, are the so-called natural occupations and where $U_{\mathbf{k}j\lambda}$, the eigenvector coefficients, form the natural orbitals.

In other words, the natural orbitals $\phi_{\mathbf{k}\lambda}(\mathbf{r}) = \langle\mathbf{r}|\mathbf{k}\lambda\rangle$ can be obtained from the unitary matrix $U_{\mathbf{k}\lambda}$:

$$ \phi_{\mathbf{k}\lambda}(\mathbf{r}) = \sum_{i=1}^{N_b} U_{\mathbf{k}i\lambda} \varphi_{\mathbf{k}i}(\mathbf{r}). \qquad (18) $$

Then all the one-body operator expectation values are readily obtained. For instance, the kinetic energy $T$ is calculated as

$$ T = \frac{1}{N_k} \sum_{\mathbf{k}\lambda} n_{\mathbf{k}\lambda} \langle \mathbf{k}\lambda | \hat{T} | \mathbf{k}\lambda \rangle \qquad (19) $$

$$ = \frac{1}{N_k} \sum_{\mathbf{k}\lambda} n_{\mathbf{k}\lambda} \sum_{ij} U_{\mathbf{k}i\lambda}^{*} U_{\mathbf{k}j\lambda} \langle \mathbf{k}i | \hat{T} | \mathbf{k}j \rangle. \qquad (20) $$

Finally, we would like to emphasize that the formal proof of the conservation of electron count [21] requires that the state range in the internal sum of $G_0$ in Eq. (2) is the same as the one used in the basis expansion in Eq. (7). This restriction is enforced in all our calculations.

Table I summarizes the numerical parameters used for the eight crystals considered in this paper. All the calculations for face-centered cubic crystals reported in this paper use four shifted **k**-point grids, as commonly used in ABINIT. The grid discretization is 4×4×4, which yields 256 **k** points in the full Brillouin zone and ten in the irreducible wedge. The calculation on hexagonal boron nitride (h-BN) uses a Γ-centered 12×12×6 **k**-point grid for exact exchange and a 4×4×2 grid for the rest. These parameters ensure smooth energy-volume curves and total energies that are converged within 0.02 hartree (see Supplemental Material [54]). We evaluate the computational effort to scale as $O(N^4)$, similar to a conventional GW calculation. However, here both the calculations of W and of Σ scale as $O(N^4)$, whereas in conventional GW calculations of band structures only the W part does.

<figure><img src="image_4.png" /><figcaption>FIG. 1. HF total energies as a function of lattice constant with different codes and pseudopotentials for bulk silicon. The equilibrium lattice constant for each calculation is given in the legend.</figcaption></figure>

## B. Adequate norm-conserving pseudopotentials

As mentioned in the previous paragraph, our implementation uses norm-conserving pseudopotentials. In the preliminary stages of our paper, we concluded that while the details of the pseudopotential are not critical when studying band structures they become of the utmost importance when investigating structural properties.

Norm-conserving pseudopotentials are designed to reproduce the electronic and energetic properties of a given mean-field approximation. For instance, using a PBE pseudopotential for a hybrid functional is not advised in principle. As no “GW-suitable” pseudopotentials exist, we have enforced the minimal requirement that the selected pseudopotential be able to reproduce HF structural properties.

In Fig. 1, we show a wide comparison among codes and techniques for bulk silicon at the HF level of theory. Silicon is chosen as a typical example. The results for two other crystals are reported in the Supplemental Material [54] with identical conclusions. The all-electron (AE) of CRYSTAL [55] with the accurate basis set designed by Heyd *et al.* [56] and the projector augmented-wave (PAW) results of VASP [39] agree very well. We consider them as the reference, since by construction the Gaussian basis set used in CRYSTAL describes all the electrons at once and since in the PAW framework, though frozen-core, the core-valence interactions are completely recalculated for each approximation.

Then we turn to the regular PBE pseudopotential obtained from the PSEUDODOJO suite [57]. This pseudopotential is highly tested and should be rather transferable as it relies on Hamann’s ONCVPSP scheme [52] that introduces several projectors per angular momentum. However, based on Fig. 1, the HF energy-volume curve departs significantly from the reference. This error is intrinsic to the pseudopotential because using it in QUANTUM ESPRESSO [58] gives the exact same result. In our opinion, the inability of the ONCVPSP pseudopotentials to reproduce HF energy-volume curves is not due to Hamann’s scheme itself, but rather due to the practical choice of large cutoff radii selected in the PSEUDODOJO
<table><caption>TABLE I. Energy cutoff for the wave functions ($E_{\text{cut}}^{\varphi}$, in hartree), energy cutoff for the plane-wave expansion of the screened Coulomb interaction $W$ ($E_{\text{cut}}^{\text{W}}$, in hartree), number of bands ($N_b$) in $G_0$ and in Eqs. (9) and (10), number of pure imaginary frequencies ($N_{\omega_W}$) used for $W(i\omega)$, and number of frequency points ($N_{\omega_{\Sigma}}$), used for the quadrature along the imaginary axis in Eq. (10).</caption><thead><tr><th></th><th>Si</th><th>C</th><th>SiC</th><th>zb-BN</th><th>AlP</th><th>AlAs</th><th>Ge</th><th>h-BN</th></tr></thead><tbody><tr><td>$E_{\text{cut}}^{\varphi}$</td><td>12</td><td>25</td><td>50</td><td>40</td><td>35</td><td>35</td><td>30</td><td>55</td></tr><tr><td>$E_{\text{cut}}^{\text{W}}$</td><td>8</td><td>15</td><td>15</td><td>15</td><td>15</td><td>10</td><td>10</td><td>25</td></tr><tr><td>$N_b$</td><td>150</td><td>175</td><td>175</td><td>175</td><td>175</td><td>175</td><td>175</td><td>1400</td></tr><tr><td>$N_{\omega_W}$</td><td>120</td><td>120</td><td>120</td><td>120</td><td>120</td><td>120</td><td>120</td><td>40</td></tr><tr><td>$N_{\omega_{\Sigma}}$</td><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td><td></td></tr></tbody></table>

natural orbitals in the gKS basis:

$$ \sum_{j=1}^{N_b} \gamma_{\mathbf{k}i j} U_{\mathbf{k}j \lambda} = n_{\mathbf{k} \lambda} U_{\mathbf{k}i \lambda}, \qquad (17) $$

where $n_{\mathbf{k}\lambda}$, the eigenvalues, are the so-called natural occupations and where $U_{\mathbf{k}j\lambda}$, the eigenvector coefficients, form the natural orbitals.

In other words, the natural orbitals $\phi_{\mathbf{k}\lambda}(\mathbf{r}) = \langle\mathbf{r}|\mathbf{k}\lambda\rangle$ can be obtained from the unitary matrix $U_{\mathbf{k}\lambda}$:

$$ \phi_{\mathbf{k}\lambda}(\mathbf{r}) = \sum_{i=1}^{N_b} U_{\mathbf{k}i\lambda} \varphi_{\mathbf{k}i}(\mathbf{r}). \qquad (18) $$

Then all the one-body operator expectation values are readily obtained. For instance, the kinetic energy $T$ is calculated as

$$ T = \frac{1}{N_k} \sum_{\mathbf{k}\lambda} n_{\mathbf{k}\lambda} \langle \mathbf{k}\lambda | \hat{T} | \mathbf{k}\lambda \rangle \qquad (19) $$

$$ = \frac{1}{N_k} \sum_{\mathbf{k}\lambda} n_{\mathbf{k}\lambda} \sum_{ij} U_{\mathbf{k}i\lambda}^{*} U_{\mathbf{k}j\lambda} \langle \mathbf{k}i | \hat{T} | \mathbf{k}j \rangle. \qquad (20) $$

Finally, we would like to emphasize that the formal proof of the conservation of electron count [21] requires that the state range in the internal sum of $G_0$ in Eq. (2) is the same as the one used in the basis expansion in Eq. (7). This restriction is enforced in all our calculations.

Table I summarizes the numerical parameters used for the eight crystals considered in this paper. All the calculations for face-centered cubic crystals reported in this paper use four shifted k-point grids, as commonly used in ABINIT. The grid discretization is 4×4×4, which yields 256 k points in the full Brillouin zone and ten in the irreducible wedge. The calculation on hexagonal boron nitride (h-BN) uses a Γ-centered 12×12×6 k-point grid for exact exchange and a 4×4×2 grid for the rest. These parameters ensure smooth energy-volume curves and total energies that are converged within 0.02 hartree (see Supplemental Material [54]). We evaluate the computational effort to scale as $O(N^4)$, similar to a conventional GW calculation. However, here both the calculations of W and of Σ scale as $O(N^4)$, whereas in conventional GW calculations of band structures only the W part does.

<figure><img src="image_4.png" /><figcaption>FIG. 1. HF total energies as a function of lattice constant with different codes and pseudopotentials for bulk silicon. The equilibrium lattice constant for each calculation is given in the legend.</figcaption></figure>

## B. Adequate norm-conserving pseudopotentials

As mentioned in the previous paragraph, our implementation uses norm-conserving pseudopotentials. In the preliminary stages of our paper, we concluded that while the details of the pseudopotential are not critical when studying band structures they become of the utmost importance when investigating structural properties.

Norm-conserving pseudopotentials are designed to reproduce the electronic and energetic properties of a given mean-field approximation. For instance, using a PBE pseudopotential for a hybrid functional is not advised in principle. As no “GW-suitable” pseudopotentials exist, we have enforced the minimal requirement that the selected pseudopotential be able to reproduce HF structural properties.

In Fig. 1, we show a wide comparison among codes and techniques for bulk silicon at the HF level of theory. Silicon is chosen as a typical example. The results for two other crystals are reported in the Supplemental Material [54] with identical conclusions. The all-electron (AE) of CRYSTAL [55] with the accurate basis set designed by Heyd *et al.* [56] and the projector augmented-wave (PAW) results of VASP [39] agree very well. We consider them as the reference, since by construction the Gaussian basis set used in CRYSTAL describes all the electrons at once and since in the PAW framework, though frozen-core, the core-valence interactions are completely recalculated for each approximation.

Then we turn to the regular PBE pseudopotential obtained from the PSEUDODOJO suite [57]. This pseudopotential is highly tested and should be rather transferable as it relies on Hamann’s ONCVPSP scheme [52] that introduces several projectors per angular momentum. However, based on Fig. 1, the HF energy-volume curve departs significantly from the reference. This error is intrinsic to the pseudopotential because using it in QUANTUM ESPRESSO [58] gives the exact same result. In our opinion, the inability of the ONCVPSP pseudopotentials to reproduce HF energy-volume curves is not due to Hamann’s scheme itself, but rather due to the practical choice of large cutoff radii selected in the PSEUDODOJO
initiative. Generating our own dedicated ONCVPSP pseudopotentials would be possible of course, but would require a significant effort. Fortunately, alternatives already exist.

In a previous work [13], one of us mitigated this problem by using pseudopotentials generated for KLI [59], which devises a local potential that simulates the nonlocal exact-exchange operator. This improves over the ONCVPSP pseudopotentials but is not quantitative enough.

Recently, in the quantum Monte Carlo (QMC) community, there has been an effort to support the design of “correlation consistent” effective core potentials (ccECPs) [60]. These norm-conserving pseudopotentials are meant to be used in combination with correlated methods beyond the usual mean-field ones. Figure 1 shows that this type of pseudopotential produces results in close agreement with AE and PAW for HF: The lattice constants match within 0.1%.

We conclude that the ccECP pseudopotentials are our preferred norm-conserving pseudopotentials to obtain quantitative results.

(i) They have been designed specifically for explicit-correlation methods and GW belongs to this family.

(ii) They are best to reproduce HF lattice constants.

The main drawback of these pseudopotentials is the high cutoff energy that is necessary to converge the total energies in plane waves. In the following, all the reported results employ ccECP pseudopotentials.

## IV. GW DENSITY MATRIX IN CRYSTALLINE SOLIDS

As summarized in Appendix A, the spin-summed natural occupations $n_{k\lambda}$ should continuously span the range from 0 to 2 at variance with regular Fermi-Dirac ground-state occupations $f_{ki}$ that are only 0 or 2.

These natural occupations for realistic crystalline solids can be compared to the momentum distribution function $n_k$ for the homogeneous electron gas in Refs. [36,61], but for the homogeneous electron gas the momentum $k$ is enough to uniquely characterize the quantum state. For solids, we need an additional quantum number $\lambda$ (similar to a band index).

In Fig. 2, we represent the natural occupations $n_{k\lambda}$ for a fixed $k$ point (1/8, 0, 0). This particular $k$ point was selected as an example: The other $k$ points produce very similar results. The PBE occupations $f_{ki}$ are shown as a reference. Then the natural occupations for $\gamma^{HF}$@PBE, the static part of the density matrix, are plotted. While their sum precisely equalizes the number of electrons $N_e$, the values can exceed 2 and be below 0. These occupation values violate the constraints of the exact density matrix. However, after adding the dynamic correlation, the $\gamma^{GW}$@PBE has all its spin-summed natural occupations between 0 and 2. Four natural orbitals have an occupation close to 1.8–1.9 and then many more (15 or so) have a nonvanishing occupation. A PBE mean-field starting point was chosen to magnify the effect. Starting from HF would yield perfectly normal natural occupations for $\gamma^{HF}$. The overall shape of the occupation is similar to the homogeneous electron gas result [36,61].

However, one notices that $N(k) = \sum_{\lambda} n_{k\lambda}$ for a given $k$ slightly deviates from the number of electrons $N_e$. This intriguing observation does not violate an exact constraint. We only proved mathematically [21] that the sum of the natural

<figure><img src="image_5.png" /><figcaption>FIG. 2. Spin-summed natural occupations obtained from PBE, $\gamma^{HF}$@PBE, and $\gamma^{GW}$@PBE density matrices in bulk silicon for $k$ point (1/8, 0, 0) in reciprocal lattice vectors. The natural occupations are ordered by descending values, from 1 to 150, which is $N_b$, the dimension of the matrix. The $x$ axis has been cut to show the first and the last values.</figcaption></figure>

occupation over the whole Brillouin zone

$$ \sum_{\mathbf{k}} \sum_{\lambda} n_{\mathbf{k}\lambda} = N_e \qquad (21) $$

is valid. Appendix B demonstrates that this variation with $k$ is possible.

As this observation can be considered surprising when compared to the usual mean-field occupations $f_{ki}$, it is insightful to monitor the sum $N(k)$ across the Brillouin zone. In Fig. 3, we report the deviation $\Delta N(k) = N(k) - N_e$ in a cut plane in the Brillouin zone. The $\Delta N(k)$ function is interpolated from a refined $\Gamma$-centered 6×6×6 with four shifts

<figure><img src="image_6.png" /><figcaption>FIG. 3. Electron count deviation $\Delta N(k)$ as a function of the wave vector $k$ in the first Brillouin zone for bulk silicon. $k = k_1\mathbf{b}_1 + k_2\mathbf{b}_2 + k_3\mathbf{b}_3$ is reported in reduced coordinates ($\mathbf{b}_i$ are the reciprocal lattice vectors). The plane $k_3 = 0$ is represented. The isoline $\Delta N(k) = 0$ is drawn with a bold black line. Special points $\Gamma$ and $X$ are marked.</figcaption></figure>
initiative. Generating our own dedicated ONCVPSP pseudopotentials would be possible of course, but would require a significant effort. Fortunately, alternatives already exist.

In a previous work [13], one of us mitigated this problem by using pseudopotentials generated for KLI [59], which devises a local potential that simulates the nonlocal exact-exchange operator. This improves over the ONCVPSP pseudopotentials but is not quantitative enough.

Recently, in the quantum Monte Carlo (QMC) community, there has been an effort to support the design of “correlation consistent” effective core potentials (ccECPs) [60]. These norm-conserving pseudopotentials are meant to be used in combination with correlated methods beyond the usual mean-field ones. Figure 1 shows that this type of pseudopotential produces results in close agreement with AE and PAW for HF: The lattice constants match within 0.1%.

We conclude that the ccECP pseudopotentials are our preferred norm-conserving pseudopotentials to obtain quantitative results.

(i) They have been designed specifically for explicit-correlation methods and GW belongs to this family.

(ii) They are best to reproduce HF lattice constants.

The main drawback of these pseudopotentials is the high cutoff energy that is necessary to converge the total energies in plane waves. In the following, all the reported results employ ccECP pseudopotentials.

## IV. GW DENSITY MATRIX IN CRYSTALLINE SOLIDS

As summarized in Appendix A, the spin-summed natural occupations $n_{k\lambda}$ should continuously span the range from 0 to 2 at variance with regular Fermi-Dirac ground-state occupations $f_{ki}$ that are only 0 or 2.

These natural occupations for realistic crystalline solids can be compared to the momentum distribution function $n_k$ for the homogeneous electron gas in Refs. [36,61], but for the homogeneous electron gas the momentum $k$ is enough to uniquely characterize the quantum state. For solids, we need an additional quantum number $\lambda$ (similar to a band index).

In Fig. 2, we represent the natural occupations $n_{k\lambda}$ for a fixed $k$ point (1/8, 0, 0). This particular $k$ point was selected as an example: The other $k$ points produce very similar results. The PBE occupations $f_{ki}$ are shown as a reference. Then the natural occupations for $\gamma^{HF}$@PBE, the static part of the density matrix, are plotted. While their sum precisely equalizes the number of electrons $N_e$, the values can exceed 2 and be below 0. These occupation values violate the constraints of the exact density matrix. However, after adding the dynamic correlation, the $\gamma^{GW}$@PBE has all its spin-summed natural occupations between 0 and 2. Four natural orbitals have an occupation close to 1.8–1.9 and then many more (15 or so) have a nonvanishing occupation. A PBE mean-field starting point was chosen to magnify the effect. Starting from HF would yield perfectly normal natural occupations for $\gamma^{HF}$. The overall shape of the occupation is similar to the homogeneous electron gas result [36,61].

However, one notices that $N(k) = \sum_{\lambda} n_{k\lambda}$ for a given $k$ slightly deviates from the number of electrons $N_e$. This intriguing observation does not violate an exact constraint. We only proved mathematically [21] that the sum of the natural

<figure><img src="image_5.png" /><figcaption>FIG. 2. Spin-summed natural occupations obtained from PBE, $\gamma^{HF}$@PBE, and $\gamma^{GW}$@PBE density matrices in bulk silicon for $k$ point (1/8, 0, 0) in reciprocal lattice vectors. The natural occupations are ordered by descending values, from 1 to 150, which is $N_b$, the dimension of the matrix. The $x$ axis has been cut to show the first and the last values.</figcaption></figure>

occupation over the whole Brillouin zone

$$ \sum_{\mathbf{k}} \sum_{\lambda} n_{k\lambda} = N_e \quad (21) $$

is valid. Appendix B demonstrates that this variation with $k$ is possible.

As this observation can be considered surprising when compared to the usual mean-field occupations $f_{ki}$, it is insightful to monitor the sum $N(k)$ across the Brillouin zone. In Fig. 3, we report the deviation $\Delta N(k) = N(k) - N_e$ in a cut plane in the Brillouin zone. The $\Delta N(k)$ function is interpolated from a refined $\Gamma$-centered 6×6×6 with four shifts

<figure><img src="image_6.png" /><figcaption>FIG. 3. Electron count deviation $\Delta N(k)$ as a function of the wave vector $k$ in the first Brillouin zone for bulk silicon. $k = k_1b_1 + k_2b_2 + k_3b_3$ is reported in reduced coordinates ($b_i$ are the reciprocal lattice vectors). The plane $k_3 = 0$ is represented. The isoline $\Delta N(k) = 0$ is drawn with a bold black line. Special points $\Gamma$ and $X$ are marked.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 4. Total energies for the water molecule in the gas phase (left-hand panel) and crystalline silicon (right-hand panel) as a function of the gKS starting point. The water results were extracted from a previous study [21]. The silicon QMC value comes from Ref. [64].</figcaption></figure>

k-point grid (864 points in the full Brillouin zone). We use the Shankland-Koelling-Wood interpolation technique [62] as implemented in ABIPY [40]. The numerical integration of $\Delta N(k)$ over the whole Brillouin zone yields $4 \times 10^{-4}$, which is very close to the expected zero.

From Fig. 3, we observe an electron transfer from the $\Gamma$ point region to the Brillouin zone edge. The weight transfer is not large (at most $\approx 0.01-0.02$), but still sizable. The electron count $N(k)$ is an observable and could be possibly measured in angle-resolved photoemission spectroscopy [63]. Note that this electron count transfer is a pure electronic correlation effect. Any static approximation of the self-energy $\Sigma$ would nullify it.

## V. STRUCTURAL PROPERTIES OF CRYSTALLINE SOLIDS WITHIN THE GW DENSITY MATRIX AND RPA

### A. Covalent-bonded crystals

In this section, we analyze the calculation of structural parameters for crystalline solids using the total energy expressions introduced in Eqs. (12) and (15). Our main question is which expression works “best” in the context of one-shot calculations, i.e., which expression best approximates a hypothetical reference scGW that is not currently available for crystalline systems. In molecules, where reference scGW calculations were produced [21], the accuracy of the $\gamma^{GW}$ total energy was demonstrated.

A way to measure the robustness of a one-shot total energy expression is to explore its sensitivity to the starting point. Here we use the PBEh($\alpha$) hybrid functional family:

$$V_{xc} = \alpha \Sigma_x + (1 - \alpha) v_x^{\text{PBE}} + v_c^{\text{PBE}}, \quad (22)$$

where the parameter $\alpha$ controls the amount of exact-exchange $\Sigma_x$.

Calculations were carried out for seven covalent crystals (Si, C, SiC, zb-BN, AlP, AlAs, and Ge). In the main text, we will mostly report silicon results. However, the complete set of results is made available as Supplemental Material [54].

In Fig. 4, we compare the total energy behavior for two different systems: Water, a small molecular system, and crystalline silicon. The results for the water molecule were extracted from Ref. [21] that was using a different

implementation based on Gaussian basis [65]. The figure reports the total energies for PBEh($\alpha$), for $E_{\text{total}}^{\gamma^{GW}}$, for $E_{\text{total}}^{\text{RPA}}$, and when available for scGW. The overall similarity between the two panels is striking: RPA is rather sensitive to the starting gKS, whereas $\gamma^{GW}$ is much less so. RPA increases with $\alpha$, whereas $\gamma^{GW}$ decreases. RPA and $\gamma^{GW}$ rejoin for large values of $\alpha$.

For water, where the scGW reference exists, the RPA and $\gamma^{GW}$ total energies give the best approximation of the full scGW total energy when they are equal. Owing to the similarity between the two panels of Fig. 4, we can reasonably anticipate that in bulk silicon, the scGW total energy will be best approximated by $\gamma^{GW}$ and by RPA with PBEh(0.75), PBEh(1.00), or HF, however with no formal proof.

If supplied with the scGW Green's function, all total energy formulas should match [14,15,19]. The above results tend to make us think that the noninteracting Green's function $G_0$ for PBEh(0.75), PBEh(1.00), or HF is close to the scGW Green's function $G$ for both the molecular and the solid-state systems.

The striking difference between the molecule and the solid in Fig. 4 is the agreement or the disagreement with respect to the gKS total energies. While for the molecule the RPA and the $\gamma^{GW}$ were much undershooting the total energy with respect to PBEh($\alpha$) (too negative correlation energies), the match is very good for crystalline silicon. The accurate quantum Monte Carlo approach reports -7.8644 hartree for silicon [64], whereas $\gamma^{GW}$ and RPA respectively give -7.8708 and -7.8692 hartree with the same ccECP pseudopotentials. This excellent agreement of the absolute total energies reminds us about the amazingly accurate scGW energies in the homogeneous electron gas [4-6] that precisely match the quantum Monte Carlo values [66].

We conclude with some reasonable confidence that the scGW total energy in bulk silicon is certainly close to the gKS, close to $E_{\text{total}}^{\gamma^{GW}}$, and close to RPA@PBEh(1.00). We also stress that RPA@PBE, which is the most commonly accepted implementation of RPA functional, underestimates noticeably the total energy. Our confidence in these results is further strengthened when considering the results for the other six crystalline systems reported in the Supplemental Material [54]. They all support the same conclusion.
<figure><img src="image_7.png" /><figcaption>FIG. 4. Total energies for the water molecule in the gas phase (left-hand panel) and crystalline silicon (right-hand panel) as a function of the gKS starting point. The water results were extracted from a previous study [21]. The silicon QMC value comes from Ref. [64].</figcaption></figure>

k-point grid (864 points in the full Brillouin zone). We use the Shankland-Koelling-Wood interpolation technique [62] as implemented in ABIPY [40]. The numerical integration of $\Delta N(k)$ over the whole Brillouin zone yields $4 \times 10^{-4}$, which is very close to the expected zero.

From Fig. 3, we observe an electron transfer from the $\Gamma$ point region to the Brillouin zone edge. The weight transfer is not large (at most $\approx 0.01-0.02$), but still sizable. The electron count $N(k)$ is an observable and could be possibly measured in angle-resolved photoemission spectroscopy [63]. Note that this electron count transfer is a pure electronic correlation effect. Any static approximation of the self-energy $\Sigma$ would nullify it.

## V. STRUCTURAL PROPERTIES OF CRYSTALLINE SOLIDS WITHIN THE GW DENSITY MATRIX AND RPA

### A. Covalent-bonded crystals

In this section, we analyze the calculation of structural parameters for crystalline solids using the total energy expressions introduced in Eqs. (12) and (15). Our main question is which expression works “best” in the context of one-shot calculations, i.e., which expression best approximates a hypothetical reference scGW that is not currently available for crystalline systems. In molecules, where reference scGW calculations were produced [21], the accuracy of the $\gamma^{GW}$ total energy was demonstrated.

A way to measure the robustness of a one-shot total energy expression is to explore its sensitivity to the starting point. Here we use the PBEh($\alpha$) hybrid functional family:

$$V_{xc} = \alpha \Sigma_x + (1 - \alpha) v_x^{\text{PBE}} + v_c^{\text{PBE}}, \quad (22)$$

where the parameter $\alpha$ controls the amount of exact-exchange $\Sigma_x$.

Calculations were carried out for seven covalent crystals (Si, C, SiC, zb-BN, AlP, AlAs, and Ge). In the main text, we will mostly report silicon results. However, the complete set of results is made available as Supplemental Material [54].

In Fig. 4, we compare the total energy behavior for two different systems: Water, a small molecular system, and crystalline silicon. The results for the water molecule were extracted from Ref. [21] that was using a different

implementation based on Gaussian basis [65]. The figure reports the total energies for PBEh($\alpha$), for $E_{\text{total}}^{\gamma^{GW}}$, for $E_{\text{total}}^{\text{RPA}}$, and when available for scGW. The overall similarity between the two panels is striking: RPA is rather sensitive to the starting gKS, whereas $\gamma^{GW}$ is much less so. RPA increases with $\alpha$, whereas $\gamma^{GW}$ decreases. RPA and $\gamma^{GW}$ rejoin for large values of $\alpha$.

For water, where the scGW reference exists, the RPA and $\gamma^{GW}$ total energies give the best approximation of the full scGW total energy when they are equal. Owing to the similarity between the two panels of Fig. 4, we can reasonably anticipate that in bulk silicon, the scGW total energy will be best approximated by $\gamma^{GW}$ and by RPA with PBEh(0.75), PBEh(1.00), or HF, however with no formal proof.

If supplied with the scGW Green's function, all total energy formulas should match [14,15,19]. The above results tend to make us think that the noninteracting Green's function $G_0$ for PBEh(0.75), PBEh(1.00), or HF is close to the scGW Green's function G for both the molecular and the solid-state systems.

The striking difference between the molecule and the solid in Fig. 4 is the agreement or the disagreement with respect to the gKS total energies. While for the molecule the RPA and the $\gamma^{GW}$ were much undershooting the total energy with respect to PBEh($\alpha$) (too negative correlation energies), the match is very good for crystalline silicon. The accurate quantum Monte Carlo approach reports -7.8644 hartree for silicon [64], whereas $\gamma^{GW}$ and RPA respectively give -7.8708 and -7.8692 hartree with the same ccECP pseudopotentials. This excellent agreement of the absolute total energies reminds us about the amazingly accurate scGW energies in the homogeneous electron gas [4-6] that precisely match the quantum Monte Carlo values [66].

We conclude with some reasonable confidence that the scGW total energy in bulk silicon is certainly close to the gKS, close to $E_{\text{total}}^{\gamma^{GW}}$, and close to RPA@PBEh(1.00). We also stress that RPA@PBE, which is the most commonly accepted implementation of RPA functional, underestimates noticeably the total energy. Our confidence in these results is further strengthened when considering the results for the other six crystalline systems reported in the Supplemental Material [54]. They all support the same conclusion.
<figure><img src="image_8.png" /><figcaption>FIG. 5. Energy-lattice constant curves for crystalline silicon for RPA functional (left-hand panel) and $\gamma^{GW}$ energy functional (central panel). Equilibrium lattices are summarized in the right-hand panel. The experimental lattice constant is given as a reference.</figcaption></figure>

However if scGW is very accurate for solids and less accurate for finite systems, the atomization energies that measure the energy gain when forming a bulk crystal as compared to the isolated atoms are likely to have a low accuracy in scGW. This should be explored in the future.

Let us now focus on the complete energy versus lattice constant curves and check the sensitivity to the starting point not only for the total energy but also for the equilibrium lattice constant *a* and for the bulk modulus *B*.

Figure 5 reports on the same scale the RPA and the $\gamma^{GW}$ total energies for different gKS starting points and in the right-hand panel the equilibrium lattice constant as a function of the gKS starting point. We can see again that the RPA total energy in the left-hand panel of Fig. 5 is much more sensitive to the starting point compared to $\gamma^{GW}$ total energy in the central panel. However, when focusing on the equilibrium lattice constant itself, the sensitivity to the starting point is much weaker: At most 0.03 bohr.

All panels in Fig. 5 support again the same conclusion: The RPA and $\gamma^{GW}$ total energies agree best when using PBEh($\alpha$) with large $\alpha$ or even when using HF. The statements drawn for silicon perfectly hold for the other six crystalline systems presented in the Supplemental Material [54].

Finally, we summarize the lattice constants and bulk moduli of the seven covalent crystals that we have studied in Table II. We focus on two gKS starting points PBEh(0.00) (i.e., standard PBE) and PBEh(0.75). While for the PBE starting point the RPA@PBE and $\gamma^{GW}$@PBE lattice constants differ by about 0.06 bohr, the RPA@PBEh(0.75) and $\gamma^{GW}$@PBEh(0.75) lattice constants always agree within 0.01 bohr. The same type of conclusion holds for the bulk modulus *B*. This is another proof that PBEh(0.75) is a good starting point to evaluate the different GW-based total energy expressions. Most probably, the properties *a* and *B* evaluated with RPA@PBEh(0.75) or $\gamma^{GW}$@PBEh(0.75) are reliable estimates to the scGW result.

In the end, we also compare to experiment. Table II shows that all the GW-based energy expressions yield structural properties in excellent agreement with respect to the experiment: A 0.1% deviation for lattice constants and 8% for bulk moduli. The different expressions and starting points have a minor influence on this. This conclusion is valid for the covalent crystals. However, it is worth considering whether this conclusion still holds for weak van der Waals interactions, which are one of the attractive features of RPA.

<table><caption>TABLE II. Lattice constants *a* (bohr), bulk moduli *B* (GPa), mean absolute error (MAE), and mean absolute percentage error (MAPE) of seven covalent crystals obtained with several total energy methods.</caption><thead><tr><th></th><th></th><th></th><th>Si</th><th>C</th><th>SiC</th><th>zb-BN</th><th>AlP</th><th>AlAs</th><th>Ge</th><th>MAE</th><th>MAPE</th></tr></thead><tbody><tr><td rowspan="2">PBEh(0.00)</td><td><i>a</i></td><td>10.407</td><td>6.774</td><td>8.310</td><td>6.861</td><td>10.471</td><td>10.847</td><td>10.947</td><td>0.101</td><td>1.00</td></tr><tr><td><i>B</i></td><td>88.44</td><td>441.76</td><td>216.68</td><td>380.18</td><td>82.72</td><td>69.98</td><td>62.24</td><td>21.961</td><td>10.5</td></tr><tr><td rowspan="2">RPA@PBEh(0.00)</td><td><i>a</i></td><td>10.282</td><td>6.815</td><td>8.257</td><td>6.867</td><td>10.308</td><td>10.725</td><td>10.769</td><td>0.011</td><td>0.12</td></tr><tr><td><i>B</i></td><td>97.08</td><td>477.44</td><td>235.25</td><td>401.35</td><td>96.78</td><td>83.57</td><td>86.10</td><td>12.123</td><td>8.7</td></tr><tr><td rowspan="2">γ<sup>GW</sup> @PBEh(0.00)</td><td><i>a</i></td><td>10.232</td><td>6.713</td><td>8.180</td><td>6.800</td><td>10.287</td><td>10.668</td><td>10.708</td><td>0.011</td><td>0.12</td></tr><tr><td><i>B</i></td><td>104.58</td><td>420.75</td><td>262.57</td><td>417.06</td><td>101.56</td><td>87.73</td><td>85.70</td><td>11.451</td><td>8.3</td></tr><tr><td rowspan="2">PBEh(0.75)</td><td><i>a</i></td><td>10.260</td><td>6.675</td><td>8.170</td><td>6.728</td><td>10.322</td><td>10.680</td><td>10.682</td><td>0.036</td><td>0.49</td></tr><tr><td><i>B</i></td><td>111.93</td><td>502.28</td><td>270.10</td><td>461.11</td><td>102.44</td><td>88.02</td><td>87.30</td><td>31.026</td><td>15.7</td></tr><tr><td rowspan="2">RPA@PBEh(0.75)</td><td><i>a</i></td><td>10.254</td><td>6.759</td><td>8.210</td><td>6.831</td><td>10.301</td><td>10.683</td><td>10.738</td><td>0.012</td><td>0.13</td></tr><tr><td><i>B</i></td><td>103.37</td><td>438.89</td><td>252.60</td><td>408.86</td><td>99.96</td><td>85.94</td><td>81.01</td><td>10.407</td><td>7.8</td></tr><tr><td rowspan="2">γ<sup>GW</sup> @PBEh(0.75)</td><td><i>a</i></td><td>10.255</td><td>6.748</td><td>8.207</td><td>6.825</td><td>10.301</td><td>10.681</td><td>10.732</td><td>0.010</td><td>0.11</td></tr><tr><td><i>B</i></td><td>103.91</td><td>436.08</td><td>254.39</td><td>410.23</td><td>100.22</td><td>86.59</td><td>82.07</td><td>11.619</td><td>8.4</td></tr><tr><td rowspan="2">Expt.</td><td><i>a</i></td><td>10.261</td><td>6.741</td><td>8.213</td><td>6.833</td><td>10.301</td><td>10.675</td><td>10.692</td><td></td><td></td></tr><tr><td><i>B</i></td><td>99</td><td>443</td><td>225</td><td>369-400</td><td>86</td><td>77</td><td>76</td><td></td><td></td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>FIG. 5. Energy-lattice constant curves for crystalline silicon for RPA functional (left-hand panel) and $\gamma^{GW}$ energy functional (central panel). Equilibrium lattices are summarized in the right-hand panel. The experimental lattice constant is given as a reference.</figcaption></figure>

However if scGW is very accurate for solids and less accurate for finite systems, the atomization energies that measure the energy gain when forming a bulk crystal as compared to the isolated atoms are likely to have a low accuracy in scGW. This should be explored in the future.

Let us now focus on the complete energy versus lattice constant curves and check the sensitivity to the starting point not only for the total energy but also for the equilibrium lattice constant *a* and for the bulk modulus *B*.

Figure 5 reports on the same scale the RPA and the $\gamma^{GW}$ total energies for different gKS starting points and in the right-hand panel the equilibrium lattice constant as a function of the gKS starting point. We can see again that the RPA total energy in the left-hand panel of Fig. 5 is much more sensitive to the starting point compared to $\gamma^{GW}$ total energy in the central panel. However, when focusing on the equilibrium lattice constant itself, the sensitivity to the starting point is much weaker: At most 0.03 bohr.

All panels in Fig. 5 support again the same conclusion: The RPA and $\gamma^{GW}$ total energies agree best when using PBEh($\alpha$) with large $\alpha$ or even when using HF. The statements drawn for silicon perfectly hold for the other six crystalline systems presented in the Supplemental Material [54].

Finally, we summarize the lattice constants and bulk moduli of the seven covalent crystals that we have studied in Table II. We focus on two gKS starting points PBEh(0.00) (i.e., standard PBE) and PBEh(0.75). While for the PBE starting point the RPA@PBE and $\gamma^{GW}$@PBE lattice constants differ by about 0.06 bohr, the RPA@PBEh(0.75) and $\gamma^{GW}$@PBEh(0.75) lattice constants always agree within 0.01 bohr. The same type of conclusion holds for the bulk modulus *B*. This is another proof that PBEh(0.75) is a good starting point to evaluate the different GW-based total energy expressions. Most probably, the properties *a* and *B* evaluated with RPA@PBEh(0.75) or $\gamma^{GW}$@PBEh(0.75) are reliable estimates to the scGW result.

In the end, we also compare to experiment. Table II shows that all the GW-based energy expressions yield structural properties in excellent agreement with respect to the experiment: A 0.1% deviation for lattice constants and 8% for bulk moduli. The different expressions and starting points have a minor influence on this. This conclusion is valid for the covalent crystals. However, it is worth considering whether this conclusion still holds for weak van der Waals interactions, which are one of the attractive features of RPA.

<table><caption>TABLE II. Lattice constants *a* (bohr), bulk moduli *B* (GPa), mean absolute error (MAE), and mean absolute percentage error (MAPE) of seven covalent crystals obtained with several total energy methods.</caption><thead><tr><th></th><th></th><th></th><th>Si</th><th>C</th><th>SiC</th><th>zb-BN</th><th>AlP</th><th>AlAs</th><th>Ge</th><th>MAE</th><th>MAPE</th></tr></thead><tbody><tr><td rowspan="2">PBEh(0.00)</td><td><i>a</i></td><td>10.407</td><td>6.774</td><td>8.310</td><td>6.861</td><td>10.471</td><td>10.847</td><td>10.947</td><td>0.101</td><td>1.00</td></tr><tr><td><i>B</i></td><td>88.44</td><td>441.76</td><td>216.68</td><td>380.18</td><td>82.72</td><td>69.98</td><td>62.24</td><td>21.961</td><td>10.5</td></tr><tr><td rowspan="2">RPA@PBEh(0.00)</td><td><i>a</i></td><td>10.282</td><td>6.815</td><td>8.257</td><td>6.867</td><td>10.308</td><td>10.725</td><td>10.769</td><td>0.011</td><td>0.12</td></tr><tr><td><i>B</i></td><td>97.08</td><td>477.44</td><td>235.25</td><td>401.35</td><td>96.78</td><td>83.57</td><td>86.10</td><td>12.123</td><td>8.7</td></tr><tr><td rowspan="2">γ<sup>GW</sup> @PBEh(0.00)</td><td><i>a</i></td><td>10.232</td><td>6.713</td><td>8.180</td><td>6.800</td><td>10.287</td><td>10.668</td><td>10.708</td><td>0.011</td><td>0.12</td></tr><tr><td><i>B</i></td><td>104.58</td><td>420.75</td><td>262.57</td><td>417.06</td><td>101.56</td><td>87.73</td><td>85.70</td><td>11.451</td><td>8.3</td></tr><tr><td rowspan="2">PBEh(0.75)</td><td><i>a</i></td><td>10.260</td><td>6.675</td><td>8.170</td><td>6.728</td><td>10.322</td><td>10.680</td><td>10.682</td><td>0.036</td><td>0.49</td></tr><tr><td><i>B</i></td><td>111.93</td><td>502.28</td><td>270.10</td><td>461.11</td><td>102.44</td><td>88.02</td><td>87.30</td><td>31.026</td><td>15.7</td></tr><tr><td rowspan="2">RPA@PBEh(0.75)</td><td><i>a</i></td><td>10.254</td><td>6.759</td><td>8.210</td><td>6.831</td><td>10.301</td><td>10.683</td><td>10.738</td><td>0.012</td><td>0.13</td></tr><tr><td><i>B</i></td><td>103.37</td><td>438.89</td><td>252.60</td><td>408.86</td><td>99.96</td><td>85.94</td><td>81.01</td><td>10.407</td><td>7.8</td></tr><tr><td rowspan="2">γ<sup>GW</sup> @PBEh(0.75)</td><td><i>a</i></td><td>10.255</td><td>6.748</td><td>8.207</td><td>6.825</td><td>10.301</td><td>10.681</td><td>10.732</td><td>0.010</td><td>0.11</td></tr><tr><td><i>B</i></td><td>103.91</td><td>436.08</td><td>254.39</td><td>410.23</td><td>100.22</td><td>86.59</td><td>82.07</td><td>11.619</td><td>8.4</td></tr><tr><td rowspan="2">Expt.</td><td><i>a</i></td><td>10.261</td><td>6.741</td><td>8.213</td><td>6.833</td><td>10.301</td><td>10.675</td><td>10.692</td><td></td><td></td></tr><tr><td><i>B</i></td><td>99</td><td>443</td><td>225</td><td>369-400</td><td>86</td><td>77</td><td>76</td><td></td><td></td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>FIG. 6. Energy as a function of the spacing between the layers for h-BN with different gKS approximations and different starting points for RPA.</figcaption></figure>

## B. Van der Waals bonded layered material

In order to test if the weak van der Waals interactions would correctly be described by scGW, we analyzed a layered material, namely the h-BN. RPA@LDA and RPA@PBE were proven to be able to describe properly the spacing between the layers [13,67].

The task of evaluating the $\gamma^{GW}$ for h-BN is beyond our current computational capabilities. Indeed, for the weak van der Waals interactions, the energy scales are so low that extremely converged calculations are required. Fortunately, based on the previous discussion, we assume that the scGW total energy is also well approximated by RPA@PBEh(0.75).

In Fig. 6, we report the energy versus spacing between the layers. Local-density approximation (LDA) is known to give the correct spacing thanks to a lucky compensation of errors [67], whereas PBE that improves the exchange over LDA does not benefit from this and yields a much too large spacing. Our RPA@LDA reproduces within 0.25 bohr the earlier estimates from Refs. [13,67]. This agreement is quite good considering the computational power difference and considering the fact that we use newly developed pseudopotentials.

Next, let us comment that the RPA@PBEh(0.75) result (see Fig. 6) shows our best approximate to scGW. The obtained lattice spacing is in very good agreement with respect to the experiment (6.18 versus 6.25 bohrs).

This example shows that tuning the starting point in RPA does not destroy the quantitative agreement with respect to the experiment. This interesting conclusion calls for further studies in the future.

## VI. CONCLUSION

The linearized GW density matrix $\gamma^{GW}$ has been introduced in realistic solid-state systems. We have carefully tested the norm-conserving pseudopotential approximation and have concluded that QMC pseudopotentials, such as ccECP [60],

are compulsory in this context for accurately determining RPA and $\gamma^{GW}$ lattice constants. On a benchmark of seven covalent crystals (Si, C, SiC, zb-BN, AlP, AlAs, and Ge), we have proven numerically that $\gamma^{GW}$ actually fulfils the exact constraints: Its natural occupation numbers range from 0 to 2 (when spin is summed) and they sum up to the correct number of electrons. In addition, the correlated nature of $\gamma^{GW}$ allows the electron occupancy to reorganize across the Brillouin zone, in strong contrast with all the mean-field approaches, where the electron count remains constant in the Brillouin zone for crystals with a band gap.

The one-shot total energy expression $E_{total}^{\gamma^{GW}}$ has been found to be superior to the usual RPA total energy expression in terms of sensitivity to the gKS starting point.

We provide strong evidence to support the assumption that $\gamma^{GW}$-based total energy is a reliable substitute for scGW, which remains unachievable at present. As a cheaper alternative, RPA can also be used, but we advocate applying it on top of gKS functionals with a large content of the exact exchange, at least 75%, such as PBEh(0.75) or HF to best approximate scGW. Our last statement disagrees with the current wisdom that recommends RPA@PBE based on comparison to experiment.

Our implementation is available in the public version of the open source code ABINIT [68].

## ACKNOWLEDGMENTS

The authors acknowledge the financial support provided by the Cross-Disciplinary Program on Numerical Simulation of the French Alternative Energies and Atomic Energy Commission (ABIDM project). This work was performed using HPC resources from GENCI-CCRT-TGCC (Grant No. 2022-096018). M.R.M. thanks S. Sharma for sharing her one-body-reduced-density-matrix expertise and F.B. thanks L. Mitas for insights about the quantum Monte Carlo references.

## APPENDIX A: THE UNIT-CELL DENSITY MATRIX IN THE NATURAL ORBITAL REPRESENTATION

### 1. The diagonal density matrix $\gamma(\mathbf{r}, \mathbf{r}')$ in $\mathbf{k}$

Imposing the Born–von Kármán periodic conditions, we can write the double Fourier expansion of $\gamma(\mathbf{r}, \mathbf{r}')$ in reciprocal space as

$$ \gamma(\mathbf{r}, \mathbf{r}') = \frac{1}{N_k \Omega} \sum_{\mathbf{k}\mathbf{k}'\mathbf{G}\mathbf{G}'} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} \gamma_{\mathbf{k}\mathbf{k}'}(\mathbf{G}, \mathbf{G}')e^{-i(\mathbf{k}'+\mathbf{G}')\cdot\mathbf{r}'}, \quad (A1) $$

where $N_k$ is the number of $\mathbf{k}$ points and $\Omega$ is the volume of the unit cell and $\mathbf{G}$ and $\mathbf{G}'$ are reciprocal lattice vectors.

By virtue of the translation invariance in the unit cell, the shift of the two space indices with $\mathbf{R} = n_1\mathbf{a}_1 + n_2\mathbf{a}_2 + n_3\mathbf{a}_3$ (with $n_1, n_2, n_3 \in \mathbb{Z}$ and $\mathbf{a}_1, \mathbf{a}_2$, and $\mathbf{a}_3$ being the primitive lattice vectors) does not change the density matrix [69]:

$$ \gamma(\mathbf{r}, \mathbf{r}') = \gamma(\mathbf{r} + \mathbf{R}, \mathbf{r}' + \mathbf{R}). \quad (A2) $$

When inserting the double Fourier transform, this implies

$$ e^{i(\mathbf{k}-\mathbf{k}')\cdot\mathbf{R}} = 1. \quad (A3) $$
<figure><img src="image_9.png" /><figcaption>FIG. 6. Energy as a function of the spacing between the layers for h-BN with different gKS approximations and different starting points for RPA.</figcaption></figure>

## B. Van der Waals bonded layered material

In order to test if the weak van der Waals interactions would correctly be described by scGW, we analyzed a layered material, namely the h-BN. RPA@LDA and RPA@PBE were proven to be able to describe properly the spacing between the layers [13,67].

The task of evaluating the $\gamma^{GW}$ for h-BN is beyond our current computational capabilities. Indeed, for the weak van der Waals interactions, the energy scales are so low that extremely converged calculations are required. Fortunately, based on the previous discussion, we assume that the scGW total energy is also well approximated by RPA@PBEh(0.75).

In Fig. 6, we report the energy versus spacing between the layers. Local-density approximation (LDA) is known to give the correct spacing thanks to a lucky compensation of errors [67], whereas PBE that improves the exchange over LDA does not benefit from this and yields a much too large spacing. Our RPA@LDA reproduces within 0.25 bohr the earlier estimates from Refs. [13,67]. This agreement is quite good considering the computational power difference and considering the fact that we use newly developed pseudopotentials.

Next, let us comment that the RPA@PBEh(0.75) result (see Fig. 6) shows our best approximate to scGW. The obtained lattice spacing is in very good agreement with respect to the experiment (6.18 versus 6.25 bohrs).

This example shows that tuning the starting point in RPA does not destroy the quantitative agreement with respect to the experiment. This interesting conclusion calls for further studies in the future.

## VI. CONCLUSION

The linearized GW density matrix $\gamma^{GW}$ has been introduced in realistic solid-state systems. We have carefully tested the norm-conserving pseudopotential approximation and have concluded that QMC pseudopotentials, such as ccECP [60],

are compulsory in this context for accurately determining RPA and $\gamma^{GW}$ lattice constants. On a benchmark of seven covalent crystals (Si, C, SiC, zb-BN, AlP, AlAs, and Ge), we have proven numerically that $\gamma^{GW}$ actually fulfils the exact constraints: Its natural occupation numbers range from 0 to 2 (when spin is summed) and they sum up to the correct number of electrons. In addition, the correlated nature of $\gamma^{GW}$ allows the electron occupancy to reorganize across the Brillouin zone, in strong contrast with all the mean-field approaches, where the electron count remains constant in the Brillouin zone for crystals with a band gap.

The one-shot total energy expression $E_{total}^{\gamma^{GW}}$ has been found to be superior to the usual RPA total energy expression in terms of sensitivity to the gKS starting point.

We provide strong evidence to support the assumption that $\gamma^{GW}$-based total energy is a reliable substitute for scGW, which remains unachievable at present. As a cheaper alternative, RPA can also be used, but we advocate applying it on top of gKS functionals with a large content of the exact exchange, at least 75%, such as PBEh(0.75) or HF to best approximate scGW. Our last statement disagrees with the current wisdom that recommends RPA@PBE based on comparison to experiment.

Our implementation is available in the public version of the open source code ABINIT [68].

## ACKNOWLEDGMENTS

The authors acknowledge the financial support provided by the Cross-Disciplinary Program on Numerical Simulation of the French Alternative Energies and Atomic Energy Commission (ABIDM project). This work was performed using HPC resources from GENCI-CCRT-TGCC (Grant No. 2022-096018). M.R.M. thanks S. Sharma for sharing her one-body-reduced-density-matrix expertise and F.B. thanks L. Mitas for insights about the quantum Monte Carlo references.

## APPENDIX A: THE UNIT-CELL DENSITY MATRIX IN THE NATURAL ORBITAL REPRESENTATION

### 1. The diagonal density matrix $\gamma(\mathbf{r}, \mathbf{r}')$ in $\mathbf{k}$

Imposing the Born-von Kármán periodic conditions, we can write the double Fourier expansion of $\gamma(\mathbf{r}, \mathbf{r}')$ in reciprocal space as

$$ \gamma(\mathbf{r}, \mathbf{r}') = \frac{1}{N_k \Omega} \sum_{\mathbf{k}\mathbf{k}'\mathbf{G}\mathbf{G}'} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} \gamma_{\mathbf{k}\mathbf{k}'}(\mathbf{G}, \mathbf{G}')e^{-i(\mathbf{k}'+\mathbf{G}')\cdot\mathbf{r}'}, \quad (A1) $$

where $N_k$ is the number of $\mathbf{k}$ points and $\Omega$ is the volume of the unit cell and $\mathbf{G}$ and $\mathbf{G}'$ are reciprocal lattice vectors.

By virtue of the translation invariance in the unit cell, the shift of the two space indices with $\mathbf{R} = n_1\mathbf{a}_1 + n_2\mathbf{a}_2 + n_3\mathbf{a}_3$ (with $n_1, n_2, n_3 \in \mathbb{Z}$ and $\mathbf{a}_1, \mathbf{a}_2$, and $\mathbf{a}_3$ being the primitive lattice vectors) does not change the density matrix [69]:

$$ \gamma(\mathbf{r}, \mathbf{r}') = \gamma(\mathbf{r} + \mathbf{R}, \mathbf{r}' + \mathbf{R}). \quad (A2) $$

When inserting the double Fourier transform, this implies

$$ e^{i(\mathbf{k}-\mathbf{k}')\cdot\mathbf{R}} = 1. \quad (A3) $$
This condition is fulfilled if and only if $(\mathbf{k} - \mathbf{k}')$ belongs to the reciprocal lattice. And as $\mathbf{k}$ and $\mathbf{k}'$ both belong to the first Brillouin zone, the only possible reciprocal lattice vector the difference can match is $\mathbf{0}$.

As a consequence, one can insert the Kronecker sign $\delta_{\mathbf{k}\mathbf{k}'}$ in Eq. (A1)

$$ \gamma(\mathbf{r}, \mathbf{r}') = \frac{1}{N_k \Omega} \sum_{\mathbf{k}GG'} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} \gamma_{\mathbf{k}}(\mathbf{G}, \mathbf{G}') e^{-i(\mathbf{k}+\mathbf{G}')\cdot\mathbf{r}'} \quad (A4) $$

and obtain the desired expression that shows that $\gamma$ is block diagonal with respect to $\mathbf{k}$.

## 2. Natural orbitals and occupations

Now for each discrete value of $\mathbf{k}$, one can diagonalize $\gamma_{\mathbf{k}}$:

$$ \sum_{\mathbf{G}'} \gamma_{\mathbf{k}}(\mathbf{G}, \mathbf{G}') \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}') = n_{\mathbf{k}\lambda} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}), \quad (A5) $$

where $\tilde{u}_{\mathbf{k}\lambda}(\mathbf{G})$ are the eigenvectors and $n_{\mathbf{k}\lambda}$ the eigenvalues, indexed with $\lambda$.

By construction,

$$ \gamma(\mathbf{r}, \mathbf{r}') = \gamma(\mathbf{r}', \mathbf{r})^* \quad (A6) $$

implies

$$ \gamma_{\mathbf{k}}(\mathbf{G}, \mathbf{G}') = \gamma_{\mathbf{k}}(\mathbf{G}', \mathbf{G})^*. \quad (A7) $$

$\gamma_{\mathbf{k}}(\mathbf{G}, \mathbf{G}')$ is a Hermitian matrix and thus the $n_{\mathbf{k}\lambda}$ are real valued and the $\tilde{u}_{\mathbf{k}\lambda}$ form a unitary matrix

Using this decomposition,

$$ \gamma_{\mathbf{k}}(\mathbf{G}, \mathbf{G}') = \sum_{\lambda} n_{\mathbf{k}\lambda} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}) \tilde{u}_{\mathbf{k}\lambda}^{*}(\mathbf{G}'). \quad (A8) $$

Inserting this in Eq. (A4), we obtain

$$ \gamma(\mathbf{r}, \mathbf{r}') = \frac{1}{N_k \Omega} \sum_{\mathbf{k}\lambda} \sum_{\mathbf{G}} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}) \times n_{\mathbf{k}\lambda} \sum_{\mathbf{G}'} \tilde{u}_{\mathbf{k}\lambda}^{*}(\mathbf{G}') e^{-i(\mathbf{k}+\mathbf{G}')\cdot\mathbf{r}'} \quad (A9) $$

Let us introduce the natural orbital in real space $\phi_{\mathbf{k}\lambda}(\mathbf{r})$ that has a Bloch wave form [70]:

$$ \phi_{\mathbf{k}\lambda}(\mathbf{r}) = \frac{1}{\sqrt{N_k \Omega}} \sum_{\mathbf{G}} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}) \quad (A10) $$

$$ = \frac{1}{\sqrt{N_k \Omega}} e^{i\mathbf{k}\cdot\mathbf{r}} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{r}). \quad (A11) $$

The final expression in real space reads

$$ \gamma(\mathbf{r}, \mathbf{r}') = \sum_{\mathbf{k}\lambda} n_{\mathbf{k}\lambda} \phi_{\mathbf{k}\lambda}(\mathbf{r}) \phi_{\mathbf{k}\lambda}^{*}(\mathbf{r}'), \quad (A12) $$

that looks extremely similar to the mean-field (gKS) expression

$$ \gamma^{\text{gKS}}(\mathbf{r}, \mathbf{r}') = \sum_{\mathbf{k}\lambda} f_{\mathbf{k}\lambda} \varphi_{\mathbf{k}\lambda}(\mathbf{r}) \varphi_{\mathbf{k}\lambda}^{*}(\mathbf{r}'), \quad (A13) $$

where $f_{\mathbf{k}\lambda}$ are the Fermi-Dirac occupations and $\varphi_{\mathbf{k}\lambda}(\mathbf{r})$ the mean-field wave functions.

However, there are subtle differences that are much more meaningful.

(1) In gKS DFT, $\varphi_{\mathbf{k}\lambda}(\mathbf{r})$ sums up to the exact electronic density, whereas the natural orbital sums up to the exact density and density matrix.

(2) For the ground state, the spin-summed $f_{\mathbf{k}\lambda}$ are constrained to be 0 or 2, whereas the spin-summed $n_{\mathbf{k}\lambda}$ continuously span the range from 0 to 2 (not proven here).

(3) Since $\int d\mathbf{r}\gamma(\mathbf{r}, \mathbf{r}')$ is normalized to the number of electrons $N_e$, both $f_{\mathbf{k}\lambda}$ and $n_{\mathbf{k}\lambda}$ sum up to $N_e$. However, for insulators, while $\sum_i f_{\mathbf{k}\lambda} = N_e$, for each $\mathbf{k}$ individually, no equivalent exists for the natural occupations $n_{\mathbf{k}\lambda}$, as we show in the Appendix B.

# APPENDIX B: NONINTEGER $N(\mathbf{k})$ VALUES FOR INDIVIDUAL $\mathbf{k}$

The noninteger values reported for the electron count, $N(\mathbf{k})$, are a consequence of the electron correlation effects. In this Appendix, we gain some insights into this result.

## 1. The many-electron wave function in crystals

The basis of Bloch waves is complete; therefore, the real-space $N$-electron wave function can be written as a linear combination of Slater determinants $[(N!)^{-1/2}|\varphi_{\mathbf{k}_\mu i}(\mathbf{r}_1) \dots \varphi_{\mathbf{k}_\nu j}(\mathbf{r}_N)|]$ built using Bloch waves:

$$ \varphi_{\mathbf{k}\lambda}(\mathbf{r}) = \frac{1}{\sqrt{N_k \Omega}} \sum_{\mathbf{G}} e^{i\mathbf{k}\cdot\mathbf{r}} u_{\mathbf{k}\lambda}(\mathbf{r}), \quad (B1) $$

which are usually the ones obtained from a mean-field method (like the ones obtained from a gKS DFT calculation). However, other basis sets can be used to build the Slater determinants, such as the Bloch waves corresponding to the natural orbitals. In this representation, the real-space (spinless) $N$-electron wave function can be written as a configuration interaction (CI) expansion

$$ \Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N) = \frac{1}{\sqrt{N! (\sqrt{N_k \Omega})^N}} \sum_{\mathbf{k}_\mu \dots \mathbf{k}_\nu} \sum_{i \in \Omega_{\mathbf{k}_\mu}} \dots \sum_{j \in \Omega_{\mathbf{k}_\nu}} C_{\mathbf{k}_\mu i \dots \mathbf{k}_\nu j} \times e^{i(\mathbf{k}_\mu \cdot \mathbf{r}_1 + \dots + \mathbf{k}_\nu \cdot \mathbf{r}_N)} u_{\mathbf{k}_\mu i}(\mathbf{r}_1) \dots u_{\mathbf{k}_\nu j}(\mathbf{r}_N), \quad (B2) $$

where the sum $\sum_{\mathbf{k}_\mu \dots \mathbf{k}_\nu}$ runs over all $\mathbf{k}$ points in the first Brillouin zone, and the $C_{\mathbf{k}_\mu i \dots \mathbf{k}_\nu j}$ are the expansion coefficients (that are adequately adjusted to ensure that $\Psi$ preserves the correct symmetries). Let us highlight that in Eq. (B2) the Hartree product of Bloch's waves contains waves that belong to different $\mathbf{k}$ points (i.e., $\Psi$ is the many-body wave function of the supercell).

The Born-von Kármán periodic conditions imposed to the many-electron wave function [71] state that $\Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N) = \Psi(\mathbf{r}_1 + \mathbf{T}, \mathbf{r}_2 + \mathbf{T}, \dots, \mathbf{r}_N + \mathbf{T})$ for $\mathbf{T} = n_1 N_1 \mathbf{a}_1 + n_2 N_2 \mathbf{a}_2 + n_3 N_3 \mathbf{a}_3$ (with $N_k = N_1 N_2 N_3$). As a consequence, $e^{i(\mathbf{k}_\mu + \dots + \mathbf{k}_\nu) \cdot \mathbf{T}} = e^{i\mathbf{k} \cdot \mathbf{T}} = 1$ and

$$ \mathbf{K} = \mathbf{k}_\mu + \dots + \mathbf{k}_\nu \quad (B3) $$

$$ = \frac{\chi_1}{N_1} \mathbf{b}_1 + \frac{\chi_2}{N_2} \mathbf{b}_2 + \frac{\chi_3}{N_3} \mathbf{b}_3, \quad (B4) $$
This condition is fulfilled if and only if $(\mathbf{k} - \mathbf{k}')$ belongs to the reciprocal lattice. And as $\mathbf{k}$ and $\mathbf{k}'$ both belong to the first Brillouin zone, the only possible reciprocal lattice vector the difference can match is $\mathbf{0}$.

As a consequence, one can insert the Kronecker sign $\delta_{\mathbf{k}\mathbf{k}'}$ in Eq. (A1)

$$ \gamma(\mathbf{r}, \mathbf{r}') = \frac{1}{N_k \Omega} \sum_{\mathbf{k}GG'} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} \gamma_k(\mathbf{G}, \mathbf{G}') e^{-i(\mathbf{k}+\mathbf{G}')\cdot\mathbf{r}'} \quad (A4) $$

and obtain the desired expression that shows that $\gamma$ is block diagonal with respect to $\mathbf{k}$.

## 2. Natural orbitals and occupations

Now for each discrete value of $\mathbf{k}$, one can diagonalize $\gamma_k$:

$$ \sum_{\mathbf{G}'} \gamma_k(\mathbf{G}, \mathbf{G}') \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}') = n_{\mathbf{k}\lambda} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}), \quad (A5) $$

where $\tilde{u}_{\mathbf{k}\lambda}(\mathbf{G})$ are the eigenvectors and $n_{\mathbf{k}\lambda}$ the eigenvalues, indexed with $\lambda$.

By construction,

$$ \gamma(\mathbf{r}, \mathbf{r}') = \gamma(\mathbf{r}', \mathbf{r})^* \quad (A6) $$

implies

$$ \gamma_k(\mathbf{G}, \mathbf{G}') = \gamma_k(\mathbf{G}', \mathbf{G})^*. \quad (A7) $$

$\gamma_k(\mathbf{G}, \mathbf{G}')$ is a Hermitian matrix and thus the $n_{\mathbf{k}\lambda}$ are real valued and the $\tilde{u}_{\mathbf{k}\lambda}$ form a unitary matrix

Using this decomposition,

$$ \gamma_k(\mathbf{G}, \mathbf{G}') = \sum_{\lambda} n_{\mathbf{k}\lambda} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}) \tilde{u}_{\mathbf{k}\lambda}^*(\mathbf{G}'). \quad (A8) $$

Inserting this in Eq. (A4), we obtain

$$ \gamma(\mathbf{r}, \mathbf{r}') = \frac{1}{N_k \Omega} \sum_{\mathbf{k}\lambda} \sum_{\mathbf{G}} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}) \times n_{\mathbf{k}\lambda} \sum_{\mathbf{G}'} \tilde{u}_{\mathbf{k}\lambda}^*(\mathbf{G}') e^{-i(\mathbf{k}+\mathbf{G}')\cdot\mathbf{r}'} \quad (A9) $$

Let us introduce the natural orbital in real space $\phi_{\mathbf{k}\lambda}(\mathbf{r})$ that has a Bloch wave form [70]:

$$ \phi_{\mathbf{k}\lambda}(\mathbf{r}) = \frac{1}{\sqrt{N_k \Omega}} \sum_{\mathbf{G}} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{G}) \quad (A10) $$

$$ = \frac{1}{\sqrt{N_k \Omega}} e^{i\mathbf{k}\cdot\mathbf{r}} \tilde{u}_{\mathbf{k}\lambda}(\mathbf{r}). \quad (A11) $$

The final expression in real space reads

$$ \gamma(\mathbf{r}, \mathbf{r}') = \sum_{\mathbf{k}\lambda} n_{\mathbf{k}\lambda} \phi_{\mathbf{k}\lambda}(\mathbf{r}) \phi_{\mathbf{k}\lambda}^*(\mathbf{r}'), \quad (A12) $$

that looks extremely similar to the mean-field (gKS) expression

$$ \gamma^{\text{gKS}}(\mathbf{r}, \mathbf{r}') = \sum_{\mathbf{k}\lambda} f_{\mathbf{k}\lambda} \varphi_{\mathbf{k}\lambda}(\mathbf{r}) \varphi_{\mathbf{k}\lambda}^*(\mathbf{r}'), \quad (A13) $$

where $f_{\mathbf{k}\lambda}$ are the Fermi-Dirac occupations and $\varphi_{\mathbf{k}\lambda}(\mathbf{r})$ the mean-field wave functions.

However, there are subtle differences that are much more meaningful.

(1) In gKS DFT, $\varphi_{\mathbf{k}\lambda}(\mathbf{r})$ sums up to the exact electronic density, whereas the natural orbital sums up to the exact density and density matrix.

(2) For the ground state, the spin-summed $f_{\mathbf{k}\lambda}$ are constrained to be 0 or 2, whereas the spin-summed $n_{\mathbf{k}\lambda}$ continuously span the range from 0 to 2 (not proven here).

(3) Since $\int d\mathbf{r}\gamma(\mathbf{r}, \mathbf{r}')$ is normalized to the number of electrons $N_e$, both $f_{\mathbf{k}\lambda}$ and $n_{\mathbf{k}\lambda}$ sum up to $N_e$. However, for insulators, while $\sum_i f_{\mathbf{k}\lambda} = N_e$, for each $\mathbf{k}$ individually, no equivalent exists for the natural occupations $n_{\mathbf{k}\lambda}$, as we show in the Appendix B.

# APPENDIX B: NONINTEGER $N(\mathbf{k})$ VALUES FOR INDIVIDUAL $\mathbf{k}$

The noninteger values reported for the electron count, $N(\mathbf{k})$, are a consequence of the electron correlation effects. In this Appendix, we gain some insights into this result.

## 1. The many-electron wave function in crystals

The basis of Bloch waves is complete; therefore, the real-space $N$-electron wave function can be written as a linear combination of Slater determinants $[(N!)^{-1/2}|\varphi_{\mathbf{k}_\mu i}(\mathbf{r}_1) \dots \varphi_{\mathbf{k}_\nu j}(\mathbf{r}_N)|]$ built using Bloch waves:

$$ \varphi_{\mathbf{k}\lambda}(\mathbf{r}) = \frac{1}{\sqrt{N_k \Omega}} \sum_{\mathbf{G}} e^{i\mathbf{k}\cdot\mathbf{r}} u_{\mathbf{k}\lambda}(\mathbf{r}), \quad (B1) $$

which are usually the ones obtained from a mean-field method (like the ones obtained from a gKS DFT calculation). However, other basis sets can be used to build the Slater determinants, such as the Bloch waves corresponding to the natural orbitals. In this representation, the real-space (spinless) $N$-electron wave function can be written as a configuration interaction (CI) expansion

$$ \Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N) = \frac{1}{\sqrt{N! (\sqrt{N_k \Omega})^N}} \sum_{\mathbf{k}_\mu \dots \mathbf{k}_\nu} \sum_{i \in \Omega_{\mathbf{k}_\mu}} \dots \sum_{j \in \Omega_{\mathbf{k}_\nu}} C_{\mathbf{k}_\mu i \dots \mathbf{k}_\nu j} \times e^{i(\mathbf{k}_\mu \cdot \mathbf{r}_1 + \dots + \mathbf{k}_\nu \cdot \mathbf{r}_N)} u_{\mathbf{k}_\mu i}(\mathbf{r}_1) \dots u_{\mathbf{k}_\nu j}(\mathbf{r}_N), \quad (B2) $$

where the sum $\sum_{\mathbf{k}_\mu \dots \mathbf{k}_\nu}$ runs over all $\mathbf{k}$ points in the first Brillouin zone, and the $C_{\mathbf{k}_\mu i \dots \mathbf{k}_\nu j}$ are the expansion coefficients (that are adequately adjusted to ensure that $\Psi$ preserves the correct symmetries). Let us highlight that in Eq. (B2) the Hartree product of Bloch's waves contains waves that belong to different $\mathbf{k}$ points (i.e., $\Psi$ is the many-body wave function of the supercell).

The Born-von Kármán periodic conditions imposed to the many-electron wave function [71] state that $\Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N) = \Psi(\mathbf{r}_1 + \mathbf{T}, \mathbf{r}_2 + \mathbf{T}, \dots, \mathbf{r}_N + \mathbf{T})$ for $\mathbf{T} = n_1 N_1 \mathbf{a}_1 + n_2 N_2 \mathbf{a}_2 + n_3 N_3 \mathbf{a}_3$ (with $N_k = N_1 N_2 N_3$). As a consequence, $e^{i(\mathbf{k}_\mu + \dots + \mathbf{k}_\nu) \cdot \mathbf{T}} = e^{i\mathbf{k} \cdot \mathbf{T}} = 1$ and

$$ \mathbf{K} = \mathbf{k}_\mu + \dots + \mathbf{k}_\nu \quad (B3) $$

$$ = \frac{\chi_1}{N_1} \mathbf{b}_1 + \frac{\chi_2}{N_2} \mathbf{b}_2 + \frac{\chi_3}{N_3} \mathbf{b}_3, \quad (B4) $$
where $\chi_n \in \mathbb{Z}$. The many-electron wave functions are eigenfunctions of the many-body translation operator $\hat{T}_R$, i.e.,

$$ \hat{T}_R \Psi_K(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N) = \Psi_K(\mathbf{r}_1 + \mathbf{R}, \mathbf{r}_2 + \mathbf{R}, \dots, \mathbf{r}_N + \mathbf{R}) \tag{B5} $$

$$ = e^{i\mathbf{K} \cdot \mathbf{R}} \Psi_K(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N), \tag{B6} $$

with $e^{i\mathbf{K} \cdot \mathbf{R}}$ eigenvalues. Since the many-body Hamiltonian $\hat{H}$ (within the Born-Oppenheimer approximation) commutes with the $\hat{T}_R$ operator [71], the solutions to the many-body Hamiltonian can be taken associated with a given $\mathbf{K}$ value (i.e., $\hat{H}\Psi_K = E\Psi_K$). Actually, the many-electron wave functions associated with $\mathbf{K} + \mathbf{G}$ (with $\mathbf{G} = \chi_1\mathbf{b}_1 + \chi_2\mathbf{b}_2 + \chi_3\mathbf{b}_3$ being a reciprocal lattice vector) also lead to the same eigenvalue $e^{i\mathbf{K} \cdot \mathbf{R}}$ upon application of the $\hat{T}_R$ operator and may contribute to the CI expansion.

## 2. The density matrix and the second-order reduced density matrix in crystals

Let us define the density matrix elements (from the many-body wave function $\Psi$) as

$$ \gamma_{ki j} = \langle \Psi | \hat{b}_{ki}^{\dagger} \hat{b}_{kj} | \Psi \rangle \tag{B7} $$

and the second-order reduced density matrix (2-RDM) elements as

$$ \Gamma_{\mathbf{k}_\mu, \mathbf{k}_\nu, \mathbf{l}} = \frac{1}{2} \langle \Psi | \hat{b}_{\mathbf{k}_\mu, \mathbf{l}}^\dagger \hat{b}_{\mathbf{k}_\nu, \mathbf{l}}^\dagger \hat{b}_{\mathbf{k}_\mu, \mathbf{m}} \hat{b}_{\mathbf{k}_\nu, \mathbf{l}} | \Psi \rangle \tag{B8} $$

whose $\mathbf{k}_n$ values fulfill the condition

$$ \mathbf{k}_{\mu} + \mathbf{k}_{\nu} - \mathbf{k}_{\tau} - \mathbf{k}_{\vartheta} + \mathbf{G} = 0, \tag{B9} $$

which ensures the correct translational symmetry of the 2-RDM, i.e.,

$$ \Gamma(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}'_1, \mathbf{r}'_2) = \Gamma(\mathbf{r}_1 + \mathbf{R}, \mathbf{r}_2 + \mathbf{R}, \mathbf{r}'_1 + \mathbf{R}, \mathbf{r}'_2 + \mathbf{R}), \tag{B10} $$

where $\Gamma(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}'_1, \mathbf{r}'_2)$ is the 2-RDM in space representation (see for example Refs. [72,73] for more details). The 2-RDM contains more information about the system than the density matrix. Indeed, the matrix elements of the density matrix can be obtained from the partial trace of the 2-RDM:

$$ \gamma_{ki j} = \frac{2}{N-1} \sum_{\mathbf{k}'} \sum_{l \in \Omega_{\mathbf{k}'}} \Gamma^{\mathbf{k} \mathbf{j} \mathbf{k}' l}_{\mathbf{k} i, \mathbf{k}' l}, \tag{B11} $$

where $\Omega_{\mathbf{k}'}$ is the subspace formed by all the one-electron wave functions sharing the same $\mathbf{k}'_\nu$ value. For completeness, let us mention that the matrix formed using the $\gamma_{ki j}$ elements is Hermitian and upon diagonalization produces the occupation numbers $n_{\mathbf{k}\lambda}$ discussed in Appendix A.

## 3. The APSG ansatz for crystals

Here we prove in the specific case of an antisymmetrized product of strongly orthogonal geminals (APSG) ansatz [74] for the many-electron wave function that electron count transfer can occur across $\mathbf{k}$ points due to electronic correlation effects. If it is true for this subclass of wave functions, then the statement also holds for the exact wave function.

The APSG ansatz for the many-electron wave function with spin, where an even number of electrons present in the system is assumed (as we employed throughout this paper), reads

$$ \Psi_{\mathbf{K}}^{\text{APSG}}(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N) = \hat{A} \prod_{P=1}^{N/2} \psi_P(\mathbf{x}_{2P-1}, \mathbf{x}_{2P}), \tag{B12} $$

where $\mathbf{x} = (\mathbf{r}, \sigma)$ is the spatial and spin coordinate with $\sigma = \alpha, \beta$ referring to the spin index, $\hat{A}$ stands for the antisymmetrizer responsible for intergeminal permutations of electron coordinates, and the geminal wave functions $\psi_P(\mathbf{x}_{2P-1}, \mathbf{x}_{2P})$ are wave functions containing one $\alpha$ and one $\beta$ electron. Because of this, the geminal wave function is a two-electron wave function; the sum of the occupation numbers for each spin channel must be

$$ \sum_{[\mathbf{k}, i \in \Omega_{\mathbf{k}}] \in P} n_{\mathbf{k} i} = 1 \tag{B13} $$

with $[\mathbf{k}, i \in \Omega_{\mathbf{k}}] \in P$ indicating that the $i$th Bloch wave belonging to the $k$th $\mathbf{k}$ point (i.e., the Bloch wave's natural orbital $\phi_{\mathbf{k}i}$) is one of the Bloch waves used in the construction of the $P$th geminal. The $P$th geminal wave function written in terms of the natural orbital Bloch waves reads as

$$ \begin{aligned} & \psi_P(\mathbf{x}_{2P-1}, \mathbf{x}_{2P}) \\ &= 2^{-1/2} \sum_{[\mathbf{k}, i \in \Omega_{\mathbf{k}}] \in P} c_{\mathbf{k}i} [\phi_{\mathbf{k}i}(\mathbf{r}_{2P-1}) \alpha_{2P-1} \phi_{\mathbf{k}i}^*(\mathbf{r}_{2P}) \beta_{2P} \\ & \qquad - \phi_{\mathbf{k}i}(\mathbf{r}_{2P}) \alpha_{2P} \phi_{\mathbf{k}i}^*(\mathbf{r}_{2P-1}) \beta_{2P-1}], \end{aligned} \tag{B14} $$

where the time-reversal symmetry is being employed to relate the degenerated Bloch waves containing electrons with opposite spin [75,76] forming a Kramers pair (i.e., $\phi_{\mathbf{k}i}$ and $\phi_{\mathbf{k}i}^*$ form a Kramers pair). Let us remark that the presence of complex-conjugated Bloch waves (natural orbitals) is related to states filled by $\beta$ spin electrons; this choice is completely arbitrary. Also, the time-reversal symmetry imposed on the many-electron wave function leads to $\Psi_{\mathbf{K}=0}$ for spin-compensated systems.

Since the $\{\psi_P\}$ geminal wave functions are built with the strong orthonormality requirement, i.e., the condition that $\forall_{P \neq Q} \int d\mathbf{x}_2 \psi_P(\mathbf{x}_1, \mathbf{x}_2) \psi_Q(\mathbf{x}'_1, \mathbf{x}_2) = \delta_{PQ}$, then the Bloch waves' natural orbitals are present in only one geminal wave function. When all the $\psi_P$ are built containing only one Kramers pair as

$$ \psi_P(\mathbf{x}_{2P-1}, \mathbf{x}_{2P}) = 2^{-1/2} [\phi_{\mathbf{k}i}(\mathbf{r}_{2P-1}) \alpha_{2P-1} \phi_{\mathbf{k}i}^*(\mathbf{r}_{2P}) \beta_{2P} \\ - \phi_{\mathbf{k}i}(\mathbf{r}_{2P}) \alpha_{2P} \phi_{\mathbf{k}i}^*(\mathbf{r}_{2P-1}) \beta_{2P-1}], \tag{B15} $$

the many-body wave function ($\Psi_{\mathbf{K}}$) defined in Eq. (B12) corresponds to a single Slater determinant; thus, the Hartree-Fock approximation is recovered, where the natural orbital basis and the so-called canonical orbitals (the mean-field ones) coincide.

The structure of the APSG wave function allows us to express the total energy in terms of the natural orbitals, the occupation numbers, and some undetermined phases. Then,
where $\chi_n \in \mathbb{Z}$. The many-electron wave functions are eigenfunctions of the many-body translation operator $\hat{T}_R$, i.e.,

$$ \hat{T}_R \Psi_K(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N) = \Psi_K(\mathbf{r}_1 + \mathbf{R}, \mathbf{r}_2 + \mathbf{R}, \dots, \mathbf{r}_N + \mathbf{R}) \tag{B5} $$

$$ = e^{i\mathbf{K} \cdot \mathbf{R}} \Psi_K(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N), \tag{B6} $$

with $e^{i\mathbf{K} \cdot \mathbf{R}}$ eigenvalues. Since the many-body Hamiltonian $\hat{H}$ (within the Born-Oppenheimer approximation) commutes with the $\hat{T}_R$ operator [71], the solutions to the many-body Hamiltonian can be taken associated with a given $\mathbf{K}$ value (i.e., $\hat{H}\Psi_K = E\Psi_K$). Actually, the many-electron wave functions associated with $\mathbf{K} + \mathbf{G}$ (with $\mathbf{G} = \chi_1\mathbf{b}_1 + \chi_2\mathbf{b}_2 + \chi_3\mathbf{b}_3$ being a reciprocal lattice vector) also lead to the same eigenvalue $e^{i\mathbf{K} \cdot \mathbf{R}}$ upon application of the $\hat{T}_R$ operator and may contribute to the CI expansion.

## 2. The density matrix and the second-order reduced density matrix in crystals

Let us define the density matrix elements (from the many-body wave function $\Psi$) as

$$ \gamma_{ki j} = \langle \Psi | \hat{b}_{ki}^{\dagger} \hat{b}_{kj} | \Psi \rangle \tag{B7} $$

and the second-order reduced density matrix (2-RDM) elements as

$$ \Gamma_{\mathbf{k}_\mu, \mathbf{k}_\nu, \mathbf{l}} = \frac{1}{2} \langle \Psi | \hat{b}_{\mathbf{k}_\mu, \mathbf{l}}^\dagger \hat{b}_{\mathbf{k}_\nu, \mathbf{l}}^\dagger \hat{b}_{\mathbf{k}_\mu, \mathbf{m}} \hat{b}_{\mathbf{k}_\nu, \mathbf{l}} | \Psi \rangle \tag{B8} $$

whose $\mathbf{k}_n$ values fulfill the condition

$$ \mathbf{k}_{\mu} + \mathbf{k}_{\nu} - \mathbf{k}_{\tau} - \mathbf{k}_{\vartheta} + \mathbf{G} = 0, \tag{B9} $$

which ensures the correct translational symmetry of the 2-RDM, i.e.,

$$ \Gamma(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}'_1, \mathbf{r}'_2) = \Gamma(\mathbf{r}_1 + \mathbf{R}, \mathbf{r}_2 + \mathbf{R}, \mathbf{r}'_1 + \mathbf{R}, \mathbf{r}'_2 + \mathbf{R}), \tag{B10} $$

where $\Gamma(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}'_1, \mathbf{r}'_2)$ is the 2-RDM in space representation (see for example Refs. [72,73] for more details). The 2-RDM contains more information about the system than the density matrix. Indeed, the matrix elements of the density matrix can be obtained from the partial trace of the 2-RDM:

$$ \gamma_{ki j} = \frac{2}{N-1} \sum_{\mathbf{k}'} \sum_{l \in \Omega_{\mathbf{k}'}} \Gamma^{\mathbf{k} \mathbf{j} \mathbf{k}' l}_{\mathbf{k} i, \mathbf{k}' l}, \tag{B11} $$

where $\Omega_{\mathbf{k}'}$ is the subspace formed by all the one-electron wave functions sharing the same $\mathbf{k}'_\nu$ value. For completeness, let us mention that the matrix formed using the $\gamma_{ki j}$ elements is Hermitian and upon diagonalization produces the occupation numbers $n_{\mathbf{k}\lambda}$ discussed in Appendix A.

## 3. The APSG ansatz for crystals

Here we prove in the specific case of an antisymmetrized product of strongly orthogonal geminals (APSG) ansatz [74] for the many-electron wave function that electron count transfer can occur across $\mathbf{k}$ points due to electronic correlation effects. If it is true for this subclass of wave functions, then the statement also holds for the exact wave function.

The APSG ansatz for the many-electron wave function with spin, where an even number of electrons present in the system is assumed (as we employed throughout this paper), reads

$$ \Psi_{\mathbf{K}}^{\text{APSG}}(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N) = \hat{A} \prod_{P=1}^{N/2} \psi_P(\mathbf{x}_{2P-1}, \mathbf{x}_{2P}), \tag{B12} $$

where $\mathbf{x} = (\mathbf{r}, \sigma)$ is the spatial and spin coordinate with $\sigma = \alpha, \beta$ referring to the spin index, $\hat{A}$ stands for the antisymmetrizer responsible for intergeminal permutations of electron coordinates, and the geminal wave functions $\psi_P(\mathbf{x}_{2P-1}, \mathbf{x}_{2P})$ are wave functions containing one $\alpha$ and one $\beta$ electron. Because of this, the geminal wave function is a two-electron wave function; the sum of the occupation numbers for each spin channel must be

$$ \sum_{[\mathbf{k}, i \in \Omega_{\mathbf{k}}] \in P} n_{\mathbf{k} i} = 1 \tag{B13} $$

with $[\mathbf{k}, i \in \Omega_{\mathbf{k}}] \in P$ indicating that the $i$th Bloch wave belonging to the $\mathbf{k}$th $\mathbf{k}$ point (i.e., the Bloch wave's natural orbital $\phi_{\mathbf{k}i}$) is one of the Bloch waves used in the construction of the $P$th geminal. The $P$th geminal wave function written in terms of the natural orbital Bloch waves reads as

$$ \begin{aligned} & \psi_P(\mathbf{x}_{2P-1}, \mathbf{x}_{2P}) \\ &= 2^{-1/2} \sum_{[\mathbf{k}, i \in \Omega_{\mathbf{k}}] \in P} c_{\mathbf{k}i} [\phi_{\mathbf{k}i}(\mathbf{r}_{2P-1}) \alpha_{2P-1} \phi_{\mathbf{k}i}^*(\mathbf{r}_{2P}) \beta_{2P} \\ & \qquad - \phi_{\mathbf{k}i}(\mathbf{r}_{2P}) \alpha_{2P} \phi_{\mathbf{k}i}^*(\mathbf{r}_{2P-1}) \beta_{2P-1}], \end{aligned} \tag{B14} $$

where the time-reversal symmetry is being employed to relate the degenerated Bloch waves containing electrons with opposite spin [75,76] forming a Kramers pair (i.e., $\phi_{\mathbf{k}i}$ and $\phi_{\mathbf{k}i}^*$ form a Kramers pair). Let us remark that the presence of complex-conjugated Bloch waves (natural orbitals) is related to states filled by $\beta$ spin electrons; this choice is completely arbitrary. Also, the time-reversal symmetry imposed on the many-electron wave function leads to $\Psi_{\mathbf{K}=0}$ for spin-compensated systems.

Since the $\{\psi_P\}$ geminal wave functions are built with the strong orthonormality requirement, i.e., the condition that $\forall_{P \neq Q} \int d\mathbf{x}_2 \psi_P(\mathbf{x}_1, \mathbf{x}_2) \psi_Q(\mathbf{x}'_1, \mathbf{x}_2) = \delta_{PQ}$, then the Bloch waves' natural orbitals are present in only one geminal wave function. When all the $\psi_P$ are built containing only one Kramers pair as

$$ \psi_P(\mathbf{x}_{2P-1}, \mathbf{x}_{2P}) = 2^{-1/2} [\phi_{\mathbf{k}i}(\mathbf{r}_{2P-1}) \alpha_{2P-1} \phi_{\mathbf{k}i}^*(\mathbf{r}_{2P}) \beta_{2P} \\ - \phi_{\mathbf{k}i}(\mathbf{r}_{2P}) \alpha_{2P} \phi_{\mathbf{k}i}^*(\mathbf{r}_{2P-1}) \beta_{2P-1}], \tag{B15} $$

the many-body wave function ($\Psi_{\mathbf{K}}$) defined in Eq. (B12) corresponds to a single Slater determinant; thus, the Hartree-Fock approximation is recovered, where the natural orbital basis and the so-called canonical orbitals (the mean-field ones) coincide.

The structure of the APSG wave function allows us to express the total energy in terms of the natural orbitals, the occupation numbers, and some undetermined phases. Then,
the total APSG energy takes the following form:

$$ 
\begin{aligned}
E^{\text{APSG}}[\{f_{\mathbf{k}_{\mu i}}, \{n_{\mathbf{k}_{\mu i}}, \{\phi_{\mathbf{k}_{\mu i}}\}\}] \\
= 2 \sum_{\mathbf{k}_{\mu}} \sum_{i \in \Omega_{\mathbf{k}_{\mu}}} n_{\mathbf{k}_{\mu} i} \langle \mathbf{k}_{\mu} i | \hat{h} | \mathbf{k}_{\mu} i \rangle \\
+ \sum_{P}^{N/2} \sum_{\substack{\mathbf{k}_{\mu}, i \in \Omega_{\mathbf{k}_{\mu}} \\ [\mathbf{k}_{\mu}, j \in \Omega_{\mathbf{k}_{\mu}}]}} \zeta_{\mathbf{k}_{\mu} i} \zeta_{\mathbf{k}_{\nu} j} \sqrt{n_{\mathbf{k}_{\mu} i} n_{\mathbf{k}_{\mu} j}} \langle \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j | \mathbf{k}_{\nu} j \mathbf{k}_{\mu} i \rangle \\
+ \sum_{P \neq Q} \sum_{\substack{\mathbf{k}_{\mu}, i \in \Omega_{\mathbf{k}_{\mu}} \\ [\mathbf{k}_{\mu}, j \in \Omega_{\mathbf{k}_{\mu}}]}} n_{\mathbf{k}_{\mu} i} n_{\mathbf{k}_{\nu} j} (2 \langle \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j | \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j \rangle \\
- \langle \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j | \mathbf{k}_{\nu} j \mathbf{k}_{\mu} i \rangle),
\end{aligned}
\quad (B16) $$

where we have employed the known condition for APSG wave functions that allows us to express the CI coefficients in terms of occupation numbers ($c_{\mathbf{k}_{\mu}}^2 = n_{\mathbf{k}_{\mu}}$ or $c_{\mathbf{k}_{\mu}} = \zeta_{\mathbf{k}_{\mu}} \sqrt{n_{\mathbf{k}_{\mu}}}$), the phases $\zeta_{\mathbf{k}_{\mu}}^2 = 1$, $\hat{h}$ refers to all one-body operators of the electronic Hamiltonian (i.e., the kinetic energy and the interaction with the external potential), and $\langle \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j | \mathbf{k}_{\tau} k_{\mathbf{x}} l \rangle$ are the usual two-electron integrals. Notice that the energy minimization procedure implies optimization of the occupation numbers, natural orbitals, and phases. Since this wave function can be entirely written in terms of the natural orbitals and occupation numbers, it has been widely used in the context of reduced density matrix functional theory to propose energy functionals [77–85].

The energy contribution in the second line of Eq. (B16) is coming from the geminal wave functions and describes intergeminal interactions (i.e., the ones among the two electrons belonging to the geminal wave function). On the other hand, energy contribution in the third line of Eq. (B16) describes intergeminal interactions, which are taken at the mean-field level (i.e., as Hartree-Fock interactions). Also, let us highlight that imposing the correct translation symmetry to $\Psi_K^{\text{APSG}}$ automatically enforces the correct symmetry in the 2-RDM elements, making the 2-RDM elements fulfill the condition presented in Eq. (B9).

Since the natural orbitals belonging to different $\Omega_k$ subspaces (e.g., $\mathbf{k}$ and $\mathbf{k}'$) can be employed in the construction of the geminal wave functions [see Eq. (B14)] the constraint given by Eq. (B9) takes the following form for the geminal's wave function:

$$ \mathbf{k}' - \mathbf{k}' - \mathbf{k} + \mathbf{k} = \mathbf{0} = \mathbf{G}, \quad (B17) $$

which ensures the correct translation symmetry in $\Psi_K^{\text{APSG}}$, the 2-RDM matrix elements, and the density matrix. Notice that the complex conjugation associated with the time-reversal symmetry of the Kramers pairs was employed.

The occupation numbers of the natural orbitals that belong to the geminal are optimized under the constraint given in Eq. (B13) during the energy minimization procedure (recalling that the CI coefficients can be written as $c_{\mathbf{k}_{\mu}} = \zeta_{\mathbf{k}_{\mu}} \sqrt{n_{\mathbf{k}_{\mu}}}$ for the $\Psi_K^{\text{APSG}}$ ansatz). Hence, the coupling of natural orbitals belonging to different $\mathbf{k}$ points is allowed; thus, a reorganization of electrons among $\mathbf{k}$ points can take place during the energy minimization procedure. Moreover, it is known [86] that the

occupation numbers are not likely to become zero; then, the reorganization of electrons among $\mathbf{k}$ points is not forbidden.

In summary, the $\Psi_K^{\text{APSG}}$ ansatz is a valid approximation to the many-electron wave function that permits us to illustrate the reasons leading to the reorganization of electrons among $\mathbf{k}$ points. Obviously, a more general valid CI expansion ansatz (or the exact full-CI expansion) could also lead to a reorganization of the electrons among $\mathbf{k}$ points since the $\Psi_K^{\text{APSG}}$ ansatz is a particular case of the exact $\Psi_K$, where the electron pairs do not interact. Then, let us conclude that the reorganization of electrons among $\mathbf{k}$ points that lead to noninteger $N_e$ values for each $\mathbf{k}$ value is purely a consequence of the electronic correlation effects. In this paper, the electronic correlation effects are captured with the $\gamma^{\text{GW}}$ approximation, which produces the reorganization of electrons among $\mathbf{k}$ points. And, this leads to the noninteger $N_e$ values obtained for each $\mathbf{k}$ point that were used to compute the $\Delta N(\mathbf{k})$ values presented in Fig. 3.

In the next section, we present an example based on the Si crystal where the reorganization of electrons among $\mathbf{k}$ points is allowed using a $\Psi_K^{\text{APSG}}$ ansatz.

### 4. Example of an allowed electronic density reorganization among the k points in the Si crystal

For a working example, let us take the Si crystal computed excluding all the core states (i.e., using a pseudopotential and retaining only eight electrons per unit cell). At the Hartree-Fock (or gKS DFT) level, eight states forming four Kramers pairs are occupied for each $\mathbf{k}$-point value. From the band structure, it is easy to recognize that the highest (in terms of energy) occupied state with an $\alpha$ electron is localized at the $\Gamma$ point $[\mathbf{k} = (0, 0, 0) = \mathbf{k}_{\Gamma}]$. On the other hand, the lowest (in terms of energy) unoccupied state for the electrons with $\alpha$ spin belongs to the $X$ point $[\mathbf{k} = (0.5, 0.5, 0) = \mathbf{k}_{X}]$. Let us label these states as $\varphi_{\mathbf{k}_{\Gamma}4}$ and $\varphi_{\mathbf{k}_{X}5}$, respectively. The energy difference between the $\varphi_{\mathbf{k}_{\Gamma}4}$ and $\varphi_{\mathbf{k}_{X}5}$ states is small (the experimental value is approximately 1 eV), which leads to the small indirect band gap obtained for this system.

In the following, let us organize in ascending order in terms of energy all the mean-field Bloch waves for the whole system (i.e., of the supercell). And, as it is usually done in the search for the optimal $\Psi_K^{\text{APSG}}$, let us write the initial guess for the APSG ansatz in terms of the mean-field Bloch waves. But, let us search for a particular $\Psi_K^{\text{APSG}}$ ansatz where all the geminal wave functions contain only one Kramers pair (i.e., are treated at the Hartree-Fock level) except for the last geminal [the $P = N/2$ in (B14)] that is built coupling the Bloch waves $\varphi_{\mathbf{k}_{\Gamma}4}$ and $\varphi_{\mathbf{k}_{X}5}$, i.e.,

$$ \begin{aligned}
\psi_{N/2}(\mathbf{x}_{N-1}, \mathbf{x}_N) ={}& 2^{-1/2}[c_{\mathbf{k}_{\Gamma}4}(\varphi_{\mathbf{k}_{\Gamma}4}(\mathbf{r}_{N-1})\alpha_{N-1}\varphi_{\mathbf{k}_{\Gamma}4}^*(\mathbf{r}_N)\beta_N \\
& + \varphi_{\mathbf{k}_{\Gamma}4}(\mathbf{r}_N)\alpha_N\varphi_{\mathbf{k}_{\Gamma}4}^*(\mathbf{r}_{N-1})\beta_{N-1}) \\
& + c_{\mathbf{k}_{X}5}(\varphi_{\mathbf{k}_{X}5}(\mathbf{r}_{N-1})\alpha_{N-1}\varphi_{\mathbf{k}_{X}5}^*(\mathbf{r}_N)\beta_N \\
& + \varphi_{\mathbf{k}_{X}5}(\mathbf{r}_N)\alpha_N\varphi_{\mathbf{k}_{X}5}^*(\mathbf{r}_{N-1})\beta_{N-1})],
\end{aligned} \quad (B18) $$

with $c_{\mathbf{k}_{\Gamma}4} = \zeta_{\mathbf{k}_{\Gamma}4} \sqrt{n_{\mathbf{k}_{\Gamma}4}}$ and $c_{\mathbf{k}_{X}5} = \zeta_{\mathbf{k}_{X}5} \sqrt{n_{\mathbf{k}_{X}5}}$ being variational parameters subject to the condition $n_{\mathbf{k}_{\Gamma}4} + n_{\mathbf{k}_{X}5} = 1$ [to fulfill the requirement presented in Eq. (B13)]. This type of geminal approach, where only two states (four considering spin) are present in the geminal wave function, is known
the total APSG energy takes the following form:

$$ 
\begin{aligned}
E^{\text{APSG}}[\{f_{\mathbf{k}_{\mu i}}, \{n_{\mathbf{k}_{\mu i}}, \{\phi_{\mathbf{k}_{\mu i}}\}\}] \\
= 2 \sum_{\mathbf{k}_{\mu}} \sum_{i \in \Omega_{\mathbf{k}_{\mu}}} n_{\mathbf{k}_{\mu} i} \langle \mathbf{k}_{\mu} i | \hat{h} | \mathbf{k}_{\mu} i \rangle \\
+ \sum_{P}^{N/2} \sum_{\substack{\mathbf{k}_{\mu}, i \in \Omega_{\mathbf{k}_{\mu}} \\ [\mathbf{k}_{\mu}, j \in \Omega_{\mathbf{k}_{\mu}}] \in P}} \zeta_{\mathbf{k}_{\mu} i} \zeta_{\mathbf{k}_{\nu} j} \sqrt{n_{\mathbf{k}_{\mu} i} n_{\mathbf{k}_{\mu} j}} \langle \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j | \mathbf{k}_{\nu} j \mathbf{k}_{\mu} i \rangle \\
+ \sum_{P \neq Q} \sum_{\substack{\mathbf{k}_{\mu}, i \in \Omega_{\mathbf{k}_{\mu}} \\ [\mathbf{k}_{\mu}, j \in \Omega_{\mathbf{k}_{\mu}}] \in P}} n_{\mathbf{k}_{\mu} i} n_{\mathbf{k}_{\nu} j} (2 \langle \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j | \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j \rangle \\
- \langle \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j | \mathbf{k}_{\nu} j \mathbf{k}_{\mu} i \rangle),
\end{aligned} 
\quad (B16) $$

where we have employed the known condition for APSG wave functions that allows us to express the CI coefficients in terms of occupation numbers ($c_{\mathbf{k}_{\mu}}^2 = n_{\mathbf{k}_{\mu}}$ or $c_{\mathbf{k}_{\mu}} = \zeta_{\mathbf{k}_{\mu}} \sqrt{n_{\mathbf{k}_{\mu}}}$), the phases $\zeta_{\mathbf{k}_{\mu}}^2 = 1$, $\hat{h}$ refers to all one-body operators of the electronic Hamiltonian (i.e., the kinetic energy and the interaction with the external potential), and $\langle \mathbf{k}_{\mu} i \mathbf{k}_{\nu} j | \mathbf{k}_{\tau} k_{\mathbf{x}} l \rangle$ are the usual two-electron integrals. Notice that the energy minimization procedure implies optimization of the occupation numbers, natural orbitals, and phases. Since this wave function can be entirely written in terms of the natural orbitals and occupation numbers, it has been widely used in the context of reduced density matrix functional theory to propose energy functionals [77–85].

The energy contribution in the second line of Eq. (B16) is coming from the geminal wave functions and describes intergeminal interactions (i.e., the ones among the two electrons belonging to the geminal wave function). On the other hand, energy contribution in the third line of Eq. (B16) describes intergeminal interactions, which are taken at the mean-field level (i.e., as Hartree-Fock interactions). Also, let us highlight that imposing the correct translation symmetry to $\Psi_K^{\text{APSG}}$ automatically enforces the correct symmetry in the 2-RDM elements, making the 2-RDM elements fulfill the condition presented in Eq. (B9).

Since the natural orbitals belonging to different $\Omega_k$ subspaces (e.g., $\mathbf{k}$ and $\mathbf{k}'$) can be employed in the construction of the geminal wave functions [see Eq. (B14)] the constraint given by Eq. (B9) takes the following form for the geminal's wave function:

$$ \mathbf{k}' - \mathbf{k}' - \mathbf{k} + \mathbf{k} = \mathbf{0} = \mathbf{G}, \quad (B17) $$

which ensures the correct translation symmetry in $\Psi_K^{\text{APSG}}$, the 2-RDM matrix elements, and the density matrix. Notice that the complex conjugation associated with the time-reversal symmetry of the Kramers pairs was employed.

The occupation numbers of the natural orbitals that belong to the geminal are optimized under the constraint given in Eq. (B13) during the energy minimization procedure (recalling that the CI coefficients can be written as $c_{\mathbf{k}_{\mu}} = \zeta_{\mathbf{k}_{\mu}} \sqrt{n_{\mathbf{k}_{\mu}}}$ for the $\Psi_K^{\text{APSG}}$ ansatz). Hence, the coupling of natural orbitals belonging to different $\mathbf{k}$ points is allowed; thus, a reorganization of electrons among $\mathbf{k}$ points can take place during the energy minimization procedure. Moreover, it is known [86] that the

occupation numbers are not likely to become zero; then, the reorganization of electrons among $\mathbf{k}$ points is not forbidden.

In summary, the $\Psi_K^{\text{APSG}}$ ansatz is a valid approximation to the many-electron wave function that permits us to illustrate the reasons leading to the reorganization of electrons among $\mathbf{k}$ points. Obviously, a more general valid CI expansion ansatz (or the exact full-CI expansion) could also lead to a reorganization of the electrons among $\mathbf{k}$ points since the $\Psi_K^{\text{APSG}}$ ansatz is a particular case of the exact $\Psi_K$, where the electron pairs do not interact. Then, let us conclude that the reorganization of electrons among $\mathbf{k}$ points that lead to noninteger $N_e$ values for each $\mathbf{k}$ value is purely a consequence of the electronic correlation effects. In this paper, the electronic correlation effects are captured with the $\gamma^{\text{GW}}$ approximation, which produces the reorganization of electrons among $\mathbf{k}$ points. And, this leads to the noninteger $N_e$ values obtained for each $\mathbf{k}$ point that were used to compute the $\Delta N(\mathbf{k})$ values presented in Fig. 3.

In the next section, we present an example based on the Si crystal where the reorganization of electrons among $\mathbf{k}$ points is allowed using a $\Psi_K^{\text{APSG}}$ ansatz.

### 4. Example of an allowed electronic density reorganization among the k points in the Si crystal

For a working example, let us take the Si crystal computed excluding all the core states (i.e., using a pseudopotential and retaining only eight electrons per unit cell). At the Hartree-Fock (or gKS DFT) level, eight states forming four Kramers pairs are occupied for each $\mathbf{k}$-point value. From the band structure, it is easy to recognize that the highest (in terms of energy) occupied state with an $\alpha$ electron is localized at the $\Gamma$ point $[\mathbf{k} = (0, 0, 0) = \mathbf{k}_{\Gamma}]$. On the other hand, the lowest (in terms of energy) unoccupied state for the electrons with $\alpha$ spin belongs to the $X$ point $[\mathbf{k} = (0.5, 0.5, 0) = \mathbf{k}_{\chi}]$. Let us label these states as $\varphi_{\mathbf{k}_{\Gamma}4}$ and $\varphi_{\mathbf{k}_{\chi}5}$, respectively. The energy difference between the $\varphi_{\mathbf{k}_{\Gamma}4}$ and $\varphi_{\mathbf{k}_{\chi}5}$ states is small (the experimental value is approximately 1 eV), which leads to the small indirect band gap obtained for this system.

In the following, let us organize in ascending order in terms of energy all the mean-field Bloch waves for the whole system (i.e., of the supercell). And, as it is usually done in the search for the optimal $\Psi_K^{\text{APSG}}$, let us write the initial guess for the APSG ansatz in terms of the mean-field Bloch waves. But, let us search for a particular $\Psi_K^{\text{APSG}}$ ansatz where all the geminal wave functions contain only one Kramers pair (i.e., are treated at the Hartree-Fock level) except for the last geminal [the $P = N/2$ in (B14)] that is built coupling the Bloch waves $\varphi_{\mathbf{k}_{\Gamma}4}$ and $\varphi_{\mathbf{k}_{\chi}5}$, i.e.,

$$ \begin{aligned} \psi_{N/2}(\mathbf{x}_{N-1}, \mathbf{x}_N) = & 2^{-1/2}[c_{\mathbf{k}_{\Gamma}4}(\varphi_{\mathbf{k}_{\Gamma}4}(\mathbf{r}_{N-1})\alpha_{N-1}\varphi_{\mathbf{k}_{\Gamma}4}^*(\mathbf{r}_N)\beta_N \\ & + \varphi_{\mathbf{k}_{\Gamma}4}(\mathbf{r}_N)\alpha_N\varphi_{\mathbf{k}_{\Gamma}4}^*(\mathbf{r}_{N-1})\beta_{N-1}) \\ & + c_{\mathbf{k}_{\chi}5}(\varphi_{\mathbf{k}_{\chi}5}(\mathbf{r}_{N-1})\alpha_{N-1}\varphi_{\mathbf{k}_{\chi}5}^*(\mathbf{r}_N)\beta_N \\ & + \varphi_{\mathbf{k}_{\chi}5}(\mathbf{r}_N)\alpha_N\varphi_{\mathbf{k}_{\chi}5}^*(\mathbf{r}_{N-1})\beta_{N-1})], \end{aligned} \quad (B18) $$

with $c_{\mathbf{k}_{\Gamma}4} = \zeta_{\mathbf{k}_{\Gamma}4} \sqrt{n_{\mathbf{k}_{\Gamma}4}}$ and $c_{\mathbf{k}_{\chi}5} = \zeta_{\mathbf{k}_{\chi}5} \sqrt{n_{\mathbf{k}_{\chi}5}}$ being variational parameters subject to the condition $n_{\mathbf{k}_{\Gamma}4} + n_{\mathbf{k}_{\chi}5} = 1$ [to fulfill the requirement presented in Eq. (B13)]. This type of geminal approach, where only two states (four considering spin) are present in the geminal wave function, is known
as a perfect pairing approach. Next, let us assume that the mean-field Bloch waves $\varphi_{k\Gamma 4}$ and $\varphi_{k\chi 5}$ coincide with the optimal natural orbitals in order to skip the orbital optimization procedure.

Next, let us focus on the energy contribution arising from the N/2 geminal to the second term in the right-hand side of the APSG energy [see Eq. (B16)]:

$$n_{k\Gamma 4}(k_{\Gamma}4k_{\Gamma}4|k_{\Gamma}4k_{\Gamma}4) + n_{k\chi 5}(k_{\chi}5k_{\chi}5|k_{\chi}5k_{\chi}5) + 2\zeta_{k\Gamma 4}\zeta_{k\chi 5}\sqrt{n_{k\Gamma 4}n_{k\chi 5}}(k_{\Gamma}4k_{\chi}5|k_{\chi}5k_{\Gamma}4). \quad (B19)$$

We set the usual approximation for (fixing) the phases (i.e., $\zeta_{k\Gamma 4}\zeta_{k\chi 5} = -1$ [77,78,80,81]) for the interaction among the states above and below the Fermi level and allow all two-electron integrals to be equal (which can occur in the extreme case when degenerated states are involved). The occupation numbers that would minimize the energy contribution are $n_{k\Gamma 4} = n_{k\chi 5} = 1/2$, illustrating that a reorganization of electrons occurs among the $\Gamma$ and $X$ k points. In the Si crystal, the mean-field Bloch waves are not completely degenerate in terms of energy and do not correspond to the optimal natural

orbitals; then, the actual optimal occupation numbers differ from 1/2. But, they also differ from the initial values at the mean-field level, where $n_{k\Gamma 4} = 1$ and $n_{k\chi 5} = 0$. Moreover, beyond the perfect pairing approach, the coupling of states to form a geminal wave function can include states belonging to other k points. Since the geminal wave functions are built with states for the $\alpha$ and the $\beta$ electron; the state for the $\alpha$ electrons is associated with a $k''$ k point while the state for the $\beta$ is related to a $-k''$ k point. The Hartree product in Eq. (B14) conserves the $K = k'' - k'' = 0$ value, which could lead to further reorganization of electrons among different k points beyond the perfect pairing approach. Thus, for example, the coupling of states belonging to $\Gamma$, $X$, $\Delta$, etc., is allowed in the Si crystal. Actually, the coupling of all k-point values in the first Brillouin zone is valid to build geminal wave functions.

Finally, let us remark that this example is based on a valid approximation to the many-electron wave function (i.e., a $\Psi_k^{\text{APSG}}$ ansatz), where we illustrate that the reorganization of electrons among k points is purely a consequence of the electronic correlation effects.

[1] W. Ku and A. G. Eguiluz, Band-Gap Problem in Semiconductors Revisited: Effects of Core States and Many-Body Self-Consistency, Phys. Rev. Lett. **89**, 126401 (2002).

[2] H. Cao, Z. Yu, P. Lu, and L.-W. Wang, Fully converged plane-wave-based self-consistent GW calculations of periodic solids, Phys. Rev. B **95**, 035139 (2017).

[3] M. Grumet, P. Liu, M. Kaltak, J. Klimeš, and G. Kresse, Beyond the quasiparticle approximation: Fully self-consistent GW calculations, Phys. Rev. B **98**, 155143 (2018).

[4] B. Holm, Total Energies from GW Calculations, Phys. Rev. Lett. **83**, 788 (1999).

[5] P. García-González and R. W. Godby, Self-consistent calculation of total energies of the electron gas using many-body perturbation theory, Phys. Rev. B **63**, 075112 (2001).

[6] A. Kutepov, S. Y. Savrasov, and G. Kotliar, Ground-state properties of simple elements from GW calculations, Phys. Rev. B **80**, 041103(R) (2009).

[7] J. F. Dobson and J. Wang, Successful Test of a Seamless van der Waals Density Functional, Phys. Rev. Lett. **82**, 2123 (1999).

[8] P. García-González and R. W. Godby, Many-body GW calculations of ground-state properties: Quasi-2D electron systems and van der Waals forces, Phys. Rev. Lett. **88**, 056406 (2002).

[9] J. Harl and G. Kresse, Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory, Phys. Rev. B **77**, 045136 (2008).

[10] D. Lu, Y. Li, D. Rocca, and G. Galli, Ab Initio Calculation of van der Waals Bonded Molecular Crystals, Phys. Rev. Lett. **102**, 206411 (2009).

[11] L. Schimka, J. Harl, A. Stroppa, A. Grueneis, M. Marsman, F. Mittendorfer, and G. Kresse, Accurate surface and adsorption energies from many-body perturbation theory, Nat. Mater. **9**, 741 (2010).

[12] S. Lebègue, J. Harl, T. Gould, J. G. Ángyán, G. Kresse, and J. F. Dobson, Cohesive Properties and Asymptotics of the Dispersion Interaction in Graphite by the Random Phase Approximation, Phys. Rev. Lett. **105**, 196401 (2010).

[13] F. Bruneval, Range-Separated Approach to the RPA Correlation Applied to the van der Waals Bond and to Diffusion of Defects, Phys. Rev. Lett. **108**, 256403 (2012).

[14] C.-O. Almbladh, U. V. Barth, and R. V. Leeuwen, Variational total energies from $\Phi$- and $\Psi$-derivable theories, Int. J. Mod. Phys. B **13**, 535 (1999).

[15] N. E. Dahlen and U. von Barth, Variational second-order Møller-Plesset theory based on the Luttinger-ward functional, J. Chem. Phys. **120**, 6826 (2004).

[16] N. E. Dahlen and U. von Barth, Variational energy functionals tested on atoms, Phys. Rev. B **69**, 195102 (2004).

[17] N. E. Dahlen, R. van Leeuwen, and U. von Barth, Variational energy functionals of the green function and of the density tested on molecules, Phys. Rev. A **73**, 012511 (2006).

[18] A. Stan, N. E. Dahlen, and R. van Leeuwen, Fully self-consistent GW calculations for atoms and molecules, Europhys. Lett. **76**, 298 (2006).

[19] M. Hellgren, F. Caruso, D. R. Rohr, X. Ren, A. Rubio, M. Scheffler, and P. Rinke, Static correlation and electron localization in molecular dimers from the self-consistent rpa and GW approximation, Phys. Rev. B **91**, 165110 (2015).

[20] F. Bruneval, Assessment of the linearized GW density matrix for molecules, J. Chem. Theory Comput. **15**, 4069 (2019).

[21] F. Bruneval, M. Rodríguez-Mayorga, P. Rinke, and M. Dvorak, Improved one-shot total energies from the linearized GW density matrix, J. Chem. Theory Comput. **17**, 2126 (2021).

[22] F. Caruso, P. Rinke, X. Ren, M. Scheffler, and A. Rubio, Unified description of ground and excited states of finite systems: The self-consistent GW approach, Phys. Rev. B **86**, 081102(R) (2012).

[23] F. Caruso, D. R. Rohr, M. Hellgren, X. Ren, P. Rinke, A. Rubio, and M. Scheffler, Bond Breaking and Bond Formation: How Electron Correlation is Captured in Many-Body Perturbation Theory and Density-Functional Theory, Phys. Rev. Lett. **110**, 146403 (2013).
as a perfect pairing approach. Next, let us assume that the mean-field Bloch waves $\phi_{k\Gamma 4}$ and $\phi_{k\chi 5}$ coincide with the optimal natural orbitals in order to skip the orbital optimization procedure.

Next, let us focus on the energy contribution arising from the N/2 geminal to the second term in the right-hand side of the APSG energy [see Eq. (B16)]:

$$n_{k\Gamma 4}(\mathbf{k}_{\Gamma}4\mathbf{k}_{\Gamma}4|\mathbf{k}_{\Gamma}4\mathbf{k}_{\Gamma}4) + n_{k\chi 5}(\mathbf{k}_{\chi}5\mathbf{k}_{\chi}5|\mathbf{k}_{\chi}5\mathbf{k}_{\chi}5) + 2\zeta_{k\Gamma 4}\zeta_{k\chi 5}\sqrt{n_{k\Gamma 4}n_{k\chi 5}}(\mathbf{k}_{\Gamma}4\mathbf{k}_{\chi}5|\mathbf{k}_{\chi}5\mathbf{k}_{\Gamma}4). \quad (B19)$$

We set the usual approximation for (fixing) the phases (i.e., $\zeta_{k\Gamma 4}\zeta_{k\chi 5} = -1$ [77,78,80,81]) for the interaction among the states above and below the Fermi level and allow all two-electron integrals to be equal (which can occur in the extreme case when degenerated states are involved). The occupation numbers that would minimize the energy contribution are $n_{k\Gamma 4} = n_{k\chi 5} = 1/2$, illustrating that a reorganization of electrons occurs among the $\Gamma$ and $X$ k points. In the Si crystal, the mean-field Bloch waves are not completely degenerate in terms of energy and do not correspond to the optimal natural

orbitals; then, the actual optimal occupation numbers differ from 1/2. But, they also differ from the initial values at the mean-field level, where $n_{k\Gamma 4} = 1$ and $n_{k\chi 5} = 0$. Moreover, beyond the perfect pairing approach, the coupling of states to form a geminal wave function can include states belonging to other k points. Since the geminal wave functions are built with states for the $\alpha$ and the $\beta$ electron; the state for the $\alpha$ electrons is associated with a $k''$ k point while the state for the $\beta$ is related to a $-k''$ k point. The Hartree product in Eq. (B14) conserves the $\mathbf{K} = \mathbf{k}'' - \mathbf{k}' = 0$ value, which could lead to further reorganization of electrons among different k points beyond the perfect pairing approach. Thus, for example, the coupling of states belonging to $\Gamma$, $X$, $\Delta$, etc., is allowed in the Si crystal. Actually, the coupling of all k-point values in the first Brillouin zone is valid to build geminal wave functions.

Finally, let us remark that this example is based on a valid approximation to the many-electron wave function (i.e., a $\Psi_{\mathbf{K}}^{\text{APSG}}$ ansatz), where we illustrate that the reorganization of electrons among k points is purely a consequence of the electronic correlation effects.

[1] W. Ku and A. G. Eguiluz, Band-Gap Problem in Semiconductors Revisited: Effects of Core States and Many-Body Self-Consistency, Phys. Rev. Lett. **89**, 126401 (2002).

[2] H. Cao, Z. Yu, P. Lu, and L.-W. Wang, Fully converged plane-wave-based self-consistent GW calculations of periodic solids, Phys. Rev. B **95**, 035139 (2017).

[3] M. Grumet, P. Liu, M. Kaltak, J. Klimeš, and G. Kresse, Beyond the quasiparticle approximation: Fully self-consistent GW calculations, Phys. Rev. B **98**, 155143 (2018).

[4] B. Holm, Total Energies from GW Calculations, Phys. Rev. Lett. **83**, 788 (1999).

[5] P. García-González and R. W. Godby, Self-consistent calculation of total energies of the electron gas using many-body perturbation theory, Phys. Rev. B **63**, 075112 (2001).

[6] A. Kutepov, S. Y. Savrasov, and G. Kotliar, Ground-state properties of simple elements from GW calculations, Phys. Rev. B **80**, 041103(R) (2009).

[7] J. F. Dobson and J. Wang, Successful Test of a Seamless van der Waals Density Functional, Phys. Rev. Lett. **82**, 2123 (1999).

[8] P. García-González and R. W. Godby, Many-body GW calculations of ground-state properties: Quasi-2D electron systems and van der Waals forces, Phys. Rev. Lett. **88**, 056406 (2002).

[9] J. Harl and G. Kresse, Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory, Phys. Rev. B **77**, 045136 (2008).

[10] D. Lu, Y. Li, D. Rocca, and G. Galli, Ab Initio Calculation of van der Waals Bonded Molecular Crystals, Phys. Rev. Lett. **102**, 206411 (2009).

[11] L. Schimka, J. Harl, A. Stroppa, A. Grueneis, M. Marsman, F. Mittendorfer, and G. Kresse, Accurate surface and adsorption energies from many-body perturbation theory, Nat. Mater. **9**, 741 (2010).

[12] S. Lebègue, J. Harl, T. Gould, J. G. Ángyán, G. Kresse, and J. F. Dobson, Cohesive Properties and Asymptotics of the Dispersion Interaction in Graphite by the Random Phase Approximation, Phys. Rev. Lett. **105**, 196401 (2010).

[13] F. Bruneval, Range-Separated Approach to the RPA Correlation Applied to the van der Waals Bond and to Diffusion of Defects, Phys. Rev. Lett. **108**, 256403 (2012).

[14] C.-O. Almbladh, U. V. Barth, and R. V. Leeuwen, Variational total energies from $\Phi$- and $\Psi$-derivable theories, Int. J. Mod. Phys. B **13**, 535 (1999).

[15] N. E. Dahlen and U. von Barth, Variational second-order Møller-Plesset theory based on the Luttinger-ward functional, J. Chem. Phys. **120**, 6826 (2004).

[16] N. E. Dahlen and U. von Barth, Variational energy functionals tested on atoms, Phys. Rev. B **69**, 195102 (2004).

[17] N. E. Dahlen, R. van Leeuwen, and U. von Barth, Variational energy functionals of the green function and of the density tested on molecules, Phys. Rev. A **73**, 012511 (2006).

[18] A. Stan, N. E. Dahlen, and R. van Leeuwen, Fully self-consistent GW calculations for atoms and molecules, Europhys. Lett. **76**, 298 (2006).

[19] M. Hellgren, F. Caruso, D. R. Rohr, X. Ren, A. Rubio, M. Scheffler, and P. Rinke, Static correlation and electron localization in molecular dimers from the self-consistent rpa and GW approximation, Phys. Rev. B **91**, 165110 (2015).

[20] F. Bruneval, Assessment of the linearized GW density matrix for molecules, J. Chem. Theory Comput. **15**, 4069 (2019).

[21] F. Bruneval, M. Rodríguez-Mayorga, P. Rinke, and M. Dvorak, Improved one-shot total energies from the linearized GW density matrix, J. Chem. Theory Comput. **17**, 2126 (2021).

[22] F. Caruso, P. Rinke, X. Ren, M. Scheffler, and A. Rubio, Unified description of ground and excited states of finite systems: The self-consistent GW approach, Phys. Rev. B **86**, 081102(R) (2012).

[23] F. Caruso, D. R. Rohr, M. Hellgren, X. Ren, P. Rinke, A. Rubio, and M. Scheffler, Bond Breaking and Bond Formation: How Electron Correlation is Captured in Many-Body Perturbation Theory and Density-Functional Theory, Phys. Rev. Lett. **110**, 146403 (2013).
[24] F. Bruneval, Improved density matrices for accurate molecular ionization potentials, *Phys. Rev. B* **99**, 041118(R) (2019).

[25] R. M. Martin, *Electronic Structure: Basic Theory and Practical Methods* (Cambridge University Press, New York, 2004), Vol. 1.

[26] D. Pines and P. Nozières, *Theory of Quantum Liquids* (Benjamin, New York, 1966).

[27] A. Klein, Perturbation theory for an infinite medium of fermions. II, *Phys. Rev. 121*, 950 (1961).

[28] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, *Phys. Rev. Lett.* **77**, 3865 (1996).

[29] C. Adamo and V. Barone, Toward reliable density functional methods without adjustable parameters: The PBE0 model, *J. Chem. Phys.* **110**, 6158 (1999).

[30] V. M. Galitskii and A. B. Migdal, Application of quantum field theory methods to the many body problem, *Sov. Phys. JETP* **34**, 96 (1958).

[31] G. Baym and L. P. Kadanoff, Conservation laws and correlation functions, *Phys. Rev. 124*, 287 (1961).

[32] A. Schindlmayr, P. García-González, and R. W. Godby, Diagrammatic self-energy approximations and the total particle number, *Phys. Rev. B* **64**, 235106 (2001).

[33] F. Caruso, Self-consistent GW approach for the unified description of ground and excited states of finite systems, Ph.D. thesis, Freie Universität Berlin, 2013.

[34] L. J. Sham and M. Schlüter, Density-Functional Theory of the Energy Gap, *Phys. Rev. Lett.* **51**, 1888 (1983).

[35] L. Hedin, New method for calculating the one-particle green's function with application to the electron-gas problem, *Phys. Rev.* **139**, A796 (1965).

[36] L. Hedin and S. Lundqvist, in *Effects of Electron-Electron and Electron-Phonon Interactions on the One-Electron States of Solids*, edited by D. T. Frederick Seitz and H. Ehrenreich, Solid State Physics Vol. 23 (Academic Press, New York, 1970), pp. 1-181.

[37] W. G. Aulbur, L. Jönsson, and J. W. Wilkins, in *Quasiparticle Calculations in Solids*, edited by H. Ehrenreich and F. Spaepen, Solid State Physics Vol. 54 (Academic Press, New York, 2000), pp. 1-218.

[38] L. Reining, The GW approximation: Content, successes and limitations, *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **8**, e1344 (2018).

[39] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, *Phys. Rev. B* **54**, 11169 (1996).

[40] X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, N. Brouwer, F. Bruneval, G. Brunin, T. Cavignac, J.-B. Charraud, W. Chen, M. Côté, S. Cottenier, J. Denier, G. Geneste et al., The Abinit project: Impact, environment and recent developments, *Comput. Phys. Commun.* **248**, 107042 (2020).

[41] J. Klimeš, M. Kaltak, E. Maggio, and G. Kresse, Singles correlation energy contributions in solids, *J. Chem. Phys.* **143**, 102816 (2015).

[42] M. Fuchs and X. Gonze, Accurate density functionals: Approaches using the adiabatic-connection fluctuation-dissipation theorem, *Phys. Rev. B* **65**, 235109 (2002).

[43] H.-V. Nguyen and G. Galli, A first-principles study of weakly bound molecules using exact exchange and the random phase approximation, *J. Chem. Phys.* **132**, 044109 (2010).

[44] J. G. Ángyán, R.-F. Liu, J. Toulouse, and G. Jansen, Correlation energy expressions from the adiabatic-connection fluctuation-dissipation theorem approach, *J. Chem. Theory Comput.* **7**, 3116 (2011).

[45] B. Ramberger, T. Schäfer, and G. Kresse, Analytic Interatomic Forces in the Random Phase Approximation, *Phys. Rev. Lett.* **118**, 106403 (2017).

[46] P. E. Blöchl, Projector augmented-wave method, *Phys. Rev. B* **50**, 17953 (1994).

[47] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, *Phys. Rev. B* **59**, 1758 (1999).

[48] M. Torrent, F. Jollet, F. Bottin, G. Zérah, and X. Gonze, Implementation of the projector augmented-wave method in the ABINIT code: Application to the study of iron under pressure, *Comput. Mater. Sci.* **42**, 337 (2008).

[49] B. Arnaud and M. Alouani, All-electron projector-augmented-wave GW approximation: Application to the electronic properties of semiconductors, *Phys. Rev. B* **62**, 4464 (2000).

[50] J. Klimeš, M. Kaltak, and G. Kresse, Predictive GW calculations using plane waves and pseudopotentials, *Phys. Rev. B* **90**, 075125 (2014).

[51] L. Kleinman and D. M. Bylander, Efficacious Form for Model Pseudopotentials, *Phys. Rev. Lett.* **48**, 1425 (1982).

[52] D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, *Phys. Rev. B* **88**, 085117 (2013).

[53] M. Kaltak and G. Kresse, Minimax isometry method: A compressive sensing approach for Matsubara summation in many-body perturbation theory, *Phys. Rev. B* **101**, 205145 (2020).

[54] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevB.108.125107 for the pseudopotential study, the convergence study, and the complete results for the seven covalently bound solids studied here.

[55] R. Dovesi, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. Bush, P. D'Arco, Y. Noël, M. Rérat, P. Carbonnière, M. Causà, S. Salustro, V. Lacivita, B. Kirtman, A. M. Ferrari, F. S. Gentile, J. Baima, M. Ferrero, R. Demichelis, and M. De La Pierre, The crystal code, 1976-2020 and beyond, a long story, *J. Chem. Phys.* **152**, 204111 (2020).

[56] J. Heyd, J. E. Peralta, G. E. Scuseria, and R. L. Martin, Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional, *J. Chem. Phys.* **123**, 174101 (2005).

[57] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, The pseudodojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, *Comput. Phys. Commun.* **226**, 39 (2018).

[58] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos et al., Quantum espresso: A modular and open-source software project for quantum simulations of materials, *J. Phys.: Condens. Matter* **21**, 395502 (2009).

[59] J. B. Krieger, Y. Li, and G. J. Iafrate, Exact relations in the optimized effective potential method employing an arbitrary $E_{xc}[\psi_{i\sigma}]$, *Phys. Lett. A* **146**, 256 (1990).

[60] M. C. Bennett, C. A. Melton, A. Annaberdiyev, G. Wang, L. Shulenburger, and L. Mitas, A new generation of effective
[24] F. Bruneval, Improved density matrices for accurate molecular ionization potentials, *Phys. Rev. B* **99**, 041118(R) (2019).

[25] R. M. Martin, *Electronic Structure: Basic Theory and Practical Methods* (Cambridge University Press, New York, 2004), Vol. 1.

[26] D. Pines and P. Nozières, *Theory of Quantum Liquids* (Benjamin, New York, 1966).

[27] A. Klein, Perturbation theory for an infinite medium of fermions. II, *Phys. Rev. 121*, 950 (1961).

[28] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, *Phys. Rev. Lett.* **77**, 3865 (1996).

[29] C. Adamo and V. Barone, Toward reliable density functional methods without adjustable parameters: The PBE0 model, *J. Chem. Phys.* **110**, 6158 (1999).

[30] V. M. Galitskii and A. B. Migdal, Application of quantum field theory methods to the many body problem, *Sov. Phys. JETP* **34**, 96 (1958).

[31] G. Baym and L. P. Kadanoff, Conservation laws and correlation functions, *Phys. Rev. 124*, 287 (1961).

[32] A. Schindlmayr, P. García-González, and R. W. Godby, Diagrammatic self-energy approximations and the total particle number, *Phys. Rev. B* **64**, 235106 (2001).

[33] F. Caruso, Self-consistent GW approach for the unified description of ground and excited states of finite systems, Ph.D. thesis, Freie Universität Berlin, 2013.

[34] L. J. Sham and M. Schlüter, Density-Functional Theory of the Energy Gap, *Phys. Rev. Lett.* **51**, 1888 (1983).

[35] L. Hedin, New method for calculating the one-particle green's function with application to the electron-gas problem, *Phys. Rev.* **139**, A796 (1965).

[36] L. Hedin and S. Lundqvist, in *Effects of Electron-Electron and Electron-Phonon Interactions on the One-Electron States of Solids*, edited by D. T. Frederick Seitz and H. Ehrenreich, Solid State Physics Vol. 23 (Academic Press, New York, 1970), pp. 1-181.

[37] W. G. Aulbur, L. Jönsson, and J. W. Wilkins, in *Quasiparticle Calculations in Solids*, edited by H. Ehrenreich and F. Spaepen, Solid State Physics Vol. 54 (Academic Press, New York, 2000), pp. 1-218.

[38] L. Reining, The GW approximation: Content, successes and limitations, *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **8**, e1344 (2018).

[39] G. Kresse and J. Furthmüller, Efficient iterative schemes for *ab initio* total-energy calculations using a plane-wave basis set, *Phys. Rev. B* **54**, 11169 (1996).

[40] X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, N. Brouwer, F. Bruneval, G. Brunin, T. Cavignac, J.-B. Charraud, W. Chen, M. Côté, S. Cottenier, J. Denier, G. Geneste et al., The Abinit project: Impact, environment and recent developments, *Comput. Phys. Commun.* **248**, 107042 (2020).

[41] J. Klimeš, M. Kaltak, E. Maggio, and G. Kresse, Singles correlation energy contributions in solids, *J. Chem. Phys.* **143**, 102816 (2015).

[42] M. Fuchs and X. Gonze, Accurate density functionals: Approaches using the adiabatic-connection fluctuation-dissipation theorem, *Phys. Rev. B* **65**, 235109 (2002).

[43] H.-V. Nguyen and G. Galli, A first-principles study of weakly bound molecules using exact exchange and the random phase approximation, *J. Chem. Phys.* **132**, 044109 (2010).

[44] J. G. Ángyán, R.-F. Liu, J. Toulouse, and G. Jansen, Correlation energy expressions from the adiabatic-connection fluctuation-dissipation theorem approach, *J. Chem. Theory Comput.* **7**, 3116 (2011).

[45] B. Ramberger, T. Schäfer, and G. Kresse, Analytic Interatomic Forces in the Random Phase Approximation, *Phys. Rev. Lett.* **118**, 106403 (2017).

[46] P. E. Blöchl, Projector augmented-wave method, *Phys. Rev. B* **50**, 17953 (1994).

[47] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, *Phys. Rev. B* **59**, 1758 (1999).

[48] M. Torrent, F. Jollet, F. Bottin, G. Zérah, and X. Gonze, Implementation of the projector augmented-wave method in the ABINIT code: Application to the study of iron under pressure, *Comput. Mater. Sci.* **42**, 337 (2008).

[49] B. Arnaud and M. Alouani, All-electron projector-augmented-wave GW approximation: Application to the electronic properties of semiconductors, *Phys. Rev. B* **62**, 4464 (2000).

[50] J. Klimeš, M. Kaltak, and G. Kresse, Predictive GW calculations using plane waves and pseudopotentials, *Phys. Rev. B* **90**, 075125 (2014).

[51] L. Kleinman and D. M. Bylander, Efficacious Form for Model Pseudopotentials, *Phys. Rev. Lett.* **48**, 1425 (1982).

[52] D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, *Phys. Rev. B* **88**, 085117 (2013).

[53] M. Kaltak and G. Kresse, Minimax isometry method: A compressive sensing approach for Matsubara summation in many-body perturbation theory, *Phys. Rev. B* **101**, 205145 (2020).

[54] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevB.108.125107 for the pseudopotential study, the convergence study, and the complete results for the seven covalently bound solids studied here.

[55] R. Dovesi, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. Bush, P. D'Arco, Y. Noël, M. Rérat, P. Carbonnière, M. Causà, S. Salustro, V. Lacivita, B. Kirtman, A. M. Ferrari, F. S. Gentile, J. Baima, M. Ferrero, R. Demichelis, and M. De La Pierre, The crystal code, 1976-2020 and beyond, a long story, *J. Chem. Phys.* **152**, 204111 (2020).

[56] J. Heyd, J. E. Peralta, G. E. Scuseria, and R. L. Martin, Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional, *J. Chem. Phys.* **123**, 174101 (2005).

[57] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, The pseudodojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, *Comput. Phys. Commun.* **226**, 39 (2018).

[58] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos et al., Quantum espresso: A modular and open-source software project for quantum simulations of materials, *J. Phys.: Condens. Matter* **21**, 395502 (2009).

[59] J. B. Krieger, Y. Li, and G. J. Iafrate, Exact relations in the optimized effective potential method employing an arbitrary $E_{xc}[\psi_{i\sigma}]$, *Phys. Lett. A* **146**, 256 (1990).

[60] M. C. Bennett, C. A. Melton, A. Annaberdiyev, G. Wang, L. Shulenburger, and L. Mitas, A new generation of effective
core potentials for correlated calculations, *J. Chem. Phys.* **147**, 224106 (2017).

[61] V. Olevano, A. Titov, M. Ladisa, K. Hämäläinen, S. Huotari, and M. Holzmann, Momentum distribution and Compton profile by the ab initio GW approximation, *Phys. Rev. B* **86**, 195123 (2012).

[62] D. D. Koelling, The band description of materials with localizing orbitals, *Int. J. Quantum Chem.* **30**, 377 (1986).

[63] A. Damascelli, Z. Hussain, and Z.-X. Shen, Angle-resolved photoemission studies of the cuprate superconductors, *Rev. Mod. Phys.* **75**, 473 (2003).

[64] A. Annaberdiyev, G. Wang, C. A. Melton, M. C. Bennett, and L. Mitas, Cohesion and excitations of diamond-structure silicon by quantum Monte Carlo: Benchmarks and control of systematic biases, *Phys. Rev. B* **103**, 205206 (2021).

[65] F. Bruneval, T. Rangel, S. M. Hamed, M. Shao, C. Yang, and J. B. Neaton, Molgw 1: Many-body perturbation theory software for atoms, molecules, and clusters, *Comput. Phys. Commun.* **208**, 149 (2016).

[66] D. M. Ceperley and B. J. Alder, Ground State of the Electron Gas by a Stochastic Method, *Phys. Rev. Lett.* **45**, 566 (1980).

[67] A. Marini, P. García-González, and A. Rubio, First-Principles Description of Correlation Effects in Layered Materials, *Phys. Rev. Lett.* **96**, 136404 (2006).

[68] http://www.abinit.org.

[69] S. Sharma, J. K. Dewhurst, N. N. Lathiotakis, and E. K. U. Gross, Reduced density matrix functional for many-electron systems, *Phys. Rev. B* **78**, 201103(R) (2008).

[70] F. Bloch, Über die quantenmechanik der elektronen in kristallgittern, *Z. Phys.* **52**, 555 (1929).

[71] A. Stoyanova, Delocalized and correlated wave functions for excited states in extended systems, Ph.D. thesis, University of Groningen, 2006.

[72] J. McClain, Q. Sun, G. K.-L. Chan, and T. C. Berkelbach, Gaussian-based coupled-cluster theory for the ground-state and band structure of solids, *J. Chem. Theory Comput.* **13**, 1209 (2017).

[73] T. Schäfer, B. Ramberger, and G. Kresse, Quartic scaling mp2 for solids: A highly parallelized algorithm in the plane wave basis, *J. Chem. Phys.* **146**, 104101 (2017).

[74] P. R. Surján, An introduction to the theory of geminals, *Correl. Localiz.* **203**, 63 (1999).

[75] H. A. Kramers, Théorie générale de la rotation paramagnétique dans les cristaux, Proc. Acad. Amst. **33** (1930).

[76] G. Aucar, H. A. Jensen, and J. Oddershede, Operator representations in Kramers bases, *Chem. Phys. Lett.* **232**, 47 (1995).

[77] P.-O. Löwdin and H. Shull, Natural orbitals in the quantum theory of two-electron systems, *Phys. Rev.* **101**, 1730 (1956).

[78] P.-O. Löwdin, Quantum theory of many-particle systems. ii. study of the ordinary Hartree-Fock approximation, *Phys. Rev.* **97**, 1490 (1955).

[79] K. Pernal, The equivalence of the piris natural orbital functional 5 (pnof5) and the antisymmetrized product of strongly orthogonal geminal theory, *Comput. Theor. Chem.* **1003**, 127 (2013).

[80] M. Piris, J. Matxain, and X. Lopez, The intrapair electron correlation in natural orbital functional theory, *J. Chem. Phys.* **139**, 234109 (2013).

[81] M. Piris, A natural orbital functional based on an explicit approach of the two-electron cumulant, *Int. J. Quantum Chem.* **113**, 620 (2013).

[82] M. Piris, Global Method for Electron Correlation, *Phys. Rev. Lett.* **119**, 063002 (2017).

[83] I. Mitxelena, M. Rodríguez-Mayorga, and M. Piris, Phase dilemma in natural orbital functional theory from the n-representability perspective, Eur. Phys. J. B **91**, 109 (2018).

[84] M. Piris, Global Natural Orbital Functional: Towards the Complete Description of the Electron Correlation, *Phys. Rev. Lett.* **127**, 233001 (2021).

[85] M. Rodríguez-Mayorga, K. Giesbertz, and L. Visscher, Relativistic reduced density matrix functional theory, *SciPost Chem.* **1**, 004 (2022).

[86] K. Giesbertz and R. van Leeuwen, Natural occupation numbers: When do they vanish? *J. Chem. Phys.* **139**, 104109 (2013).
core potentials for correlated calculations, *J. Chem. Phys.* **147**, 224106 (2017).

[61] V. Olevano, A. Titov, M. Ladisa, K. Hämäläinen, S. Huotari, and M. Holzmann, Momentum distribution and Compton profile by the ab initio GW approximation, *Phys. Rev. B* **86**, 195123 (2012).

[62] D. D. Koelling, The band description of materials with localizing orbitals, *Int. J. Quantum Chem.* **30**, 377 (1986).

[63] A. Damascelli, Z. Hussain, and Z.-X. Shen, Angle-resolved photoemission studies of the cuprate superconductors, *Rev. Mod. Phys.* **75**, 473 (2003).

[64] A. Annaberdiyev, G. Wang, C. A. Melton, M. C. Bennett, and L. Mitas, Cohesion and excitations of diamond-structure silicon by quantum Monte Carlo: Benchmarks and control of systematic biases, *Phys. Rev. B* **103**, 205206 (2021).

[65] F. Bruneval, T. Rangel, S. M. Hamed, M. Shao, C. Yang, and J. B. Neaton, Molgw 1: Many-body perturbation theory software for atoms, molecules, and clusters, *Comput. Phys. Commun.* **208**, 149 (2016).

[66] D. M. Ceperley and B. J. Alder, Ground State of the Electron Gas by a Stochastic Method, *Phys. Rev. Lett.* **45**, 566 (1980).

[67] A. Marini, P. García-González, and A. Rubio, First-Principles Description of Correlation Effects in Layered Materials, *Phys. Rev. Lett.* **96**, 136404 (2006).

[68] http://www.abinit.org.

[69] S. Sharma, J. K. Dewhurst, N. N. Lathiotakis, and E. K. U. Gross, Reduced density matrix functional for many-electron systems, *Phys. Rev. B* **78**, 201103(R) (2008).

[70] F. Bloch, Über die quantenmechanik der elektronen in kristallgittern, *Z. Phys.* **52**, 555 (1929).

[71] A. Stoyanova, Delocalized and correlated wave functions for excited states in extended systems, Ph.D. thesis, University of Groningen, 2006.

[72] J. McClain, Q. Sun, G. K.-L. Chan, and T. C. Berkelbach, Gaussian-based coupled-cluster theory for the ground-state and band structure of solids, *J. Chem. Theory Comput.* **13**, 1209 (2017).

[73] T. Schäfer, B. Ramberger, and G. Kresse, Quartic scaling mp2 for solids: A highly parallelized algorithm in the plane wave basis, *J. Chem. Phys.* **146**, 104101 (2017).

[74] P. R. Surján, An introduction to the theory of geminals, *Correl. Localiz.* **203**, 63 (1999).

[75] H. A. Kramers, Théorie générale de la rotation paramagnétique dans les cristaux, Proc. Acad. Amst. **33** (1930).

[76] G. Aucar, H. A. Jensen, and J. Oddershede, Operator representations in Kramers bases, *Chem. Phys. Lett.* **232**, 47 (1995).

[77] P.-O. Löwdin and H. Shull, Natural orbitals in the quantum theory of two-electron systems, *Phys. Rev.* **101**, 1730 (1956).

[78] P.-O. Löwdin, Quantum theory of many-particle systems. ii. study of the ordinary Hartree-Fock approximation, *Phys. Rev.* **97**, 1490 (1955).

[79] K. Pernal, The equivalence of the piris natural orbital functional 5 (pnof5) and the antisymmetrized product of strongly orthogonal geminal theory, *Comput. Theor. Chem.* **1003**, 127 (2013).

[80] M. Piris, J. Matxain, and X. Lopez, The intrapair electron correlation in natural orbital functional theory, *J. Chem. Phys.* **139**, 234109 (2013).

[81] M. Piris, A natural orbital functional based on an explicit approach of the two-electron cumulant, *Int. J. Quantum Chem.* **113**, 620 (2013).

[82] M. Piris, Global Method for Electron Correlation, *Phys. Rev. Lett.* **119**, 063002 (2017).

[83] I. Mitxelena, M. Rodríguez-Mayorga, and M. Piris, Phase dilemma in natural orbital functional theory from the n-representability perspective, Eur. Phys. J. B **91**, 109 (2018).

[84] M. Piris, Global Natural Orbital Functional: Towards the Complete Description of the Electron Correlation, *Phys. Rev. Lett.* **127**, 233001 (2021).

[85] M. Rodríguez-Mayorga, K. Giesbertz, and L. Visscher, Relativistic reduced density matrix functional theory, *SciPost Chem.* **1**, 004 (2022).

[86] K. Giesbertz and R. van Leeuwen, Natural occupation numbers: When do they vanish? *J. Chem. Phys.* **139**, 104109 (2013).