python lmp.extract_variable for atom variable

Hi lammps users !

In python coupling, is it possible to get the vector of an atom variable ?
This is not working for me :

  vector = lmp.extract_variable(variable,0,0)

where variable is the atom variable name in lammps script. Maybe with a loop ?
An other related question would be how to impose an atom variable value to lammps given a vector in python ?
something like

  lmp.command("variable Ux atom %g" % vector)

However, this would work for scalar, but not vectors ?..
Some suggestions ?

Thanks to all,
joris

You need to look at the functions in lib/library.cpp.
They are all wrapped by the Python wrapper.
The example python/example/simple.py does what
you are asking via the gather/scatter functions.

You also need to think about whether you are running
in parallel or serial. You also need to realize that
reading is easier than writing. What you see from
the Python side is a copy of what is inside LAMMPS.
So you need a method that copies the Python values
back into the LAMMPS internals. That is what
the scatter method does.

Steve

ok, I will try scatter.
scatter in parallel is not possible ?
When you say "reading is easier than writing", does it mean that writing affect badly the computationnal effort ?

Hi all,

As only lammps_scatter_atoms is available in python coupling, I wrote a routine lammps_scatter_variable in library.cpp that scatter a vector data into a atom-type variable in lammps. I also add the routine definition to library.h and add it to lammps.py. I built again the whole shared lammps library (make clean-all, rm *.so, make makeshlib, make -f Makefile.shlib openmpi).
However, when I call it in python I get the error:

  Traceback (most recent call last):
   File "plot.py", line 80, in <module>
     lmp.scatter_variable(vari,"nve_group",value3)
   File "/home/joris/LIGGGHTS/LIGGGHTS-PUBLIC/python/lammps.py", line 171, in scatter_variable
     self.lib.lammps_scatter_variable(self.lmp,name,group,data)
   File "/usr/lib/python2.7/ctypes/__init__.py", line 378, in __getattr__
     func = self.__getitem__(name)
   File "/usr/lib/python2.7/ctypes/__init__.py", line 383, in __getitem__
     func = self._FuncPtr((name_or_ordinal, self))
AttributeError: /home/joris/LIGGGHTS/LIGGGHTS-PUBLIC/src/liblmp_openmpi.so: undefined symbol: lammps_scatter_variable

*Somehow like if it didn't find the new function I created. Is there an other modify in order to take the new function into account ?**

Ok I saw my error : void lammps_scatter_variable(void *ptr, char *name, char *group, *double* *data).

Anyway, it is not working, the variable atom values seems not to be changed. I guess that lammps is evaluating again the variable at each time step so that no change can be saved permanently ...

Basically, I want to add a constant force to each atoms from python vector (command addforce in lammps), then compute Lammps for a few time step (with the same force), then collect some atom properties to compute new forces in python and start again.

Any ideas ?

Thanks,
j

Gather and scatter work in parallel - they take contributions
from all procs and gather it to one proc, or scatter it out
to all the procs. Extract_atom gets you each proc's atomic
data. If you want to change it and write that back into each proc's data,
I think you would need another function to do the write
from Python to C internal. Note that you don't have this
problem using the lib from C/C++. If you extract the
atom pointer, you can read/write directly into the internals
of LAMMPS.

Steve

I don't think it makes any sense to write to a LAMMPS
variable, atom-style or otherwise. LAMMPS creates
the values in a variable, when it needs to evaluate the
variable. It won't look at any values you have put there.

I think you want to modify/overwrite the per-atom force
which is an atom property.

Steve

Ok, I realize now that it is not a good way of doing.

To overwrite a per atom force, lmp.scatter_atoms() should thus do the job from python ?

My problem is that if I scatter the force at the begining and then run lammps for a few time step,
  the added force won't remain during the several time step...

Coupling would work only if I make a python loop over one lammps time step and then scatter the force again?...

If you want to modulate the forces on atoms every step,
from an external program, then I don't see how you
have any choice but to make your external program
do something every step. A few ideas:

(1) fix addforce allows for use of a LAMMPS variable
(atom or equal style) which gives a great deal of flexibility
in modifying forces every steps. Are you sure you can't
do what you want with that command?

(2) If you do use the Python wrapper, then be sure to
use the run pre no post no options so that you can
do successive one-step runs with minimal overhead.

(3) There is a fix external command which allows LAMMPS
to call out to an external program every timestep. It could
possibly be modified/enhanced to do something like you
want with forces However, I've never thought about using that with
Python.

(4) There is probably a way via Numpy to get a Python NumPy array
to directly point to a LAMMPS internal array, like per-atom forces.
That would be the lowest overhead way to modify a LAMMPS internal
from Python. Of course you would have to worry about LAMMPS
reallocating the array, which it can do occasionally.

Steve

Ok, thanks for the suggestions.

Maybe I should have explain better my model.

Basically, I wanted to compute variables that depend on their previous step value
(to generate time-correlated auto-regressive random process) in order to use them as an added atom force (drag force from external field).

