Introduction

The isolation of graphene in 2004 can be regarded as a milestone in materials science that initiated the research field of atomically thin 2D materials1. Compared to their 3D counterparts, 2D materials have a higher surface-to-volume ratio, making them ideal candidates for catalysts and sensors2,3. Due to the confinement of electrons, holes, phonons, and photons in the 2D plane, the electronic, thermal, and optical properties of 2D materials present unusual features not found in their 3D counterparts4,5,6. For instance, their electronic structure – especially band gaps – can be easily adjusted by acting on the vertical quantum confinement through, e.g., the number of atomic layers, or external perturbations, such as an external electric field, and strain7,8. The sensitivity to strain, i.e., to structural details, implies that 2D materials also exhibit strong electron-phonon coupling8. In addition, exciton binding energies are significantly larger than in 3D materials, and they can be tuned by the dielectric environment, e.g., by encapsulation or deposition on substrates9,10,11. All these characteristics make them outstanding components in novel applications for electronics and optoelectronics12,13,14,15,16,17.

For a deep understanding of 2D materials, an accurate description of their band-structure is a must. Many-body perturbation theory within the GW approach has become the state-of-the-art for ab initio electronic-structure calculations of materials. In this sense, many studies have employed GW to investigate the electronic properties of 2D materials18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65. Surprisingly, as illustrated in Fig. 1 for monolayer MoS2, they show a wide dispersion in the fundamental band gap. The same has been found for a number of 2D materials that have been extensively studied in the last years. Results for MoS2, MoSe2, MoTe2, WS2, WSe2, BN, and phosphorene are summarized in the Supplementary Information. In the extreme cases of MoS2, WS2, WSe2, and BN, the calculated band gaps are scattered between 2.31–2.97, 2.43–3.19, 1.70–2.89, and 6.00–7.74 eV, respectively; in the worst case, the deviation (ratio between largest and smallest values) is as much as 61%. Moreover, for some materials, such as for MoS2, MoTe2, WS2, WSe2, and BN, not even the gap character is uniquely obtained—being direct or indirect, depending on the details of the calculation.

Many factors contribute to this confusing and unsatisfactory situation:

  • In various works, different geometries have been adopted. In this context, it must be said that the properties of 2D materials are highly sensitive to structural parameters18,19,66. Small changes in the lattice constant a already have a large impact on the energy gap, as seen in Fig. 1 and Supplementary Tables 4–10. Moreover, often the lattice parameter alone is not sufficient to unambiguously determine the structure of a 2D material. For instance, phosphorene is characterized by four structural parameters; transition metal dichalcogenides require, besides a, the distance between two chalcogens (for MoS2, dSS as depicted in Fig. 2). These “other” structural parameters have a notable effect on the electronic properties as well66,67. Unfortunately, in several studies, only the lattice parameter is reported, which prevents not only a fair comparison between published results but also reproducibility.

  • A second reason can be attributed to the various ways of performing G0W0 calculations. First, there is the well-known starting problem68,69,70,71,72,73,74,75,76. Then, especially for 2D materials, G0W0 energy gaps converge very slowly with respect to the vacuum thickness and the number of k-points. Even slabs with a vacuum layer of 60 Å together with a 33 × 33 × 1 k-grid have been shown to be insufficient to obtain fully converged results28. However, by truncating the Coulomb potential, convergence can be achieved with a reasonable amount of vacuum49. The number of k-points can be drastically reduced by an analytic treatment of the q = 0 singularity of the dielectric screening or by using nonuniform k-grids24,38.

  • Last but not least, also spin-orbit coupling (SOC) plays an important role in many cases. Besides decreasing the size of the fundamental gap mainly through a splitting of the valence band77, in some 2D materials and for certain geometries and methods, disregarding or including this effect may change its character from indirect to direct or vice versa66.

    Fig. 1: Literature results for the G0W0 energy gap of MoS2 at the K point of the Brillouin zone as a function of the in-plane lattice parameter a.
    figure 1

    They are obtained with and without SOC (filled and open symbols, respectively), as well as with and without truncation of the Coulomb potential (triangles and circles, respectively), using LDA, PBE, and HSE as starting points (green, orange and blue, respectively). Supplementary Table 4 summarizes these data.

    Full size image
    Fig. 2: Top view (left) and side view (middle) of the MoS2 slab geometry, determined by the in-plane lattice constant a and the distance between sulfur atoms, dSS.
    figure 2

    The Mo-S bond length dMoS and the S-Mo-S angle θ are shown as well. L is the unit-cell size along the out-of-plane direction z, including a vacuum layer. High-symmetry points and paths used in the band-structure plots (Fig. 4) are highlighted in the BZ (right panel).

    Full size image

