  There has been a lot of excitement around the observation of
superconductivity in twisted bilayer graphene, associated to flat bands close
to the Fermi level. Such correlated electronic states also occur in multilayer
rhombohedral stacked graphene (RG), which has been receiving increasing
attention in the last years. In both natural and artificial samples however,
multilayer stacked Bernal graphene (BG) occurs more frequently, making it
desirable to determine what is their relative stability and under which
conditions RG might be favored. Here, we study the energetics of BG and RG in
bulk and also multilayer stacked graphene using first-principles calculations.
It is shown that the electronic temperature, not accounted for in previous
studies, plays a crucial role in determining which phase is preferred. We also
show that the low energy states at room temperature consist of BG, RG and mixed
BG-RG systems with a particular type of interface. Energies of all stacking
sequences (SSs) are calculated for N = 12 layers, and an Ising model is used to
fit them, which can be used for larger N as well. In this way, the ordering of
low energy SSs can be determined and analyzed in terms of a few parameters. Our
work clarifies inconsistent results in the literature, and sets the basis to
studying the effect of external factors on the stability of multilayer graphene
systems in first principles calculations.
