Please find attached bandstructure calculated with atk that match the paper.
The detail missing in your calculation was the conduction band starting 3.5 eV above the gamma point, while the occupied bands are well described. The work function of graphene is 4.5 eV, thus, at 3.5 eV the electron is almost not bound and is mainly in the vacuum region. Thus, to get a description of this band you need to describe the electron wavefunction in vacuum.
Using planewaves the wavefunctions are everywhere, and vacuum bands are automatically described by planewaves. In a LCAO method as ATK, the basis functions are located at the atom positions, this makes the methods more efficient than plane-waves, however, they only describe accurately the electron bands where the electron is moving in the vicinity of the atoms (which is the case for the occupied bands, but not for bands high above the fermi level).
To describe such bands with a LCAO method you need to put basis functions in the vacuum, socalled ghost atoms, and this is what is done in the attached calculation. (in the attached calculations I just put a lot to make sure everything is converged, and you probably only need a fraction of the basis functions I use).
As for the reference paper, eventhough the rest of the calculations they present are done with a LCAO method, the band structure calculation of graphene (which is mainly an illustration) must be done in a different way as the rest of the calculations, I guess they obtained that using a plane-wave method. I am contacting the authors to clarify this, and I will let you know once I get their reply.
However, it is important to note that these vacuum states are of no importance for calculating most properties of graphene, and the LCAO basis set of graphene without the ghost atoms can describe very accurately most graphene properties.
