Passing back force using extract_atom() for python wrapper


I am trying to use python to run a simple LAMMPS simulation where I modify the force on each individual atom in a unique way. However when I try to modify the force using both extract_atom() and gather_atom()/scatter() neither one appears to modify the force. I am fairly sure it is not working because I set the force to 0 just to test it and am still having atoms move as if forces were applied. Below is a snippet of what I am trying

tmp_force = lmp.extract_atom(“f”,3)
for idx in range(n_atoms):
for jdx in range(3):
tmp_force[idx][jdx] = 0

Run simulation a single step

lmp.command(“fix simul all nve”)
lmp.command(“run 1”)

And the results are

simulation step:  0

 Atom:  1
force:  -0.9984898859777808 
 position:  0.0
force:  0.0015101140222192285 
 position:  0.0
force:  0.0 
 position:  0.0

simulation step:  1

 Atom:  1
force:  -0.9982779080583614 
 position:  -1.248112357472226e-05
force:  0.0016166350850116125 
 position:  1.887642527774036e-08
force:  0.0 
 position:  0.0

where I have only included the first atom as well as the initial forces it feels

To my understanding extract_atom() is essentially a pointer to the actual lammps variable and changing the value in the python script should change it in LAMMPS as well. I get the same results if I try using the gather_atom/scatter technique as well.

Any help would be appreciated, thanks

There is a conceptual problem with what you are trying to do. You are making changes to the forces outside the MD loop(s), so they will be lost during the MD stepping. With “run 1” the force arrays are wiped out and forces recomputed twice: once during the setup phase and a second time, when doing the one time step.

If you want to make modifications to the forces you have to either store those in an atom style variable and then use either fix add force or fix set force with that variable, of you have to put your python code into a function and call it from inside the MD loop via fix python/invoke.

Please also note, that arrays that are made accessible with extract_atom() are of the dimension nlocal and not natoms (except when you have only one MPI rank when nlocal == natoms)