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 ( H h h , H h M , H h X , R h h , R h M , and R h X ) of MX₂ (M = Mo, W; X = S, Se) bilayers and compare them with the dispersion-corrected density-functional theory (DFT) functionals Perdew-Burke-Ernzerhof (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 of less than 2 meV/atom for interlayer interaction energies. We find interlayer distances obtained by PBE+D3(BJ) as being too short, with an impact on the electronic structure, resulting in the incorrect prediction of the band gap character in some cases. A further result of this work is the significant lowering of the interlayer energy and increasing of the interlayer distance in the high-energy stacking configurations. These stackings can be accessible via shear strain and promote exfoliation.