In this manuscript, we address all these issues and provide a benchmark data set of density functional theory (DFT) and G0W0 calculations, taking monolayer MoS2 as a representative 2D material. Due to its unique properties, it can be considered the most important 2D material after graphene. MoS2 exhibits high electron mobility14,78; moderate SOC that can be exploited in spin- and valleytronics79,80,81,82,83,84,85; a direct fundamental band gap with intermediately strong exciton binding, which is suitable for (opto)electronic devices operating at room temperature14,29,31,38,78,86. For these reasons, there are many experimental and theoretical works in the literature that investigate MoS2, allowing for a better comparison with our results.

We employ the linearized augmented planewave + local orbital (in short LAPW+LO) method as implemented in the exciting code. LAPW+LO is known to achieve ultimate precision for solving the Kohn-Sham (KS) equations of DFT87 and high-level GW results88. Besides the local and semilocal DFT functionals LDA89 and PBE90,91 respectively, we include HSE0692,93,94 both for geometry optimization and as a starting point for G0W0. So far, HSE06 has not often been used for such calculations of MoS220,22,29,33, and to the best of our knowledge, neither a Coulomb truncation nor an adequate treatment of the singularity at q = 0 was applied. For brevity, hereafter, we will refer to HSE06 as HSE. In our G0W0 calculations, we truncate the Coulomb potential28,49,95, and apply a special analytical treatment for the q = 0 singularity24. Moreover, we investigate the role of SOC at all levels. We carefully evaluate the impact of all these elements and conclude what leads to the most reliable electronic structure of this important material. Besides a detailed analysis of energy gaps, we address effective masses and spin-orbit splittings.

Results

Ground-state geometries

The geometry of MoS2, depicted in Fig. 2, is determined by the in-plane lattice parameter a and the distance between sulfur atoms, dSS. In Table 1, we list these structural parameters as obtained with LDA, PBE, and HSE, and include the Mo-S bond length dMoS and the angle θ between Mo and S atoms as well. As expected, LDA underestimates the lattice spacing, PBE slightly overestimates it, and HSE shows the best performance with respect to experiment. All three exchange-correlation (xc) functionals underestimate the S-S bond length, PBE being closest to its measured counterpart. Comparison with computed literature data reveals good agreement.

Table 1 Equilibrium geometry of MoS2 obtained with LDA, PBE, and HSE, compared with literature values
Full size table

Electronic structure

Table 2 summarizes the energy gaps obtained with different functionals for the different geometries. We consider here the direct gap at the K point (Eg(KK)) as well as the indirect gaps between Γ and K (Eg(ΓK)) and between K and T (Eg(KT)). For the definition of the T point, see Fig. 2. For each geometry and methodology, the bold font highlights the fundamental gap. For the calculations that include SOC, Fig. 3 displays the energy gaps given in Table 2 with respect to the lattice parameter. In the DFT calculations, regardless of the calculation method, Eg(KT) (squares) shows a weak dependence on the geometry. The fundamental gap obtained with LDA, PBE, and HSE, is always direct at K. In G0W0, the fundamental gap is Eg(KT) for the structures with smaller lattice parameter and Eg(KK) for larger lattice constants.

Table 2 Energy gaps (in eV) obtained for the 18 different cases considered
Full size table
Fig. 3: Calculated energy gaps of MoS2 as a function of the lattice parameter a.
figure 3

The values for KK, ΓK, and KT are represented as circles, triangles, and squares, respectively. Dotted (solid) lines stand for DFT (G0W0) results; all include SOC. Note that these values only indirectly reflect the S-S distance dSS.

Full size image

As to be expected and also observed in ref. 38, for a fixed geometry, the energy gaps obtained by LDA and PBE are quite similar, with the largest difference being 0.02 eV. The two functionals also agree on the location of the valence-band maximum (VBM) and the conduction-band minimum (CBm). For both, the band gap is direct if SOC is included, and indirect otherwise for the PBE and HSE geometries. In contrast, HSE gives a direct gap, regardless of whether SOC is considered.

Also G0W0@LDA and G0W0@PBE are very close to each other, the largest difference being 0.04 eV. This can be attributed to the similarity between LDA and PBE when the same geometry is adopted. However, the similarity between G0W0@LDA and G0W0@PBE results is material dependent, as observed in other works96,97,98,99.

