1 Introduction↩︎

GW calculations [1]–[3] have become a standard method for calculating electron addition and removal energies of molecules [4]–[7], two-dimensional materials [8]–[10] and bulk solids [11]–[13]. Recent advancements of the \(GW\) method span a broad spectrum, including the application to deep core excitations [14]–[22], relativistic GW schemes [23]–[28] and vertex corrections [29]–[38]. These developments have firmly established GW as a powerful and versatile approach within the domain of many-body perturbation theory.

Despite the methodological maturity of \(GW\), computational challenges remain, particularly when applied to atomically thin materials. One of the primary limitations arises from the use of plane-wave basis sets in systems with vacuum, such as molecules or low-dimensional systems. The need to represent the vacuum leads to large plane-wave basis set size and significant computational cost. These constraints motivate the development of alternative basis representations. One compelling solution involves atom-centered basis functions, which are localized and naturally adapted to such geometries. Atom-centered basis functions are already the standard for \(GW\) implementations targeting molecules [3]–[7]. Also, several implementations of \(GW\) with atom-centered basis functions and periodic boundary conditions have been reported [39]–[47].

The first pioneering \(GW\) calculations in semiconductors were carried by Strinati et al. [39], [40], who studied the band structure of diamond using a local-orbital formulation of the \(GW\) calculation. In the 90ies, Rohlfing et al. [41], [42] studied bulk semiconductors and a silicon surface finding good agreement to plane-wave based \(GW\) implementations in the order of 0.1 eV or better. The more recent periodic \(GW\) implementations with atom-centered basis functions [43]–[46] focus on 3D crystals and report similar precision. The \(GW\) implementations [39]–[46] use a formulation of \(GW\) in frequency and rely on \(k\)-point sampling to account for the periodic boundary conditions.

An alternative is the \(GW\) space-time method [48], [49] which has been originally formulated using plane waves and real-space representations as well as time and frequency representations. The \(GW\) space-time method can be also reformulated using an atom-centered basis set [50]–[53] which eliminates the use of plane-waves and real-space grids and allows for low-scaling \(GW\) calculations on large molecules. The \(GW\) space-time approach in an atom-centered basis can be also combined with periodic boundary conditions, as we have demonstrated in our previous work [47] where we employed a \(\Gamma\)-only approach for the density response function and the self-energy, while relying on dense \(k\)-point sampling for the screened Coulomb interaction \(W\) to account for the divergence of \(W\) at the \(\Gamma\)-point. This implementation enabled the study of a twisted transition-metal dichalcogenide heterobilayers with almost 1000 atoms in the unit cell. The drawback of the \(\Gamma\)-only implementation [47] is that large unit cells are required where the density response needs to vanish on the length scale of the unit cell.

In this work, we address the limitation of the \(\Gamma\)-point only approach [47], which requires large unit cells, by introducing a lattice summation over neighbor cells for both the density response and self-energy, as originally proposed in the \(GW\) space-time method [48]. This extension enables the accurate and efficient treatment of crystals with small unit cells. The resulting \(GW\) algorithm is particularly well-suited for low-dimensional materials, as the number of neighboring cells required in the lattice sums is significantly reduced compared to 3D bulk crystals. Furthermore, the Gaussian basis set is efficient in simulations involving large vacuum regions, a common requirement when modeling low-dimensional systems.

Our implementation also supports the inclusion of relativistic effects via spin-orbit coupling (SOC) from Gaussian dual-space pseudopotentials [54]–[56] and a perturbative correction to the quasiparticle energies. We focus on \(GW\) band structure calculations of atomically thin transition metal dichalcogenides (TMDCs) MoS\(_2\), MoSe\(_2\), WS\(_2\), WSe\(_2\) and we benchmark our \(GW\) band structures against state-of-the-art plane-wave-based \(GW\) implementations in BerkeleyGW [57] and VASP [12]. We demonstrate that our \(GW\) algorithm yields accurate and converged quasiparticle band structures across multiple convergence parameters, including basis set size, \(k\)-point sampling, the summation of neighbor cells, and the time- and frequency-mesh. We also discuss the computational efficiency of our \(GW\) band structure algorithm, which enables a \(GW\) band structure calculation of an atomically thin material on a laptop (Intel Xeon i5, 192 GB RAM) within roughly a day.

2 GW space-time algorithm↩︎

Many efficient \(GW\) algorithms [49], [51], [52] build on the \(GW\) space-time method [48]. In order to introduce the basic idea of the \(GW\) space-time method, we use a generic formulation in this section for non-periodic systems projecting all quantities on real-space grids. It is important to note that this formulation differs from the original \(GW\) space-time method [48] where some quantities were calculated using a plane-wave basis set.

In this work, we employ the \(G_0W_0\) scheme which starts from a self-consistent Kohn-Sham density functional theory (KS-DFT) calculation [58], \[[h_0 (\mathbf{r}) + v^\text{xc} (\mathbf{r})] \psi_n (\mathbf{r}) = \varepsilon_n^{\operatorname{DFT}} \psi_n (\mathbf{r})\,. \label{e1}\tag{1}\] \(h_0\) contains the kinetic energy, the Hartree potential and the external potential, while the exchange-correlation potential \(v^\text{xc}\) contains all electron-electron interactions beyond Hartree. \(\psi_n(\mathbf{r})\) is the KS orbital \(n\) and \(\varepsilon_n^{\operatorname{DFT}}\) the associated KS eigenvalue. The terms \(G_0\) and \(W_0\) indicate that the Green’s function \(G\) and the screened Coulomb interaction \(W\) are both computed using KS orbitals and KS eigenvalues, i.e., self-consistent updates of \(G\) and \(W\) from Green’s function theory are omitted in \(G_0W_0\).

KS orbitals and eigenvalues are used to calculate the single-particle Green’s function in imaginary time, \[\label{green95spacetime} G (\mathbf{r}, \mathbf{r}', i \tau) = \begin{cases} \;\;\;\, i \sum\limits_i^{\operatorname{occ}} \psi_i (\mathbf{r}) \psi_i^* (\mathbf{r}') e^{ -|(\varepsilon_i^\text{DFT}-\varepsilon_\text{F}) \tau|}, & \tau < 0\;,\\[1em] - i \sum\limits_a^{\operatorname{empty}} \psi_a^* (\mathbf{r}) \psi_a (\mathbf{r}') e^{ -|(\varepsilon_a^\text{DFT}-\varepsilon_\text{F}) \tau|}, & \tau > 0\;, \end{cases}\tag{2}\] where the sum over the index \(i\) runs over all occupied KS orbitals and the sum over the index \(a\) over all virtual, i.e., empty KS orbitals. \(\varepsilon_\text{F}\) is the Fermi level. The irreducible polarizability follows, \[\chi (\mathbf{r}, \mathbf{r}', i \tau) = - iG (\mathbf{r}, \mathbf{r}', i \tau) G (\mathbf{r}, \mathbf{r}', - i \tau)\;, \label{e3}\tag{3}\] which is then transformed to imaginary frequency, \[\begin{align} \chi (\mathbf{r}, \mathbf{r}', i\omega) = i \int\limits_{-\infty}^\infty e^{-i\omega\tau}\; \chi (\mathbf{r}, \mathbf{r}',i\tau)\;\mathrm{d}\tau\,. \end{align}\] This transform can be understood as Laplace transform followed by analytic continuation to the imaginary axis and effectively takes the form of a Fourier transform [59]. Next, the dielectric function \(\epsilon\) can be calculated in imaginary frequency from the irreducible polarizability, \[\epsilon (\mathbf{r}, \mathbf{r}', i \omega) = \delta (\mathbf{r}- \mathbf{r}') - \int \mathrm{d}\mathbf{r}'' v (\mathbf{r}, \mathbf{r}'') \,\chi (\mathbf{r}'', \mathbf{r}', i \omega)\;, \label{eps}\tag{4}\] using the Dirac delta function \(\delta(\mathbf{r})\) and the Coulomb interaction \(v (\mathbf{r}, \mathbf{r}'){=}1/|\mathbf{r}{-}\mathbf{r}'|\). The screened Coulomb interaction can be computed from the inverse dielectric function, \[W (\mathbf{r}, \mathbf{r}', i \omega) = \int \mathrm{d}\mathbf{r}''\, \epsilon^{- 1} (\mathbf{r}, \mathbf{r}'', i \omega) \,v (\mathbf{r}'', \mathbf{r}')\;. \label{W}\tag{5}\] It is convenient in \(GW\) implementations to split the screened interaction \(W\) into the bare Coulomb interaction \(v\) and the correlation part \(W^\text{c}\), \[\begin{align} W^\text{c} (\mathbf{r}, \mathbf{r}', i \omega) = W (\mathbf{r}, \mathbf{r}', i \omega) - v(\mathbf{r}, \mathbf{r}')\,. \end{align}\] In the space-time method, \(W^\text{c}\) is required in imaginary time, \[\begin{align} W^\text{c} (\mathbf{r}, \mathbf{r}', i\tau) = \frac{i}{2\pi} \int\limits_{-\infty}^\infty e^{i\omega\tau}\; W^\text{c} (\mathbf{r}, \mathbf{r}',i\omega)\;\mathrm{d}\omega\,, \end{align}\] and the correlation self-energy is given as product of the Green’s function and the screened Coulomb interaction, \[\label{sigma-ls} \Sigma^\text{c} (\mathbf{r}, \mathbf{r}', i \tau) = i\,G (\mathbf{r}, \mathbf{r}', i \tau) W^\text{c} (\mathbf{r}, \mathbf{r}', i \tau)\;.\tag{6}\] The self-energy is then transformed to imaginary frequency, \[\begin{align} \Sigma^\text{c} (\mathbf{r}, \mathbf{r}', i\omega) = i \int\limits_{-\infty}^\infty e^{-i\omega\tau}\; \Sigma^\text{c} (\mathbf{r}, \mathbf{r}',i\tau)\;\mathrm{d}\tau\,,\label{Sigmactime} \end{align}\tag{7}\] and we calculate its \((n,n)\)-diagonal element, \[\begin{align} \Sigma_n^\text{c} (i\omega) &= \langle \psi_n| \Sigma^\text{c} (i\omega)|\psi_n\rangle \nonumber \\&= \int \mathrm{d}\mathbf{r}\,\mathrm{d}\mathbf{r}' \,\psi_n^*(\mathbf{r})\, \Sigma^\text{c} (\mathbf{r}, \mathbf{r}', i\omega)\,\psi_n(\mathbf{r}')\,. \end{align}\] The self-energy is then analytically continued to real frequency [3], [60].

Focusing on the \(G_0W_0\) method already introduced before, we use KS orbitals to approximate the QP wavefunctions and we compute \(G\) and \(W\) only once using KS orbitals and KS eigenvalues from Eqs. 2  – 5 . The QP energies \(\varepsilon_n^{G_0 W_0}\) can finally be calculated by solving the QP equation, \[\label{energy-it} \varepsilon_n^{G_0 W_0} = \varepsilon_n^{\operatorname{DFT}} + \operatorname{Re} \Sigma_n^\text{c} (\varepsilon_n^{G_0 W_0}) + \Sigma_n^\text{x} - v_n^{\operatorname{xc}}\,,\tag{8}\] where \(\Sigma_n^\text{x}\) and \(v_n^{\operatorname{xc}}\) are the \(n,n\)-diagonal elements of the exchange self-energy and the exchange-correlation potential.

3 KS-DFT with periodic boundary conditions and Gaussian basis functions↩︎

We use KS-DFT with periodic boundary conditions [61]–[63], i.e., \(h_0 (\mathbf{r})\) and \(v^\text{xc} (\mathbf{r})\) from Eq. 1 are lattice periodic, \[\begin{align} h_0 (\mathbf{r}+\mathbf{R}) = h_0(\mathbf{r}) \;,v^\text{xc} (\mathbf{r}+\mathbf{R}) = v^\text{xc} (\mathbf{r})\,, \end{align}\] for every lattice vector \[\begin{align} \mathbf{R}{=}\sum_{j=1}^{d} n_j\,\mathbf{a}_j\,, \end{align}\] where \(d\) is the dimension, \(n_j\) are integers and \(\mathbf{a}_j\) the primitive vectors of the lattice.

Bloch’s theorem [64] states that the solutions \(\psi_{n\mathbf{k}}(\mathbf{r})\) of the Kohn-Sham equations \[[h_0 (\mathbf{r}) + v^\text{xc} (\mathbf{r})] \psi_{n\mathbf{k}} (\mathbf{r}) = \varepsilon_{n\mathbf{k}} ^{\operatorname{DFT}} \psi_{n\mathbf{k}} (\mathbf{r}) \label{e14}\tag{9}\] with lattice-periodic \(h_0(\mathbf{r})\) and \(v^\text{xc} (\mathbf{r})\) are Bloch functions \[\begin{align} \psi_{n\mathbf{k}} (\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}}u_{n\mathbf{k}}(\mathbf{r}) \label{e16} \end{align}\tag{10}\] where \(\mathbf{k}\) is the crystal momentum in the first Brillouin zone (BZ) and \(u_{n\mathbf{k}}(\mathbf{r})\) is a lattice-periodic function, \(u_{n\mathbf{k}}(\mathbf{r}){=}u_{n\mathbf{k}}(\mathbf{r}{+}\mathbf{R})\). The eigenvalues \(\varepsilon_{n\mathbf{k}}^\text{DFT}\) of band \(n\) and crystal momentum \(\mathbf{k}\) of the first Brillouin zone are the DFT bandstructure.

The requirement 10 on Bloch functions \(\psi_{n\mathbf{k}}(\mathbf{r})\) can be fulfilled by the basis expansion \[\begin{align} \psi_{n\mathbf{k}}(\mathbf{r}) = \sum_\mu C_{ \mu n}(\mathbf{k})\sum_{\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}\phi^\mathbf{R}_\mu(\mathbf{r})\,,\label{kpointsbf} \end{align}\tag{11}\] where \(C_{\mu n}(\mathbf{k})\) are molecular orbital (MO) coefficients and \(\phi_\mu^\mathbf{R}(\mathbf{r})\) are atom-centered Gaussian-type basis functions centered on an atom in the cell with lattice vector \(\mathbf{R}\). For computing the MO coefficients \(C_{\mu n}(\mathbf{k})\), one inserts Eq. 11 into Eq. 9 , multiplies with an atom-centered Gaussian function \(\phi_\nu^\mathbf{0}(\mathbf{r})\) in the unit cell \(\mathbf{0}\), and integrates over the whole real space, which gives \[\begin{align} \sum_\nu h_{\mu\nu} (\mathbf{k})\, C_{ \nu n}(\mathbf{k}) = \sum_\nu S_{\mu\nu}(\mathbf{k}) \, C_{ \nu n}(\mathbf{k}) \, \varepsilon_{n\mathbf{k}} ^{\operatorname{DFT}}\,,\label{e18a} \end{align}\tag{12}\] with the Kohn-Sham matrix \[\begin{align} \begin{aligned} & h_{\mu\nu} (\mathbf{k})= \sum_\mathbf{R}\,e^{i\mathbf{k}\cdot\mathbf{R}}\; h_{\mu\nu}^\mathbf{R}\;,\label{e18} \\[0.3em] &h_{\mu\nu}^\mathbf{R}= \int d\mathbf{r}\;\phi_\mu^\mathbf{0}(\mathbf{r}) \left[h_0(\mathbf{r})+v^\text{xc}(\mathbf{r})\right]\phi_\nu^\mathbf{R}(\mathbf{r}) \,, \end{aligned} \end{align}\tag{13}\] and the overlap matrix \[\begin{align} &S_{\mu\nu} (\mathbf{k}) = \sum_\mathbf{R}\,e^{i\mathbf{k}\cdot\mathbf{R}} \; S_{\mu\nu}^\mathbf{R}\;,S_{\mu\nu}^\mathbf{R}= \int d\mathbf{r}\;\phi_\mu^\mathbf{0}(\mathbf{r}) \,\phi_\nu^\mathbf{R}(\mathbf{r})\,. \label{e19} \end{align}\tag{14}\] Note that the sums over lattice vectors \(\mathbf{R}\) in Eqs. 13 and 14 can be restricted to \(\mathbf{R}\) with small \(|\mathbf{R}|\) because the atom-centered Gaussian function \(\phi_\mu^\mathbf{R}(\mathbf{r}){\equiv}\phi_\mu^\mathbf{0}(\mathbf{r}{-}\mathbf{R})\) quickly decays for large \(|\mathbf{r}{-}\mathbf{R}|\).

