The structures obtained from ReaxFF calculations were verified by comparing the highsymmetry bilayer stackings to structures fully optimized using FHI-vibes and the density functional theory (DFT) with the Perdew-Burke-Ernzerhof (PBE) functional augmented with the non-local many-body dispersion (MBD) corrections . The intermediate basis set of tier 1 and a converged 12 × 12 × 1 k-grid were used, and a maximum force of 0.02 eV Å -1 was set as the threshold.