I started using AIREBO and ran into the same problems as previous users have had with it crashing.
I’ve been able to trace the bug in pair_airebo.cpp to the Lennard-Jones calculation of forces that occurs in the bondorderLJ() function, where just before the kinetic energy first spikes, I see the “f2” and “f3” variables going very high. This occurs because “sin321” is calculated to be 0.0007, which is large enough to get through the “if” statement in line 2682 (since it’s greater than TOL=1.0E-9) of the November 17 2016 version of the source code (attached), but small enough when squared to mess up the forces. Richard and Axel also reported this issue on Github, but I guess their fix of putting in that tolerance just wasn’t strict enough. When I change TOL to 1.0E-3, the program runs without crashing.
As further evidence that it’s the Lennard-Jones calculation that contains the bug, when I set the LJ_flag of AIREBO to 0, I also found that the program runs without crashing. Further, the AIREBO/Morse style also runs without crashing.
I’m still concerned about 2 things.
1: Other bugs might be present in the code since Steve Stuart mentioned that there are “several bugs in the force calculation.” If this is just one of several bugs, and catching this one allows the program to run, but it runs incorrectly, it could make matters worse since it leads to false confidence in the code. Steve, do you happen to know in which parts of the code the known errors you mentioned are supposed to be lurking?
2: When TOL is set so low, perhaps calculations are skipped more regularly than they should be. Perhaps a separate tolerance should be set for line 2682 than the other spots? Is it okay for the “if” loop that line 2682 leads into to be skipped every once in a while?
Any thoughts? If other AIREBO users that were getting this error can check whether this fixes their simulations too, that could be helpful.
pair_airebo.cpp (137 KB)