The Kohn-Sham matrix \(h_{\mu\nu} (\mathbf{k})\) depends on the electron density \(n(\mathbf{r})\), \[\begin{align} n(\mathbf{r}) &= \sum_i^\text{occ} \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\; |\psi_{i\mathbf{k}}(\mathbf{r})|^2 \,,\label{e21a} \end{align}\tag{15}\] where we integrate over the crystal momentum \(\mathbf{k}\) in the BZ and \(\Omega_\text{BZ}\) denotes the volume of the BZ. To obtain an efficient algorithm for computing \(n(\mathbf{r})\), we use Eq. 11 to arrive at \[\begin{align} n(\mathbf{r}) = \sum_{\mu\nu} \sum_{\mathbf{R}_1\mathbf{R}_2} D_{\mu\nu}^{\mathbf{R}_2-\mathbf{R}_1} \;\phi_\mu^{\mathbf{R}_1}(\mathbf{r})\,\phi_\nu^{\mathbf{R}_2}(\mathbf{r}) \label{e22a} \end{align}\tag{16}\] using the density matrix \[\begin{align} &D_{\mu\nu}^{\mathbf{R}} = \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\; e^{-i\mathbf{k}\cdot\mathbf{R}}\;D_{\mu\nu}(\mathbf{k})\,, \tag{17} \\[0.5em] & D_{\mu\nu}(\mathbf{k}) = \sum_n^\text{occ} C_{\mu n} (\mathbf{k})\,C_{\nu n}^*(\mathbf{k})\,.\tag{18} \end{align}\] The integration over the BZ in Eq. 17 is executed in every self-consistent field cycle of the KS-DFT calculation using a discrete \(N_1{\times}N_2{\times}N_3\) Monkhorst-Pack \(k\)-point mesh \(\{\mathbf{k}_\ell\}\) [65] for two-dimensional periodicity: \[\begin{align} &D_{\mu\nu}^{\mathbf{R}} \simeq \frac{1}{N_1N_2N_3} \sum_{\mathbf{k}_\ell}^\text{BZ} e^{-i\mathbf{k}_\ell\cdot\mathbf{R}} \;\sum_n^\text{occ} C_{\mu n}(\mathbf{k}_\ell)\,C_{\nu n}^*(\mathbf{k}_\ell)\,.\label{e24} \end{align}\tag{19}\] For a periodic direction \(j\), we choose \(N_j\) as even integers, which leads to a \(k\)-mesh that excludes the \(\Gamma\)-point, \[\begin{align} \mathbf{k}_\ell = \sum_{j=1}^{d} \frac{ \ell_j}{2N_j}\,\mathbf{b}_j\;, \label{kmesh} \end{align}\tag{20}\] where we define \(\ell{=}\{\ell_j\}_{j=1}^{d}\) and each \(\ell_j\) takes one of the following odd integers \[\begin{align} \ell_j\in \{\pm\,1,\, \pm\,3, \;\ldots\;,\, \pm\, (N_j-1)\} \,.\label{ell} \end{align}\tag{21}\] \(\mathbf{b}_j\) are the primitive translation vectors of the reciprocal lattice that fulfill \(\mathbf{a}_{j_1}{\cdot}\mathbf{b}_{j_2}{=}2\pi\delta_{j_1j_2}\).

Note that the density matrix \(D_{\mu\nu}^{\mathbf{R}}\) computed from Eq. 19 features an erroneous periodicity in a superlattice with primitive translation vectors \(\{\mathbf{T}_j\}_{j=1}^d\), \[\begin{align} D_{\mu\nu}^{\mathbf{R}+\mathbf{T}} = D_{\mu\nu}^{\mathbf{R}} \;,\mathbf{T}= \sum_{j=1}^d t_j \mathbf{T}_j \;,\mathbf{T}_j = 2 N_j\, \mathbf{a}_j\;, \end{align}\] where \(\{t_j\}_{j=1}^d\) are integers. To avoid issues in a practical KS-DFT calculation, we restrict the lattice vectors \(\mathbf{R}\) in Eq. 19 to a single supercell (SC) of the \(\mathbf{T}\)-superlattice, \[\begin{align} \mathbf{R}\in \text{SC}\Leftrightarrow\mathbf{R}= \sum_{j=1}^d n_j\,\mathbf{a}_j \;, n_j\in \{0\,,\,\pm\,1\,,\,\ldots\,,\,\pm\,N_j\}\,.\label{e28} \end{align}\tag{22}\] This restricts the lattice summation in Eq. 16 to \({\mathbf{R}_1}{-}{\mathbf{R}_2}{\in}\text{SC}\). In the limit of a dense \(k\)-point mesh, i.e., large \(N_1,N_2,N_3\), the SC 22 is large and the lattice summation 16 converges, because the overlap \(\phi_\mu^{\mathbf{R}_1}(\mathbf{r})\,\phi_\nu^{\mathbf{R}_2}(\mathbf{r})\) quickly decays for large \(|\mathbf{R}_1{-}\mathbf{R}_2|\).

Starting from the eigenvalues \(\varepsilon_{n\mathbf{k}} ^{\operatorname{DFT}}\) from Eq. 12 , we add spin-orbit coupling (SOC) \(V^{\operatorname{SOC}}\) from Hartwigsen-Goedecker-Hutter (HGH) pseudopotentials [54]–[56], see details in Appendix 17, to obtain the spin-orbit perturbed Hamiltonian: \[\begin{align} h^{\operatorname{DFT+SOC}}_{n\sigma,\,n'\sigma'}(\mathbf{k}) = \delta_{nn'}\,\delta_{\sigma\sigma'}\,\varepsilon_{n\mathbf{k}}^{\operatorname{DFT}} + V_{nn',\sigma\sigma'~}^{\operatorname{SOC}}(\mathbf{k})\,. \end{align}\] We diagonalize the Hamiltonian with SOC to obtain the DFT band structure with SOC: \[\begin{align} \sum_{n'\sigma'} h^{\operatorname{DFT+SOC}}_{n\sigma,\,n'\sigma'}(\mathbf{k})\; C_{n'\sigma'}^{(j)}(\mathbf{k}) = \varepsilon^{\operatorname{DFT+SOC}}_{j\mathbf{k}} C_{n \sigma }^{(j)}(\mathbf{k})\,.\label{e31} \end{align}\tag{23}\]

4 GW space-time method with Gaussian basis functions and periodic boundary conditions↩︎

In this section, we reformulate the \(GW\) space-time method shown in Sec. 2 in an atomic-orbital basis with periodic boundary conditions. The starting point is a DFT calculation with periodic boundary conditions in an atomic-orbital basis (Sec. 3). In the main text, we only present the working equations that are implemented in the algorithm. We give a detailed derivation of these equations in the Appendix.

As shown in Appendix 11, the Green’s function in imaginary time is given by [49] \[\begin{align} \begin{aligned} G_{\mu \nu}(i\tau,\mathbf{k}) =& \; \theta(\tau) \sum_a^\text{empty} C_{ \mu a}(\mathbf{k})C_{ \nu a}^*(\mathbf{k})e^{-(\varepsilon_{a\mathbf{k}}-\varepsilon_\text{F})\tau} \\ & - \theta(-\tau) \sum_i^\text{occ} C_{ \mu i}(\mathbf{k})C_{ \nu i}^*(\mathbf{k})e^{-(\varepsilon_{i\mathbf{k}}-\varepsilon_\text{F})\tau}\,. \end{aligned} \label{e29} \end{align}\tag{24}\] Note that Eq. 24 is the analogue to Eq. 18 for the density matrix. As further shown in Appendix 11, the Green’s function is transformed to real space via an integration over the BZ, \[\begin{align} G^\mathbf{R}_{\mu \nu}(i\tau) &= \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\;e^{-i\mathbf{k}\cdot\mathbf{R}}\,G_{\mu \nu}(i\tau,\mathbf{k}) \\ &{\color{black} \simeq \frac{1}{N_1N_2N_3}\sum_{\mathbf{k}_\ell}^\text{BZ}e^{-i\mathbf{k}_\ell\cdot\mathbf{R}}\,G_{\mu \nu}(i\tau,\mathbf{k}_\ell)} \,.\label{e21} \end{align}\tag{25}\] The \(\mathbf{k}_\ell\)-mesh used in Eq. 25 is the odd Monkhorst-Pack mesh from the DFT calculation used for the density matrix, Eqs. 19  – 21 .

