Pair_lj_cut forces

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)

Your observation is correct. The whole point of what is done when applying the force in the compute function is to avoid 1 division, and (more importantly) 1 square root.

To explain: while you compute the force between two particles as a “radial force” (i.e. based on the distance), you have to then “project” that force on the x-, y-, and z-component of the positions of the distance vector. The main reason that people use a 12-6 Lennard-Jones potential (and not a 9-6 one) is that you can formulate it such that you don’t need to use the square root operation (starting from not comparing “r” against “rcut”, but rather “r2" against "rcut2”. Similarly for projecting forces on the distance vector you need to compute delx/r, dely/r and delz/r and multiply it with the radial force which would have to be in terms of 1/r**5 and 1/r**11 (due to rules for derivatives of polynomials). Thus what is done is to “move” the 1/r term from the application of the force to the force vector already to the step of computing the radial force component.

If you believe this is something that would need to be explained in more detail as comments in the source code and/or the documentation, you are welcome to submit a pull request implementing such change(s) or submit a feature request issue or submit a change to the programmer guide documentation to the LAMMPS github repository.

About the name of factor_lj. This is for historical reasons (as for a lot of LAMMPS, keep in mind that this is software that has been around for over 25 years in some form or another). Initially, there were only Lennard-Jones and Coulomb interactions. But factor_lj (as well as input script parameters like evdwl or special_bonds lj) apply to all non-coulomb interactions. In the case of the yukawa potential that is a bit ambivalent, since this is actually a damped coulomb interaction, but since it is not computed based on the per-atom charge parameter Atom::q[i] it counts toward the non-coulomb energy.

The purpose of factor_lj (and factor_coul) is to apply “exclusions” for molecular force fields (like CHARMM, Amber, GROMOS, OPLS, etc). In those cases the non-bonded interactions between pairs of atoms must not be computed if they are connected with a bond (1-2 pair) or the neighbor of an atom that is bonded (1-3 pair) or one atom further away (1-4 pair), those corresponding to topological bond, angles, or dihedrals. The default settings are have all those factors (factor_lj and factor_coul) set to 0.0 and in that case the pairs are removed from the neighbor list entirely. Otherwise the property whether a pair is to be considered as a 1-2, 1-3, or 1-4 pair is encoded into the 2 upper bits of the neighbor atom index (you’ll get 0, 1, 2, or 3 from → sbmask(j), and then strip off those bits with &= NEIGMASK). This is particularly important for long-range coulomb, where you have to subtract the full coulomb and not just the damped interaction computed in the non-bonded term since the kspace solver will always be based on the full interaction and thus not consider the exclusion.

Thank you very much for the confirmation and more detailed explanation, you also managed to clear up a few things I didn’t even ask! Perhaps we are thinking about it differently, but wouldn’t the radial force be in terms of 1/r**7 and 1/r**13? I believe this would get you the 1/r**8 terms and 1/r**14 terms present in fpair after you move the 1/r term from the normalized vector.

For factor_lj, I’m still not quite clear whether it is a potential-dependent parameter (which I would then need to define somewhere), or if it now strictly used for potential-independent bookkeeping and exclusions. I suspect it is the latter.

(edit: I’ll also make a pull request with comments on that file, but this may take a little bit.)

How exclusions are applied and what scaling factors to use are integral part of a force field parameterization for a model that includes explicit bonds. They are not used for models without bonds, but it is good practice to honor them always just in case somebody will use the potential in combination with bonded interactions. Those factors are set with the special_bonds lj or special_bonds lj/coul command. special_bonds command — LAMMPS documentation

You are correct. I overlooked the multiplication with r2inv in the fpair = ... line.