As we can't do that directly in lammps, I thought about python coupling.

As the external fluctuating field is changing slowly (long time correlation), my idea was to compute it from python,
add constant forces to atom via a scatter type function, and then run 1000 steps of Lammps with this constant forcing.
This would be repeated several time.

I think you could do all of this within LAMMPS if you wrote a fix.
Comments below.

Steve

Ok, thanks for the suggestions.

Maybe I should have explain better my model.

Basically, I wanted to compute variables that depend on their previous step
value
(to generate time-correlated auto-regressive random process) in order to use
them as an added atom force (drag force from external field).

A fix can easily store some atom-based value from a previous step.

As we can't do that directly in lammps, I thought about python coupling.

As the external fluctuating field is changing slowly (long time
correlation), my idea was to compute it from python,
add constant forces to atom via a scatter type function, and then run 1000
steps of Lammps with this constant forcing.

Assuming you can write a variable formula that represents the external
field (space and time dependent), you can do this with fix addforce using
the variable.

Yeah, that would be great. My competence in C++ are a bit weak, but I'll give a try.
As you're really helpfull on the forum, it could work!

Tell me if I'm wrong, there would be 3 fixes to be done:

1/ be able to access last time step atom value. This should be saved in memory somewhere at the end of the time step.
I should write a function like: /save atom variable/ that is called at the end of the time step. The variable should be then available from any variable like function. But how is managed the memory in lammps? can we acces any variable from any other function ? What is the function that loops over time steps ? Where to do the call to this "save" function ?

2/ Modify the variable.cpp code to include the possibility of a call to a saved variable. Maybe define a new class like s_variable (similar to v_variable or c_variable) that refer to this kind of saved variable ?

3/ (optionnal) Modify addforce.cpp so that the forces are evaluated only every few time step (to save cpu time) when the forcing is slow compared to lammps time step.

Thanks for the comments and for your help Steve,
Joris

None of these are necessary. The fix can save the last timestep
value, as well as compute the current timestep value and the
resulting force. It can make that force available as an output,
which fix_addforce can access (thru a variable). Fix addforce
doesn't need to be modified. Just set the accessed value to
0 or whatever you want it to be on each timestep.
You shouldn't have to modify any other part of LAMMPS to
do what you want. Just write the new fix.

Steve

(4) There is probably a way via Numpy to get a Python NumPy array
to directly point to a LAMMPS internal array, like per-atom forces.
That would be the lowest overhead way to modify a LAMMPS internal
from Python. Of course you would have to worry about LAMMPS
reallocating the array, which it can do occasionally.

Playing with python/example/simple.py, I realized that what I said
above already just works. You can use extract_fix to get a pointer
to an internal per-atom vector or array stored by the fix, or extract_atom to do
likewise for the atom coords or forces. I thought by the time
this gets back to Python, that the values pointed to were copies
of what is stored inside LAMMPS, but this is not the case.
The Python variable is accessing the internal values directly.
So this means you can overwrite the internal values directly from the
Python script. E.g.

x = lmp.extract_atom("x",3)
x[0][0] = 10.0

Steve

Ok, I think I got it, I wrote a slightly different extract_fix for property/atom vector fix, and now I can get the array value from python.
As you said, a change is applied directly into lammps memory ! Here I replace the initial 0.001 value to 1 :

   x =lmp.extract_property("Pforce")
   val=x[0][0]
   print("x=%g" %val)
   x[0][0]=1.0

   x =lmp.extract_property("Pforce")
   val=x[0][0]
   print("x=%g" %val)

Output :

x=0.001
x=1

Great !! Now I will be able to move on more interesting stuff :wink:
Thanks a lot for helping Steve,

Joris

What does extract_property() do, and how do you
invoke it externally every timestep?

Steve

My extract_property function do exactly the same as extract_atoms() but for a property/atom array field as defined in Liggghts for the drag force, or for the Pforce.cpp file that I sent you.

You're right, the post_force fix has nothing to do with python, but this "fix Pforce", together with a post_force function, allows me to add a force every lammps time steps without re-evaluate it at each time step. It is some kind of constant forcing on the atoms. It is an atom property that is used for the array of forces. I modify it every 1000 time steps via python. I could not do it with addforce, because addforce evaluate at each time step the input variable, while all the computation of variable is done in the python script and transferred every few time steps to lammps.

It seems to work nicely on a single processor, but today, I tried to upgrade it to mpi computations. I add to create a gather_property and scatter_property function similar to gather_atoms and scatter_atoms to deal with several processor. I got no errors, but it seems that the property/atom style has an other way of indexing atoms, so that the forces are a bit mixed up, and I get a different result than with single processor even if the map keyword is used. I'll try to fix that tomorrow... But if you have some suggestions, they would be welcomed !!

Cheers,
Joris

I just extended fix external to make more options for
how external driver programs interact with LAMMPS while
it is running. This may be useful for what you are trying to do.

Steve