I am working to add a new polynomial potential to LAMMPS, and I have been using the
pair_lj_cut.cpp as templates/references. For the most part everything has gone smoothly or been pretty clear (though I haven’t actually compiled yet).
However, when I was looking at the
single() methods of
PairLJCut I glanced at the force calculation lines, but did not quite see what I expected. My expectation was the
forcelj variable would be equal to the negative first derivative of the Lennard-Jones potential, let’s call this quantity
What I find (I have checked this by hand an in a Mathematica script, though I could still be wrong) is:
f_analytic != forcelj, f_analytic != fpair
f_analytic = (1/r)*forcelj, f_analytic = r*fpair.
fpair are the same as the analytical result. Now I do see that the
fpair is equal to
factor_lj * forcelj* r2inv (
r2inv=1/(r^2)), so my suspicion is that the combination of these factors and multiplication by the displacement vector between the two atoms gives a vector whose magnitude is equal to the analytical result. (Which would imply
factor_lj=1 and the vector which the
fpair scalar multiplies is of length
r, which all seems reasonable.)
My first question is “have I understood this correctly?” (I ignored the inner and outer regions/radii in the discussion above, but I think that’s safe in this case.)
If the answer to that is “yes”, I would like to follow up with “could these facts be made more clear?” If I have understood the process correctly, I suspect the function is written the way it is to minimize the number of floating point operations (in this case you’d be saving the process of normalizing). However, I can imagine the case where somebody (like me) who wants make a new potential sees the words “forcelj” and/or “fpair” and just writes in whatever the analytical expression is (thus giving them incorrect results). This error seems far too easy to make, mainly due to the naming of the variables and lack of comments/documentation. Perhaps there is a part of the documentation somewhere that describes this naunce, but I could not find one.
I am also still a bit confused about what exactly
factor_lj is, why it shows up in non-L-J potentials (Yukawa for example), and whether it is changed on a per-potential basis (and where would this change be made)? These are secondary to (but perhaps somewhat entangled with) my main questions.
(edited for formatting)