How do I get pe/atom values into a numpy array using PyLammps?
I tried this in my python script (import PyLammps)
L.compute(" eng all pe/atom “)
L.dump(" dumpall all custom “,10," Ddumps/d.* id type x y z c_eng “)
L.run(0)
engs = np.array([L.eval(“c_eng[”+str(i+1)+"]") for i in range(natoms)])
This works only sometimes. That is, it worked the first time, but then when
I tried a similar call later in the program it failed. Specifically:
L.minimize(" 1.e-6 1.e-8 10000 10000 “)
engs = np.array([L.eval(“c_eng[”+str(i+1)+"]") for i in range(natoms)])
gives
ERROR: Compute used in variable between runs is not current (…/variable.cpp:1480)
Last command: $(c_eng[1]
First, it seems like using L.eval is an awkward construction.
Second, why does it fail the second call???
There is probably a better way.
Using LAMMPS 12 Dec 2018. Perhaps I need to update?
…Murray
PyLammps doesn’t have a convenient function for what you are looking for, but the the parent type lammps, there is extract_compute.
in general, PyLammps is more for simple and interactive demos. if you want to do “serious” programming, it is almost always better to use the more low-level lammps class. but you can mix this, since you have access to the lammps class through the member variable “lmp”
so to get a numpy array of the data of the compute, you would do:
engs = L.lmp.extract_compute(“eng”,1,1)
some more comments:
per-atom data is distributed across processors. so its length is “nlocal”. the equality “nlocal == natoms” only holds true for the special case of doing a serial run.
you get the “compute between runs is not current” error, because compute pe/atom is not really a computation, but copying data, that was accumulated during the last force compute. for that to work, the force compute has to be told, that the potential energy needs to be tallied (and in per-atom arrays) as it normally is not done to speed up the run. and to trigger the tallying, you need something that “consumes” the data like your dump command. as soon as anything is changed or there is no consumer, you cannot access the per-atom compute anymore.
the workaround for that, is to have a “cached copy”. one way to do this would be to define a fix ave/atom 1 1 1 c_eng, which would be a “consumer” of the pe/atom compute and trigger its computation in every force computation and keep a local copy inside the fix, which will remain available in between runs.
this last issue is a general LAMMPS “feature” and has nothing to do with PyLammps (or the LAMMPS version).
HTH,
axel.