Similarly to the D3 correction, the PBE+MBD and PBE+D3(BJ) models yield an unphysical minimum at 3.6 Å. We also expect in this case that the AIMD trajectory will not remain stable. Conversely, PBE+TS, PBE+D2 and SCAN+rVV10 predict the AA' structure as clear global minimum. The minimum computed with PBE+TS is slightly more pronounced. Hence, we selected the PBE+TS model for our AIMD simulations. As expected, the AA'-type stacking between two adjacent layers was preserved during the AIMD simulations with the PBE+TS model.