London Dispersion‐Corrected Density Functionals Applied to van der Waals Stacked Layered Materials: Validation of Structure, Energy, and Electronic Properties

Most density functionals lack to correctly account for long‐range London dispersion interactions, and numerous a posteriori correction schemes have been proposed in recent years. In van der Waals structures, the interlayer distance controls the proximity effect on the electronic structure, and the interlayer interaction energy indicates the possibility to mechanically exfoliate a layered material. For upcoming twisted van der Waals heterostructures, a reliable but efficient and scalable theoretical scheme to correctly predict the interlayer distance is required. Therefore, the performance of a series of popular London dispersion corrections combined with computationally affordable density functionals is validated. As reference data, the experimental interlayer distance of layered bulk materials is used, and corresponding interlayer interaction energies are calculated using the random phase approximation. We demonstrate that the SCAN‐rVV10 and PBE‐rVV10L functionals predict interlayer interaction energies and interlayer distances of the studied layered systems within the range of the defined error limits of 10 meV per atom and 0.12 Å, respectively. Semi‐empirical and empirical dispersion‐corrected functionals show significantly larger error bars, with PBE+dDsC performing best with comparable quality of geometries, but with higher interlayer interaction energy error limits of ≈20 meV per atom.


Introduction
With the isolation and investigation of single-layer graphene [1] it was demonstrated that 2D crystals possess intriguing properties that are not present in their layered bulk counterparts. The most widely studied 2D crystals are graphene, hexagonal boron nitride (hBN), and transition metal dichalcogenides (TMDCs). [2][3][4] As layered bulk crystals, they exhibit strong bonds within the individual layers, but only weakly attractive interlayer interactions, known as London dispersion interactions. [5,6] The latter dominate the interlayer attraction also in van der Waals (vdW) heterostructures, [7,8] such as MoSe 2 /WSe 2 , [9] hBN/TBG [10] (twisted bilayer graphene), and WSe 2 /graphene. [11] The latter examples strikingly showed that the interlayer interaction has a severe influence on the electronic properties of the stacked materials. Layered materials are characterized by their interlayer distance and interlayer interaction energy. Both are controlled by the London dispersion interactions (LDIs) between the layers. Even though there may also be some additional contributions to the interlayer interaction, those are minor and have a rather small impact (for details, Björkman et al.'s study [12] ). To precisely predict the electronic structure of vdW heterostructures, the accurate description of their structure is required, which in turn needs a high-quality potential energy surface. While the precise value of the interlayer interaction energy is not necessarily required for the prediction of the electronic structure of vdW heterostructures, this quantity is very important as it determines if it is possible to mechanically exfoliate layered structures, or to intercalate them with molecules or ions.
Density-Functional Theory (DFT) [13,14] is one of the most widely used methods for electronic structure calculations in the several disciplines of physics, chemistry, and materials science due to its generally high predictive power, robustness, and computational affordability. [15,16] However, DFT calculations strongly depend on exchange-correlation (XC) functional choice. Despite their good performance in dense periodic systems and isolated molecules, present DFT functionals fail to account for static correlation, including LDIs. [12] High-level methods like the Random Phase Approximation (RPA) [17][18][19] within the Adiabatic www.advancedsciencenews.com www.advtheorysimul.com Connection Fluctuation Dissipation Theorem (ACFDT) [20][21][22] and Quantum Monte Carlo (QMC) [23] naturally include static correlation, but they are not practical for complex vdW heterostructures due to their immense computational demand. Therefore, it would be extremely beneficial to identify a reliable and computationally affordable DFT functional that properly captures the LDIs.
The incorporation of LDIs into DFT broadly follows two schemes. [24] The first one adds a force-field like (semi-)empirical term to the electronic energy of the employed density-functional. This term can be pairwise, like for example in the methods proposed by Tkatchenko and Scheffler (TS), [25] in Grimme's DFT-Ds, [26][27][28] and in Corminboeuf's dDsC. [29,30] The leading correction term, the C 6 term, scales with the interatomic distance R IJ as R IJ -6 . In addition to C 6 , the Grimme corrections D3 and D4 take also into account C 8 and the C 9 (three-body ATM) terms, the latter being optional in D3. To compensate for the short-range divergence of the LDI correction terms, damping functions are used. In the density-dependent dispersion-correction (dDsC) method employs a unique damping function and the coefficients of the attractive C 6 term are based on Becke and Johnson's exchange-dipole-moment (XDM) [34][35][36] formalism. Another popular additive approach is Tkatchenko's Many-Body-Dispersion (MBD) [31][32][33] correction. We employ here the most recent MBD-NL variant.
The second scheme directly includes nonlocal electron interactions into the exchange-correlation functional. Examples for these vdW functionals include Dion's vdW-DF, [37] and the revised Vydrov-van Voorhis (rVV10) [38] extension to GGA's and meta-GGAs, including a variant optimized for layered materials (PBE-rVV10L). [39] All these methods are computationally much more affordable compared to the more rigorous approaches beyond DFT that natively incorporate LDIs. However, the methods of the first scheme are more approachable due to their simple implementation and negligible computational cost. An additional, less popular way to correct for LDIs is the use of dispersion-corrected atom-centered pseudopotentials (DCACP), [40,41] which are not in the scope of this study.
In this study, we evaluate the performance of the nine London dispersion corrections for DFT, eight of them applied to the PBE functional (MBD-NL, [31] TS, D3(BJ)+ATM [27,42,43] (ATM = Axilrod-Teller-Muto three-body term), D3(BJ), D3(0)+ATM, D4, [26] dDsC, rVV10L), and SCAN-rVV10. [44] We quantify their performance on layered materials in terms of interlayer distances and interlayer interaction energies. As these methods produce different optimized geometries, they result in sizeable variations in the electronic structure, including bandgap values and even characteristics, and band curvatures for semiconducting 2D materials. We found that SCAN-rVV10 and PBE-rVV10L functionals give excellent predictions in terms of structures and interlayer energies. Among the empirically corrected functionals, PBE+dDsC performs best with interlayer distances comparable to the vdW functionals, but with a lower predictive power for interaction energies. While all above-mentioned dispersion correction schemes give good geometries and are therefore suitable for structure predictions in vdW heterostructures, the vdW functionals studied in this work outperform the (semi-)empirically corrected PBE ones in particular for interlayer energies.
We further note that the (semi)empirical London dispersion treatments typically systematically under-or overestimate the interlayer distances. While pairwise PBE+TS, PBE+D3(BJ)+ATM, PBE+D3(BJ), and PBE+D4 methods tend to predict shorter interlayer distances, PBE+MBD-NL, PBE+dDsC, PBE+D3(0)+ATM, and PBE-rVV10L methods generally overestimate them (Figure 1). The best performing method is SCAN-rVV10, which is closest to the experimental reference values and does not show outliers. PBE-rVV10L performs somewhat worse, with too large mean. Among the empirical London dispersion corrections, PBE+D4 and PBE+dDsC show the best performance. We note that most methods predict the interlayer distance of PtSe 2 bulk at significantly lower values compared to the experimental reference (see Table S16, Supporting Information).
For interlayer energies, again both vdW functionals show excellent performance with no deviation of the mean value, and with a scatter of less than 10 meV per atom. Among the empirical methods, dDsC and PBE+D3(0)+ATM are closest to the reference, MBD-NL underbinds on average by 20 meV per atom for the bulk systems. The other methods show significant overbinding, with 46 meV per atom for TS, and large scatterings. The general overestimation of the interlayer interaction energy of PBE+TS, D3(BJ)+ATM and PBE+D3(BJ) has already been reported by Otero-de-la-Roza et al. [45] It is interesting to note that the interaction energy overestimation increases with the mass of the chalcogen atoms of the materials for Grimme's corrections (see Tables S1-S15, Supporting Information).
The performance of higher-order terms in the Grimme D3 and D4 methods, i.e., the attractive C 8 term and the repulsive three-body ATM term, strongly depends on the damping function (see Figure S20, Supporting Information), where zero damping outperforms BJ. For short interatomic distances, damping functions reduce the LDI to zero or nearly constant values. For strongly compressed layered materials, one could expect that this effect limits the interlayer London dispersion attraction. However, this effect is compensated by the numerous attractive atomic interlayer contributions, and the London dispersion energy increases for shortened interlayer distances except for the dDsC method ( Figure S21, Supporting Information). The emerging interlayer Pauli repulsion becomes, however, dominant for such compressed systems.
The energy values of bilayer and bulk differ in mean and scatter by a factor of ≈2, but behave otherwise very similar. This is  easily understood as in the bulk there is one interstitial space per crystal layer, while the two layers of the bilayer system share one interstitial space.
One of the most convenient systems to exemplify the performance patterns of the interlayer distance and interlayer interaction energy per method is 2H-MoS 2 . These values, together with electronic characteristics that will be discussed below, are tabulated in Table 2. To understand the impact of interlayer interaction on the intralayer geometries, also the in-plane lattice constant a and especially transition metal-chalcogen (M-X) bond length d M−X is insignificant compared to other systems, especially with the ones with heavier chalcogen atoms.
In addition to the few other systems, as stated above, MoS 2 system gives us a clear picture to see the correlation between the interlayer distances d, hole effective masses m * h at direction [100] and electronic band gaps. While former can be seen on the Table 4 and Figure 2c, the latter can be seen in Figure 2b and Table S1 (Supporting Information). Table 3 shows the geometric and some electronic characteristics of bulk systems calculated with the SCAN-rVV10 functional,  the density functional of highest predictive power employed in this work. While SCAN-rVV10 predicts interlayer distances with absolute error less than 0.10 Å for the vast majority of the bulk systems, except for hBN and 1T-PtSe 2 , where errors are slightly larger (≈0.11 Å). hBN is also one of the few systems with 1T-TiS 2 that SCAN-rVV10 predicts its interlayer interaction energy with larger than absolute error of 10 meV per atom (≈14 meV per atom and ≈11 meV per atom, respectively). Geometric and electronic characteristics of bilayer systems optimized using SCAN-rVV10 are tabulated in Table 4. Opposed to their bulk counterparts loss of symmetry yields different M-X bond lengths in the bilayers, namely those pointing to the interstitial space (d in M−X ) and the others, pointing toward vacuum (d out M−X ). The difference between these two M-X bond types is largest for the semi-metallic 1T bilayers, where d out M−X about 0.1 Å longer than in the bulk. We further note a larger interlayer separation for bilayer systems compared to their bulk counterparts (see Tables S1-S18, Supporting Information).
It is well-known that the proximity effect has significant impact on the electronic structure of TMDCs and other layered materials. [48][49][50][51] To have a consistent picture of the impact of the structural parameters on the electronic properties, we calculated the band structure using the hybrid Heyd-Scuseria-Ernzerhof hybrid functional (HSE06) functional (see Methods for details) and obtained the characteristic values such as band gap, its character, and band curvatures at the valence band maximum (VBM) and conduction band minimum (CBM), expressed by the ef-fective electron (m * e ) and hole (m * h ) masses. As one would expect, we observe for most semiconducting bilayer and bulk systems that structures with larger interlayer distance show wider electronic band gaps (see Figure 2b for 2H-MoS 2 , 2H-MoTe 2 , and 2H-WSe 2 and Supporting Information for the remaining methods). However, for many 1T systems, especially for bilayers, the predicted electronic band gaps depend only slightly on the interlayer distance. For bulk HfSe 2 , TiS 2 and ZrSe 2 , and bilayer HfS 2 , HfSe 2 , ZrS 2 , ZrSe 2 , and TiS 2 , the different methods vary in interlayer distance by up to 0.25 Å, which is associated with a modest electronic band gap change of up to and 0.1 eV. However, contrary to the other semiconductor systems, interlayer distance and band gap are not in correlation. Detailed tabulated electronic band gaps are listed in Tables 3 and 4 for SCAN-rVV10 bulk and bilayer systems and in the Supporting Information for the remaining methods. The most striking impact that an inappropriately chosen structure optimization method might have is the difference in band gap nature: For many of the 2H semiconducting structures, there are various high-energy valence band pockets that are very close in energy, and the methods employed here do not agree on the location of the VBM for MoSe 2 bilayer and MoTe 2 bulk. The latter case is demonstrated in Figure 2a.
Even though the intralayer geometry is only slightly affected by the choice of London dispersion treatment, due to the nearly isoenergetic valleys at the valence band top the different functionals predict different electronic band gap characteristics for WSe 2 and WTe 2 monolayers (see Figures S6 and S7, Supporting Information). The band dispersions, on the other hand, differ only little and for the Γ→M path (or [100] direction) they depend linearly on the interlayer distance, as shown in Figure 2b, while they show even less impact on the band curvature for [120] and [110] lattice directions. Effective masses for semiconducting systems are tabulated in the Supporting Information (Tables S1-S6).

Conclusion
Interlayer distances and interlayer interaction energies are two essential parameters which are controlled by LDIs. We have shown that various correction schemes for LDI are capable to el-evate the PBE functional to predict interlayer distances and interlayer interaction energies with a small error margin. The highest predicted power is demonstrated by the vdW functionals PBE-rVV10L and SCAN-rVV10, whereas (semi)empirically corrected functionals show larger error margins, most notable for the interlayer interaction energy. Among the latter ones, we note the excellent performance of PBE+dDsC for geometries.
Proper choice of the LDI correction within the DFT framework becomes even more important for electronic band structure calculation due to their implications on electronic properties and possible experimental studies for the various electronic applications. We have shown that the electronic band gaps and electron-hole effective masses of the semiconducting materials are severely (up to 21% for E gap and 40% for m * h [100]) affected by interlayer distances. We conclude that it is possible to utilize DFT-based method for the prediction of geometries of vdW heterostructures. If computer time allows, SCAN-rVV10 is the method of choice, otherwise PBE+dDsC predicts geometries with sufficient accuracy to allow for further investigation of the electronic properties.

Experimental Section
Concept: In this study, interlayer distances were defined as the distance along the z-axis between two transition metals of the layers for TMDCs, carbon atoms for graphite/graphene, and boron atoms for hBN. The interlayer interaction energies, on the other hand, were determined using the following equation Total Energy − n * E Free−standing Layer Total Energy (1) where n is the number of the layers in the unit cell.
Reference Data: Reference interlayer distances, in-plane lattice constants, and M-X bond lengths of the bulk layered systems were taken from Springer Materials database [46] and Crystallography Open Database (COD). [47] For reference, interlayer interaction energies at the RPA level were calculated out using the Vienna Ab Initio Simulation Package (VASP) [52][53][54][55] code, as explained in the section below. To obtain reference interaction energies, the RPA total energies calculated on SCAN-rVV10 optimized bulk, bilayer, and monolayer geometries were chosen.
For the methods (a)-(f), the all-electron Fritz Haber Institute ab initio materials simulations (FHI-aims) [56][57][58][59] package was employed with numeric atom-centered orbitals (NAO). For the methods (g)-(i), as well as the RPA calculations, the Vienna Ab Initio Simulation Package (VASP) [52][53][54][55] was employed using projector-augmented waves (PAW). Grimme's corrections were made available through the Atomic Simulation Environment (ASE). [60] For relaxations, -centered k-grids of 16 16 6/1 (bulk/2D) for Graphite/Graphene, 9 9 4/1 for h-BN, 12 12 4/1 for 2H systems and 15 15 8/1 for 1T systems were chosen. In (a)-(d), an automatic adjustment of the Self Consistent Field (SCF) settings after three steps was set on tight tier one NAOs (default 2010). In (e)-(g), the basis set plane wave cutoff was set to 1.5 times the default given in the PAW library with Gaussian smearing of 0.05 eV and total energy convergence criterion of 10 −6 eV. In both cases, full structural relaxation including lattice vectors and atomic positions considering space-group symmetries were undertaken as implemented in FHI-vibes [61] (a-b) and the ASE (methods c-i) till forces and stresses were minimized below a threshold of 0.01 eV Å -1 . The vacuum space of monolayers and bilayers was fixed at 18 Å.
For the RPA Hartree Fock (E EXX ) energy calculations, denser -centered k-grids of 21 21 15/2 were used. For the RPA correlation energies (E C ) on the other hand, same k-grids of the relaxation calculations were used. PBE orbitals were used for the RPA calculations.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.