Energy conservation problem regrading REBO potential in LAMMPS

Hi Axel,
Thanks for your reply. I fully agree that the energy conservation cannot be exact and it can be reduced by decrease the timestep used in a simulation.

The reason why I guess there may exist some small bugs in REBO potential is quite simple: in my experience, if there is no bug in a force field, there should be no obvious global drift in the total energy. Increasing the timestep (not too much) will only increase the magnitude of the fluctuation of the total energy.

To prove this point, I use our own codes (our group has a home-made REBO potential) to simulate the same configuration with timestep equal to 1 fs, the total energy fluctuates without an obvious global drift. As the timestep decreases to 0.5 fs, the magnitude of the fluctuation also reduces. But none of them has an obvious global drift.
However, as you showed in the figure, there always exists an obvious global drift for all the timesteps used (1 fs, 0.5 fs, and 0.25 fs), although the magnitude of the drift decreases with the reducing the timestep.

Best,
Ouyang

Hi Axel,
Thanks for your reply. I fully agree that the energy conservation cannot
be exact and it can be reduced by decrease the timestep used in a
simulation.

The reason why I guess there may exist some small bugs in REBO potential
is quite simple: in my experience, if there is no bug in a force field,
there should be no obvious global drift in the total energy. Increasing the
timestep (not too much) will only increase the magnitude of the fluctuation
of the total energy.

To prove this point, I use our own codes (our group has a home-made REBO
potential) to simulate the same configuration with timestep equal to 1 fs,
the total energy fluctuates without an obvious global drift. As the
timestep decreases to 0.5 fs, the magnitude of the fluctuation also
reduces. But none of them has an obvious global drift.

​since you do have an alternate REBO implementation and i don't, you are in
a much better position than me to debug and identify where the any
unphysical energy is created in the LAMMPS implementation. ​

axel.

Hi Axel,
I’d like to identify the bug, but to be honest, I’m not very familiar with the codes, especially the codes of REBO potential in LAMMPS. So it will take a long time for me to find the bug. It will be great if someone who is familiar with the codes of REBO potential in LAMMPS can help to have a look. I will also read the codes of REBO potential in LAMMPS.

Today, I did another simulation for the nanoribbons by deleting all the Hydrogens on it. I found the energy is conserved when timestep is 1 fs. So it seems the bug is related to the Hydrogen. I hope this information is helpful.

BTW: it seems there is a problem to send the message to [email protected]…239…forge.net.

Best,
Ouyang

Hi Axel,
I'd like to identify the bug, but to be honest, I'm not very familiar with
the codes, especially the codes of REBO potential in LAMMPS. So it will
take a long time for me to find the bug. It will be great if someone who is
familiar with the codes of REBO potential in LAMMPS can help to have a
look. I will also read the codes of REBO potential in LAMMPS.

Today, I did another simulation for the nanoribbons by deleting all the
Hydrogens on it. I found the energy is conserved when timestep is 1 fs. So
it seems the bug is related to the Hydrogen. I hope this information is
helpful.

BTW: it seems there is a problem to send the message to
[email protected]

​there should not be any more.​ this week the mailing list was suspended
for several days while sourceforge was moving their servers to a new data
center.
most reasonable mail transfer agent programs should have tried again for
that period and no mail lost.

axel.

Hi Axel,
This is an update of regarding the energy conservation of REBO potential.
Since you fixed several bugs in the middle of 2017, I’m wondering if this is a newly introduced bug in REBO potential. To verify this, I download an older version of LAMMPS (10Aug2015) and use it to run the same input file with the same configuration. I found there exist similar global drift also in this version of LAMMPS. I compared the source codes of REBO potential between 10Aug2015 and 5Feb2018, I saw you made a lot of changes. However, this bug seems to be an old one. I’m surprised that people never noticed it before.

As you suggested, I use a smaller timestep (0.25 fs) for NVE simulation and I found the energy start to deviate at ~700 ps. I ran the simulation for 2 ns and the global drift of total energy is around 0.1 eV. The initial total energy is around -1306.25 eV and the fluctuation of potential energy and kinetic energy is around 1 eV. Considering there are only 204 atoms in the nanoribbons, the energy drift is not small.

In my experience, if there is no bug in the force field, the fluctuation of total energy should be at least two orders of magnitude smaller than that of potential energy and kinetic energy. When I delete all the Hydrogen atom in the nanoribbons, this condition is satisfied even for larger timestep (1 fs).
Therefore, I think the bugs of REBO potential only exist in the C-H interaction or H-H interaction.

I still didn’t test the LJ term and torsion term in AIREBO potential, I will update to you when I have these results.

Best,
Ouyang

Hi Ouyang,
I’m interesting about the problem you proposed for REBO potential since it’s important to all the researchers who use this potential in LAMMPS. So I test your input file with your configuration in the latest version of LAMMPS 5Feb2018. Unfortunately, I found my results are similar to yours. Again, I use a relative small timestep (0.25 fs).

Hi Axel and Steve,
Did you check Ouyang’s input file before? Do you think there indeed exist an old bug for C-H or H-H interactions in REBO potential?
I’m not sure how big is the influence of this bug on the simulation results, but I think it’s worth to fix this bug if you can confirm that Ouyang is correct since REBO potential is frequently used.

Best,
Huang

Hi Ouyang,
I'm interesting about the problem you proposed for REBO potential since
it's important to all the researchers who use this potential in LAMMPS. So
I test your input file with your configuration in the latest version of
LAMMPS 5Feb2018. Unfortunately, I found my results are similar to yours.
Again, I use a relative small timestep (0.25 fs).