For a given geometry, the locations of the VBM and the CBm are independent of the starting point, with the only exception being the HSE geometry when SOC is included. When comparing the three geometries, we encounter three different scenarios. First, for the LDA geometry, the fundamental gap changes from a direct KS gap at K to an indirect QP gap (between K and T), independent of the starting point. Second, for the PBE and HSE geometries, when SOC is disregarded, the indirect gap ΓK obtained withII D LDA and PBE becomes direct and located at K upon applying G0W0. Third, for the HSE geometry, and SOC being included we observe an indirect band gap for G0W0@LDA while it is direct for PBE and HSE as starting points. This can be understood in terms of the small differences between the KK and KT gaps, ΔKT, which are 0.01 eV, 0.05 eV, and 0.15 eV for G0W0@LDA, G0W0@PBE and G0W0@HSE, respectively. Including SOC, splits the conduction band state at T (K) by ~0.07 eV (~3 meV), decreasing ΔKT by ~0.03 eV. This is enough to make ΔKT negative for G0W0@LDA, but not for G0W0@PBE, and G0W0@HSE. A more detailed discussion about ΔKT can be found in Section “Discussion of energy gaps: comparison with experiment”.

Figure 4 shows the band structures obtained for the HSE geometry including SOC. As expected, the differences between G0W0@HSE and HSE bands are less pronounced than those between G0W0@LDA (G0W0@PBE) and LDA (PBE) bands. For all three starting points, the G0W0 corrections are not uniform over all k-points, i.e., a simple scissors approximation is, strictly speaking, not applicable. We will explore this in a more quantitative fashion further below. The SOC splitting in the valence band is zero at the Γ point and increases toward the K point. The SOC effect on the conduction bands is much less pronounced. These observations are in agreement with other theoretical25,29,100,101,102,103 and experimental works86,104.

Fig. 4: Band structures including SOC obtained for the HSE geometry.
figure 4

The parameters ΔKT and ΔΓK are discussed in Section “Discussion of energy gaps: comparison with experiment”.

Full size image

The impact of the self-energy correction on the energy gaps, \(\Delta {{{{\rm{E}}}}}_{{{{\rm{g}}}}}={{{{{\rm{E}}}}}_{{{{\rm{g}}}}}}^{{G}_{0}{W}_{0}}-{{{{\rm{E}}}}}_{{{{\rm{g}}}}}^{{{{\rm{DFT}}}}}\), is shown in Fig. 5 for the case when SOC is included. Clearly, ΔEg is more significant for (semi)local DFT starting points than for HSE. Interestingly, for a given starting point, ΔEg hardly depends on the geometry. For LDA and PBE as the starting points, the ranges of ΔEg(KK), ΔEg(ΓK), and ΔEg(KT) are 0.84–0.90, 1.0–1.1, and 0.61–0.65 eV, respectively. The dependence is even weaker for G0W0@HSE with values of 0.64–0.66, 0.74–0.78, and 0.39–0.40 eV, respectively. Very similar results are observed when SOC is disregarded.

Fig. 5: G0W0 self-energy correction, ΔEg, to the band gaps obtained for different starting points (LDA green, PBE orange, HSE blue) and geometries, including SOC effects.
figure 5

Left, middle, and right panels refer to ΔEg evaluated at the K point, between Γ and K, and between K and T, respectively.

Full size image

Spin-orbit splittings and effective masses

In Table 3, we report for the HSE structure the spin-orbit splitting Δval (Δcond) at the K point for the highest occupied (lowest unoccupied) band. Effective electron (hole) masses me (mh) calculated at the K point along different directions are shown as well. We observe that neither the spin-orbit splittings nor the effective masses are very sensitive to the geometry (see Supplementary Table 11 for more details).

Table 3 Spin-orbit splittings in valence (Δval) and conduction (Δcond) band (in meV) and effective hole (mh) and electron (me) masses, in units of the free electron mass m0, at the K point along different directions, obtained for the HSE geometry
Full size table

Δval exhibits a very narrow spread among all the methods employed here (DFT and G0W0), i.e., a range of 143–149 meV. These values are in excellent agreement with the experimental counterparts of Δval = 130–160 meV86,102,104,105,106,107,108. The value for the conduction band, Δcond, is 3 meV for LDA, PBE, G0W0@LDA, and G0W0@PBE; it is only slightly higher for HSE and G0W0@HSE, namely 4 meV. Again, there is excellent agreement with the measured value of Δcond = 4.3 ± 0.1 meV109.