For the computation of the density response \(\chi\), we introduce the Resolution of Identity (RI) [66], where products of Gaussian basis functions are expanded over an auxiliary RI basis set \(\{\varphi_P^\mathbf{P}\}\): \[\begin{align} \phi_\mu^{\mathbf{R}}(\mathbf{r})\,\phi_\nu^{\mathbf{T}}(\mathbf{r}) \simeq \sum_{P\mathbf{P}} B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}}\varphi_P^\mathbf{P}(\mathbf{r})\,,\label{eRI} \end{align}\tag{26}\] where the projection coefficients are defined as \[\begin{align} B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}} =\sum_{Q\mathbf{Q}} (\mu\mathbf{R}\,\nu\mathbf{T}\,|\,Q\mathbf{Q})\,(M^{-1})_{PQ}^{\mathbf{P}-\mathbf{Q}}\,,\label{epRI} \end{align}\tag{27}\] with the three-centre integrals \[\begin{align} &(\mu\mathbf{R}\,\nu\mathbf{T}\,|\,P\mathbf{P}) = {\int} d\mathbf{r}\,d\mathbf{r}'\, \phi_\mu^{\mathbf{R}}(\mathbf{r})\, \phi_\nu^{\mathbf{T}}(\mathbf{r})\, V_{r_\text{c}}(\mathbf{r}-\mathbf{r}')\, \varphi_P^{\mathbf{P}}(\mathbf{r}') \label{e5c} \end{align}\tag{28}\] that can be understood as projection of the product \(\phi_\mu^{\mathbf{R}_1} \phi_\nu^{\mathbf{R}_2}\) onto the RI basis function \(\varphi_P^\mathbf{R}\) under a metric, which we have chosen to be the truncated Coulomb operator (RI-tCm) [67]–[69], \[\begin{align} & V_{r_\text{c}}(\mathbf{r}-\mathbf{r}')=\left\{ \begin{array}{cl} 1/|\mathbf{r}-\mathbf{r}'| & \text{if } |\mathbf{r}-\mathbf{r}'|\le r_\text{c}\,, \\[0.5em] 0 & \text{else\,.} \end{array} \right. \label{e6c} \end{align}\tag{29}\] The inverse metric matrix \((M^{-1})_{PQ}^{\mathbf{P}-\mathbf{Q}}\) in Eq. 27 arises because of the non-orthogonality of the auxiliary RI basis set and is computed from \[\begin{align} M_{PQ}(\mathbf{k}) =& \sum_\mathbf{R}e^{i\mathbf{k}\cdot\mathbf{R}}\int d\mathbf{r}\,d\mathbf{r}'\, \varphi_P^\mathbf{0}(\mathbf{r}) \, V_{r_\text{c}}(\mathbf{r}-\mathbf{r}')\,\varphi_Q^\mathbf{R}(\mathbf{r}')\,, \label{e11a} \end{align}\tag{30}\] a subsequent matrix inversion and transformation to real space, as shown in detail in Appendix [app:RI] which contains a derivation of Eq. 26 /27 .

The cutoff radius \(r_\text{c}\) from Eq. 29 is typically chosen in the order of a few Ångstroms [47], [68], [69]. In the limit \(r_\text{c}{\rightarrow}0\), RI-tCm is equivalent to the RI with the overlap metric which suffers from a slow convergence with respect to the RI basis set size [66], [67]. In the limit of \(r_\text{c}{\rightarrow}\infty\), RI-tCM is equivalent to the RI with the Coulomb metric where the convergence with the RI basis set size is fast. The drawback of the Coulomb metric is that tensor elements \((\mu\mathbf{R}\,\nu\mathbf{T}\,|\,P\mathbf{P})\) from Eq. 28 with the Coulomb operator only decay polynomially for \(|\mathbf{R}{-}\mathbf{P}|{\rightarrow}\infty\) and \(|\mathbf{T}{-}\mathbf{P}|{\rightarrow}\infty\) which prohibits the application of RI with the Coulomb metric in this algorithm. When truncating the Coulomb operator as in Eq. 29 at finite \(r_\text{c}\) and when using Gaussian basis functions, the tensor \((\mu\mathbf{R}\,\nu\mathbf{T}\,|\,P\mathbf{P})\) decays like a Gaussian for \(|\mathbf{R}{-}\mathbf{P}|{\rightarrow}\infty\) and \(|\mathbf{T}{-}\mathbf{P}|{\rightarrow}\infty\). It has been shown that RI-tCm converges quickly with the size of the auxiliary basis \(\{\varphi_P^\mathbf{R}\}\). [67], [68]

The density response \(\chi_{PQ}^\mathbf{R}(i\tau){=}\braket{\varphi^\mathbf{0}_P|\chi(i\tau)|\varphi^\mathbf{R}_Q}\) in imaginary time in the RI basis \(\{\varphi^\mathbf{R}_P\}\) can be obtained as (see Appendix 11 for a derivation and Appendix 12 for parallel implementation) \[\begin{align} \chi_{PQ}^\mathbf{R}(i\tau) =& \sum_{\lambda\mathbf{R}_1\nu\mathbf{R}_2} \bigg[\sum_{\mu } \sum_{\mathbf{S}_1}^\text{SC} (\mu\mathbf{R}_1{-}\mathbf{S}_1\,\nu\mathbf{R}_2\,|\,P\mathbf{0}) \; G^{\mathbf{S}_1}_{ \lambda\mu}(-i\tau)\bigg] \nonumber \\ & \times \bigg[\sum_{\sigma}\sum_{\mathbf{S}_2}^\text{SC} (\lambda\mathbf{R}_1\,\sigma\mathbf{R}_2{-}\mathbf{S}_2\,|\,Q\mathbf{R})\; G^{\mathbf{S}_2}_{\nu \sigma}(i\tau)\bigg]\,, \label{chiT} \end{align}\tag{31}\] where \(\mathbf{R}\) is a lattice vector inside the SC [Eq. 22 ], \(\mathbf{R}_1,\mathbf{R}_2,\mathbf{S}_1,\mathbf{S}_2\) are lattice vectors. Note that we only sum over cells \(\mathbf{R}_1\) and \(\mathbf{R}_2\) with small \(|\mathbf{R}_1|,|\mathbf{R}_2|\) because of the sparsity of the three-centre integrals 28 .

Following the \(GW\) space-time method, [47], [48] we transform the polarizability from real space and time to the Brillouin zone and frequency \[\begin{align} \chi_{PQ}(\mathbf{k},i\omega) =\sum_\mathbf{R}^\text{SC} \int d\tau \, \cos(\omega\tau) \,e^{i\mathbf{k}\cdot\mathbf{R}}\; \chi_{PQ}^\mathbf{R}(i\tau)\,.\label{chitrafo} \end{align}\tag{32}\] Note that as \(\chi^\mathbf{R}\) computed from Eq. 31 decays exponentially for large \(\mathbf{R}\) in gapped systems, so one can perform the Fourier transformation 32 to an arbitrary \(\mathbf{k}\). We execute the imaginary-time integration numerically using minimax grids [49], [70], [71].

As next step, we calculate the dielectric function [47] \[\begin{align} {\boldsymbol{\epsilon}}(\mathbf{k},i\omega) = \mathbf{Id}- \mathbf{V}^{0.5}(\mathbf{k}) \mathbf{M}^{-1}(\mathbf{k})\boldsymbol{\chi}(\mathbf{k},i\omega) \mathbf{M}^{-1}(\mathbf{k})\mathbf{V}^{0.5}(\mathbf{k}) \label{e22} \end{align}\tag{33}\] where \(\mathbf{Id}\) is the identity matrix and the truncated Coulomb matrix \(\mathbf{M}(\mathbf{k})\) defined in Eq. 30 appears due to the RI-tCm. We use Tikhonov regularization [72] for the RI expansion to prevent linear dependencies of fit coefficients, as discussed in the supporting information of Ref. [47]. The regularization leads to the modified matrix inversion, \[\begin{align} \mathbf{M}^{-1}(\mathbf{k}) = \big(\mathbf{M}(\mathbf{k}){+}\alpha\mathbf{Id}\big)^{-1}\,,\label{e12} \end{align}\tag{34}\] where \(\alpha\) is the regularization parameter.

\(\mathbf{V}^{0.5}(\mathbf{k})\) in Eq. 33 is the matrix square root of the bare Coulomb interaction \(\mathbf{V}(\mathbf{k})\) [44], [73]–[75] \[\begin{align} V_{PQ}(\mathbf{k}) &= \sum_\mathbf{R}e^{i\mathbf{k}\cdot\mathbf{R}}{\int} d\mathbf{r}\,d\mathbf{r}'\,\varphi_P^\mathbf{0}(\mathbf{r})\,\frac{ 1}{|\mathbf{r}- \mathbf{r}'|}\,\varphi_Q^\mathbf{R}(\mathbf{r}')\,.\label{Vper} \end{align}\tag{35}\] Details on the lattice summation over \(\mathbf{R}\) are given in Appendix 13. We obtain the correlation part of the screened interaction \(W^\text{c}(i\omega) {=}(\epsilon^{-1}(i\omega){-}1)V\) as \[\begin{align} \mathbf{W}^\text{c}(\mathbf{k},i\omega) = \mathbf{V}^{0.5}(\mathbf{k})({\boldsymbol{\epsilon}}^{-1}(\mathbf{k},i\omega) -\mathbf{Id}) \mathbf{V}^{0.5}(\mathbf{k}) \label{Wk} \end{align}\tag{36}\] and transform it to real space \(W_{PQ}^{\text{c},\mathbf{R}}{=}\braket{\varphi^\mathbf{0}_P|W^\text{c}|\varphi^\mathbf{R}_Q}\), \[\begin{align} W_{PQ}^{\text{c},\mathbf{R}} (i\omega) = \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\; e^{-i\mathbf{k}\cdot\mathbf{R}}\,W^\text{c}_{PQ}(\mathbf{k},i\omega)\,.\label{trafoWtoR} \end{align}\tag{37}\] Special care is required for the BZ integral as \(W_{PQ}^\text{c}(\mathbf{k},i\omega)\) diverges at the \(\Gamma\)-point with \(1/k\) for two-dimensional materials if \(\varphi_P\) and \(\varphi_Q\) are s-type basis functions [43]–[45]. We evaluate \(W_{PQ}^\text{c}(\mathbf{k},i\omega)\) using two Monkhorst-Pack meshes [65]: \(\{\mathbf{k}_\ell\}\) has \(4N_j\) \(k\)-points in periodic directions \(j\) and \(\{\mathbf{q}_\ell\}\) has \(8N_j\) \(k\)-points in periodic directions \(j\). The number of \(k\)-points in the \(\{\mathbf{k}_\ell\}\) mesh and \(\{\mathbf{q}_\ell\}\) mesh is thus \[\begin{align} N_k = \prod_{j=1}^d 4N_j\;,N_q = \prod_{j=1}^d 8N_j\;. \end{align}\] We extrapolate the BZ integration 37 with the inverse square root of the number of \(k\)-points [44]. Reformulating Eq. 37 , the \(k\)-extrapolated screened Coulomb interaction becomes \[\begin{align} W_{PQ}^{\text{c},\mathbf{R}}(i\omega) &= \sum_{\ell} v_\ell\,e^{-i\mathbf{q}_\ell\cdot\mathbf{R}} \,W_{PQ}^\text{c}(\mathbf{q}_\ell,i\omega) \nonumber \\ & - \sum_{\ell} w_\ell\,e^{-i\mathbf{k}_\ell\cdot\mathbf{R}} \,W_{PQ}^\text{c}(\mathbf{k}_\ell,i\omega)\,, \label{e16a} \end{align}\tag{38}\] where the extrapolation is incorporated into the integration weights: \[\begin{align} v_\ell = \frac{1}{(1-\sqrt{N_k/N_q })\,N_q}\;\;,w_\ell = \frac{1}{({\sqrt{N_q/N_k }}-1)\,N_k}\;. \label{e17b} \end{align}\tag{39}\]

We transform \(W_{PQ}^{\text{c},\mathbf{R}}(i\omega)\) back to \(k\)-points, \[\begin{align} \hat{W}^\text{c}_{PQ}(\mathbf{k},i\omega) = \sum_\mathbf{R}^\text{SC}e^{i\mathbf{k}\cdot\mathbf{R}}\; W_{PQ}^{\text{c},\mathbf{R}}(i\omega) \,, \end{align}\] where we only sum over the cells \(\mathbf{R}{\in}\text{SC}\), see Eq. 22 , to prevent for the divergence of \(\hat{W}^\text{c}(\mathbf{k},i\omega)\) at \(\mathbf{k}{=}0\). We incorporate the RI metric, \[\begin{align} \tilde{\mathbf{W}}^\text{c}(\mathbf{k},i\omega) = \mathbf{M}^{-1}(\mathbf{k})\, \hat{\mathbf{W}}^\text{c}(\mathbf{k},i\omega)\,\mathbf{M}^{-1}(\mathbf{k})\,, \end{align}\] and transform to real space using the \(k\)-point mesh from DFT, Eqs. 19 and 20 \[\begin{align} \tilde{W}_{PQ}^{\text{c},\mathbf{R}} (i\omega) &= \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\; e^{-i\mathbf{k}\cdot\mathbf{R}}\,\tilde{W}^\text{c}_{PQ}(\mathbf{k},i\omega) \\ &{\color{black} \simeq \frac{1}{N_1N_2N_3}\sum_{\mathbf{k}}^{\textrm{BZ}}e^{-i\mathbf{k}_\ell\cdot\mathbf{R}}\,\tilde{W}^\text{c}_{PQ}(\mathbf{k}_\ell,i\omega)\,.}\label{e49} \end{align}\tag{40}\]

Following the \(GW\) space-time method [48], we transform \(W^\text{c}(i\omega)\) to imaginary time using minimax grids [49], [70], [71]. This completes the ingredients for the \(GW\) self-energy \(\Sigma^\text{c}(i\tau){=}iG(i\tau)W^\text{c}(i\tau)\). We calculate the self-energy \(\Sigma^{\text{c},\mathbf{R}}_{\lambda\sigma}(i\tau){=}\braket{\phi_\lambda^\mathbf{0}|\Sigma^\text{c}(i\tau)|\phi_\sigma^\mathbf{R}}\) in real space (derivation in Appendix 16), \[\begin{align} \Sigma^{\text{c},\mathbf{R}}_{\lambda\sigma}(i\tau) &= i \sum_{P\nu} \sum_{\mathbf{R}_1 }^\text{SC} \sum_{ \mathbf{S}_1}^\text{SC} \left[ \sum_\mu \sum_{ \mathbf{S}_2}^\text{SC}(\lambda\mathbf{0}\;\mu\mathbf{S}_1{-}\mathbf{S}_2 \,|\,P\mathbf{R}_1) \; G^{\mathbf{S}_2}_{\mu\nu}(i\tau)\right] \nonumber \\[0.2em] &\times \left[\sum_Q\sum_{\mathbf{R}_2 }^\text{SC} ( \sigma\mathbf{R}\;\nu\mathbf{S}_1 \,|\,Q\mathbf{R}_1{-}\mathbf{R}_2) \;\tilde{W}_{QP}^{\text{c},\mathbf{R}_2}(i\tau)\right] \label{SigmaR} \end{align}\tag{41}\] and in the Kohn-Sham basis \(\psi_{n\mathbf{k}}(\mathbf{r})\), \[\begin{align} &\Sigma_{\lambda\sigma}^\text{c}(\mathbf{k},i\tau) = \sum_\mathbf{R}^\text{SC} e^{i\mathbf{k}\cdot\mathbf{R}}\; \Sigma_{\lambda\sigma}^{\text{c},\mathbf{R}}(i\tau) \,,\tag{42} \\[0.3em] &\Sigma_{n\mathbf{k}}^\text{c}(i\tau) = \sum_{\lambda\sigma} C_{\lambda n}^*(\mathbf{k}) \, \Sigma^\text{c}_{\lambda\sigma}(\mathbf{k},i\tau)\,C_{\sigma n}(\mathbf{k})\;. \tag{43} \end{align}\]

A major advantage of the present \(GW\) algorithm is that the self-energy is computed in real space, Eq. 41 , involving only a modest number of neighbor cells \(\mathbf{R}_1,\mathbf{R}_2,\mathbf{S}_1,\mathbf{S}_2\) to be included. Instead, when using \(k\)-point sampling instead of lattice summation to compute the self-energy, special care needs to be taken due to the diverging \(W(\mathbf{k})\) at the \(\Gamma\)-point, which requires special correction schemes [76]. In the present algorithm, the divergence of \(W(\mathbf{k})\) at the \(\Gamma\)-point is taken into account via BZ integration of \(W(\mathbf{k})\) in Eq. 37 with suitable \(k\)-point grids, Eq. 38 .

We calculate the exchange self-energy similarly to previous work [69], [77], \[\begin{align} \Sigma^{\text{x},\mathbf{R}}_{\lambda\sigma}(i\tau) &= i \sum_{P\nu} \sum_{\mathbf{R}_1 }^\text{SC} \sum_{ \mathbf{S}_1}^\text{SC} \left[ \sum_\mu \sum_{ \mathbf{S}_2}^\text{SC}(\lambda\mathbf{0}\;\mu\mathbf{S}_1{-}\mathbf{S}_2 \,|\,P\mathbf{R}_1) \; D^{\mathbf{S}_2}_{\mu\nu}\right] \nonumber \\[0.2em] &\times \left[\sum_Q\sum_{\mathbf{R}_2 }^\text{SC} ( \sigma\mathbf{R}\;\nu\mathbf{S}_1 \,|\,Q\mathbf{R}_1{-}\mathbf{R}_2) \; \tilde{V}_{QP}^{\text{tr},\mathbf{R}_2} \right] \,, \label{SigmaxR} \end{align}\tag{44}\] where \(D^{\mathbf{S}}_{\mu\nu}\) is the density matrix 17 from KS-DFT and \[\begin{align} &\tilde{\mathbf{V}}^{\text{tr},\mathbf{R}} = \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\;e^{-i\mathbf{k}\cdot\mathbf{R}}\; \mathbf{M}^{-1}(\mathbf{k})\, \mathbf{V}^\text{tr}(\mathbf{k})\,\mathbf{M}^{-1}(\mathbf{k})\,, \\[0.5em] & V_{PQ}^\text{tr}(\mathbf{k}) = \sum_\mathbf{R}e^{i\mathbf{k}\cdot\mathbf{R}} \int d\mathbf{r}\,d\mathbf{r}'\, \varphi_P^\mathbf{0}(\mathbf{r}) \, V_{r_\text{SC}}(\mathbf{r},\mathbf{r}')\,\varphi_Q^\mathbf{R}(\mathbf{r}') \end{align}\] is the truncated Coulomb matrix with truncation radius \(r_\text{SC}{=}\min_j N_j|\mathbf{a}_j|/2\), i.e., half of the shortest primitive translation vector of the supercell SC 22  [77], [78]. Note that we restrict all lattice vector differences appearing in Eq. 41 and 44 to the SC 22 , i.e., \(\mathbf{S}_1{-}\mathbf{S}_2,\mathbf{R}_1{-}\mathbf{R}_2{\in}\text{SC}\).

We obtain the quasiparticle energies \(\varepsilon_{n\mathbf{k}}^{G_0W_0}\) by replacing the DFT exchange-correlation contribution \(v^\text{xc}_{n\mathbf{k}}\) with the self-energy, \[\begin{align} \varepsilon_{n\mathbf{k}}^{G_0W_0} = \varepsilon_{n\mathbf{k}}^\text{DFT}+ \Sigma^\text{x}_{n\mathbf{k}} + \text{Re}\,\Sigma^\text{c}_{n\mathbf{k}}(\varepsilon_{n\mathbf{k}}^{G_0W_0}) -v^\text{xc}_{n\mathbf{k}} \,.\label{qpeq} \end{align}\tag{45}\]

Identically to the DFT band structure 23 , we add the spin-orbit potential \(V^{\operatorname{SOC}}\) to the eigenenergies to obtain the spin-orbit corrected Hamiltonian (the details are given in Appendix 17): \[\begin{align} h^{{G_0W_0}\text{+SOC}}_{n\sigma,\,n'\sigma'}(\mathbf{k}) = \delta_{nn'}\,\delta_{\sigma\sigma'}\,\varepsilon_{n\mathbf{k}}^{G_0W_0} + V_{nn',\sigma\sigma'~}^{\operatorname{SOC}}(\mathbf{k})\,. \end{align}\] Diagonalization leads to the \(GW\) band structure with SOC: \[\begin{align} \sum_{n'\sigma'} h^{{G_0W_0}\text{+SOC}}_{n\sigma,\,n'\sigma'}(\mathbf{k})\; C_{n'\sigma'}^{(j)}(\mathbf{k}) = \varepsilon^{{G_0W_0}\text{+SOC}}_{j\mathbf{k}} C_{n \sigma }^{(j)}(\mathbf{k})\,. \label{e58} \end{align}\tag{46}\]

The primary difference between our algorithm and the \(\Gamma\)-point-only \(GW\) implementation using Gaussian basis functions reported in Ref. [47] lies in the treatment of periodic boundary conditions. In our approach, both the irreducible density response function \(\chi\) and the self-energy \(\Sigma\) are computed via explicit lattice summations [Eqs. 31 and 41 ] over all unit cells. In contrast, Ref. [47] restricts the lattice sum to the nearest-neighbor cell, effectively circumventing the need for a full lattice summation. This approach is only exact in the limit of large unit cells, and the presented \(GW\) algorithm enables the treatment of small unit cells.

5 Computational Details↩︎

5.1 Overview↩︎

There are several numerical approaches, also involving approximations, to compute \(GW\) band structures. Commonly employed approximations are the use of pseudopotentials to exclude core electrons from the computation, and the use of plasmon-pole models to simplify the frequency dependence in \(GW\). In this work, we employ pseudopotentials, details given in the following, but we avoid plasmon-pole models by treating the full frequency (and time) dependence.

The evaluation of intermediate quantities in the \(GW\) method, which depend on the real-space coordinates \(\mathbf{r}\) and \(\mathbf{r}'\), requires the use of a basis set. Real-space grids, while conceptually straightforward, typically involve a large number of grid points and are therefore computationally less efficient. Common alternatives include plane-wave and atomic-orbital basis sets. When using identical pseudopotentials—or when the influence of the pseudopotential is negligible—different basis sets yield the same \(GW\) band structure, provided the basis is sufficiently large to ensure convergence. In practice, agreement between different \(GW\) implementations is typically within 0.1 eV when the same pseudopotentials are used across codes [41], [79]. A comprehensive benchmark study comparing \(GW\) results across different basis sets (e.g., plane-waves vs.atom-centered orbitals) for a large and diverse set of materials—on the order of 100 solids—with agreement at the 10 meV level or better has yet to be conducted, and the very strong sensitivity of the band gaps with respect to parameters such as the lattice constant [80] makes comparison with already available theoretical literature fairly difficult. To assess the numerical precision of our \(GW\) implementation (Sec. 4), we carry out illustrative tests on four selected reference materials, namely monolayer MoS\(_\text{2}\), MoS\(\text{e}_\text{2}\), WS\(_\text{2}\) and WS\(\text{e}_\text{2}\). We use two well-established \(GW\) codes, BerkeleyGW [57] and VASP [12], to compute reference band gaps and band structures for comparison.

5.2 Atomic geometries of MoS\(_\text{2}\), MoS\(\text{e}_\text{2}\), WS\(_\text{2}\), WS\(\text{e}_\text{2}\)↩︎

For our benchmark calculations, we employ monolayer transition metal dichalcogenides MoS\(_\text{2}\), MoS\(\text{e}_\text{2}\), WS\(_\text{2}\), and WS\(\text{e}_\text{2}\). These materials are non-magnetic and stable, and they attract widespread interest thanks to a rare combination of properties: they are atomically thin, have a direct bandgap and strong spin–orbit coupling which make them ideal for both fundamental studies and emerging applications in electronics, spintronics, optoelectronics and energy harvesting [81]. We take atomic geometries and lattice parameters for these materials from the C2DB database  [9].

5.3 GW space-time calculations (CP2K)↩︎

We have implemented the low-scaling \(GW\) space-time algorithm presented in this work in the CP2K software package [82]–[84]. CP2K employs a Gaussian basis set for representing KS orbitals [Eq. 11 ] and a plane-waves basis set for the electron density to evaluate the Hartree potential via Ewald summation [82], [83]. We use Gaussian dual-space pseudopotentials [54]. In the DFT calculation, we employ the PBE exchange-correlation functional [85]. The plane-wave cutoff for the electron density is set to \(500\) Ry. This value was converged beforehand on the DFT-level results, as the plane-wave grid is not used for the \(GW\) part of the calculations.

In the \(GW\) algorithm, we employ a minimax time-frequency grid [49], [70], [71]. We compute two- and three-centre integrals over Gaussians using analytical schemes [86], [87]. The self-energy is analytically continued from imaginary frequency to the real frequency using a Padé model [6], [60], [88] with 16 parameters. Unless otherwise noted, we employ a cutoff radius \(r_c {=}7\) Å, an RI regularization \(\alpha{=}10^{-2}\) and a box height for the 2D materials of 15 Å (for computing the Fourier transform in the Hartree potential). As already stated, the extrapolation of the \(k\)-integration of \(W\) [Eq. 39 ] relies on two \(k\)-meshes, one 4 times denser than the equivalent DFT \(k\)-mesh along each direction and one 8 times denser (e.g. for a \(32{\times}32\) DFT \(k\)-mesh, the corresponding coarse W \(k\)-mesh is \(128{\times}128\) and the dense \(k\)-mesh is \(256{\times}256\)). The remaining computational parameters include the number of minimax time and frequency points \(N_{\tau/\omega}\) [49], [70], [71], the DFT \(k\)-mesh [Eqs. 15 , 17 , 24 ; the DFT \(k\)-mesh also defines the SC supercell 22 ], and the filter threshold for sparse matrix-tensor operations in Eqs. 31 , 41 and 44 . In this paper, we will use two sets of parameters: tight settings (\(N_{\tau/\omega}{=}30\), DFT \(k\)-mesh: \(32{\times}32\), filter: \(10^{-11}\)) which correspond to a reference set of parameters that we have defined through extensive testing in order to yield well-converged \(GW\) band gaps, and a set of light settings (\(N_{\tau/\omega}{=}10\), DFT \(k\)-mesh: \(24{\times}24\), filter: \(10^{-6}\)) that have sufficiently reduced memory and computational costs in order to be used for laptop calculations while still giving decently accurate results (quantitative values are given in Sec. 6 and Sec. 8). Note that the memory bottleneck comes from the parallelization strategy, as one can see in Appendix 12, which implies the storage of all the three-centre integrals on each parallel group. In theory, one could also recompute these at each loop of the code, which would in practice prohibitively increase the computational time.

We use the single-, double- and triple-zeta MOLOPT basis sets [89] for expanding the KS orbitals, Eq. 11 . These basis sets have been optimized for the total energy of the ground state so that they might exhibit a slow convergence behaviour for excited states, and therefore for band gaps at the \(GW\) level. This motivated us to also use augmented single-, double- and triple-zeta Gaussian basis sets (aug-SZV-MOLOPT, aug-DZVP-MOLOPT and aug-TZVP-MOLOPT). We have created those by augmenting Gaussian SZV-MOLOPT, DZVP-MOLOPT and TZVP-MOLOPT bases [89] of S, Se, Mo, W with an additional \(s\), \(p\), \(d\) function, an additional \(f\) function (for Mo, W and aug-DZVP-MOLOPT and aug-TZVP-MOLOPT of S, Se) and an additional \(g\) function (all aug-TZVP-MOLOPT and aug-DZVP-MOLOPT of Mo, W). We have optimized the contraction coefficients of the additional functions by optimizing the lowest five \(GW\)+Bethe-Salpeter excitation energies [90] of a molecular set [89]. We will report DFT and \(GW\) with the original, i.e., non-augmented MOLOPT basis sets and the augmented MOLOPT basis sets.

For the density response 31 , dielectric function 33 and screened Coulomb interaction 36 , an auxiliary Gaussian RI basis set is required. There is no general approach for constructing optimally sized RI basis sets, as their design is closely tied to the chosen AO basis. In the case of the widely used cc-pVNZ Dunning basis sets [91], for example, the corresponding RI basis sets introduced in Ref. [92] provide a single, fixed RI basis for each AO basis. As a result, convergence with respect to the size of the RI basis is not commonly investigated. It is also possible to generate RI basis sets on the fly [93], although this usually leads to fairly large numbers of basis functions. In our case, we have optimized the RI basis set \(\{\varphi_P\}\) by matching the RI-MP2 correlation energy [94] of single atoms to the MP2 correlation energy. Unless stated otherwise, all reference calculations in this paper were carried using an RI basis with a relative RI-MP2 correlation energy difference of \(10^{-3}\) with respect to the corresponding MP2 correlation energy.

We incorporate SOC via parameters from HGH Gaussian dual-space pseudopotentials [54]–[56], see Eqs. 23 , 46 and Appendix 17 for details. We compute the SOC for states in a window of 40 eV, so 20 eV below the valence band maximum to 20 eV above the conduction band minimum, in order to avoid possible numerical instabilities with bands far from the gap. For the case of WSe\(_2\), we chose a window to 20 eV for aug-DZVP-MOLOPT and aug-TZVP-MOLOPT, as discussed in Appendix 17.

All inputs and outputs of the calculations are openly available, see data and code availability statement.

5.4 Reference GW calculations with a plane-wave-based algorithm (BerkeleyGW)↩︎

Quantum Espresso.

We performed the Quantum Espresso (QE) [95] DFT calculations employing the PBE exchange-correlation functional [85]. Fully relativistic norm-conserving pseudopotentials were used, as provided by the PseudoDojo database [96]. A plane-wave energy cut-off of 100 Ry was applied, and the self-consistent charge density was converged on a \(30 \times 30 \times 1\) \(\mathbf{k}\)-grid with a total energy convergence threshold of \(10^{-9}\) Ry.

BerkeleyGW.

Using the QE DFT energies and states, we performed for each material, a one-shot \(GW\) calculation (\(G_0W_0\)) using the BerkeleyGW package [57], [97]. We considered the full spinor implementation of BerkeleyGW [98], which incorporates SOC non-perturbatively. The dielectric matrix was computed with a dielectric cut-off of 25 Ry, considering a total of 3999 occupied and empty bands on a \(12 \times 12 \times 1\) uniform \(\mathbf{k}\)-grid. For completeness, we computed the quasi-particle band-gap using the generalized plasmon-pole model of Hybertsen-Louie [97] (see Appendix 18) and the full-frequency evaluation of the self-energy. In the explicit full frequency calculation, we used the contour-deformation method with the Adler-Wiser formula. We employed a frequency step and broadening of 0.25 eV, using 15 frequency points along the imaginary axis within the contour deformation approach. A low-frequency cutoff of 20.0 eV was set to restrict the real-axis integration range. To accelerate convergence in the vicinity of \(|\mathbf{q}| \rightarrow 0\), a non-uniform neck subsampling approach was employed [99] and the spurious interactions between periodic replicas in the perpendicular direction to the surface were removed with a Coulomb interaction truncation scheme [100]. The full frequency dependence of self-energy was evaluated with a frequency step of \(0.2\) eV in the frequency range \([-2.0, 2.0]\) eV and centering the frequency grid around each mean-field quasi-particle energy. For the band structures, the quasi-particle energies computed on the coarse \(\mathbf{k}\)-grid were interpolated along high-symmetry \(\mathbf{k}\)-paths.

5.5 Reference GW calculations with the projector-augmented-wave scheme (VASP)↩︎

To obtain the quasi-particle band gap, we carry out a single-shot \(G_0W_0\) calculation using the Vienna Ab initio Simulation Package (VASP) [101], [102]. In this process, the initial KS wavefunctions and energy levels, derived from a preceding DFT calculation also performed with VASP, are used. Core electrons are treated using the \(GW\) adaptation of the projector-augmented-wave (PAW) pseudopotentials [103], [104]. For all TMDCs, we apply PAWs constructed on the Perdew–Burke–Ernzerhof (PBE) functional [85] with an energy cutoff of 500 eV. The considered valence electron configurations are \(4s^24p^65s^14d^5\) for Mo, \(3s^23p^4\) for S, \(5p^66s^25d^4\) for W, \(4s^24p^4\) for Se. A total of 384 bands, comprising both occupied and unoccupied states, are taken into account, along with a uniform \(\mathbf{k}\)-grid of \(18 \times 18 \times 1\) for MoS\(_2\) and WS\(_2\), and of \(15 \times 15 \times 1\) for MoSe\(_2\) and WSe\(_2\), ensuring smooth convergence in the vicinity of \(|\mathbf{q}| \rightarrow 0\). Moreover, the dielectric tensor is broadened with a Lorentzian of 0.01 eV in all cases except for MoSe\(_2\), where a broadening of 0.1 eV is applied. Finally, the self-energy is calculated using a full-frequency implementation with 100 points along the imaginary frequency grid, as provided in VASP.

6 DFT band structure of \(\text{MoS}_\text{2}\), \(\text{MoSe}_\text{2}\), WS\(_\text{2}\), WS\(\text{e}_\text{2}\)↩︎

We begin by analyzing the DFT band structures and band gaps of monolayer \(\text{MoS}_2\), \(\text{MoSe}_2\), \(\text{WS}_2\), and \(\text{WSe}_2\), with particular focus on the agreement between the three numerical approaches described in Sec. 5, namely Gaussian basis sets versus plane waves, and different treatments of core electrons via pseudopotentials. This comparison is important because \(G_0W_0\) band structures are computed on top of the underlying Kohn-Sham DFT results (see, for example, Eq. 45 ). Since discrepancies in the \(GW\) band structures computed from different numerical approaches are expected to be larger than those at the DFT level, close agreement among the DFT band structures is a necessary prerequisite for a reliable \(GW\) benchmark.

We show the DFT band structure computed with the PBE exchange-correlation functional [85] and SOC in Fig. 1, for the aug-SZV-MOLOPT basis set (light settings) and the aug-DZVP-MOLOPT (tight settings), and the band structure computed from QE (with 8 empty bands). We observe excellent agreement between the aug-DZVP-MOLOPT calculation and the QE calculation, with a difference of 18 meV on average between their respective PBE+SOC bandgaps (see Table [t1]), and also with VASP results with a difference of 10 meV. The agreement of the PBE+SOC gap computed with the small aug-SZV-MOLOPT basis is also good, the average deviation is 29 meV to QE and 23 meV to VASP (Table [t1]). This shows a good agreement of DFT band structures between each code, which validates their use as a starting point for a \(GW\) benchmark. Table [t1] also reports the DFT band gaps without SOC and the SOC splitting at the K-point, demonstrating that our approach yields SOC splittings in good agreement with plane-wave reference calculations.

For comparison, we also computed the band gaps using the original, non-augmented MOLOPT basis sets [89] (see Table [t1]). At the DFT level, the results show good agreement with those from augmented basis sets—except for SZV-MOLOPT—indicating that augmentation is not strictly required for fast convergence of the DFT band gap with respect to basis set size.

To assess the accuracy of the SOC implementation, we report in Table [t1] the spin-orbit splitting of the valence band maximum at the K-point. This splitting is crucial as it is involved in determining the (optical) energy difference between A and B excitons, observed in reflectance and photoluminescence spectra [105], [106]. It thus plays a crucial role for valley selective optical excitations. We observe good agreement between our perturbative SOC implementation and the fully relativistic implementation in QE, with an average difference of 17 meV between the aug-DZVP-MOLOPT and QE results, as an example. The difference is more significant for the selenium-based TMDs, with 32 meV on average whereas the difference is 2 meV for the sulfur-based TMDs. The results are similar for the comparison with the VASP SOC-splitting. This finding validates our perturbative SOC treatment from HGH pseudopotentials [Eqs. 23 , 46 and Appendix 17].

We demonstrate the impact of SOC on the band structure by giving in Appendix 17 the comparison between the PBE and PBE+SOC bandstructures for all monolayers (Fig. 10). One can therefore see that the SOC lifts the spin degeneracy and therefore splits the band structure, especially at the K-point as we have already discussed, This leads to a reduction of the band gap with respect to the calculation without SOC, especially in the case of WS\(_2\) and WSe\(_2\) where the band gap is lowered by 0.3 eV (Table [t1]).

7 GW band gap of WS\(\text{e}_\text{2}\): Convergence study↩︎

In this section, we analyze the convergence of the \(GW\) bandgap of monolayer WSe\(_2\) with the numerical parameters summarized in Sec. 5.3. We focus on the direct band gap of WSe\(_2\) at K, calculated using \(G_0W_0\)@PBE including SOC. In Appendix 20, we provide a similar benchmark for MoS\(_2\).

Using our tight convergence settings (Sec. 5.3), we obtain a \(G_0W_0\)@PBE+SOC band gap of 1.95 eV with the aug-SZV-MOLOPT basis and 1.93 eV with the aug-DZVP-MOLOPT basis (Table [t1]). A larger augmented basis set is not used due to current memory limitations in our implementation. Fig. 2 shows how the band gap changes with various computational parameters. In Fig. 2a, we test different sizes of the time-frequency grid. We find that the gap changes by up to 12 meV in \(G_0W_0\)@PBE (blue trace), when using 14, 20, 26, 30 and 34 points. These variations can arise by poles in the self-energy [112], which is an unphysical feature of \(G_0W_0\)@PBE and eigenvalue self-consistency in \(G\) will cure this deficiency [111], [112]. We apply Hedin’s shift [1], [3], [107]–[110] to approximate eigenvalue self-consistency in \(G\), which reduces the variation between 14 and 34 time-frequency points to 3 meV (red curve in Fig. 1a).

The band gap is also converged with respect to the \(k\)-point mesh from DFT (Eq. 20 , Fig. 2c), changing by less than 10 meV between a \(24{\times}24\) and \(48{\times}48\) \(k\)-point meshes. The filter threshold for three-centre integrals 28 decides about removing small three-centre integrals from the calculation. Thus, decreasing the filter threshold increases the numerical precision and we demonstrate in Fig. 2e that a filter threshold of \(10^{-11}\) is sufficient for sub-meV convergence.

Note that the number of cells included in the lattice sums 31 , 41 and 44 depend on both the filter threshold of three-centre integrals (Fig. 2e) and the supercell "SC" defined in Eq. 22 . The size of the supercell "SC" is determined by the DFT \(k\)-mesh [Eq. 22 ]. We thus demonstrate convergence w.r.t. the number of cells in the lattice sums 31 , 41 and 44 by demonstrating convergence with the \(k\)-mesh (Fig. 2c) and the filter threshold (Fig. 2e). Similarly, we show the fast convergence of the lattice sum 35 of the Coulomb matrix in Fig. 2i (the size factor of the lattice summation is defined in Appendix 13). Also, the box height (Fig. 2g) can be well-converged.

The computational parameters with the most impact on the computation time are the \(k\)-mesh (Fig. 2d) and the basis set (Fig. 2l). When applying the present \(GW\) algorithm to other materials, we recommend to employ safe numerical parameters for the time-frequency integration (30 points), filter (10\(^{-12}\)) and box height (15 Å), and to focus convergence tests on the \(k\)-mesh and the basis set.

We also study the convergence of the \(GW\) gap with respect to the RI basis size in Fig. 3. Without regularization \((\alpha{=}0\)), the band gap increases with RI basis size (Fig. 3a). We assign this issue to linear dependencies in the RI basis set, which leads to numerical instabilities in the RI basis expansion, see Appendix [app:RI], Eq. 48 . In a nutshell, numerical instabilities arise when two spatially close diffuse \(s\)-type RI function can partially compensate each other, leading to large expansion coefficients \(B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}}\) in Eq. 48 with alternating sign.

To fix this issue, we use RI regularization [47] via \(\alpha{=}10^{-3}\) and \(\alpha{=}10^{-2}\) in Eq. 34 . RI regularization ensures numerical stability, see Fig. 3b,c: For an RI basis set size of a single WSe\(_2\) unit between 130 and 187 and a cutoff radius \(r_\text{c}{\in}\{5\,\text{\AA},7\,\text{\AA},9\) Å}, the \(GW\) gap is identical within 27 meV.

8 GW band structure of \(\text{MoS}_\text{2}\), \(\text{MoSe}_\text{2}\), WS\(_\text{2}\), WS\(\text{e}_\text{2}\)↩︎

We employ our \(GW\) algorithm to compute the band structure of \(\text{MoS}_\text{2}\), \(\text{MoSe}_\text{2}\), WS\(_\text{2}\) and WS\(\text{e}_\text{2}\) using two sets of parameters: an aug-SZV-MOLOPT basis with light settings and an aug-DZVP-MOLOPT basis with tight settings, as shown in Fig. 4. For comparison, we also include band structures computed using a reference plane-wave-based implementation (BerkeleyGW) with 8 empty bands. The overall qualitative agreement is good; for quantitative assessment, we focus on the band gap at the K-point (see Table [t1]). With the aug-DZVP-MOLOPT basis and tight settings, the \(GW\) band gaps agree with the plane-wave results on average within 50 meV. The agreement varies across the four materials, and is better for materials containing sulfur (20 meV) than materials containing selenium (75 meV). Note that this discrepancy is also observed between the plane-wave codes (Table [t1]), as the average deviation is 10 meV for the materials containing sulfur and 30 meV for the materials containing selenium. For the aug-SZV-MOLOPT basis and light settings, the average agreement is 70 meV.

This shows that one can achieve decent numerical precision of \(GW\) band structures using the relatively small aug-SZV-MOLOPT basis set. The observed material dependence of the \(GW\) band gaps can be rationalized by several factors. First, some differences are already present at the DFT level with gaps that differ by 5 meV for MoS₂ and up to 30 meV for WSe₂ (Table [t1]). The fact that, already at the DFT level with SOC, the differences between the codes are relatively minor in MoS₂ explains, at least in a naïve way, the better agreement among codes after applying the \(GW\) correction for this material compared to WSe₂. Second, the different implementations of SOC in each code (and their corresponding basis sets) can lead to slight differences in the computed dielectric screening, which then propagate to the final \(GW\) band gaps. Finally, as noted above, deviations are stronger in selenium-based TMDCs compared to sulfur-based ones. Selenium is more polarizable than sulfur and has more spatially extended valence orbitals, making selenium-based TMDCs more sensitive to how high-energy and local-field effects are represented.

The average deviations of \(GW\) band gaps between the numerical approaches are summarized in Fig. 5.

The discrepancies may be attributed to several factors, including the use of different pseudopotentials, the limited size of the aug-DZVP-MOLOPT basis compared to high plane-wave cutoffs, and sensitivity to convergence parameters in frequency integration or dielectric matrix evaluation. Time-frequency resolution may also contribute to residual differences, which could be mitigated using Hedin’s shift [112]. This is also demonstrated in Appendix 19, where we report some available results in the literature for the band gaps of the TMDCs at the \(G_0W_0\)@PBE+SOC. As we have stated earlier, one should note however that the comparison with the already existing calculations is fairly difficult as most of these use different lattice constants, which has a huge impact on the band gap [80].

Comparison to experimental measurements is inherently challenging, as the band gap is highly sensitive to external influences that are difficult to control, such as substrate screening effects [113] and strain [80]. Reported experimental band gaps for the four materials range from 1.9 eV to 2.4 eV [113]–[116], which aligns with our \(GW\) results (Table [t1]). However, achieving a precise one-to-one correspondence for each material remains elusive at this stage, especially given that many experimentally reported band gaps are measured in the presence of a substrate, which is not considered in this work, as it lies outside the scope of this manuscript.

For completeness, we also performed benchmark full-spinor calculations with BerkeleyGW using the generalized plasmon pole model [97], see Table [t2] in Appendix 18. This approximation yields bandgaps with significant deviations - up to \(0.25\) eV - when compared with the full frequency calculation of the quasi-particle self-energy. These discrepancies sharply contrast with the excellent agreement of below \(50\) meV between the different codes, as reported in Fig. 5.

To demonstrate the computational efficiency of our \(GW\) algorithm, we carried out \(GW\) band structure calculations with settings "aug-SZV-MOLOPT, light" on a laptop with an Intel i5 processor (14 cores) and 192 GB RAM. The computation time on this hardware is between 20 and 45 hours, depending on the 2D material. To study the scalability of our implementation, we ran a computational performance test on an HPC system, using the same parameters as for the laptop calculations in order for the timings to be comparable. The results show that on this HPC system using 1024 cores, the same calculations completed in approximately 23 minutes (Fig. 6). These results demonstrate that our method enables efficient and scalable \(GW\) band structure calculations of 2D materials.

9 Conclusion and outlook↩︎

We have developed a \(GW\) space-time algorithm for periodic systems using Gaussian basis functions with spin-orbit coupling, enabling accurate quasiparticle band structure calculations for atomically thin materials. Our implementation, available in the open-source code CP2K, achieves high accuracy: for monolayer MoS\(_2\), MoSe\(_2\), WS\(_2\), and WSe\(_2\), \(GW\) band gaps agree on average within 50 meV to plane-wave \(GW\) band gaps from BerkeleyGW and VASP. Using lighter computational settings, the average agreement is 70 meV, and a complete \(GW\) band structure calculation can be performed on a laptop (Intel i5, 192GB RAM) in approximately 24 hours, or in less than 30 minutes on 1024 HPC cores. Future extensions of this work will combine the present \(GW\) algorithm with real-space grids for \(G\) and \(W\) [52], [117], [118]. We expect this to significantly reduce the number of required lattice summations, further lowering the computational cost. The Gaussian basis sets investigated in this work yield converged band gaps and band structures within \(\sim\) 100 meV. Developing Gaussian basis sets that enable convergence of periodic \(GW\) calculations ideally within 10 meV is subject of ongoing work.

Data and Code availability↩︎

Inputs and outputs of all calculations reported in this work are available in a Github repository [119]. The algorithm presented in this work is available in the open-source package CP2K [82]–[84].

Appendix↩︎

10 Periodic RI↩︎

  For computing \(\chi_{PQ}^\mathbf{R}(i\tau)\) from Eq. 31 , the three-centre integrals \((\mu\mathbf{R}\,\nu\mathbf{T}|P\mathbf{R})\) are appearing. Moreover, in Eq. 30 the metric matrix \(M_{PQ}(\mathbf{k})\) is appearing. \((\mu\mathbf{R}\,\nu\mathbf{T}|P\mathbf{R})\) and \(M_{PQ}(\mathbf{k})\) are stemming from the application of the resolution of the identity (RI) [44], [47], [66]. In this Appendix, we give details on RI with periodic boundary conditions.

In quantum chemistry, using Gaussian basis functions without periodic boundary conditions, RI is expressed as \[\begin{align} \phi_\mu(\mathbf{r})\,\phi_\nu(\mathbf{r}) \simeq \sum_P B_P^{\mu\nu}\varphi_P(\mathbf{r})\,,\label{ec1} \end{align}\tag{47}\] where the product of two Gaussian functions \(\phi_\mu(\mathbf{r})\) and \(\phi_\nu(\mathbf{r})\) is approximated as a linear combination of auxiliary Gaussian functions \(\varphi_P(\mathbf{r})\) with coefficients \(B_P^{\mu\nu}\). \(B_P^{\mu\nu}\) are determined by fitting procedures like least-squares minimization [66] to get a good approximation in Eq. 47 .

In the periodic case, the basis functions \(\phi_\mu^{\mathbf{R}}(\mathbf{r})\) and \(\phi^{\mathbf{T}}_\nu(\mathbf{r})\) can be located in cell \(\mathbf{R}\) and \(\mathbf{T}\). The RI expansion then uses auxiliary Gaussians \(\varphi_P^\mathbf{P}(\mathbf{r})\) in any cell \(\mathbf{P}\), \[\begin{align} \phi_\mu^{\mathbf{R}}(\mathbf{r})\,\phi_\nu^{\mathbf{T}}(\mathbf{r}) \simeq \sum_{P\mathbf{P}} B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}}\varphi_P^\mathbf{P}(\mathbf{r})\,.\label{ec2} \end{align}\tag{48}\] We show in this Appendix that the expansion coefficients \(B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}}\) are determined by least-squares minimization under a metric \(m(\mathbf{r})\) [66] as \[\begin{align} B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}} =\sum_{Q\mathbf{Q}} (\mu\mathbf{R}\,\nu\mathbf{T}\,|\,Q\mathbf{Q})_m\,(M^{-1})_{PQ}^{\mathbf{P}-\mathbf{Q}}\,,\label{ec3} \end{align}\tag{49}\] where \((\mu\mathbf{R}\,\nu\mathbf{T}\,|\,Q\mathbf{Q})_m\) is the three-centre integral of the metric \(m(\mathbf{r})\) (\(\mathbf{Q}\) is a lattice vector), \[\begin{align} (\mu\mathbf{R}\,\nu\mathbf{T}&\,|\,Q\mathbf{Q})_m = {\int} d\mathbf{r}\,d\mathbf{r}'\; \phi_\mu^{\mathbf{R}}(\mathbf{r})\, \phi_\nu^{\mathbf{T}}(\mathbf{r})\, m(\mathbf{r}-\mathbf{r}') \, \varphi_Q^{\mathbf{Q}}(\mathbf{r}') \,, \end{align}\] and \[\begin{align} (M^{-1})_{PQ}^{\mathbf{R}} = \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\; e^{-i\mathbf{k}\cdot\mathbf{R}}\,M^{-1}_{PQ}(\mathbf{k}) \label{ec4} \end{align}\tag{50}\] where \(M^{-1}_{PQ}(\mathbf{k})\) are matrix elements of the inverse of the metric matrix \[\begin{align} M_{PQ}(\mathbf{k}) =& \sum_\mathbf{R}e^{i\mathbf{k}\cdot\mathbf{R}}{\int} d\mathbf{r}\,d\mathbf{r}'\, \varphi_P^\mathbf{0}(\mathbf{r}) \, m(\mathbf{r}-\mathbf{r}') \,\varphi_Q^\mathbf{R}(\mathbf{r}')\,. \label{ec5} \end{align}\tag{51}\] Note that we regularize \(\mathbf{M}^{-1}(\mathbf{k})\), Eq. 34 , to prevent for linear dependencies and thus large expansion coefficients \(B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}}\) with alternating sign in the RI expansion 48 .

