questions on angles_class2.cpp


I’d like to write a new angle style that has some similarities to the class2 angle style. I’ve been looking at the code from “angles_class2.cpp” in order to get an idea of how I should write the new angle style. I have a few questions about how the forces are tabulated for the bond angle term E_ba. Specifically, I don’t understand why the term ‘aa11’ is defined as:

aa11 = aa1 * c/rsq1

and not

aa11 = aa1*(1/r1)

Likewise, I don’t understand why ‘aa12’ is defined as:

aa12 = -aa1 / (r1 * r2)

Actually, I don’t even understand what the ‘aa12’ is doing there at all- I would have thought that the ‘aa11’ term would be sufficient to differentiate the theta term in E_ba w.r.t. x . It looks as if there’s an issue related to dealing with the specific symmetry of this particular potential that I’ve neglected. Anyone out there who has experience with the code for this angle style want to answer this? Apologies if this question has already been answered on the mail list, I searched but didn’t find anything.


PS This code comes from version LAMMPS-28May12

It's been a very long time since I looked at that code.
You can try asking Paul Saxe - psaxe at