If this line fixes the problem, the issue is caused by too infrequent neighbour list rebuiliding. To solve the problem you can increase frequency of updating the neighbour list (neigh_modify command), increase skin size (neighbor), decrease timestep, or use a fix dt/reset.
You are using a quite low value for a Tdamp parameter. This can be a root of the problem too.
Your description of a problem with densities predicted by a LAMMPS implementation of AIREBO is too vague to comment. Please provide the numerical values of densities and an input script (along with data files with simulated systems) used for the calculations.
We tried your suggestion (neigh_modify every 1 delay 0 check no) and it
did not entirely solve the problem. The simulation still crashes (although
at a 10 K higher energy and the overall potential energy progressively
goes to lower values than it was in the former simulation).
Note that we already had tried this comment without success when trying to
fix other crashes.
The timestep used (0.2 fs) should not cause any problem as my advisor (JM
Leyssale) has performed MD simulations on similar systems up to very high
temperatures for many years using his own serial code without experiencing
any crashes. By the way, in the litterature, timesteps for MD simulations
with AIREBO are always between 0.2 and 0.5 fs. Actually, a timestep of 0.2
fs is normally even suitable to perform Born Oppenheimer MD simulations
with full DFT energy for such systems, so it should not be a problem.
Regarding the damping parameter of the thermostat, we agree that it is
low, but it is almost mandatory for tempertature ramp simulations with
finite ramp rate. Again, there are hundreds of papers reporting on
simulations achieved with similar parameters.
Also, we noticed a recent post from Carla de Tomas Andres, reporting on
similar issues as those we observed. Furthermore, we actually experienced
the exact same energy/temperature spikes as those they report on when
quenching a C/H liquid with low H content. We report here on an extreme
case, when we have high H content, for which we can even not stabilize the
system at room temperature without crash. This is not normal with AIREBO!
We really believe that there is a problem here with the AIREBO
implementation in LAMMPS as these anomalous behaviors never occur with our
serial code.
Since the problem we get on the equations of state of alkanes is probably
linked to the one reported here we will not expand on it right now but
rather focus on solving the main issue.
This is not normal with AIREBO!
We really believe that there is a problem here with the AIREBO
implementation in LAMMPS as these anomalous behaviors never occur with our
serial code.
It sound like there is a problem, as you’ve checked all
the other things people have mentioned. I’m CCng Steve
Stuart, as his group looked into the accuracy of the LAMMPS
AIREBO a few years ago and carefully checked things. They
didn’t find any remaining issues (some were fixed at that time).
Maybe he will have additional ideas beyond these.
It sounds like you have 3 implementations: LAMMPS, yours, the
Stuart group code.
a) Are you sure these are all running the same version of AIREBO?
The LAMMPS version is an older AIREBO from an early paper (2000,
see the pair airebo doc page).
Could the issue be that yours or the Stuart group are using an
improved AIREBO? Ditto
for the AIREBO param files that the various codes are using?
b) Can you reproduce the problem with a small system and a short
simulation? One way to make it short is to write a restart file
right before things blow up. Then see if you start from the restart
file if things still quickly blow up. This is important to be able
to do, b/c there is unlikely that two implementation will track each other
for 50K steps (your output), even w/out a bug, due to numerical
round-off. So that makes it hard to debug.
c) At the time that LAMMPS blows up, can you write a snapshot,
e.g. step 55049 in your output above. Then read it back in and
run a 0 step simulation that just compute energy/pressure and
dumps the forces. Do the same with yours and Stuart’s code.
If there are significant differences (ideally the LAMMPS one
will produce the NaNs you see on step 55050), and the
others do not, then this should be easy to debug. If you
can do something similar for a small system (as in (b)) then
even better.
Thanks for posting careful messages about the issue,
Steve