Precise structure and energy of group 6 transition metal dichalcogenide homo- and heterobilayers in high-symmetry configurations

Two-dimensional group 6 transition metal dichalcogenide (2D TMDC) bilayers show various high-symmetry stacking configurations, which have also been observed in extended domains formed in their twisted homo- and heterobilayers. The interlayer energy varies for these stacking configurations, and the energy differences determine the relative size of the stacking domains. Therefore, the precise prediction of the composition- and stacking-dependent interlayer energy is crucial to model the domain structure of 2D TMDCs in their twisted bilayer homo- and heterostructures. For the validation of approximate methods that are necessary to tackle these systems encompassing thousands of atoms precise reference data is still lacking. Here, we employ the Random Phase Approximation (RPA) on previously validated SCAN-rVV10 geometries to obtain interaction energies of state-of-the-art accuracy on the six high-symmetry stacking configurations (Hhh, HhM, HhX, Rhh, RhM, and RhX) of MX2 (M=Mo, W, X=S, Se) bilayers and compare them with dispersion-corrected DFT functionals PBE-D3(BJ), PBE-rVV10L, and SCAN-rVV10. We identify SCAN-rVV10 as most reliable DFT variant with an average deviation of 1.2 meV/atom in relative energies from the RPA reference, and a root mean squared error (RMSE) of less than 2 meV/atom for interlayer interaction energies. We find interlayer distances obtained by PBE-D3(BJ) too short, with an impact on the electronic structure, which results in the incorrect prediction of the band gap character in some cases. A further result of this work is the significant lowering of interlayer energy and increasing of interlayer distance in the high-energy stacking configurations. These stackings can be accessible via shear strain and promote exfoliation.