[lammps-users] lammps.numpy.get_neighlist

Hi

I’m trying to understand lammps.numpy.get_neighlist() in lammps.py. In order to do this, I took the ‘dam breakage’ sph example in …\LAMMPS 64-bit 29Oct2020\Examples\USER\sph\water_collapse

I only deleted the line ‘run ${nrun}’ in the water_collapse.lmp file and I ran the following python script

from lammps import lammps

import numpy as np

lmp = lammps()
infile = “water_collapse.lmp”
lmp.file(infile)
lmp.command(“run 100”)
ida = lmp.numpy.extract_atom(“id”)
id_corner=np.where(ida == 6001)[0][0]
print("particle id fluid bottom/left corner = ", id_corner)

nl0=lmp.numpy.get_neighlist(0)
nl1=lmp.numpy.get_neighlist(1)
nl2=lmp.numpy.get_neighlist(2)
n0= lmp.numpy.get_neighlist_element_neighbors(0, id_corner)
n1= lmp.numpy.get_neighlist_element_neighbors(1, id_corner)
n2= lmp.numpy.get_neighlist_element_neighbors(2, id_corner)
print(nl0) # what list is this?
print(nl1) #taiwater list: half list
print(nl2) #rhosum list: full list
print(n0)
print(n1)
print(n2)

particle id fluid bottom/left corner = 977

Neighbor List (9702 atoms) --> this is neighlist(0)
Neighbor List (15702 atoms)--> this is neighlist(1) 
Neighbor List (15702 atoms)--> this is neighlist(2)
(2766, array([[2619], --> this is from neighlist(0,977), but why says 2766? 
       [2621],
       [2763],
       [2618],
       [2765],
       [2620],
       [2767],
       [2622],
       [2769],
       [2901],
       [2764],
       [2903],
       [2905],
       [2768],
       [2907],
       [2902],
       [2904],
       [2906],
       [3025],
       [3027]]))
(977, array([[1153],  --> this is from neighlist(1,977)
      [ 978],
       [1155],
       [ 979],
       [1147],
       [1148],
       [1149],
       [1151],
       [1152],
       [1365],
       [1154],
       [1367],
       [1156],
       [1364],
       [1366]]))
(977, array([[   1], --> this is from neighlist(2,977) 
       [   2],
       [   3],
       [   4],
       [   5],
       [   6],
       [   7],
       [   8],
       [   9],
       [  10],
       [  11],
       [  12],
       [  13],
       [  14],
       [ 973],
       [ 974],
       [ 975],
       [ 976],
       [1153],
       [ 978],
       [1155],
       [ 979],
       [1147],
       [1148],
       [1149],
       [1151],
       [1152],
       [1365],
       [1154],
       [1367],
       [1156],
       [1364],
       [1366]]))

By looking at the source code I can see that pair_sph_taitwater.cpp uses half neighbour list and pair_sph_rhosum.cpp uses the full list. Therefore I assume that neighlist(1) refers to taitwater and neighlist(2) rhosum. I checked all the elements of neighlist(1) and neighlist(2) and it makes sense. 

But what is neighlist(0)? Besides, despite in **get_neighlist_element_neighbors** I ask for element 977, it says that the element is 2766, which is a different atom.
On a similar subject, if, in the python script, I add 

n_rhosum=lmp.find_pair_neighlist('sph/rhosum', request=0)
n_taitwater=lmp.find_pair_neighlist('sph/taitwater', request=1)
print(n_rhosum) #why 0, shouldn't be 2 since the rhosum list is nl2?
print(n_taitwater)

I get
n_rhosum=0

n_taitwater=1

but why n_rhosum=0 shouldn't be n_rhosum=2 since the rhosum list is neighlist(2)? Am I missing something? I'm not sure I understand the request=1 part, but if I use n_rhosum=lmp.find_pair_neighlist('sph/rhosum', request=2), I get n_rhosum=-1

I’m sure I’m misunderstanding something; but can anybody point me in the right direction?

Best

Alessio

OutlookEmoji-15801352865415db36086-22f7-4766-8791-c2d7630e8f89.png

Hi

I’m trying to understand lammps.numpy.get_neighlist() in lammps.py. In order to do this, I took the ‘dam breakage’ sph example in …\LAMMPS 64-bit 29Oct2020\Examples\USER\sph\water_collapse

you are making your life complicated by picking one of the more complex examples. before going over your script code you should just look at the neighbor list info as printed by LAMMPS with the normal input.

Neighbor list info …
update every 5 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 0.039
ghost atom cutoff = 0.039
binsize = 0.0195, bins = 206 411 1
3 neighbor lists, perpetual/occasional/extra = 3 0 0
(1) pair sph/rhosum, perpetual, skip from (3)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(2) pair sph/taitwater, perpetual, half/full from (3)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(3) neighbor class addition, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/2d
bin: standard

