  Multilayer graphene with rhombohedral and Bernal stacking are supposed to be
metallic, as predicted by density functional theory calculations using
semi-local functionals. However recent angular resolved photoemission and
transport data have questioned this point of view. In particular, rhombohedral
flakes are suggested to be magnetic insulators. Bernal flakes composed of an
even number of layers are insulating, while those composed of an odd number of
layers are pseudogapped. Here, by systematically benchmarking with plane waves
codes, we develop very accurate all-electron Gaussian basis sets for graphene
multilayers. We find that, in agreement with our previous calculations,
rhombohedral stacked multilayer graphene are gapped for and magnetic. However,
the valence band curvature and the details of the electronic structure depend
crucially on the basis set. Only substantially extended basis sets are able to
correctly reproduce the effective mass of the valence band top at the K point,
while the popular POB-TZVP basis set leads to a severe overestimation. In the
case of Bernal stacking, we show that exact exchange gaps the flakes composed
by four layers and opens pseudogaps for N = 3, 6, 7, 8. However, the gap or
pseudogap size and its behaviour as a function of thickness are not compatible
with experimental data. Moreover, hybrid functionals lead to a metallic
solution for 5 layers and a magnetic ground state for 5, 6 and 8 layers.
Magnetism is very weak with practically no effect on the electronic structure
and the magnetic moments are mostly concentrated in the central layers. Our
hybrid functional calculations on trilayer Bernal graphene multilayers are in
excellent agreement with non-magnetic GW calculations. For thicker multilayers,
our calculations are a benchmark for manybody theoretical modeling of the low
energy electronic structure.