In the following, we outline our proof of Eq. 49 . As in non-periodic RI [66], we define the residual \(\mathfrak{R}\) of Eq. 48 , \[\begin{align} \mathfrak{R}(\mathbf{r}) = \phi_\mu^{\mathbf{R}}(\mathbf{r})\,\phi_\nu^{\mathbf{T}}(\mathbf{r}) {-}\sum_{P\mathbf{P}} B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}}\varphi_P^\mathbf{P}(\mathbf{r})\,.\label{ec7} \end{align}\tag{52}\] Now, we vary the expansion coefficients \(B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}}\) in Eq. 52 to minimize the repulsion of \(\mathfrak{R}\) with itself in the metric \(m\), \[\begin{align} (\mathfrak{R}|\mathfrak{R})_m ={\int} d\mathbf{r}\,d\mathbf{r}'\; \mathfrak{R}(\mathbf{r})\; m(\mathbf{r}-\mathbf{r}') \;\mathfrak{R}(\mathbf{r}')\geq 0\;\rightarrow\;\text{min}\,. \end{align}\] In the ideal case, we have \(\mathfrak{R}{=}0\) yielding zero repulsion of \(\mathfrak{R}\) with itself. In the general case, we are looking for a minimum of \((\mathfrak{R}|\mathfrak{R})_m\), i.e. we take \(\partial (\mathfrak{R}|\mathfrak{R})_m /\partial B_{P\mathbf{P}}^{\mu\mathbf{R}\nu\mathbf{T}}{=}0\) which gives \[\begin{align} -2 (\mu\mathbf{R}\,\nu\mathbf{T}\,|\,P\mathbf{P})_m +2\sum_{Q\mathbf{Q}} B_{Q\mathbf{Q}}^{\mu\mathbf{R}\nu\mathbf{T}} (Q\mathbf{Q}|P\mathbf{P})_m = 0 \end{align}\] and thus \[\begin{align} \sum_{Q\mathbf{Q}} B_{Q\mathbf{Q}}^{\mu\mathbf{R}\nu\mathbf{T}} M_{PQ}^{\mathbf{P}-\mathbf{Q}} = (\mu\mathbf{R}\,\nu\mathbf{T}\,|\,P\mathbf{P})_m\,.\label{ec10} \end{align}\tag{53}\] We insert \[\begin{align} M_{PQ}^{\mathbf{P}-\mathbf{Q}} = \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\; e^{-i\mathbf{k}\cdot(\mathbf{P}-\mathbf{Q})} \,M_{PQ}(\mathbf{k}) \end{align}\] into Eq. 53 , we multiply with \(e^{i\mathbf{q}\cdot \mathbf{P}}\) and we sum over all lattice vectors \(\mathbf{P}\): \[\begin{align} \sum_{Q\mathbf{Q}} B_{Q\mathbf{Q}}^{\mu\mathbf{R}\nu\mathbf{T}} \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\; e^{ i\mathbf{k}\cdot\mathbf{Q}}& \,M_{PQ}(\mathbf{k}) \sum_\mathbf{P}e^{i(\mathbf{q}-\mathbf{k})\cdot\mathbf{P}}\nonumber \\ & = \sum_{\mathbf{P}} e^{ i\mathbf{q}\cdot\mathbf{P}} (\mu\mathbf{R}\,\nu\mathbf{T}\,|\,P\mathbf{P})_m\,. \end{align}\] We use \(\sum\limits_\mathbf{P}e^{i(\mathbf{q}-\mathbf{k})\cdot\mathbf{P}}{=}\Omega_\text{BZ}\, \delta(\mathbf{q}{-}\mathbf{k})\), \(\int\limits_\text{BZ}d\mathbf{k}\, f(\mathbf{k})\delta(\mathbf{q}{-}\mathbf{k}){=}f(\mathbf{q})\) and we multiply with the matrix \(\mathbf{M}^{-1}(\mathbf{q})\) to obtain \[\begin{align} \sum_{\mathbf{Q}} B_{Q\mathbf{Q}}^{\mu\mathbf{R}\nu\mathbf{T}} e^{i\mathbf{q}\cdot\mathbf{Q}} = \sum_{P\mathbf{P}} e^{ i\mathbf{q}\cdot\mathbf{P}} (\mu\mathbf{R}\,\nu\mathbf{T}\,|\,P\mathbf{P})_m\,M^{-1}_{PQ}(\mathbf{q})\,. \end{align}\] We then multiply with \(e^{-i\mathbf{q}\cdot\mathbf{Q}'}\), we integrate over the BZ (\(\mathbf{q}\)) and we use Eq. 50 as well as \[\begin{align} \int\limits_\text{BZ} \frac{d\mathbf{q}}{\Omega_\text{BZ}}\,e^{i\mathbf{q}(\mathbf{Q}-\mathbf{Q}')} = \delta_{\mathbf{Q}\mathbf{Q}'} \end{align}\] (\(\mathbf{Q}\) and \(\mathbf{Q}'\) are lattice vectors; \(\delta_{\mathbf{Q}\mathbf{Q}'}\) is the Kronecker-\(\delta\)) to obtain Eq. 49 .

