Not sure this rather physical rather than lammps-specific technical question is allowed to be asked here, so feel free to ignore or delete it.
For a bilayer graphene system, I have plotted the forces that I obtain directly from lammps using dump ID group-ID style N file fx fy fz. If I plot these forces, I correctly and as expected get a map where the values of each atom are smaller than my force tolerance criterion given for the minimization.
Now, when instead, I calculate the potential energy profile using compute pe all pe/atom
followed by dump ID group-ID style N file c_pe, after which I take the negative gradient of this profile to get the forces, I get non-zero forces on the atoms.
What am I missing, physically or computationally, to match these two outcomes?
On the one hand, I do expect forces to be nearly zero as this is the definition of the stopping criterion, but at the same time, I also expect some atoms to feel the need to move away from the energetically less favorable stacking configurations that exist in twisted bilayer graphene, hence a non-zero force vector.
Dump files output floating point numbers with only 6-7 digits precision (using the printf() style %g format). That can lead to significant differences to values computed based on the native precision of double precision floating point numbers (15-16). You can use dump_modify to change this to, e.g. %.16g).
Another possible reason for differences may be inconsistencies between energy and force in the force styles. This can be tested internally in LAMMPS using the fix numdiff command — LAMMPS documentation
I’m not sure the mathematics quite works out. Consider a simulation box with two identical LJ particles. At almost all distances (within the cutoff) the force between them will be non-zero – but I think LAMMPS will always print that they both have the same potential energy using compute pe/atom.
The force is the gradient of the energy with respect to changes in position of one particle, not the gradient of the energy between one particle’s position and another particle’s position. In any case the per-atom energy generally depends on your convention for dividing up interaction energies (of the three atoms in an angle interaction, how much energy does each atom own? What about the four atoms in a dihedral?).
Apologies if you had a very different bit of math in mind.
Indeed, I got confused with the force calculated using the gradient using the changes in the position of one particle (as done by lammps using numdiff for instance) and the forces I calculated using the position of one atom and another atom… the latter observable is the one i’m interested in to explain some of my physics.
Using the numdiff suggestion from Alex I indeed got zero forces similar to the ones obtained with fx fy fz values, so no error in the way the code was calculating the forces.