Factor of 0.5 in AIREBO force calculation

LAMMPS LAMMPS Mailing List Mirror
1

In the AIREBO potential that is currently in LAMMPS (or at least the version of LAMMPS that I am using from mid summer), there is a 0.5 in “pij forces”, “dwik forces” and “PIJ forces” calculations in the subroutine “bondorder”. It is not clear to me where this 0.5 contribution comes from. The -0.5 from the derivative of the bondorder term (db/dr) is included in the variable “tmp”. My student and I have both done the math by hand and do not get a 0.5 factor in these force calculations.

I know there has been some discussion on this mail list about a factor of 2 between the Brenner REBO paper and the Stuart AIREBO paper for some coefficients. Is this factor of 0.5 in certain force calculations an artifact of this previous discussion? Or is it due to double counting in the summation over neighbors k in some way.

Thank you in advance for your response,

-Doug

2

Not clear if you are asking this about REBO or AIREBO,
since both use the same lo-level code. I'm guessing
Steve Stuart (CCd) can answer this, as they have looked
at the LAMMPS AIREBO in detail.

Steve

3

Doug,

Very quick answer: this factor of 0.5 cancels an extra factor of 2 that appears elsewhere.

When the bond order is calcualated, it is carried as 2*bij. This extra factor of 2 is canceled by calculating only 0.5 of various terms that the bond order multiplies. This carries over into the force calculations as well.

(This is definitely confusing, and was carried over when the LAMMPS implementation was ported from a prior code. But it is correct, even if opaque and inconvenient.)

-Steve Stuart

4

Steve and Steve-

Thank you for the quick responses. Can you be a bit more clear regarding how (or where) the bond order term is calculated as 2*bij in the code?

We did some test calculations both in LAMMPS and by hand (as we were learning how to use the REBO and AIREBO implementations in LAMMPS) and got expected values for the bond order term and the energies for sample dimer, trimer and cubic structures in the material system we are studying. If the bond order is carried as 2*bij, where is it halved before it is used to compute energies in LAMMPS?

Thanks again,

-Doug