11 Derivation of Eq. 31 for computing \(\chi^\mathbf{R}_{PQ}(i\tau)\)↩︎

For deriving the irreducible response \(\chi^\mathbf{R}_{PQ}(i\tau)\) for a periodic system, Eq. 31 , we start from Eq. 3 and 2 for computing \(\chi\) for a molecule. For a molecule, we sum over the molecular quantum number \(i\) and \(a\) for occupied and empty MOs in Eq. 2 . For a periodic system, the quantum numbers of the one-electronic states are \(i\mathbf{k}\) and \(a\mathbf{k}\) labeling an occupied and empty band at \(k\)-point \(\mathbf{k}\) in the BZ. We thus replace and abbreviate \[\begin{align} \sum_i^\text{occ}\;&\rightarrow \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\;\;\sum_i^\text{occ} \simeq\sum_{i\mathbf{k}}^{\text{occ}} \\[0.5em] \sum_a^\text{empty}\;&\rightarrow \int\limits_\text{BZ} \frac{d\mathbf{k}}{\Omega_\text{BZ}}\;\sum_a^\text{empty} \;\;\simeq \; \sum_{a\mathbf{k}}^{\text{empty}} \end{align}\] so that we get from Eq. 3 /2 for \(\tau{>}0\) omitting the imaginary prefactor: \[\begin{align} \chi (\mathbf{r}, \mathbf{r}', i \tau) =&\; - \sum_{i\mathbf{k}}^{\text{occ}}\psi_{i\mathbf{k}} ^*(\mathbf{r}) \psi_{i\mathbf{k}} (\mathbf{r}') e^{-|(\varepsilon_{i\mathbf{k}}^\text{DFT}-\varepsilon_\text{F}) \tau|}\nonumber \\ &\times \sum_{a\mathbf{k}}^{\text{empty}}\psi_{a\mathbf{k}}(\mathbf{r})\psi_{a\mathbf{k}}^*(\mathbf{r}') e^{ -|(\varepsilon_{a\mathbf{k}}^\text{DFT}-\varepsilon_\text{F}) \tau|} \\[0.5em] \overset{\eqref{kpointsbf},\eqref{e29}}{=} \sum_{\lambda\mathbf{R}_1\mu\mathbf{S}_1}& \left[ \sum_{\mathbf{k}} e^{i\mathbf{k}\cdot(\mathbf{R}_1-\mathbf{S}_1)} G_{\lambda\mu}(-i\tau,\mathbf{k}) \right] \phi_\mu^{\mathbf{S}_1}(\mathbf{r})\, \phi_\lambda^{\mathbf{R}_1}(\mathbf{r}')\nonumber \\ \times \sum_{\sigma\mathbf{S}_2\nu\mathbf{R}_2}& \left[ \sum_{\mathbf{q}} e^{i\mathbf{q}\cdot(\mathbf{R}_2-\mathbf{S}_2)} G_{\nu\sigma}(i\tau,\mathbf{q}) \right] \phi_\nu^{\mathbf{R}_2}(\mathbf{r})\, \phi_\sigma^{\mathbf{S}_2}(\mathbf{r}') \\[1em] \overset{\eqref{e23}}{=}\;\sum_{\lambda\mathbf{R}_1\mu\mathbf{S}_1}& G_{\lambda\mu}^{\mathbf{R}_1-\mathbf{S}_1}(-i\tau ) \; \phi_\mu^{\mathbf{S}_1}(\mathbf{r})\, \phi_\lambda^{\mathbf{R}_1}(\mathbf{r}')\nonumber \\ \times \sum_{\sigma\mathbf{S}_2\nu\mathbf{R}_2}& G_{\nu\sigma}^{\mathbf{R}_2-\mathbf{S}_2}(i\tau ) \; \phi_\nu^{\mathbf{R}_2}(\mathbf{r})\, \phi_\sigma^{\mathbf{S}_2}(\mathbf{r}')\,.\label{a5} \end{align}\tag{54}\] We obtain the matrix element \(\chi_{PQ}^\mathbf{R}(i\tau)\) from the projection of \(\chi(\mathbf{r},\mathbf{r}',i\tau)\) in the RI basis \(\{\varphi_P^\mathbf{R}(\mathbf{r})\}\), incorporating the RI metric which is the truncated Coulomb operator 29 : \[\begin{align} \chi_{PQ}^\mathbf{R}(i\tau) =& \braket{\varphi^\mathbf{0}_P|\chi(i\tau)|\varphi^\mathbf{R}_Q} \\[0.5em] =& \int d\mathbf{r}_1 \,d\mathbf{r}_2\,d\mathbf{r}_3\,d\mathbf{r}_4\; \varphi_P^\mathbf{0}(\mathbf{r}_1)\; V_{r_\text{c}}(\mathbf{r}_1,\mathbf{r}_2)\;\nonumber \\&\times \chi(\mathbf{r}_2,\mathbf{r}_3,i\tau)\; V_{r_\text{c}}(\mathbf{r}_3,\mathbf{r}_4)\; \varphi_Q^\mathbf{R}(\mathbf{r}_4) \\[1em] =& \sum_{\lambda\mathbf{R}_1}\sum_{\mu\mathbf{S}_1} \sum_{\sigma\mathbf{S}_2}\sum_{\nu\mathbf{R}_2} G_{ \lambda\mu}^{\mathbf{R}_1-\mathbf{S}_1}(-i\tau ) \; G_{\nu\sigma}^{\mathbf{R}_2-\mathbf{S}_2}(i\tau )\nonumber \\[0.5em] &\times (\mu\mathbf{S}_1\,\nu\mathbf{R}_2\, |\,P\,\mathbf{0})\; (\lambda\mathbf{R}_1\,\sigma \mathbf{S}_2 \,|\,Q\mathbf{R})\,. \end{align}\] We replace \(\mathbf{S}_1{\rightarrow}\,\mathbf{R}_1{-}\,\mathbf{S}_1\) and \(\mathbf{S}_2{\rightarrow}\,\mathbf{R}_2{-}\,\mathbf{S}_2\) such that Eq. 31 follows: \[\begin{align} \chi_{PQ}^\mathbf{R}(i\tau)& = \sum_{\lambda\mathbf{R}_1}\sum_{\mu\mathbf{S}_1} \sum_{\sigma\mathbf{S}_2}\sum_{\nu\mathbf{R}_2} G_{ \lambda\mu}^{ \mathbf{S}_1}(-i\tau ) \; G_{\nu\sigma}^{ \mathbf{S}_2}(i\tau )\nonumber \\[0.5em] &\times (\mu{\mathbf{R}_1}{-}\mathbf{S}_1\,\nu\mathbf{R}_2\, |\,P\,\mathbf{0})\; (\lambda\mathbf{R}_1\,\sigma \mathbf{R}_2{-}\mathbf{S}_2 \,|\,Q\mathbf{R})\,. \end{align}\]

