Questions about REBO/AIREBO implementation

LAMMPS LAMMPS Mailing List Mirror
1

Hi all,

I am writing an implementation of REBO from the Brenner 2002 paper, and am comparing against LAMMPS as a test. During the implementation, a number of questions have come up that may be able to be answered by people who better understand the history of REBO and AIREBO.

Fortunately, I’ve managed to find the answers to at least two of my concerns in advance (the factors of 1/2 in piRC: https://lammps.sandia.gov/threads/msg20308.html, and the changes in PCC(1, 1) and PCC(0, 2): https://github.com/lammps/lammps/pull/986). But many questions remain:

(1). REBO and AIREBO define different fitting parameters for the g splines, particularly for hydrogen atoms. AIREBO also adds terms to lambda_ijk (Stuart, A14). In both cases, LAMMPS uses the AIREBO values for both REBO and AIREBO. Is this done on purpose?

(2). In REBO, the coordination numbers used in the g, T, and piRC splines are NiX = sum of all weights from site i to atoms of type X (e.g. NiC = 3 in graphite). AIREBO instead uses NijX, which does not count the bond for the current interaction (e.g. NijC = 2 in graphite). Yet both papers use the same reference values for them just about everywhere (e.g. both define the transition range for the g spline as N = 3.2 to 3.7).

LAMMPS always uses NijX regardless of the potential. Given how the two formulations are irreconcilable, I assume this was simply an error in Brenner (2002)?

(3). There’s a bunch of madness surrounding FCC (a.k.a. piCC) even once the scale factor of 2 in Brenner is taken into account, including one thing that appears to be wrong in both papers as well as LAMMPS.

In Brenner (2002) Table 4, the following values appear for the FCC spline: (comments from the table)

FCC(2, 2, 2) 0.02200000 Benzene
FCC(2, 2, 3) 0.03970587 Average
FCC(2, 2, 4) 0.03308822 from
FCC(2, 2, 5) 0.02647058 difference
FCC(2, 2, 6) 0.01985293 F(2, 2, 2)
FCC(2, 2, 7) 0.01323529 to
FCC(2, 2, 8) 0.00661764 difference
FCC(2, 2, 9) 0.00000000 F(2, 2, 9)

dFCC(2, 2, 4–8)/dk = -0.006618

FCC(1, 2, 3) = -0.0200000000 C6H5
FCC(1, 2, 4) = -0.0233778774 Average from
FCC(1, 2, 5) = -0.0267557548 F(1,2,3) to F(1,2,6)
FCC(1, 2, 6–9) = -0.030133632 Graphite vacancy
dFCC(1, 2, 4–5)/dk = -0.020044

The AIREBO paper and LAMMPS both replicate these identically. (fixing the factor of 2, of course)

(3a): In the first group, the comment contradicts the values. Only elements from F(2,2,3)-F(2,2,9) follow a line equation; F(2,2,2) is not part of that line. Because the range of k used in dFCC/dk is consistent with a line from (2,2,3)-(2,2,9), I assume the the comment is wrong and that F(2,2,3) was obtained from a fitting species?

(3b): In the second group, F(1,2,3)-F(1,2,6) does follow a line equation, but the value written for dFCC/dk is not even close to the actual slope of the line. (-0.0033779 expected versus -0.020044 written). Has anyone brought this up before?

Thank you for your time,
Michael

2

Hi Michael,

Just a reminder that the AIREBO/REBO implementation in LAMMPS is the first generation REBO whereas the Brenner 2002 formulism is the 2nd generation REBO. Therefore many of the differences you noted below.

For 2nd generation REBO, I believe a number of faculty staff have a copy of the stand-alone code: Susan Sinnott, Judith Harrison, J. David Schall, Doug Spearout, etc so you don’t have to re-invent the wheel.

3

…this does not at all seem to be true? (nor would it account for any of the things I’ve pointed out)

The LAMMPS implementation has the Tij spline, and it switches the g spline based on coordination number like the second-generation rebo. Also, in the first link I posted, the author of the LAMMPS code himself says that the implementation is based off of the second-generation REBO code written by Brenner.

Michael

4

You are right, my mistake. Haven’t worked on AIREBO/REBO for too long.

Ray

5

Continuing with the implementation effort, I found a few more things:

(4) Brenner’s TCC covers T(2,2,1) and T(2,2,9), while Stuart’s TCC covers T(2,2,1) and T(2,2,2-9). As a result, 7 elements which are zero in Brenner’s paper are nonzero in Stuart’s. Since this isn’t brought up explicitly in Stuart’s paper, I imagine that the one in Brenner’s paper is an error…

(5) The circumstances surrounding lambda_ijk are much stranger than I thought.

Brenner says “lambda_HHC and lambda_HCH were fitted to remove spurious wells,” but I don’t see any value defined for them; all I can find is:

lambda_HHH = 4

Meanwhile, Stuart’s lambdas, if you work them out, are:

lambda_HHH = 4(r_ij - r_ik)
lambda_HHC = 4(r_ij - r_ik + rho_CH - rho_HH)
lambda_HCH = 4(r_ij - r_ik + rho_HH - rho_CH)

lambda_HCC = 4(r_ij - r_ik)

lambda_CXY = 0 for all X, Y

Not only are these functions of r_ij, but they also have units of length! However, lambda_ijk appears inside of an exponential, so it should be dimensionless. One could handwave this by saying that the 4 is supposed to represent a quantity of 4*(1/Angstrom), but that value seems too… convenient?

LAMMPS faithfully replicates Stuart’s lambdas. I’m not sure what lammps does internally if you tell it to use different units…but if it uses the new units throughout all computations, then there is no parameter you can adjust that will fix the scaling of lambda_ijk for use with other units of length.

6

Even more to round it off.

(6) I wasn’t sure until now, but there are indeed differences between Stuart’s and Brenner’s piRC.

An incomplete list:

  • Stuart’s piRC_CC has a block of elements at (0, 0<=y<=3, 3<=z<=9) (and by symmetry at (0<=x<=3, 0, 3<=z<=9)) that are assigned a value of 0.0049586079, which have a value of 0 in Brenner’s paper.
  • The differences in piRC_CH are huge! Of all the values printed for this function in table 9 (Brenner) and table IX (Stuart), only one of them actually matches, and that’s the value fitted to C6H6 in F(0,2,5-9). Some of the values differ by even an order of magnitude!

LAMMPS uses Stuart’s values in both REBO and AIREBO.

(7) It seems to me there is a bug in LAMMPS where it attempts to satisfy the F(i>3,j,k)=F(3,j,k) and F(i,j>3,k)=F(i,3,k) conditions by copying the values from the knots at coordinate 3 to the knots at coordinate 4:

// make top end of piCC flat instead of zero
i = 4;
for (j = 0; j < 4; j++){
for (k = 1; k < 11; k++){
piCCf[i][j][k] = piCCf[i-1][j][k];
}
}

But this is not fully accurate; it should copy the perpendicular derivatives piCCdfdy and piCCdfdz as well. Such perpendicular derivatives are similarly lost as values are copied from z=9 to z=10, but in that case it doesn’t matter because the z coordinate is clipped to 9 and so the curves for 9<z<10 are never generated.

…which leads to an even better fix: It seems to me that if the piRCSpline truncated the point to (3,3,9) rather than (4,4,9), then there would be no need to copy any of these knot values or solve for the curves at 3<x<4 and 3<y<4 in the first place, and bugs like this couldn’t happen! piCC, piCH, piHH ought to be double[3][3][9][64], and the only conceivable reason I can even imagine for them to be anything else is because LAMMPS still has code that dutifully reads these arrays from the potential. (and even then, this is only to be overwritten later… this is just plain confusing to anyone reading the code and shouldn’t be there!)

Michael