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.