The effective hole mass, mh, exhibits minor variations, not larger than 0.04m0, when going from KΓ to KM. This is in line with other calculations27,110. The measured value for freestanding MoS2 is mh = (0.43 ± 0.02)m0111 and, apart from LDA and PBE, all theoretical results show excellent agreement. The electron mass, me, is more isotropic than mh. For LDA, it differs by at most 0.02m0 between KΓ and KM. Our values are in line with other calculated results18,21,27,34,101,110,112. The measured counterpart of (0.67 ± 0.08)m0113 is significantly larger than the calculated value reported here and in other theoretical works18,21,27,34,101,110,112. As discussed in ref. 34, the difference could originate from the heavy doping of the measured sample, which may introduce metallic screening.

Discussion of energy gaps: comparison with experiment

The experimental band gap for free-standing MoS2, determined by photocurrent spectroscopy, is 2.5 eV86. In order to compare with experiment, it is important to account for the zero-point renormalization energy of 75 meV114. This means that the theoretical value computed without this correction should be 2.575 ≅ 2.6 eV to match its measured counterpart. For our discussion here, we consider the HSE and PBE geometries, which are closer to experiment than that obtained by LDA. For the following assessment, we refer to the values in Table 2. For the structure optimized with HSE, G0W0 performed on LDA, PBE, and HSE as starting points gives Eg(KK) of 2.52, 2.52, and 2.76 eV, respectively, i.e., G0W0@LDA and G0W0@PBE underestimate the measured value by about 0.08 eV, whereas G0W0@HSE overestimates it by 0.16 eV. However, even though G0W0@LDA agrees better with experiment than G0W0@HSE, it erroneously predicts an indirect band gap Eg(KT) which is 0.02 eV smaller than Eg(KK). G0W0@PBE shows the best agreement with experiment and also predicts the gap to be direct. Interestingly, considering the PBE geometry, as done in several works19,21,22,23,24,28,31,32,33,35,40,41,42,43,112, G0W0@HSE, giving a direct band gap of Eg(KK) = 2.68 eV, agrees best with experiment, the deviation being 0.08 eV only. With LDA and PBE as starting points, the calculated G0W0 band gap is direct, but 0.16 and 0.15 eV, respectively, below the experimental value.

Other relevant aspects of the band-structure concern relative energy differences, in particular ΔKT (introduced above) as well as the maximum energy at the Γ point wrt the VBM at the K point, ΔΓK = Eg(KK) − Eg(ΓK) (see Fig. 4). Experimentally, ΔKT is expected to be ≳ 60 meV101,113 and ΔΓK ≈ 140 meV101,115. Taking our calculations with SOC at the HSE geometry, (see Table 4), the G0W0@HSE value of 0.12 eV reproduces ΔKT best. On the other hand, HSE satisfies ΔΓK best with a value of 0.12 eV (0.02 eV smaller than in experiment), while the values obtained with the other methods differ from experiment by 0.09 eV (G0W0@HSE), −0.09 eV (LDA and PBE), 0.10 eV (G0W0@PBE), and 0.11 eV eV (G0W0@LDA), respectively. At the PBE geometry, G0W0@HSE and G0W0@PBE give the same value for ΔΓK. With an overestimation of 0.07 eV, it is closer to experiment than the value at the HSE geometry. Again, HSE is the only starting point for which G0W0 predicts ΔKT in agreement with experiment.

Table 4 ΔKT and ΔΓK from G0W0 calculations including SOC, compared to measured values
Full size table

In summary, considering the band gap as well as the energy differences ΔKT and ΔΓK, we conclude that at the PBE geometry, G0W0@HSE including SOC shows the best overall agreement with experimental data. Also for the HSE structure, HSE is the best starting point, with results that are overall only slightly worse. Overall, G0W0@HSE at the HSE geometry can be considered more appropriate, since only one xc functional is needed for providing decent results for both, the geometry and the electronic properties and thus the most consistent picture. Also for other materials, HSE has been found to be a superior G0W0 starting point76,116,117,118 compared to LDA and PBE. For such materials with intermediate band gaps,71,72,119,120,121,122,123,124, this functional better justifies the perturbative self-energy correction68,73,116. Figure 5 confirms this for MoS2.

Discussion of energy gaps: comparison with theoretical works