12 Parallel implementation of Eq. 31 for computing \(\chi^\mathbf{R}_{PQ}(i\tau)\)↩︎

An efficient, low-memory parallel implementation of Eq. 31 is a key for routinely executing the present algorithm. We sketch our parallel implementation in the algorithm from Fig. 7: The first main idea of the algorithm is to divide the \(N\) MPI processes in groups with \(n\) MPI processes per group, i.e., the number of groups is \(N/n\). We store all three-centre matrix elements \((\mu\mathbf{R}\,\nu\mathbf{S}\,|\,P\mathbf{0})\) in every subgroup, i.e., the memory available to the subgroup needs to be sufficiently large to fit all \((\mu\mathbf{R}\,\nu\mathbf{S}\,|\,P\mathbf{0})\) elements. This is automatically guaranteed as the program sets the group size \(n\) such that \[\begin{align} n \times \text{avail.~mem.~per MPI proc.} > \text{mem.~of }(\mu\mathbf{R}\,\nu\mathbf{S}\,|\,P\mathbf{0})\,. \end{align}\] We execute all tensor operations from Eq. 31 group-locally to avoid communication between all MPI processes. The second main idea is to distribute the computations of Eq. 31 on the different groups by distributing \(\mathbf{R}_1{-}\mathbf{R}_2\,{=}{:}\,\Delta \mathbf{R}\) to different subgroups. It turns out that with this distribution, communication between different groups is avoided; only the group-local result for \(\chi^\mathbf{R}_{PQ}(i\tau)\) needs to be summed up, see last line in the algorithm from Fig. 7. We employ the set "3c" of lattice vectors \(\mathbf{R}\) which are the lattice vectors \(\mathbf{R}\) with large three-centre integral \((\mu\mathbf{R}\,\nu\mathbf{S}\,|\,P\mathbf{0})\), i.e., \[\begin{align} \label{Frobcut} \mathbf{R}\in\text{3c}\Leftrightarrow\text{there is \mathbf{S} such that } \text{F}[(\mu\mathbf{R}\,\nu\mathbf{S}\,|\,P\mathbf{0})] > \delta\,, \end{align}\tag{55}\] where F denotes the Frobenius norm, \[\begin{align} \label{Frobnorm} \text{F}[(\mu\mathbf{R}\,\nu\mathbf{S}\,|\,P\mathbf{0})] := \sqrt{\sum_{\mu\nu P} |(\mu\mathbf{R}\,\nu\mathbf{S}\,|\,P\mathbf{0})|^2}\;, \end{align}\tag{56}\] and \(\delta\) is the filter threshold. All other quantities of the algorithm from Fig. 7 are defined in the main text. The three tensor operations are executed by the sparse-tensor library dbt [120].

13 Lattice sum of the Coulomb matrix element 35↩︎

