Energy conservation problem regrading REBO potential in LAMMPS

Hi all,
I'm using LAMMPS to simulate graphene nanoribbons recently. Since there
are some bugs in AIREBO potential reported by several groups (see
Pair AIREBO bug reported by several independent groups · Issue #59 · lammps/lammps · GitHub) and the LAMMPS developers
have fixed all the bugs after 19 May 2017. I'm using the latest version of
LAMMPS (5 Feb 2018). But when I'm testing the energy conservation for
graphene nanoribbons using the REBO potential, I found there is a global
drift in total energy.

The simulation details are as follows:
Configuration: graphene nanoribbon (GNR), the edge atoms are passivated
with Hydrogen atoms.
Potential: REBO potential.
Timestep: I choose two timesteps, 1 fs, and 0.5 fs.
The input file, the initial configuration, and the log files are attached.

Besides these tests, I also did simulations for monolayer graphene with
periodic boundary conditions and the energy drift is much smaller than that
for GNR. So I guess there are still some bugs in the REBO potential that
are not fixed. I hope my tests can provide some information.
I checked my input file carefully and I didn't find errors since it is
very simple. Can anybody else who cares about the energy conservation of
REBO potential check and run my input file on another computer to ensure
the energy is not conserved for this configuration or point out my mistakes?

​it doesn't look to me as if there is a bug. there will *always* be a
small(!) drift of the total energy​. this has multiple reasons. the two
most important ones are:
1) discretization of the numerical integration of the differential
equations (of motion). this error gets smaller with smaller time steps, and
gets larger the faster atoms move, i.e. the smaller their mass is.
2) discretization of numbers through using floating-point math. this error
gets larger the more floating point operations are required, i.e. the
smaller the time step is.

so energy conservation should improve with reducing the time step, however,
there is a limit where the error due to the use of floating point math
starts to matter. of course, there also is the "cost" in CPU time, that you
don't want to become too large. the optimum choice of time step depends on
the fastest motion, and how much you are willing to trade off energy
conservation against cost of computation.

in my personal experience, for a system containing hydrogen atoms and plain
NVE time integration without any thermostats, 1fs is far too large and even
0.5fs is pushing it. have you tried 0.25fs? that is what is required, e.g.
for flexible water molecules, to get good energy conservation. the fact,
that you see much improved energy conservation without the hydrogen atoms
and the attached plot comparing total energies for runs over 1ns with your
data for 1fs and 0.5fs and a run using a 0.25fs time step are confirming my
claim.

​axel.​