By employing coupled-cluster calculations including singles and doubles125,126, Pulkin et al. obtained energy gaps of Eg(ΓK) = 2.93 eV and Eg(KK) = 3.00 eV with an error bar of ± 0.05 eV112 for the PBE geometry of ref. 103; SOC was not included. For our PBE geometry and also omitting SOC, the G0W0@HSE results are the ones closest to these values, with Eg(ΓK) differing by 0.04 eV and Eg(KK) by 0.24 eV.

For a fair comparison with other published G0W0 values with LDA and PBE as starting points, we restrict ourselves here to results obtained by using a Coulomb truncation in combination with either a special treatment of the q = 0 singularity or a nonuniform k-grid sampling, since these methods ensure well converged gaps. In ref. 24, disregarding SOC and adopting the PBE geometry (a = 3.184 Å, dSS = 3.127 Å), a direct band gap of 2.54 eV was reported for G0W0@PBE which is very close to ours (Eg(KK) = 2.52 eV), i.e., differing by only 0.02 eV. Including SOC and the thus slightly changed PBE geometry (a = 3.18 Å, dSS = 3.13 Å), a G0W0@LDA value of 2.48 eV was obtained in ref. 23; at basically the same geometry (differences in the order of 10−3 Å), our results of Eg(KK) = 2.44 eV is only 0.04 eV smaller.

For a lattice parameter of 3.15 Å, Rasmussen et al. calculated a G0W0@PBE band gap of 2.64 eV without SOC24. For the same lattice constant, but including SOC, Qiu et al. reported a G0W0@LDA band gap of38 of Eg(KK) = 2.59 eV with the plasmon-pole model and Eg(KK) = 2.45 eV with the contour deformation method. In our case, the structure optimized with HSE (a = 3.160 Å) is closest to a = 3.15 Å. For this structure, without including SOC, we compute a G0W0@PBE band gap of Eg(KK) = 2.60 eV, which agrees quite well with the one by ref. 24, differing by less than 0.04 eV. When we include SOC, we obtain Eg(KK) = 2.52 eV with G0W0@LDA, although at this geometry, we obtain an indirect gap that is 18 meV smaller than Eg(KK). This is in good agreement with ref. 38, with a difference of 0.07 meV only.

For the experimental geometry and neglecting SOC, ref. 28 reported values of Eg(KT) = 2.58 eV and Eg(KK) = 2.77 eV for G0W0@LDA. In our case, at the HSE geometry, we obtain Eg(KT) = 2.61 eV and Eg(KK) = 2.60 eV. As the HSE geometry is close to experiment, we may attribute the discrepancies mainly to the different underlying KS states. Indeed, at the LDA level, the energy gaps in ref. 28 are Eg(KK) = 1.77 eV and Eg(ΓK) = 1.83 eV28, while ours are Eg(KK) = 1.73 eV and Eg(ΓK) = 1.71 eV. The values for ΔEg, however, compare fairly well (ΔEg(KK) = 1.00 eV, ΔEg(KT) = 0.6–0.7 eV in ref. 28, compared to ΔEg(KK) = 0.87 eV, ΔEg(KT) = 0.63 eV in the present work).

When it comes to G0W0@HSE, there are only a few results for MoS2 in the literature, neither obtained with Coulomb truncation nor by any special treatment of the q = 0 singularity. For MoS2, these two aspects lead to opposite effects, competing with each other when converging band gaps with respect to the vacuum size and the number of k-points24: Neglecting them, band gaps increase when the vacuum layer is enlarged, whereas denser k-grids make them decrease. Hence, due to fortunate error cancellation, an insufficient vacuum length combined with a coarse k-grid may lead to G0W0 band gaps that agree well with those obtained in a highly converged situation24,29. In ref. 33, using the PBE geometry and taking SOC into account, a KK gap of 2.66 eV was reported. With 15 Å of vacuum, a 6 × 6 × 1 k-grid, and adopting the PBE geometry, in ref. 22, band gaps of 2.05 and 2.82 eV at the HSE and G0W0 levels, respectively, have been obtained. The HSE band gap agrees quite well with ours (2.04 eV), whereas our G0W0@HSE06 gap is 0.14 eV smaller. The band gap of 2.72 eV reported in ref. 29 is based on the experimental lattice parameter of 3.16 Å, a 12 × 12 × 1 k-points grid, and a vacuum layer of 17 Å, and includes SOC effects. The authors state to have chosen these settings to take advantage of error cancellation in the band gap29, and, surprisingly, our band gap of 2.76 eV obtained with G0W0@HSE at the HSE geometry (a = 3.160 Å) agrees quite well.

