Total force from pair style

Hello,

The total force on each atom in my system comes from 2 fixes and the lj/cut pair_style.
If I dump fx,fy then it will include all three sources, how can I only dump the force that comes from pair_style?

I think compute pair/local can dump all the forces but between each two pair and then I will have to sum them to get what I wanted, is there any other direct way?

Thank you in advance

That would be hugely inefficient.

It is a bit of a complex issue. What you need is to get access to the forces after the pair style and before any fixes that operate in the “post_force” step. The direct approach would be to have a fix that is defined before the other fixes make a copy of the forces and then make those available. This could be done with custom C++ code that could be adapted from existing fixes that do similar things, but that can be tedious.

After some thinking, I found a way to do this without having to write C++ code, only python.
This uses the LAMMPS python module and the PYTHON package which provides the fix python/invoke command.

The general strategy:

  • use fix property atom to define a custom per-atom array with 3 elements (no ghost atoms needed)
  • define a python function that uses the LAMMPS python module to extract pointers to the force array and the custom property array and copy the force data to the custom property
  • use fix python/invoke to run this python function at every step during post_force processing

Here is a modified version of the melt example that does just this:

# 3d Lennard-Jones melt

units		lj
atom_style	atomic

lattice		fcc 0.8442
region		box block 0 10 0 10 0 10
create_box	1 box
create_atoms	1 box
mass		1 1.0

velocity	all create 3.0 87287 loop geom

pair_style	lj/cut 2.5
pair_coeff	1 1 1.0 1.0 2.5

neighbor	0.3 bin
neigh_modify	every 20 delay 0 check no

fix 0 all property/atom d2_pairforce 3 ghost no

python post_force_callback here """
from lammps import lammps
from lammps.constants import LAMMPS_DOUBLE_2D
import numpy as np

def post_force_callback(handle, vflag):
    lmp = lammps(ptr=handle)
    nlocal = lmp.extract_setting('nlocal')
    forces = lmp.numpy.extract_atom('f', LAMMPS_DOUBLE_2D, nlocal, 3)
    storage = lmp.numpy.extract_atom('d2_pairforce', LAMMPS_DOUBLE_2D, nlocal, 3)
    for i in range(0, nlocal):
        storage[i][0] = forces[i][0]
        storage[i][1] = forces[i][1]
        storage[i][2] = forces[i][2]
"""

fix		1 all nve
fix             2 all python/invoke 1 post_force post_force_callback
#fix             3 all addforce 0.1 0.2 0.3
dump		id all custom 50 dump.melt id x y z fx fy fz d2_pairforce[1] d2_pairforce[2] d2_pairforce[3]

thermo		50
run		250

You can see the copy of the forces in the dump file. If you uncomment the fix addforce line, you will see that the pairforce values remain unchanged while the fx, fy, fz properties have the forces added from the fix.

There is just one little glitch in this. Fix python/invoke is not used during the “setup” phase, so the forces for step 0 are not copied. This behavior is unexpected and can be considered a bug.
A fix for that would be the following patch to the LAMMPS sources:

  diff --git a/src/PYTHON/fix_python_invoke.cpp b/src/PYTHON/fix_python_invoke.cpp
  index 786c48568d..281379cda4 100644
  --- a/src/PYTHON/fix_python_invoke.cpp
  +++ b/src/PYTHON/fix_python_invoke.cpp
  @@ -105,6 +105,13 @@ void FixPythonInvoke::end_of_step()
   
   /* ---------------------------------------------------------------------- */
   
  +void FixPythonInvoke::setup(int vflag)
  +{
  +  if (selected_callback == POST_FORCE) post_force(vflag);
  +}
  +
  +/* ---------------------------------------------------------------------- */
  +
   void FixPythonInvoke::post_force(int vflag)
   {
     if (update->ntimestep % nevery != 0) return;
  diff --git a/src/PYTHON/fix_python_invoke.h b/src/PYTHON/fix_python_invoke.h
  index 12f463501f..09382e5780 100644
  --- a/src/PYTHON/fix_python_invoke.h
  +++ b/src/PYTHON/fix_python_invoke.h
  @@ -30,6 +30,7 @@ class FixPythonInvoke : public Fix {
     FixPythonInvoke(class LAMMPS *, int, char **);
     ~FixPythonInvoke() override;
     int setmask() override;
  +  void setup(int) override;
     void end_of_step() override;
     void post_force(int) override;
   

What about fix store/force?

Absolutely. I forgot about that one. This is exactly the C++ code required.
LAMMPS has too many features… :frowning:

At least, by doing this python hack, I found a bug in fix python/invoke. :smiley: