Why is the band gap of graphene "opened" in this VASP calculation?

Question

I used VASP to do the band structure calculation on Graphene. Originally I would expect that the Dirac cone can be observed in the K point. If we watched qualitatively, we can see this behavior. But when I zoom in the region where the conduction band touches with the valences band, a bandgap can be observed as shown in the figure. That's weird because I've already adjusted my K point meshing to be an odd number and increase the number of meshing, it doesn't help. The upper band was still lifted.

Answer 1

The reason why graphene has Dirac points at K and $-$K is because of a combination of time reversal and inversion symmetries. Therefore, if you impose these symmetries in your VASP calculation, there is no reason why there should be a gap at K. Even if you do not impose these symmetries explicitly, you should be able to get the degeneracy to very high accuracy (better than what you are showing) if you do a well converged calculation in terms of the number of $\mathbf{k}$-points that you are using to sample the Brillouin zone.

So what could be going wrong? From the partial information that you provide, my best guess is that: (i) you don't include the symmetries in the calculation, (ii) the $\mathbf{k}$-point grid you are using for the self-consistent part of the calculation is not fine enough, and in particular it does not explicitly include the K point, and (iii) when you then do the non-self-consistent calculation to obtain the band structure along the high-symmetry lines, then you have very few points along the lines and the code simply naively interpolates between them to generate the plot.

So what would I suggest you do? The first thing is to impose the relevant symmetries. The second is to make sure that the $\mathbf{k}$-point you are using in the self-consistent calculations is large enough, and in particular that it explicitly includes the K point, and (iii) that in the non-self-consistent calculation you include a large number of $\mathbf{k}$-points along the high symmetry points near K.

All of this is based on my interpretation of the partial information you provide.

The discussion above is somewhat different if you included spin-orbit coupling. In that case, a gap develops at the K point in graphene, an observation that led to the field of topological insulators. However, I don't think that this explains your results, because spin-orbit coupling in graphene is very small, on the $\mu$eV scale, orders of magnitude smaller than what you are seeing. For this reason spin-orbit coupling can be neglected when doing calculations for graphene, so I imagine that this is what you are doing.