Closing remarks

We have employed the LAPW+LO method to provide a set of benchmark G0W0 calculations of the electronic structure of two-dimensional MoS2. We have addressed the impact of geometry, SOC, and DFT starting point on the energy gaps, spin-orbit splittings, and effective masses. We find that the self-energy corrections to the band gaps hardly depend on the adopted geometry. As could be expected, employing LDA and PBE as starting points does not make a significant difference when the same structure is used. The best agreement with experimental results is achieved by G0W0@HSE at either the HSE or PBE geometry, considering SOC. The spin-splittings obtained with all methods agree well with experimental results. This also holds true for the effective hole mass, using either HSE or G0W0 on top of any of the considered starting points (LDA, PBE, and HSE). In line with other theoretical works, we highlight the importance of a Coulomb truncation and an adequate treatment of the Coulomb singularity around q = 0 as being fundamental for high-quality calculations. Our findings are expected to be valid for other two-dimensional materials as well.

Methods

Linearized augmented planewave + local orbital methods

The full-potential all-electron code exciting127 implements the family of LAPW+LO methods. In this framework, the unit cell is divided into non-overlapping muffin-tin (MT) spheres, centered at the atomic positions, and the interstitial region in between the spheres. exciting treats all electrons in a calculation by distinguishing between core and valence/semi-core states. For core electrons, assumed to be confined inside the respective MT sphere, the KS potential is employed to solve the Dirac equation which captures relativistic effects including SOC. The KS wavefunctions \(\left\vert {\Psi }_{n{{{\bf{k}}}}}\right\rangle\) for valence and semicore states, characterized by band index n and wavevector k, are expanded in terms of augmented planewaves, \(\left\vert {\phi }_{{{{\bf{G}}}}+{{{\bf{k}}}}}\right\rangle\), and local orbitals, \(\vert {\phi }_{\gamma }\rangle\), as

$$\left\vert {\Psi }_{n{{{\bf{k}}}}}\right\rangle =\mathop{\sum}\limits_{{{{\bf{G}}}}}{C}_{{{{\bf{G}}}}n,{{{\bf{k}}}}}\left\vert {\phi }_{{{{\bf{G}}}}+{{{\bf{k}}}}}\right\rangle +\mathop{\sum}\limits_{\gamma }{C}_{\gamma n,{{{\bf{k}}}}}\vert {\phi }_{\gamma }\rangle .$$
(1)

\(\left\vert {\phi }_{{{{\bf{G}}}}+{{{\bf{k}}}}}\right\rangle\) are constructed by augmenting each planewave with reciprocal lattice vector G, living in the interstitial region, by a linear combination of atomic-like functions inside the MT spheres. In contrast, the LOs \(\vert {\phi }_{\gamma }\rangle\) are non-zero only inside a specific MT sphere. They are used for reducing the linearization error127,128, for the description of semicore states, as well as for improving the basis set for unoccupied states87,88. The quality of the basis set can be systematically improved by increasing the number of augmented planewaves (controlled in exciting by the dimensionless parameter rgkmax) and by introducing more LOs87,127. With all these features, exciting can be regarded as a reference code not only for solving the KS equation129, where it is capable of reaching micro-Hartree precision87, but also for G0W0 calculations88.

G 0 W 0 approximation

In the G0W0 approximation, one takes a set of KS eigenvalues {εnk} and eigenfunctions {Ψnk} as a reference, and evaluates first-order quasi-particle (QP) corrections to the KS eigenvalues in first-order perturbation theory as

$${\varepsilon }_{n{{{\bf{k}}}}}^{{{{\rm{QP}}}}}={\varepsilon }_{n{{{\bf{k}}}}}+{Z}_{n{{{\bf{k}}}}}\langle {\Psi }_{n{{{\bf{k}}}}}| \Sigma ({\varepsilon }_{n{{{\bf{k}}}}})-{V}_{{{{\rm{xc}}}}}| {\Psi }_{n{{{\bf{k}}}}}\rangle ,$$
(2)

where Znk is the renormalization factor, Vxc is the xc potential, and Σ is the self-energy. The latter is given as the convolution

$$\Sigma ({{{\bf{r}}}},{{{{\bf{r}}}}}^{{\prime} },\omega )=\frac{{{{\rm{i}}}}}{2\pi }\int\,G({{{\bf{r}}}},{{{{\bf{r}}}}}^{{\prime} },\omega +{\omega }^{{\prime} })W({{{\bf{r}}}},{{{{\bf{r}}}}}^{{\prime} },{\omega }^{{\prime} })\,d{\omega }^{{\prime} },$$
(3)

