Hi, all,

I am working to add a new polynomial potential to LAMMPS, and I have been using the `pair_lj_cut.h`

and `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 `compute()`

and `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 `f_analytic`

.

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`

but instead

`f_analytic = (1/r)*forcelj, f_analytic = r*fpair`

.

So neither `forcelj`

nor `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.

Thank you.

Best,

Gabriel

(edited for formatting)