  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 interactions 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.