with G being the single-particle Green function and W the screened Coulomb potential.

In this non-selfconsistent method, the quality of the QP eigenvalues may depend critically on the starting point. In many cases, LDA and PBE have been proven to be good starting points for G0W0, leading to QP energies that agree well with experiments73,130,131. However, e.g., for materials containing d electrons, such as Mo, hybrid functionals, like HSE, usually provide an improved reference for QP corrections compared to semilocal functionals71,72,119,120,121,122,123,124. Here, we evaluate the quality of each of these three as a starting point.

Coulomb truncation

In calculations with periodic boundary conditions, for treating 2D systems, a sufficient amount of vacuum is required to avoid spurious interaction between the replica along the out-of-plane direction. Local and semilocal density functionals have a (nonphysical) asymptotic decay much faster than 1/r, facilitating convergence of unoccupied states –and thus KS gaps– with respect to the vacuum size. In G0W0, the 1/r decay of the self-energy complicates this task. Specifically for MoS2, even a vacuum layer with 60 Å thickness turned out not being sufficient to obtain a fully converged band gap28,49. Truncating the Coulomb potential along the out-of-plane direction z, however, leads to well-converged G0W0 band gaps with a considerably smaller vacuum size28,49.

Here, we truncate the Coulomb potential with a step function along z. Setting the cutoff length to L/2, where L is the size of the supercell along z (Fig. 2), the truncated Coulomb potential can be written in a planewave basis as95:

$${v}_{{{{\bf{G}}}}{{{{\bf{G}}}}}^{{\prime} }}({{{\bf{q}}}})={\delta }_{{{{\bf{G}}}}{{{{\bf{G}}}}}^{{\prime} }}\frac{4\pi }{{Q}^{2}}\left[1-{{{{\rm{e}}}}}^{-{Q}_{xy}L/2}\cos ({G}_{z}L/2)\right],$$
(4)

where Q = q + G, and q is a vector in the first Brillouin zone (BZ).

Treatment of the q = 0 singularity

On the down-side, truncating the Coulomb interaction slows down the convergence in terms of k-points because of the non-smooth behavior of the dielectric function around the singularity at q = 023,24,41. To bypass this problem, we follow an analytical treatment of Wc, the correlation part of W, close to the singularity as proposed in ref. 24. Without this treatment, the correlation part of the self-energy Σc(ω) can be written as132

$$\begin{array}{l}\langle {\Psi }_{n{{{\bf{k}}}}}| {\Sigma }^{{{{\rm{c}}}}}(\omega )| {\Psi }_{n{{{\bf{k}}}}}\rangle =\frac{{{{\rm{i}}}}}{2\pi }\mathop{\sum}\limits_{mij}\int\nolimits_{-\infty }^{\infty }d{\omega }^{{\prime} }\frac{1}{{N}_{{{{\bf{q}}}}}}\mathop{\sum}\limits_{{{{\bf{q}}}}}\frac{1}{\omega +{\omega }^{{\prime} }-{\tilde{\epsilon }}_{m{{{\bf{k-q}}}}}}\\\qquad\qquad\qquad\qquad\;\;{\left[{M}_{nm}^{i}({{{\bf{k}}}},{{{\bf{q}}}})\right]}^{* }{M}_{nm}^{j}({{{\bf{k}}}},{{{\bf{q}}}}){W}_{ij}^{{{{\rm{c}}}}}({{{\bf{q}}}},{\omega }^{{\prime} }),\end{array}$$
(5)

where \({\tilde{\epsilon }}_{n{{{\bf{k}}}}}={\epsilon }_{n{{{\bf{k}}}}}+{{{\rm{i}}}}\,\eta \,{{{\rm{sign}}}}({E}_{{{{\rm{F}}}}}-{\epsilon }_{n{{{\bf{k}}}}})\) and EF the Fermi energy. \({M}_{nm}^{i}({{{\bf{k}}}},{{{\bf{q}}}})\) are the expansion coefficients of the mixed-product basis, an auxiliary basis to represent products of KS wavefunctions. Like LAPWs and LOs, they have distinct characteristics in the MT spheres and interstitial region132,133. To treat the q = 0 case separately, the corresponding term in Eq. (5) is replaced by