Hi Axel and Steve,
Did you check Ouyang's input file before? Do you think there indeed exist
an old bug for C-H or H-H interactions in REBO potential?
I'm not sure how big is the influence of this bug on the simulation
results, but I think it's worth to fix this bug if you can confirm that
Ouyang is correct since REBO potential is frequently used.

​identifying a problem and being able (or having the time and means) to
fix​ it are two different things. while it may be desirable to do
something, that does not automatically mean, that you are entitled to have
this addressed by the core LAMMPS developers.

please keep in mind, that most styles in LAMMPS are contributed code (so is
AIREBO/REBO) and that most of the recent corrections and bugfixes were
contributed by multiple people, that also invested the time to meticulously
check and review the code and its results. people like steve and me may
identify and resolve (obvious) syntax errors and correct programming style
problems or performance issues, but we cannot correct issues with the
science. we simply don't have the time and even if we had - given the
nature and size of the LAMMPS project - our time and effort is much better
invested in overall maintenance and improvement of the package, in
integrating pending contributions, figuring out better/faster ways to run
LAMMPS calculations and so on. in short, stuff that *many* LAMMPS users
benefit from, not just a subsection of them.

thus fixing issues like this one requires people that have a vested
personal interest (like you) and the means to compare the results computed
by LAMMPS step-by-step and across each piece of the force field against a
suitable reference code (which neither steve nor i have).

axel.

Hi Axel,
I see. Thanks for your detailed explanation.
I would like to invest time to identify the bug, unfortunately, I also don’t have a reference code.

Hi Ouyang
Since you have a home-made code at hand, could you take some time to identify the bug?
It will be great if anyone has a reference code (like the people discussing the issue in https://github.com/lammps/lammps/issues/59 ) can help to have a look at the code since they have more experience.

Best,
Huang

Hi Huang,
Thanks for your tests!

I agree with you that it will be much faster if someone who already has a lot of experience can have a look at the code since the possible bugs should be only related to the interactions with Hydrogen. However, everyone is busy with their own stuff, they may not have time or enough motivation to fix this small bug.

I really want to identify this bug and fix it. Actually I’m reading the codes of AIREBO potential in LAMMPS now. But to be honest, I’m not familiar with the codes and REBO potential is a quite complicated many-body potential. So it will take a long time for me to find this bug and fix it. But I will continue and I will update if I succeed.

Best,
Ouyang

How do you know that there's a bug in the code?

Your energy isn't conserved when using hydrogens on a nanoribbon. The AIREBO parameters that you're using may not be transferable to that scenario.

Drew Rohskopf
Atomistic Simulation & Energy Group
Massachusetts Institute of Technology
(404)403-0313

In my experience, if the energy and force are calculated consistently (i.e., no bug), the energy conservation is independent of the values of the parameters used in the force field. We can understand this point by taking LJ potential as an example. People use different parameters (epsilon and sigma) for different systems, although the parameters, as you said, may not be transferable to some scenario, it will only influence the “correctness” of the energy and force, but the total energy should be conserved in NVE simulations.

Best,
Ouyang

That makes sense, but in this scenario you have hydrogens (very small mass) and bad parameters can cause unstable vibrations with extremely high frequencies due to this small mass. If they vibrate too fast, the time step may not be sufficient to resolve their motion. I've seen this happen before; it may be worth testing by increasing the hydrogen mass to see if energy is still drifting.

Drew Rohskopf
Atomistic Simulation & Energy Group
Massachusetts Institute of Technology
(404)403-0313

Thanks for your suggestion.

I did NVE simulation with the same configuration by setting the mass of
Hydrogen equal to that of Carbon. I run the simulation for 500 ps and the
timestep is 1 fs. In this case, the energy drift is reduced but still
there (the absolute value of total energy is ~1238 eV and the energy drift
is ~0.15 eV, the fluctuation of the kinetic energy is ~ 3 eV). Considering
there are only 204 atoms in the nanoribbons, the energy drift is not small.

I also did NVE simulation with the same timestep by deleting all the
hydrogens, and in this case, the total energy is conserved (there is no
global drift and the fluctuation of the total energy is two orders of
magnitude smaller than that of potential energy and kinetic energy).

Thus, I think there may exist some discrepancies between energy and force
when calculating the interaction with Hydrogens.

Best,
Ouyang

That makes sense, but in this scenario you have hydrogens (very small
mass) and bad parameters can cause unstable vibrations with extremely high
frequencies due to this small mass. If they vibrate too fast, the time step
may not be sufficient to resolve their motion. I've seen this happen
before; it may be worth testing by increasing the hydrogen mass to see if
energy is still drifting.

​drew,

when setting up the same system (that is including regular mass hydrogens)
with CHARMM parameters and topology, there is much better energy
conservation. in all cases.

this is the most nasty kind of problem, when the result is only a little
bit off. please keep in mind, that energy conservation is a higher order
effect.

axel. ​

A better way to address this issue is to contact the person who
solved the previous issues with AIREBO in LAMMPS.
I don’t recall his name and am on travel, but Axel will know
who it was (from Aachen as I recall). He had a version
of the KIM AIREBO (from the S Stuart group) and could run careful tests
comparing both codes on simple inputs to high precision. He found a couple
code differences and fixed them in LAMMPS. After staring
at the 2 codes, my recollections is he thought there were no other differences.

If he runs his test harness on your input, he might
find another issue, or he might find that both codes
give the same answer.

Steve