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.