AIREBO: Discrepancies in its Lennard-Jones force calculation component

Dear all,

I've recently done a comparison of the OpenKIM AIREBO model and the LAMMPS AIREBO code and found a number of discrepancies in the force calculation of the AIREBO LJ interaction.
In the past, there have been a number of reports of issues with AIREBO by a number of groups [1,2,3], which this discussion is probably related to.

The code from OpenKIM seems to be a version of the serial code developed by the Stuart group.
In all cases I believe that these discrepancies are indeed bugs in LAMMPS' AIREBO implementation (analytically taking the derivatives of the potential to see which variant should be right).
I've collected all the fixes that are necessary in the AIREBO code in LAMMPS in a pull request [4].
For details, see the individual commit messages---each commit addresses exactly one issue.

To list the issues:
1. In bondorderLJ(), Etmp is calculated twice, so it's value is two time the expected value.
2. In bondorderLJ(), the derivatives of the cosines in p_ij do not take into account that rij is scaled.
3. In bondorderLJ(), the sum omega term does not take into account that rij is scaled.
4. In FLJ(), the LJ term is modified, and the formula contributing to the derivative of the modified LJ term mixes up format and signs.

In addition, I also found a discrepancy in the pi^rc_CC spline for Nconj > 8.
However, I am not sure about this one; the original papers can IMO be interpreted such that any of the splines is correct.
This discrepancy also probably has a much lower impact on force accuracy, so it is not as important.

It would be great if someone else could try my patch to see if there are any issues, and possibly run it against other AIREBO implementations to verify its efficacy.

Thanks a lot,

PS. I tried to conform to the Github workflow that you have, but please forgive me if I overlooked something. I am more than happy to get pointers how to improve this in the future.

PPS. Steps to reproduce my observations:
1. Grab the current LAMMPS version
2. Enable Manybody and OpenKIM
3. Install OpenKIM API with AIREBO model included [5]
4. Build LAMMPS, once with patch, once without
5. Run provided input files [6] with lammps w/o patch
6. vimdiff dump.lammps
     -> see the deviations, especially in timestep 8-10
7. Run provided input files with lammps w/ patch
8. vimdiff dump.lammps
      -> note that the deviations are noticeably reduced

The input file is a reduced (to 50 atoms/10 timesteps) example from a bigger verification suite.
The scripts rerun a trajectory that exercises the offending code regions, so each trial really is based on the same atom positions.
For your convenience, I've also included the resulting dump files.

5: You might (depending on your system) have to add "perror.f" to the source line in the src/models/*AIREBO*/Makefile

Thanks Markus - this sounds very helpful. We’ve never had

a copy of source code for any other AIREBO implementation to

compare to. One recent issue that was flagged was for

co-linear atoms (3 I believe) in the LJ dihedral calculation,

which can cause the angle term to blow up (reportedly).

Steve Stuart (Clemson) said in their AIREBO code they

taper the interaction as the angle -> 0.0 to avoid this issue.

Do you see anything similar in the KIM version? It doesn’t

seem related to anything in your list of changes …


Dear Steve,

The issue of collinear atoms likely is not related to the changes I'm suggesting.
However, they might (this is speculative though) be part of the reason why some people's simulations even reach states where collinearity becomes an issue.

Regardless: The code to switch off the dihedral term for collinear atoms exists both in the KIM code and the LAMMPS code.
Very broadly speaking, the two codes seem to be related in at least some capacity (many variable names have the same names/code is structured similarly).

The part of the code that is supposed to switch off the dihedral for collinear atoms in LAMMPS is the part that looks like
     tspijl = Sp2(costh, thmin, thmax, &dtsijl) or
     tspjik = Sp2(costh, thmin, thmax, &dtsjik)
Later, we multiply by the dihedral term with (1-tspijl) * (1-tspijl).
Because thmin = -1, thmax = -0.995, once the cosine is at -0.995 or lower (which corresponds to an angle between 174 deg and 196 deg), the interaction thus should switch off.
When I started looking at the code this was quite confusing, since there is no mention of such a thing in the paper.

Now there is one change in the patch set (namely de4af6f) that might have an impact here: The cosine calculation used for the switching is boguous in bondorderLJ(), see for example line 2699.
Note that in the current implementation of that function, rijmag != r32mag, so calculating costh = 1/2(r23mag^2+rikmag^2-rjkmag^2)/rijmag/rikmag can't possibly result in something useful---see how in the numerator we use r23mag, but in the denominator rijmag?
Whether this was the cause of that issue or not is hard to tell, but it is possible.

It might be wise to not skip iterations based on the sine of the angle and some tolerance, but rather directly check whether tspijl == 1.0 etc.
This might even be cheaper/cost neutral, since computing the sine costs us a sqrt, whereas the cutoff is a polynomial---not sure here though.

Best, Markus

Released these changes/bug-fixes in a 7Mar17 patch today.

Thanks Markus for the careful analysis!