We compute the lattice sum 35 , \[\begin{align} V_{PQ}(\mathbf{k}) &= \sum_\mathbf{R}e^{i\mathbf{k}\cdot\mathbf{R}}\, (\varphi_P^\mathbf{0}| \varphi_Q^\mathbf{R})\;, \label{ed1} \\ (\varphi_P^\mathbf{0}| \varphi_Q^\mathbf{R}) &= {\int} d\mathbf{r}\,d\mathbf{r}'\,\varphi_P^\mathbf{0}(\mathbf{r})\,\frac{ 1}{|\mathbf{r}- \mathbf{r}'|}\,\varphi_Q^\mathbf{R}(\mathbf{r}')\,, \end{align}\tag{57}\] by explicit summation of \(\sum_\mathbf{R}\). Several schemes for executing this lattice sum have been described, see Ref. [121], Appendix of Ref. [73] and references therein.

We first note that the lattice summation in Eq. 57 is divergent if \(\varphi_P^\mathbf{0}\) and \(\varphi_Q^\mathbf{R}\) are \(s\)-type basis functions in the limit \(\mathbf{k}{\rightarrow}0\). This is because \((\varphi_P^\mathbf{0}| \varphi_Q^\mathbf{R}){\sim}|\mathbf{R}|^{-1}\) and the lattice sum \(\sum\limits_\mathbf{R}(\varphi_P^\mathbf{0}| \varphi_Q^\mathbf{R}){\sim}\sum\limits_\mathbf{R}|\mathbf{R}|^{-1}\) (for \(\mathbf{k}{=}0\)) can be approximated as an integral \(\int_{\mathbb{R}^d} d\mathbf{R}\,|\mathbf{R}|^{-1}{\sim}\int_0^\infty dR \,R^{d-1} \,R^{-1}\) which diverges for dimensionality \(d{=}1,2,3\) (1d-molecular chains, 2d-materials or surfaces, 3d-bulk solids). More precisely, the Coulomb matrix element \(V_{PQ}(\mathbf{k})\) diverges at \(\mathbf{k}{\rightarrow}0\) as \(V_{PQ}(\mathbf{k}){\sim}|\mathbf{k}|^{-(d-1)}\) for \(s\)-functions \(P,Q\) for \(d{=}2,3\) (for \(d{=}1\) the divergence is logarithmic). [76] The whole \(GW\) algorithm is not divergent because this \(1/k^{d-1}\) divergence is integrable in the Brillouin zone, see as an example Eq. 37 : \(\int_\text{BZ}d\mathbf{k}\,1/k^{d-1}{\sim}\int_0^{k_\text{max}} dk\,k^{d-1} / k^{d-1}\), where \(k_\text{max}\) is determined by the Brillouin zone size.

The requirements for convergence of the lattice sum 57 have been extensively discussed in the literature [73], [121] and references therein. We reproduce the arguments here to make the discussion self-contained. We consider the lattice sum \[\begin{align} V_{P\rho}(\mathbf{k}) = \sum_\mathbf{T}e^{i\mathbf{k}\cdot\mathbf{T}}\, (\varphi_P^\mathbf{0}| \rho^\mathbf{T}) \label{ed3} \end{align}\tag{58}\] where \(\mathbf{T}\) are lattice vectors and \(\rho^\mathbf{T}(\mathbf{r}){=}\rho(\mathbf{r}{-}\mathbf{T})\) where \(\rho(\mathbf{r})\) is a function of \(\mathbf{r}\) which decays exponentially or faster for large \(|\mathbf{r}|\). The function \(\varphi_P^\mathbf{0}(\mathbf{r})\) is assigned to cell \(\mathbf{0}\) and decays exponentially or faster for large \(|\mathbf{r}|\).

The lattice sum 58 is absolutely convergent for dimensionality \(d{=}1,2,3\) if \[\begin{align} k_\text{min} + l_\text{min} \ge d \end{align}\] where \(2^{k_\text{min}}\) and \(2^{l_\text{min}}\) are the lowest nonvanishing dipole moments of \(\varphi^\mathbf{0}_P\) and \(\rho^\mathbf{T}\), respectively [121]. In our case, \(\varphi^\mathbf{0}_P\) can be an \(s\)-type basis function, i.e., \(k_\text{min}{=}0\), so \[\begin{align} l_\text{min} \ge d\label{ed6a} \end{align}\tag{59}\] guarantees absolute convergence of the lattice sum 58 for all Gaussian basis functions \(\varphi^\mathbf{0}_P\). Absolute convergence implies that the order of \(\mathbf{T}\) when executing the lattice sum Eq. 58 is irrelevant for the result.

Eq. 59 means it is required for absolute convergence of lattice sum 58 for \(d{=}1,2,3\) that the function \(\rho(\mathbf{r})\) (and thus \(\rho^\mathbf{T}(\mathbf{r})\)) has zero monopole moment, i.e., \[\begin{align} \int\limits_{\mathbb{R}^3} d\mathbf{r}\;\rho^\mathbf{T}(\mathbf{r}) = \int\limits_{\mathbb{R}^3} d\mathbf{r}\;\rho(\mathbf{r}) = 0\,.\label{ed6b} \end{align}\tag{60}\] Additionally, for \(d{=}2,3\), all dipole moments of \(\rho^\mathbf{T}\) need to vanish, \[\begin{align} \int\limits_{\mathbb{R}^3} d\mathbf{r}\;r_\alpha\,\rho^\mathbf{T}(\mathbf{r}) = \int\limits_{\mathbb{R}^3} d\mathbf{r}\;r_\alpha\,\rho(\mathbf{r}) = 0\,.\label{ed6c} \end{align}\tag{61}\] Additionally, for \(d{=}3\), all quadrupole moments of \(\rho^\mathbf{T}\) need to vanish, \[\begin{align} \int\limits_{\mathbb{R}^3} d\mathbf{r}\;r_\alpha\,r_\beta\,\rho^\mathbf{T}(\mathbf{r}) = \int\limits_{\mathbb{R}^3} d\mathbf{r}\;r_\alpha\,r_\beta\,\rho(\mathbf{r}) = 0\,,\label{ed6d} \end{align}\tag{62}\] i.e., Eq. 62 needs to hold for all combinations of \(\alpha,\beta{=}1,2,3\).

We apply this theorem to derive an absolutely convergent expression for the lattice sum \(V_{PQ}(\mathbf{k})\) [Eq. 35 , reproduced in Eq. 57 ], here illustrated for the case \(d{=}2\). An absolutely convergent lattice summation is essential, as it guarantees that the summation result is independent of the summation order. This is particularly important in numerical computations, where the lattice sum must be truncated, thereby imposing a specific summation order.

We start the derivation of an absolutely convergent lattice summation by defining a \((2N_1){\times}(2N_2)\) supercell "SC" as sketched in Fig. 8. Here, \(N_1\) and \(N_2\) correspond to the \(N_1{\times}N_2\) Monkhorst-Pack \(k\)-point mesh [65] used for the BZ integration 37 of \(W\). A lattice vector \(\mathbf{R}\) belongs to SC if \[\begin{align} \mathbf{R}\in \text{SC}\Leftrightarrow\mathbf{R}= \sum_{j=1}^d n_j\,\mathbf{a}_j \;, n_j\in \{0, 1,\;\ldots\,,\, 2N_j-1\}\,.\label{ed5} \end{align}\tag{63}\] It is then easy to show that the sum of the phase factors \(e^{i\mathbf{k}\cdot\mathbf{R}}\) over the SC lattice vectors vanishes, \[\begin{align} \sum_{\mathbf{R}\in\text{SC}}e^{i\mathbf{k}\cdot\mathbf{R}} = 0 \label{ed6} \end{align}\tag{64}\] if \(\mathbf{k}\) is contained in the \({N_1}{\times}N_2\) Monkhorst-Pack \(k\)-point mesh in case \(N_1\) and \(N_2\) are even integers (see Appendix 14; in case \(N_1,N_2\) are both odd, the \(k\)-mesh contains the \(\Gamma\)-point \(\mathbf{k}{=}0\) which violates Eq. 64 ).

Eq. 64 motivates us to reorder the lattice sum 57 using the superlattice (SL) with lattice vectors \(\mathbf{T}\), \[\begin{align} \mathbf{T}\in \text{SL}\Leftrightarrow\mathbf{T}= 2 N_1\, t_1\, \mathbf{a}_1 + 2 N_2\, t_2\, \mathbf{a}_2 \;, \label{ed7} \end{align}\tag{65}\] where each \(t_1, t_2\) is integer. As illustrated in Fig. 8, we carry out the infinite lattice sum 57 by an infinite sum over the superlattice SL: \[\begin{align} V_{PQ}(\mathbf{k}) &= \sum_{\mathbf{T}\in\text{SL}} \, \sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot(\mathbf{R}+\mathbf{T})}\, (\varphi_P^\mathbf{0}| \varphi_Q^{\mathbf{R}+\mathbf{T}}) \tag{66} \\[0.5em] &\equiv \sum_{\mathbf{T}\in\text{SL}} e^{i\mathbf{k}\cdot \mathbf{T}} \,(\varphi_P^\mathbf{0}|\rho^\mathbf{T}) \tag{67} \end{align}\] where \[\begin{align} \rho^\mathbf{T}(\mathbf{r})= \sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \mathbf{R}}\, \varphi_Q^{\mathbf{R}+\mathbf{T}}(\mathbf{r})\,. \end{align}\] The benefit of lattice sum 67 is that \(\rho^\mathbf{T}\) has a zero monopole moment, \[\begin{align} \int\limits_{\mathbb{R}^3} d\mathbf{r}\;\rho^\mathbf{T}(\mathbf{r}) & = \sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \mathbf{R}} \int\limits_{\mathbb{R}^3} d\mathbf{r}\; \varphi_Q^{\mathbf{R}+\mathbf{T}}(\mathbf{r})\nonumber \\[0.5em] &= \sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \mathbf{R}} \int\limits_{\mathbb{R}^3} d\mathbf{r}\; \varphi_Q^{\mathbf{0}}(\mathbf{r}) \overset{\eqref{ed6}}{=} 0\,, \label{ed11} \end{align}\tag{68}\] and zero dipole moments, \[\begin{align} & \int\limits_{\mathbb{R}^3} d\mathbf{r}\;r_\alpha\,\rho^\mathbf{T}(\mathbf{r}) = \sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \mathbf{R}} \int\limits_{\mathbb{R}^3} d\mathbf{r}\;r_\alpha\, \varphi_Q^{\mathbf{R}+\mathbf{T}}(\mathbf{r})\nonumber \\[0.5em] &= \sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \mathbf{R}} \int\limits_{\mathbb{R}^3} d\mathbf{r}\;(r_\alpha-R_\alpha-T_\alpha)\, \varphi_Q^{\mathbf{0}}(\mathbf{r}-\mathbf{R}-\mathbf{T})\nonumber \\ & + \sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \mathbf{R}} \,R_\alpha \int\limits_{\mathbb{R}^3} d\mathbf{r}\; \varphi_Q^{\mathbf{R}+\mathbf{T}}(\mathbf{r}) + T_\alpha \sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \mathbf{R}} \int\limits_{\mathbb{R}^3} d\mathbf{r}\; \varphi_Q^{\mathbf{R}+\mathbf{T}}(\mathbf{r}) \nonumber \\[0.5em] & \overset{\eqref{ed6}}{=} 0 \,+\,\sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \mathbf{R}} \,R_\alpha \int\limits_{\mathbb{R}^3} d\mathbf{r}\; \varphi_Q^{\mathbf{R}+\mathbf{T}}(\mathbf{r}) \,+\, 0 = 0\,, \label{ed11a} \end{align}\tag{69}\] where the last equality stems from: \[\begin{align} &\sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \mathbf{R}} \,R_\alpha=\sum_{(R_\alpha,R_\beta)\in\text{SC}} e^{i(k_\alpha R_\alpha+k_\beta R_\beta)}\,R_\alpha \nonumber\\[0.5em] &=\sum_{R_\alpha} e^{ik_\alpha R_\alpha}\,R_\alpha\;\cdot\;\sum_{R_\beta} e^{ik_\beta R_\beta} = 0\,, \end{align}\] where we used \[\begin{align} \sum_{R_\beta} e^{ik_\beta R_\beta} = 0\label{ed17} \end{align}\tag{70}\] as shown in Appendix 14, Eq. 73 . In general, one can easily extend this result and see that for a \(d\)-dimensional supercell, all moments of order \(k{<}d\) vanish if Eq. 64 is verified. As such, it is a sufficient condition for the suppression of the quadrupole \((k{=}2)\) moment for \(d{=}3\), which is condition Eq. 62 , and therefore also guarantees the absolute convergence of the three-dimensional lattice sum  58 . Comparing Eqs. 68 and 69 to Eqs. 61 and 62 implies that the lattice sum over \(\mathbf{T}\) in Eq. 66 is absolutely convergent, i.e., the result of the lattice sum 66 /67 over \(\mathbf{T}\) is independent of the actual order of summation and we thus can truncate the summation after a finite number of superlattice vectors \(\mathbf{T}\).

In practice, we check the convergence of the \(\mathbf{T}\) sum by restricting \(\mathbf{T}\) from Eq. 65 in lattice sum 66 to \[\begin{align} \mathbf{T}_{t_1,\,t_2}= 2 N_1\, t_1\, \mathbf{a}_1 + 2 N_2\, t_2\, \mathbf{a}_2\;,t_j\in \{0,\pm1,\ldots,\pm M\} \end{align}\] leading to the lattice sum \[\begin{align} V_{PQ}(\mathbf{k}) &= \sum_{t_1,\,t_2 = -M}^{M} \; \sum_{\mathbf{R}\in\text{SC}} e^{i\mathbf{k}\cdot \left(\mathbf{R}+ \mathbf{T}_{t_1,\,t_2}\right)}\; (\varphi_P^\mathbf{0}| \varphi_Q^{\mathbf{R}+ \mathbf{T}_{t_1,\,t_2}}) \label{ed13} \end{align}\tag{71}\] for \(\mathbf{k}\) from the even \(N_1{\times}N_2\) Monkhorst-Pack \(k\)-point mesh. Eq. 71 is implemented in the code and can be converged when increasing \(M\). We define the size factor \(\Delta\) for Fig. 2 as: \[\label{lattsacle} M=2^{\Delta-1}\tag{72}\] We show in Appendix 15 that the lattice sum 71 reproduces the Coulomb integrals used in plane-wave \(GW\) algorithms.

14 Proof of Eq. 64 and 70↩︎

The \(N_1{\times}N_2\) Monkhorst-Pack \(k\)-points (\(N_1, N_2\) even) are \[\begin{align} \mathbf{k}_\ell = \sum_{j=1}^{2} \frac{ \ell_j}{2N_j}\,\mathbf{b}_j\;, \end{align}\] where we define \(\ell{=}(\ell_1,\ell_2)\) and \(\ell_j\) takes as value one of the following odd integers \[\begin{align} \ell_j\in \{\pm\,1,\, \pm\,3, \;\ldots\;,\, \pm\, (N_j-1)\} \,. \end{align}\] \(\mathbf{b}_j\) are primitive translation vectors of the reciprocal lattice with \(\mathbf{a}_{j_1}{\cdot}\mathbf{b}_{j_2}{=}2\pi\delta_{j_1j_2}\). Then \[\begin{align} \sum_{\mathbf{R}\in\text{SC}}e^{i\mathbf{k}\cdot\mathbf{R}} &\overset{\eqref{ed5}}{=} \sum_{n_1=0}^{2N_1-1} \;\sum_{n_2=0}^{2N_2-1}{\exp}\left[i \left( \frac{\ell_1\mathbf{b}_1}{2N_1} + \frac{\ell_2\mathbf{b}_1 }{2N_2} \right)\cdot \left(n_1{\mathbf{a}_1}{+}\,n_2\mathbf{a}_2\right)\right] \nonumber \\[0.5em] &= \sum_{n_1=0}^{2N_1-1} \;\sum_{n_2=0}^{2N_2-1} \exp\left[2\pi i \left( \frac{\ell_1n_1}{2N_1} + \frac{\ell_2n_2}{2N_2} \right)\right] \nonumber \\[0.5em] &= \sum_{n_1=0}^{2N_1-1}\left( \exp \frac{\pi i \ell_1}{N_1} \right)^{n_1} \;\cdot\;\sum_{n_2=0}^{2N_2-1} \left( \exp \frac{\pi i \ell_2}{N_2} \right)^{n_2}\nonumber \\[0.5em] &= \frac{1-\left( \exp \frac{\pi i \ell_1}{N_1} \right)^{2N_1}}{1- \exp \frac{\pi i \ell_1}{N_1} } \cdot \frac{1-\left( \exp \frac{\pi i \ell_2}{N_2} \right)^{2N_2}}{1- \exp \frac{\pi i \ell_2}{N_2} }\nonumber \\[0.5em] &= \frac{1- \exp \left( 2\pi i \ell_1 \right) }{1- \exp \frac{\pi i \ell_1}{N_1} } \cdot \frac{1- \exp \left(2\pi i \ell_2\right) }{1- \exp \frac{\pi i \ell_2}{N_2} } = 0\,. \end{align}\] In the derivation, we have used \(\sum_{n=0}^{N-1}q^n{=}(1{-}q^N)/(1{-} q)\). In particular, this implies 70 as for each \(\beta=1,2\) direction \(k_\beta\) of any of the \(\mathbf{k}_\ell\) with \(\ell=(\ell_\alpha,\ell_\beta)\), we have: \[\begin{align} \sum_{R_\beta} e^{ik_\beta R_\beta}&= \sum_{n_\beta=0}^{2N_\beta-1}{\exp}\left[i \left( \frac{\ell_\beta\mathbf{b}_\beta}{2N_\beta} \right)\cdot \left(n_\beta{\mathbf{a}_\beta}\right)\right] \nonumber \\ &= \frac{1- \exp \left( 2\pi i \ell_\beta \right) }{1- \exp \frac{\pi i \ell_\beta}{N_\beta} }=0 \label{ae4} \end{align}\tag{73}\]

15 Proof that absolutely convergent lattice sum 71 reproduces plane-wave GW Coulomb integrals↩︎

