AIREBO energy conservation during an NVE MD simulation

Hello LAMMPS community,

I am currently simulating hydrogen adsorption on graphene and have run
into an energy conservation issue with the AIREBO potential (I'm using
the 15Aug14 version). Decreasing the simulation time-step (all the way
down to 0.0001 fs) had no impact, which makes me think the interatomic
potential is causing the problem.

I attached a simplistic input script that reproduces the issue in
which a single hydrogen atom is fired at a graphene substrate and
chemically adsorbs onto the surface. When I run the same simulation
with the REBO potential, which is AIREBO with the LJ and torsion terms
turned off, the energy discrepancy is much less.

Looking at the trajectory of the runs and the total energy vs time, it
seems the loss of energy conservation occurs right when the REBO C-H
"bond" formation occurs. For the AIREBO case, the energy increases by
over 0.12 eV. Comparing the REBO vs AIREBO runs, perhaps the LJ term
is not switching off smoothly. I know a few years ago the AIREBO
potential was updated by the Stuart group to fix a lot of issues to
bring the LAMMPS version in agreement with their orginal code, maybe a
bug was missed.

Does anyone have any advice or ideas how I can pinpoint and resolve
this lack of energy conservation in my simulations with the AIREBO
potential?

Thanks to everyone for your time,

-CBL

MD_run.in (458 Bytes)

energy_plot.png

graphene_opt (53 KB)

Hello LAMMPS community,

I am currently simulating hydrogen adsorption on graphene and have run
into an energy conservation issue with the AIREBO potential (I'm using
the 15Aug14 version). Decreasing the simulation time-step (all the way
down to 0.0001 fs) had no impact, which makes me think the interatomic
potential is causing the problem.

I attached a simplistic input script that reproduces the issue in
which a single hydrogen atom is fired at a graphene substrate and
chemically adsorbs onto the surface. When I run the same simulation
with the REBO potential, which is AIREBO with the LJ and torsion terms
turned off, the energy discrepancy is much less.

Looking at the trajectory of the runs and the total energy vs time, it
seems the loss of energy conservation occurs right when the REBO C-H
"bond" formation occurs. For the AIREBO case, the energy increases by
over 0.12 eV. Comparing the REBO vs AIREBO runs, perhaps the LJ term
is not switching off smoothly. I know a few years ago the AIREBO
potential was updated by the Stuart group to fix a lot of issues to
bring the LAMMPS version in agreement with their orginal code, maybe a
bug was missed.

Does anyone have any advice or ideas how I can pinpoint and resolve
this lack of energy conservation in my simulations with the AIREBO
potential?

why should there be energy conservation, when a chemical reaction happens?

I was under the impression that the total energy should be conserved
for a chemical reaction when using a classical conservative potential
(AIREBO) in the NVE ensemble, as long as the switching functions
ensure there are no discontinuities in the energy and forces. So the
total energy remains constant, it just transforms between KE and PE.

From your question, I take it that this is not the case here?

-CBL

I was under the impression that the total energy should be conserved
for a chemical reaction when using a classical conservative potential
(AIREBO) in the NVE ensemble, as long as the switching functions
ensure there are no discontinuities in the energy and forces. So the
total energy remains constant, it just transforms between KE and PE.

From your question, I take it that this is not the case here?

I don't know. you would have to contact the folks that maintain the code.

it was just a blind guess, since you transition from a pair-wise
additive interaction (LJ) to a manybody (angles, dihedrals)
interaction.

axel.

As far as I’m aware there are no known discrepancies
between the LAMMPS version of AIREBO and that
of the Stuart group.

Steve