trouble with python interface

I seem to be having trouble using the python interface to do something.

I want to be able to displace an atom and then query the forces that are generated.

So I create a ctype thing using gather_atoms
x0ctype = lmp.gather_atoms("x",1,3)

["lmp" is my instance of lammps]

As I understand it, this array now has x,y,z for the first atom, then x,y,z for the second, and so forth

Then I displace an atoms by modifying x0ctype
For example,
x0ctype[0] += h

then I use scatter_atoms to put the entire array back into LAMMPS' memory locations:
lmp.scatter_atoms("x",1,3,x0ctype)

[I have set
atom_modify map hash
which I understand is necessary for this to work.]

Then I execute a single step

lmp.command("run 0 post no")

To query the forces and make an array which corresponds to the atoms as contained in x0ctype,
I do this:

    fx = lmp.extract_variable("fx","all",1)
    fy = lmp.extract_variable("fy","all",1)
    fz = lmp.extract_variable("fz","all",1)
    fs = np.array([ [fx[i],fy[i],fz[i]] for i in range(len(fz)) ]).flatten()

(np is "numpy")

In the input to lammps, I have

variable fx atom fx
variable fy atom fy
variable fz atom fz

to define the variables.

This does not seem to be working. The forces are not as expected.
I seem to be doing something fundamentally wrong.

Any suggestions would be appreciated.

If works for me if I add these lines to python/examples/simple.py

(at the end)

f = lmp.extract_atom(“f”,3)
print f[0][0]

fx = lmp.extract_variable(“fx”,“all”,1)
print fx[0]

and these lines to python/examples/in.simple

variable fx atom fx
dump 1 all custom 1 tmp.dump id type fx

In particular the fx on the 1st atom in the dump file

(for step 21), is the same as what Python prints out.

Note however, that by default LAMMPS sorts the
atoms spatially (for performance). This happens at
the beginning and during each run. Which means
the internal force array (and the internal force variable)
that you are accessing, will not typically list atoms in order from

1 to N, even on step 0. You can see this in the dump file,
after the 2nd run.

You can prevent that by using atom_modify sort 0 0.0.

Or you could reorder the force array by atom ID on
the Python side, once you have grabbed it.

Steve

I should add that the library interface functions for gather/scatter atoms

pull atom properties out of LAMMPS into a Python array
and reorder them for you. You get the full atom list, even in parallel.

The extract_variable function pulls the values out of LAMMPS
and copies them into a Python vector, but does no reordering.

You are also only getting one processor’s values (for its atoms).

The extract_atom function simply returns a pointer to the

atom values stored (by one processor) within LAMMPS. So
if you overwrite a value in Python, it is changed inside LAMMPS.

So you would not want to reorder that Python vector in place.

Steve