The stacking orders in layered hexagonal boron nitride bulk and bilayers are studied using high-level ab initio theory (local second-order Moller-Plesset perturbation theory, LMP2). Our results show that both electrostatic and London dispersion inter
actions are responsible for interlayer distance and stacking order, with AA being the most stable one. The minimum energy sliding path includes only the AA high-symmetry stacking, and the energy barrier is 3.4 meV per atom for the bilayer. State-of-the-art Density-functionals with and without London dispersion correction fail to correctly describe the interlayer energies with the exception of PBEsol that agrees very well with our LMP2 results and experiment.