What you can see is that there are 3 neighbor list requests, #1 and #2 are requested by the corresponding pair styles, #3 by the neighbor class itself.
request #3 triggest a regular recurring neighbor list build of a full neighbor list for atoms without bond topology.
request #2 triggers a half-from-full neighbor list build with newton_pair set to on. This will post process the neighbor list from #3 which is much faster than a full build from scratch. request #1 triggerst a skip-list build from #3, where the entries of the complete list are selected by atom type according to which pair type combinations the hybrid substyle is active for. again, this is much faster than a full build from scratch, which is the reason why the 3rd request was added, so that only one full build for the superset of neighbor lists is done and then the other two are created via post-processing of the full list.

I only deleted the line ‘run ${nrun}’ in the water_collapse.lmp file and I ran the following python script

from lammps import lammps

import numpy as np

lmp = lammps()
infile = “water_collapse.lmp”
lmp.file(infile)
lmp.command(“run 100”)
ida = lmp.numpy.extract_atom(“id”)
id_corner=np.where(ida == 6001)[0][0]
print("particle id fluid bottom/left corner = ", id_corner)

nl0=lmp.numpy.get_neighlist(0)
nl1=lmp.numpy.get_neighlist(1)
nl2=lmp.numpy.get_neighlist(2)
n0= lmp.numpy.get_neighlist_element_neighbors(0, id_corner)
n1= lmp.numpy.get_neighlist_element_neighbors(1, id_corner)
n2= lmp.numpy.get_neighlist_element_neighbors(2, id_corner)
print(nl0) # what list is this?
print(nl1) #taiwater list: half list
print(nl2) #rhosum list: full list

this is not correct. as you can see from the output, the order is different.
nl0 is the full neighbor list for sph/rhosum
it must be this list, since sph/rhosum is only applied to atoms of type 1 and not type 2. thus there are fewer entries (i.e. for 9702 out of 15702 atoms).

nl2 is the full list for all atoms from which the two other lists are built through post-processing.

particle id fluid bottom/left corner = 977

Neighbor List (9702 atoms) --> this is neighlist(0)
Neighbor List (15702 atoms)--> this is neighlist(1) 
Neighbor List (15702 atoms)--> this is neighlist(2)
(2766, array([[2619], --> this is from neighlist(0,977), but why says 2766? 

the “element” index is not the local atom index, but an index from 0 to the size of the neighbor list class.
rather than doing the lookup again, you can loop over the neighbor list with:

nl0=lmp.numpy.get_neighlist(idx_rhosum)
for i in range(0,nl0.size):
idx, nlist = nl0.get(i)
print("id = “,idx,” list = ",nlist)

the output will then be the local atom index in “idx” for which nlist is the corresponding neighbor list.
as it happens, local index 977 is the first atom in the list (since it is the first atom of type 1 in the data file).
mind you this order can change at every reneighboring.


By looking at the source code I can see that pair_sph_taitwater.cpp uses half neighbour list and pair_sph_rhosum.cpp uses the full list. Therefore I assume that neighlist(1) refers to taitwater and neighlist(2) rhosum. I checked all the elements of neighlist(1) and neighlist(2) and it makes sense. 

no it doesn’t. see my explanation from before.


But what is neighlist(0)? Besides, despite in **get_neighlist_element_neighbors** I ask for element 977, it says that the element is 2766, which is a different atom.
On a similar subject, if, in the python script, I add 

n_rhosum=lmp.find_pair_neighlist('sph/rhosum', request=0)
n_taitwater=lmp.find_pair_neighlist('sph/taitwater', request=1)
print(n_rhosum) #why 0, shouldn't be 2 since the rhosum list is nl2?

the request argument is actually not what it should be. this needs to be corrected in the LAMMPS code.
it should be the copy of the neighbor list request ID, which can be set by styles if the style requests multiple neighbor lists with different settings (some pair styles list meam/c and fixes do that). this will be corrected in a future LAMMPS patch release. right now the request should be the index of the list of neighbor lists in the LAMMPS output minus one.

so if you want to find the neighbors of atom 977 you have to check against the first item of the tuple returned by the “get()” method of the neighbor list class.


print(n_taitwater)

I get
n_rhosum=0

n_taitwater=1

but why n_rhosum=0 shouldn't be n_rhosum=2 since the rhosum list is neighlist(2)? Am I missing something? I'm not sure I understand the request=1 part, but if I use n_rhosum=lmp.find_pair_neighlist('sph/rhosum', request=2), I get n_rhosum=-1

that is the correct behavior.

I’m sure I’m misunderstanding something; but can anybody point me in the right direction?

it will hopefully be better documented the next patch release.

axel.

Hi Alex

thank you very much: it makes sense now

Best

Alessio

OutlookEmoji-15801352865415db36086-22f7-4766-8791-c2d7630e8f89.png