I am hoping to extract the pairwise distances between two atom types (which I have setup using a group command for those atom types, and the relevant compute command, labelled “d0”), within Python.
If I call
lmp.extract_compute("d0", 2, 1)
(cstyle=2 <=> local, ctype=1 <=> vector) I can get access to the ctype pointer to the vector, which can be iterated through to extract values. Is there a way of identifying the length of the underlying cvector directly so that I can convert this directly into a NumPy array. Or does one need to access the size of each neighbour list for each atom involved? I noticed calling
lmp.numpy.extract_compute("d0", 2, 1)
results in a NULL pointer access (with particular reference to
nrows = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_VECTOR) on line 170 of the numpy_wrapper.py script).
Thanks in advance for any suggestions/comments/criticisms of my approach.
Using LAMMPS (29 Sep 2021 - Update 3), on MacOS-Monterey (Intel Core i9).
Which commands exactly did you issue for that?
lmp.command("group CC type 2 3")
lmp.command("compute d0 CC pair/local dist")
Thanks. I will have to take a closer look at this, but that will take some time.
One recommendation: please prefer the symbolic constants over numbers since those may change in case of errors in the documentation or implementation.
@tcnicholas after looking at the details of the implementation, this looks like a combination of incorrect usage and a bug.
- when querying local computes you should (with your version of LAMMPS) always check for
LMP_SIZE_COLS. In case you have a local vector the number of columns will be 0. This is how local computes work internally and that was copied over to the library interface. However, it makes sense to have
LMP_SIZE_VECTOR act like an alias to
LMP_SIZE_ROWS, but that needs to be done at the C++ level in the
- the numpy wrapper is incorrectly treating the LMP_TYPE_GLOBAL and LMP_TYPE_LOCAL case the same. Instead it should treat this separately and then for local computes decide the NumPy array dimensions based on whether there are zero columns or not. The type property that would distinguish between scalar, vector, or array will be ignored (since there cannot be a local scalar, and local arrays and vectors are described with the same internal variables).
If you apply the following changes, recompile, and reinstall the python module, things should work as expected. This change will be included into the next patch release -which will also become the next stable release - later this week. Because of the new stable release there will be no updates and backports to the 29 Sep 2021 release.
Please let me know ASAP, if there are still issues.
diff --git a/python/lammps/numpy_wrapper.py b/python/lammps/numpy_wrapper.py
index 20ec85d001..73ee7bbe5e 100644
@@ -165,7 +165,7 @@ class numpy_wrapper:
value = self.lmp.extract_compute(cid, cstyle, ctype)
- if cstyle in (LMP_STYLE_GLOBAL, LMP_STYLE_LOCAL):
+ if cstyle == LMP_STYLE_GLOBAL:
if ctype == LMP_TYPE_VECTOR:
nrows = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_VECTOR)
return self.darray(value, nrows)
@@ -173,6 +173,13 @@ class numpy_wrapper:
nrows = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_ROWS)
ncols = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_COLS)
return self.darray(value, nrows, ncols)
+ elif cstyle == LMP_STYLE_LOCAL:
+ nrows = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_ROWS)
+ ncols = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_COLS)
+ if ncols == 0:
+ return self.darray(value, nrows)
+ return self.darray(value, nrows, ncols)
elif cstyle == LMP_STYLE_ATOM:
if ctype == LMP_TYPE_VECTOR:
nlocal = self.lmp.extract_global("nlocal")
diff --git a/src/library.cpp b/src/library.cpp
index dbe665828e..bbe7b84324 100644
@@ -1738,6 +1738,7 @@ void *lammps_extract_compute(void *handle, const char *id, int style, int type)
if (type == LMP_TYPE_SCALAR) return (void *) &compute->size_local_rows; /* for backward compatibility */
if (type == LMP_TYPE_VECTOR) return (void *) compute->vector_local;
if (type == LMP_TYPE_ARRAY) return (void *) compute->array_local;
+ if (type == LMP_SIZE_VECTOR) return (void *) &compute->size_local_rows; /* alias for LMP_SIZE_ROWS */
if (type == LMP_SIZE_ROWS) return (void *) &compute->size_local_rows;
if (type == LMP_SIZE_COLS) return (void *) &compute->size_local_cols;