Neighbors of neighbors in a compute

Dear lammps users,

I am writing a compute that needs to access for each atom the neighbors of its neighbors, running on parallel processors. The algorithm goes like this:

  1. For each atom i, loop through its neighbors j.
  2. For each neighbor j, loop through its neighbors k.
  3. For each atom i, store a list of k neighbors of neighbors.
  4. Calculate distances between some i’s and k’s.

The problem I’m having is that when a neighbor j is a ghost atom, I do not have access to its neighbor list. How can I get access to the neighbors of ghosts? I have seen that there is an option in the NeighRequest class to set ghost = 1 if I want neighbors of ghosts but using that option I haven’t been able to find the neighbors of the ghost in the ilist array. Where are they stored, and what information about the neighbors of ghosts does this option store? I need ID’s and coordinates for all neighbors as well as the neighbors of neighbors.

Thank you,
Abdullah

This is a required list for compute CNA. How does that piece of code do it?

You’re right, and I have looked at the CNA code. My apologies, I actually left something out in my question by accident. What I actually need is also the neighbor list of k atoms. So the full algorithm would look like:

  1. For each atom i, loop through its neighbors j.

  2. For each neighbor j, loop through its neighbors k.

  3. For each atom i, store a list of k neighbors of neighbors.

  4. Calculate distances between some i’s and k’s.

  5. Loop though neighbors of k and Identify common neighbors between i and k.

The problem here arises when if j is a ghost, and then k is a neighbor of j but not a neighbor of i (and might not even be a ghost of i), how can I access the neighbors of k?

You can request a neighbor list that includes neighbors of ghost atoms

(within limitations, ghost atoms do not extend infinitely far).

Do this grep from the src dir:

grep “ghost = 1” *.cpp /.cpp | grep neighbor

and you will see some pair styles that use neighbors of neighbors.

Steve

Thank you for that suggestion Steve. I have tested my compute using the chairebo potential which uses the ghostneigh = 1 option and my compute finally works and is doing what it needs to do. Am I correct in understanding that the “ghost = 1” option used in a compute only requests the information from the neighbor list, but the neighbor list will not actually compile neighbors of ghosts unless the potential being used has the “ghostneigh = 1” option?

The problem now of course is that I don’t want to use chairebo potential, I would like to use tersoff and its parameterizations, EA most importantly. I saw that the KOKKOS package has a tersoff potential that uses “ghostneigh = 1” option. Can this potential be used for the other parameterizations of tersoff like SiC_Erhart-Albe or SiC_1994 for example? Is it possible to build kokkos_mpi as a library? I have been having issues trying to build it.

Alternatively, if I simply modify the tersoff potential file to include the “ghostneigh = 1” option, would that cause any issues downstream in the tersoff potential code?

Thank you,

Abdullah

~WRD000.jpg

Comments below.

Steve

~WRD000.jpg

Thank you again for your time Steve. So in that case there must be something else I’m missing to explain why my compute works with airebo potential but not with tersoff, all else being unchanged. For trouble shooting purposes I have the compute outputting values for inum and gnum and for some reason when I use the tersoff potential, my compute outputs gnum = 0 meaning no neighbor lists for ghosts were generated, but when all I do is switch to airebo potential I get gnum = 2995 for a 64 particle system. inum = 64 in all cases (I am testing on a single processor currently). Do you have any ideas why this might be?

Or perhaps the real question should be, why does my compute get gnum = 0 even though I have set in my compute’s init() function the following line: neighbor->requests[irequest]->ghost = 1;

Any thoughts or insight would be highly appreciated. Thank you,
Abdullah

~WRD000.jpg

Before the LAMMPS run starts it prints info on all the neighbor lists

it is using. What does it say about yours?

Steve

~WRD000.jpg

Here is what it says:
Neighbor list info …
2 neighbor list requests
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.46
ghost atom cutoff = 6.46
binsize = 3.23, bins = 4 4 4

and after a run 0, it says:

Nlocal: 64 ave 64 max 64 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 557 ave 557 max 557 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 2944 ave 2944 max 2944 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 2944
Ave neighs/atom = 46
Neighbor list builds = 0
Dangerous builds = 0

Also it may help to know that I don’t have my compute command from the very beginning. This run is part of a coupled C++ program that calls lammps as a library. I run a few MD steps without this compute, then I call the compute, call dump outputting the compute, run 0 to calculate & dump compute values, undump, uncompute, run more md steps and repeat.

I think I’ve narrowed down my problem a little bit, it seems that my compute does not trigger a new neighbor list build (or at least not the full_bin_ghost() method). It is using whatever neighbor list was there in the last step. Although this is desirable, it would only work in this case if the already existing neighbor list was listing neighbors of ghosts as well, which for tersoff runs, it does not. This explains why my compute worked when I tried it with airebo potential which does use neighbors of ghosts. In my compute init() method I am using the following settings:

int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->occasional = 1;
neighbor->requests[irequest]->ghost = 1;

and in my compute_peratom() method I call for:

invoked_peratom = update->ntimestep;
neighbor->build_one(list);

before I start looping through any atoms or neighbor lists for processing. Again I thank you and highly appreciate your support with this issue,

Abdullah

~WRD000.jpg

You’re using an old version of LAMMPS. The current version

(for a while) prints out a lot info about individual neighbor

lists, which would be helpful in your case to diagnose if

the neigh lists are setup correctly.

Steve

~WRD000.jpg