There is a report that Lenard Jones potential does not capture the interactions between graphite layers properly [Phys. Rev. B 71, 235415 (2005)].
Since AIREBO is using LJ potential to capture the interaction between layers, is it appropriate to use this AIREBO potential to model graphite?
What is the best potential that we can use to model the graphite with LAMMPS?
The 2005 PRB paper you cited stated that a pure two-body LJ is not
sufficient in distinguishing the different stackings of graphite,
which I believe is valid. But AIREBO is not a pure LJ potential - LJ
is only a part of the formalism and it turns on/off based on different
criteria. If you take a look at the 2000 AIREBO paper (J. Chem. Phys.
112, 6472 (2000)), you will see that Stuart and co-authors discussed
the ability of AIREBO to distinguish different graphite stacking
sequences. Therefore, it should fine and appropriate to model
graphite with AIREBO.
As to your second question, you should ask "what is the best empirical
potential to use for modeling graphite?" And the answer should be in
the Literature. Once a best potential is identified, see if it is
supported by LAMMPS.