In this Appendix, we prove that the Coulomb matrix computed from the lattice sum from Eq. 71 gives the same Coulomb matrix used in plane-wave \(GW\). We start with the inverse Fourier transform of the Coulomb interaction, \[\begin{align} {\int\limits_{\mathbb{R}^3}} \frac{d\mathbf{G}}{(2\pi)^3}&\, \frac{4\pi}{G^2} \,e^{-i\mathbf{G}(\mathbf{r}-\mathbf{r}')} = \int\limits_0^\infty \frac{dG}{2\pi^2} \int\limits_0^{2\pi}d\varphi \int\limits_0^\pi d\theta \sin\theta\,e^{-iG|\mathbf{r}-\mathbf{r}'|\cos\theta} \nonumber \\[0.1em] &= \int\limits_0^\infty \frac{dG}{\pi} \frac{1}{iG|\mathbf{r}-\mathbf{r}'|} \left(e^{iG|\mathbf{r}-\mathbf{r}'| }-e^{-iG|\mathbf{r}-\mathbf{r}'| } \right)\nonumber \\[0.1em] &= \frac{2}{\pi|\mathbf{r}-\mathbf{r}'|}\,\int\limits_0^\infty dG\; \frac{\sin (G|\mathbf{r}-\mathbf{r}'|)}{G} = \frac{1}{|\mathbf{r}-\mathbf{r}'|}\,.\label{af1} \end{align}\tag{74}\] For 2D systems, the Coulomb operator is truncated in the non-periodic direction to eliminate artificial image interactions [100] and all following steps are analogous.

We assume an absolutely convergent (a.c.) lattice sum for the Coulomb matrix, like Eq. 71 , which we abbreviate as \[\begin{align} V_{PQ}(\mathbf{k}) &= \sum_\mathbf{R}^\text{a.c.} e^{i\mathbf{k}\cdot\mathbf{R}}\, (\varphi_P^\mathbf{0}| \varphi_Q^\mathbf{R}) \,.\label{af2a} \end{align}\tag{75}\] The absolutely convergent lattice summation allows us to interchange the lattice summation and the integration, \[\begin{align} V_{PQ}(\mathbf{k}) &= (\varphi_P^\mathbf{0}| \sum_\mathbf{R}^\text{a.c.} e^{i\mathbf{k}\cdot\mathbf{R}}\, \varphi_Q^\mathbf{R})\,.\label{af2} \end{align}\tag{76}\] Note that if the lattice sum is not absolutely convergent, we cannot conclude Eq. 76 from Eq. 75 and the following derivation does not hold.

For connecting the lattice sum 76 to a plane-wave expression, we insert the Fourier transform of the Coulomb interaction 74 into Eq. 76 : \[\begin{align} &V_{PQ}(\mathbf{k}) = \int\limits_{\mathbb{R}^3} \int\limits_{\mathbb{R}^3}d\mathbf{r}\,d\mathbf{r}'\,\varphi_P^\mathbf{0}(\mathbf{r})\,\frac{ 1}{|\mathbf{r}- \mathbf{r}'|}\,\sum_\mathbf{R}^\text{a.c.} e^{i\mathbf{k}\cdot\mathbf{R}}\,\varphi_Q^\mathbf{R}(\mathbf{r}')\nonumber \\[0.2em] &\overset{\eqref{af1}}{=} \int\limits_{\mathbb{R}^3} \frac{d\mathbf{G}}{(2\pi)^3}\,\frac{4\pi}{G^2} \int\limits_{\mathbb{R}^3} \int\limits_{\mathbb{R}^3} d\mathbf{r}\,d\mathbf{r}'\,\varphi_P^\mathbf{0}(\mathbf{r})\,e^{-i\mathbf{G}(\mathbf{r}-\mathbf{r}')}\,\sum_\mathbf{R}^\text{a.c.} e^{i\mathbf{k}\cdot\mathbf{R}}\,\varphi_Q^\mathbf{R}(\mathbf{r}')\nonumber \\[0.35em] &\overset{\eqref{af5}}{=} \int\limits_{\mathbb{R}^3} \frac{d\mathbf{G}}{(2\pi)^3}\,\frac{4\pi}{G^2} \,\varphi_P^*(\mathbf{G})\sum_\mathbf{R}^\text{a.c.} \,e^{i\mathbf{k}\cdot\mathbf{R}} {\int\limits_{\mathbb{R}^3} } \,d\mathbf{r}'\,e^{i\mathbf{G}\mathbf{r}'}\,\varphi_Q^\mathbf{0}(\mathbf{r}'-\mathbf{R})\nonumber \\[0.35em] &= \int\limits_{\mathbb{R}^3} \frac{d\mathbf{G}}{(2\pi)^3}\,\frac{4\pi}{G^2} \,\varphi_P^*(\mathbf{G})\sum_\mathbf{R}^\text{a.c.} \,e^{i\mathbf{k}\cdot\mathbf{R}} {\int\limits_{\mathbb{R}^3} } \,d\mathbf{r}'\,e^{i\mathbf{G}(\mathbf{r}'+\mathbf{R})}\,\varphi_Q (\mathbf{r}' )\nonumber \\[0.35em] &\overset{\eqref{af5}}{=} \int\limits_{\mathbb{R}^3} \frac{d\mathbf{G}}{(2\pi)^3}\,\frac{4\pi}{G^2} \,\varphi_P^*(\mathbf{G})\sum_\mathbf{R}^\text{a.c.} \,e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{R}} \varphi_Q (\mathbf{G})\nonumber \\[0.35em] &\overset{\eqref{af6}}{=} \int\limits_{\mathbb{R}^3} \frac{d\tilde{\mathbf{G}}}{(2\pi)^3}\,\frac{4\pi}{\tilde{G}^2} \,\varphi_P^*(\tilde{\mathbf{G}}) \varphi_Q (\tilde{\mathbf{G}}) \sum_{\mathbf{G}}^\text{rlv} \Omega_\text{BZ}\,\delta(\mathbf{G}-(\tilde{\mathbf{G}}+\mathbf{k}))\nonumber \\[0.5em] &\overset{\eqref{af7}}{=} \sum_{\mathbf{G}}^\text{rlv} \frac{1}{\Omega}\,\frac{4\pi}{ |\mathbf{G}+\mathbf{k}|^2} \,\varphi_P^*(\mathbf{G}+\mathbf{k})\, \varphi_Q ( \mathbf{G}+\mathbf{k})\,.\label{af4} \end{align}\tag{77}\] In the course of deriving Eq. 77 , we have used the Fourier transforms of the Gaussian basis functions \[\begin{align} \varphi_P(\mathbf{G}+\mathbf{k}) &= \int\limits_{\mathbb{R}^3} d\mathbf{r}\; e^{i(\mathbf{G}+\mathbf{k})\cdot\mathbf{r}}\,\varphi_P^\mathbf{0}(\mathbf{r}) \,,\label{af5} \end{align}\tag{78}\] the identity \[\begin{align} \sum_\mathbf{R}\,e^{i\tilde{\mathbf{G}}\cdot\mathbf{R}} = \sum_{\mathbf{G}}^\text{rlv} \Omega_\text{BZ}\,\delta(\mathbf{G}-\tilde{\mathbf{G}}) \label{af6} \end{align}\tag{79}\] (using the abbreviation rlv for a reciprocal lattice vector \(\mathbf{G}\); the vector \(\tilde{\mathbf{G}}\) is arbitrary), and that the Brillouin zone volume \(\Omega_\text{BZ}\) is connected to the unit cell volume \(\Omega\) via \[\begin{align} \Omega_\text{BZ} = \frac{(2\pi)^3}{\Omega}\label{af7}\,. \end{align}\tag{80}\]

Note that Eq. 77 has been used in Ref. [122] to compute the Coulomb matrix element \(V_{PQ}(\mathbf{k})\) of Gaussians \(P,Q\).

Eq. 77 is closely related to the Coulomb matrix in a plane-wave basis, given by \[\begin{align} V_{\mathbf{G}\mathbf{G}'}(\mathbf{k}) = \frac{4\pi}{|\mathbf{k}+\mathbf{G}|^2}\,\delta_{\mathbf{G}\mathbf{G}'} \,, \label{ed14a} \end{align}\tag{81}\] which arises in plane-wave \(GW\) implementations. Eq. 81 and Eq. 77 are connected via a basis transformation from the plane-wave basis to the Gaussian basis. This shows that our absolutely convergent lattice summation scheme 71 yields Coulomb matrix elements fully consistent with 81 , ensuring compatibility with plane-wave-based \(GW\) implementations.

16 Derivation of Eq. 41 for computing \(\Sigma^{\text{c},\mathbf{R}}_{\lambda\sigma}(i\tau)\)↩︎

For the derivation of Eq. 41 , we review \[\begin{align} \Sigma_{n\mathbf{k}}(i\tau) &= \braket{\psi_{n\mathbf{k}}|\Sigma(i\tau)|\psi_{n\mathbf{k}}} \\[0.5em]& = \int\limits_\text{cell} d\mathbf{r}\int\limits_{\mathbb{R}^3}d\mathbf{r}'\, \psi^*_{n\mathbf{k}} (\mathbf{r})\, \Sigma(\mathbf{r},\mathbf{r}',i\tau)\, \psi_{n\mathbf{k}} (\mathbf{r}') \\[0.5em] & = \sum_{\mu\nu} C_{\mu n}^*(\mathbf{k}) \, \Sigma_{\mu\nu}(\mathbf{k},i\tau)\,C_{\nu n}(\mathbf{k}) \,. \end{align}\] To arrive at the last line, which is Eq. 43 , we have used the basis expansion 11 and the \(\mathbf{k}\leftrightarrow\mathbf{R}\) transformation 42 , \[\begin{align} &\Sigma_{\mu\nu}(\mathbf{k},i\tau) = \sum_\mathbf{R}^\text{SC} e^{i\mathbf{k}\cdot\mathbf{R}}\; \Sigma_{\mu\nu}^{\mathbf{R}}(i\tau) \,. \end{align}\] We further have \[\begin{align} &\Sigma_{\mu\nu}^{\mathbf{R}}(i\tau) = \int d\mathbf{r}\,d\mathbf{r}'\, \phi_\mu^\mathbf{0}(\mathbf{r})\, G(\mathbf{r},\mathbf{r}',i\tau) \, W(\mathbf{r},\mathbf{r}',i\tau) \, \phi_\nu^\mathbf{R}(\mathbf{r}')\nonumber \\ &\overset{\eqref{a5}}{=}\sum_{\lambda \mathbf{S}_1,\sigma\mathbf{S}_2}G_{\lambda\sigma}^{\mathbf{S}_2-\mathbf{S}_1} \int d\mathbf{r}\,d\mathbf{r}'\, \phi_\mu^\mathbf{0}(\mathbf{r}) \phi_\lambda^{\mathbf{S}_1}(\mathbf{r}) W(\mathbf{r},\mathbf{r}',i\tau)\phi_\sigma^{\mathbf{S}_2}(\mathbf{r}') \phi_\nu^\mathbf{R}(\mathbf{r}')\,.\label{eqe1} \end{align}\tag{82}\] Inserting periodic RI 48 /49 for the products \(\phi_\mu^\mathbf{0}(\mathbf{r}) \phi_\lambda^{\mathbf{S}_1}(\mathbf{r})\) and \(\phi_\sigma^{\mathbf{S}_2}(\mathbf{r}') \phi_\nu^\mathbf{R}(\mathbf{r}')\) into Eq. 82 leads to Eq. 41 .

17 Spin-orbit coupling from HGH pseudopotentials↩︎

We employ spin-orbit coupling (SOC) from Hartwigsen-Goedecker-Hutter (HGH) pseudopotentials [54], [56], \[\begin{align} \hat{V}^{\operatorname{SOC}}(\mathbf{r},\mathbf{r}') &= \sum_l \Delta V^\text{SO}_l(\mathbf{r},\mathbf{r}')\, \frac{\hbar}{2}\,\mathbf{L}\cdot\hat{\boldsymbol{\sigma}} \\ \Delta V^\text{SO}_l(\mathbf{r},\mathbf{r}')&=\sum_{i,j=1}^3\sum_{m=-l}^l Y_{lm}(\theta,\varphi)\,p^l_i(r)\,k^l_{ij}\,p^l_j(r')Y_{lm}^*(\theta',\varphi') \label{e2} \end{align}\tag{83}\] where \(\hat{V}^{\operatorname{SOC}}(\mathbf{r},\mathbf{r}')\) is the non-local SOC part to the HGH pseudopotential, \(\mathbf{L}{=}{-}i\hbar\,\mathbf{r}{\times}\nabla_\mathbf{r}\) the angular momentum, \(\hat{\boldsymbol{\sigma}}\) the Pauli matrices, \(Y_{lm}\) the spherical harmonics, \(\mathbf{r}{=}(r,\theta,\varphi)\) are spherical coordinates, \(p^l_i(r)\) are tabulated Gaussian functions and \(k_{ij}^l\) are tabulated SOC parameters. We compute SOC matrix elements in the Gaussian basis as \[\begin{align} V_{\mu\nu,\sigma\sigma'}^{\operatorname{SOC}}(\mathbf{k}) &= \sum_\mathbf{R}e^{i\mathbf{k}\cdot\mathbf{R}}{\int} d\mathbf{r}\,d\mathbf{r}'\, \phi_\mu^\mathbf{0}(\mathbf{r}) \braket{\sigma|\hat{V}^{\operatorname{SOC}}(\mathbf{r},\mathbf{r}')|\sigma'} \phi_\nu^\mathbf{R}(\mathbf{r}') \end{align}\] where \(\sigma,\sigma'{\in}\{\uparrow,\downarrow\}\) is the spin quantum number along the \(z\)-quantization axis. We compute the SOC matrix elements in the Bloch basis, \[\begin{align} V_{nn',\sigma\sigma'}^{\operatorname{SOC}}(\mathbf{k}) &= \sum_{\mu\nu} [C_{\mu n}^\sigma(\mathbf{k})]^* \; V_{\mu\nu,\sigma\sigma'}^{\operatorname{SOC}}(\mathbf{k})\; C_{\nu n'}^{\sigma'}(\mathbf{k})\,. \end{align}\]

In principle, the SOC matrix \(V_{nn',\sigma\sigma'}^{\operatorname{SOC}}(\mathbf{k})\) should be computed between all states of the Bloch basis. However, we have observed numerical instabilities originating from nonvanishing pseudopotential overlaps between different atoms. We found that this issue can be circumvented by restricting the correction to an energy window \(E_{\operatorname{W}}\) around the valence band maximum \(\varepsilon_{\operatorname{VBM}}\) and the conduction band minimum \(\varepsilon_{\operatorname{CBM}}\), so that we only include \(V_{nn',\sigma\sigma'}^{\operatorname{SOC}}(\mathbf{k})\) for bands \(n,n'\) at \(\mathbf{k}\) that obey \[\label{SOC:enerwin} \varepsilon_{n\mathbf{k}}^{G_0W_0}, \varepsilon_{n'\mathbf{k}}^{G_0W_0} \in [\varepsilon_{\operatorname{VBM}}-E_{\operatorname{W}}/2,\varepsilon_{\operatorname{CBM}}+E_{\operatorname{W}}/2]\,.\tag{84}\] We provide in Fig. 9 the evolution of the K-point splitting in PBE+SOC calculation with respect to this energy window \(E_\text{W}\). We obtain a stable SOC splitting within 1-2 meV for \(E_\text{W}{=}10\) eV and \(E_\text{W}{=}20\) eV, validating our approach 84 . For WSe\(_2\) with aug-DZVP-MOLOPT and aug-TZVP-MOLOPT (Fig. 9h), we observe deviations of more than 10 meV when increasing the energy window to \(E_\text{W}{=}40\) eV and \(E_\text{W}{=}100\) eV, so that we chose a window of 20 eV for these calculations, and 40 eV for the other ones for the SOC calculations reported in the main text.

We build the single-particle Hamiltonian with SOC, \[\begin{align} h^{{G_0W_0}\text{+SOC}}_{n\sigma,\,n'\sigma'}(\mathbf{k}) = \delta_{nn'}\,\delta_{\sigma\sigma'}\,\varepsilon_{n\mathbf{k}}^{G_0W_0} + V_{nn',\sigma\sigma'}^{\operatorname{SOC}}(\mathbf{k})\,. \end{align}\] We diagonalize \(\mathbf{h}^{{G_0W_0}\text{+SOC}}(\mathbf{k})\) to obtain the band structure \(\varepsilon^{{G_0W_0}\text{+SOC}}_{j\mathbf{k}}\) and coefficients \(C_{n \sigma }^{(j)}(\mathbf{k})\) in a perturbative manner, \[\begin{align} \sum_{n'\sigma'} h^{{G_0W_0}\text{+SOC}}_{n\sigma,\,n'\sigma'}(\mathbf{k})\; C_{n'\sigma'}^{(j)}(\mathbf{k}) = \varepsilon^{{G_0W_0}\text{+SOC}}_{j\mathbf{k}} C_{n \sigma }^{(j)}(\mathbf{k})\,. \end{align}\]

To demonstrate the effects of this correction, we provide in the following the PBE band structures of all monolayers with and without SOC (Fig. 10). We observe that SOC splits the top valence bands, predominantly at the K-point, in agreement with the literature [9].

18 GW band gaps with plasmon-pole model (BerkeleyGW)↩︎

In this Appendix, we present additional data of the band gaps for the TMDCs benchmarked in this work using BerkeleyGW (Table [t2]). We compare the generalized plasmon pole model (GPP) [97] with the full frequency implentation (data of the full-frequency calculation is reproduced from Table [t1]). Some available results in the literature are given in Table [t395GPP].

[H]
  \fontsize{10}{12}\selectfont
 \caption{{$G_0W_0$@PBE+SOC bandgap (in eV) using the plane-wave code BerkeleyGW \cite{Hybertsen1986, Deslippe2012, Barker2022} of monolayer $\text{MoS}_\text{2}$, $\text{MoSe}_\text{2}$, WS$_\text{2}$ and WS$\text{e}_\text{2}$. The bandgap is extracted at the K point.}}
{
\renewcommand{\arraystretch}{1.2}
    \begin{tabular}{lcccc}
    \hline
    Calculation & MoS$_2$ & MoSe$_2$ & WS$_2$ & WSe$_2$ 
    \\
    \hline
BerkeleyGW - GPP & 2.53  & 2.12 & 2.53 & 2.13 \\
BerkeleyGW - Full frequency  & 2.28 & 1.98 & 2.36 & 2.05 
      \\
 \hline
    \end{tabular}
}
\label{t2}

19 Numerical results on the TMD band gaps in the literature↩︎

This section contains a comparison of the results of the paper for the band gaps with the available literature on the topic in Table [t395FF] for full frequency calculations and Table [t395GPP] for plasmon-pole calculations, staying at the \(G_0W_0\)@PBE+SOC level. One can easily see a wide range of values, which can be ascribed to the different values for the lattice constants between different calculations, but also to differences in the pseudopotential and SOC calculation method used.

20 Benchmark study of MoS\(_2\)↩︎

In this section, we add a benchmark study of the convergence of the MoS\(_2\) band gap, in Fig. 11. As one can see, the convergence behaviour with respect to all parameters is similar to the one observed for the case of WSe\(_2\) in Fig. 2, although the k-point convergence is slightly more unstable (while still being converged within 10 meV)

References↩︎