$$\frac{1}{\omega +{\omega }^{{\prime} }-{\tilde{\epsilon }}_{m{{{\bf{k}}}}}}{\left[{M}_{nm}^{i}({{{\bf{k}}}},{{{\bf{0}}}})\right]}^{* }{M}_{nm}^{j}({{{\bf{k}}}},{{{\bf{0}}}})\frac{1}{{\Omega }_{0}}{\int}_{{\Omega }_{0}}{{{\rm{d}}}}{{{\bf{q}}}}\,{W}_{ij}^{{{{\rm{c}}}}}({{{\bf{q}}}},{\omega }^{{\prime} }),$$
(6)

where Ω0 is a small region around q = 0. Analytical expressions for \({W}_{ij}^{{{{\rm{c}}}}}({{{\bf{q}}}},{\omega }^{{\prime} })\) in the limit q → 024 are then employed to calculate the integral in Eq. (6).

Spin-orbit coupling

In this study, SOC is treated via the second-variational (SV) scheme134. In LDA and PBE calculations, the conventional SV approach is utilized127,135,136,137: First, the scalar-relativistic problem, i.e., omitting SOC, is solved. A subset of the resulting eigenvectors is then used as basis set for addressing the full problem. The number of eigenvectors is a convergence parameter. In this work, the SOC term is evaluated with the zero-order regular approximation (ZORA)138,139.

A ground-state calculation with HSE is performed via a nested loop140: In the outer loop, the nonlocal exchange is computed for a subset of KS wavefunctions, using a mixed-product basis. The inner loop solves the generalized KS equations self-consistently, where only the local part of the effective potential is updated in each step. Within this inner loop, the SV scheme is applied self-consistently to incorporate SOC. The corresponding term, evaluated within the ZORA, is based on PBE, which is justified by the minimal contribution due to the gradient of the nonlocal potential141,142,143.

G0W0 calculations with SOC are performed in two steps on top of ground-state calculations that include SOC. First, the QP energies are computed as explained in Section “G0W0 approximation”, using the scalar-relativistic KS eigenvalues and eigenvectors. In the second step, the obtained QP energies are used, together with the SV KS eigenvectors, to evaluate SOC through the diagonalization of the SV Hamiltonian. This approach is sufficient in the case of MoS2 since SOC does not cause band inversion133,143.

Computational details

In our calculations, we employ the all-electron full-potential code exciting127. The only exception is for obtaining the HSE equilibrium geometry, where FHI-aims144,145 is used, since so far exciting lacks geometry relaxation with hybrid functionals. Even though exciting and FHI-aims implement very different basis sets to expand the KS wavefunctions, the two codes have been shown to be among those with the best mutual agreement129. Moreover, a comparison of energy gaps and geometry relaxations for MoS2 confirms this finding (see Supplementary Information - Section I).

In all calculations, the unit-cell size L along the out-of-plane direction (Fig. 2) is set to 30 bohr. Different flavors of xc functionals are applied, namely LDA, PBE, and HSE. In the latter, we use the typical parameters94, i.e., a mixing factor of α = 0.25 and a screening range of ω = 0.11 bohr−1. To determine the respective equilibrium geometries, we relax the atomic positions until the total force on each ion is smaller than 10 μHa bohr−1. For these geometries, the electronic structure is calculated with all three functionals, with and without SOC, giving rise to a set of 18 calculations. These calculations are followed by G0W0 calculations, taking the respective DFT solutions as starting points.

The dimensionless parameter rgkmax that controls the exciting basis-set size is set to 8. In the LDA and PBE calculations, we use a 30 × 30 × 1 k-grid. In HSE and G0W0 calculations, we employ 400 empty states and an 18 × 18 × 1 k-grid. In G0W0, the correlation part of the self-energy is evaluated with 32 frequency points along the imaginary axis, and then analytic continuation to the real axis is carried out by means of Pade’s approximant. For the bare Coulomb potential, we use a 2D cutoff95 combined with a special treatment of the q = 0 singularity as described in Section “Treatment of the q = 0 singularity”. Furthermore, we carefully determine the minimal set of LOs, sufficient to converge at least the lowest 400 KS states. This is achieved with 2 and 6 LOs for sulfur s and p states, respectively, as well as 3, 6, and 10 LOs for molybdenum s, p, and d states, respectively. Overall, we estimate a numerical precision of 50–100 meV in the energy gaps obtained with our G0W0 calculations. To determine effective masses, we use parabolic fits within a range of 0.05 Å−1 around the VBM and the CBm at the K point of the BZ. Depending on the respective case, K can host either global or local extrema.