Neighbor List correlation function - issue storing neighbor list

Dear Lammps developers,

I am creating a compute style to calculate the neighbor list correlation function (NLCF) following the paper:

Rabani, E., Gezelter, J. D., & Berne, B. J. (1997). Calculating the hopping rate for self-diffusion on rough potential energy surfaces: Cage correlations. The Journal of chemical physics, 107(17), 6867-6876.

Our research group already has a python script that works as a post-processing tool; It uses the coordinates of the atoms that are in the dump files to create a neighbor list for each atom based on a cutoff (smaller than the force cutoff). However, We think it would be more convenient to have this feature as a compute style.

We choose to use ‘compute_orientorder_atom’ as a starting point for ‘compute_nlcf’ since this compute style builds neighbor lists based on a cutoff and iterates through every atom to perform the computation.

One of the main challenges we have faced during the development of compute_nlcf is that in order to compute the NLCF at any time step, we need to retrieve the neighbor list of each atom at time = 0; meaning that we have to store it in memory.

We have seen that one method to store data, such as atom coordinates, is to use fix_store. So far we think that fix_store only supports saving data within a ‘n by m’ matrix (please correct me if I am wrong). We would like to use a fix_store or some other type of LAMMPS implemented method to save our neighbor lists at time = 0, but a key issue is that every neighbor list can ultimately be a different size.

For now, we have implemented this storing process using a multi-dimensional dynamic array int** that is a member variable of the compute_nlcf class. When compute_peratom() is invoked we use a flag variable to determine if time = 0. If it is, then for each atom ‘i’ and its neighbor list ‘nearest[ ]’, we dynamically allocate enough memory (size of nearest[ ]) and store this neighbor list in the int** member variable.

Please, if there is any more efficient method already implemented that we can use to store the neighbor list at time = 0, please let me know.

I am using LAMMPS version 22 August 2018.

Thank you.

-S

Dear Lammps developers,

I am creating a compute style to calculate the neighbor list correlation function (NLCF) following the paper:

Rabani, E., Gezelter, J. D., & Berne, B. J. (1997). Calculating the hopping rate for self-diffusion on rough potential energy surfaces: Cage correlations. The Journal of chemical physics, 107(17), 6867-6876.

Our research group already has a python script that works as a post-processing tool; It uses the coordinates of the atoms that are in the dump files to create a neighbor list for each atom based on a cutoff (smaller than the force cutoff). However, We think it would be more convenient to have this feature as a compute style.

We choose to use ‘compute_orientorder_atom’ as a starting point for ‘compute_nlcf’ since this compute style builds neighbor lists based on a cutoff and iterates through every atom to perform the computation.

One of the main challenges we have faced during the development of compute_nlcf is that in order to compute the NLCF at any time step, we need to retrieve the neighbor list of each atom at time = 0; meaning that we have to store it in memory.

We have seen that one method to store data, such as atom coordinates, is to use fix_store. So far we think that fix_store only supports saving data within a ‘n by m’ matrix (please correct me if I am wrong). We would like to use a fix_store or some other type of LAMMPS implemented method to save our neighbor lists at time = 0, but a key issue is that every neighbor list can ultimately be a different size.

For now, we have implemented this storing process using a multi-dimensional dynamic array int** that is a member variable of the compute_nlcf class. When compute_peratom() is invoked we use a flag variable to determine if time = 0. If it is, then for each atom ‘i’ and its neighbor list ‘nearest[ ]’, we dynamically allocate enough memory (size of nearest[ ]) and store this neighbor list in the int** member variable.

Please, if there is any more efficient method already implemented that we can use to store the neighbor list at time = 0, please let me know.

unless i misunderstand what you are doing, i see a more serious problem here:

the indices stored in the neighbor list are “local” atom indices for “local” and “ghost” atoms. those change as atoms move from subdomain to subdomain or are sorted for better performance, and thus provide no consistent identification of atoms. if you want to store the identity of neighbors, you would have to store the atom “tags” of the neighbors and then use atom->map() to convert the atom index back to a “local” index (if actually still present in the same subdomain), possibly looping over “atom->sametag[]” to identify the most suitable periodic image, if the box is small and the cutoff large.

You can look at the src/PERI package and its fix_peri_neigh.cpp/h class.
In peridynamics (PD) the initial set of bonds is like a time=0 neighbor list.
Each atom carries around its list of bonds forever, until a bond breaks.

As Axel said, you will have to store the true ID of the neighbors, not
the local index. And you will have to account for the fact
that at some point in time, the time=0 neighbor may have moved
far enough away that it can no longer be found by the central atom.

Steve

Thank you Alex and Steve for your answers. Now we are using tag[i] to store the “real” id of the atoms. However I have a question: We can use atom-> map to convert the real ids to local ids as Axel suggested and then we do the computations. However, I believe that we can also convert the locals to real using tag[i]. Is there a specific reason that you suggest to do it in that way? I really appreciate your help.

Thank you.
-Spencer

Thank you Alex and Steve for your answers. Now we are using tag[i] to store the “real” id of the atoms. However I have a question: We can use atom-> map to convert the real ids to local ids as Axel suggested and then we do the computations. However, I believe that we can also convert the locals to real using tag[i]. Is there a specific reason that you suggest to do it in that way?

these are the inverse operation of each other. which one is appropriate depends on how you structure your algorithm. i assumed, you would be looping over all “tags” rather than local ids, and then you need the map function to determined which index has the tag (and don’t forget the “sametag” member, in case multiple periodic copies exist of the same atom as ghosts, and you have to take that into account as well).

axel.

Thank you Axel and Steve again,

I’ve continued working on the implementation of our compute (compute_nlcf), as was explained in our previous email, and also have started working on a fix (fix_storeNL) that will perform the creation and storing of the neighbor list of every atom at t=0. For the implementation of the fix_storeNL I am using the ‘fix_peri_neigh’ class as a starting point. First, to understand the PERI package, I’ve ran some examples to get more familiar with its functionality as it relates to our problem. The compute_nlcf and fix_storeNL work in this way: When the compute_nlcf is first called, it creates and initializes the fix_storeNL within its constructor. At t=0 the ‘setup’ function within fix_storeNL creates the NL for each atom and stores within some member variable (equivalent to ‘tagint **partner’ in fix_peri_neigh). Using the ‘compute_peratom’ function within compute_nlcf, I am able to retrieve and print out at different timesteps the NL that was created and stored by fix_storeNL at t=0. We want the NL created at t=0 to be unchanged as the simulation runs (contrary to PERI package where the pair style updates the NL when bonds break).

The code works well in serial, but it hangs in parallel when the atoms begin migrating to different processors. I’ve been tracking the problem for a while and I think it happens in either the unpack_exchange or pack_exchange functions (same name as in the PERI package). I have compared the behavior of these functions in fix_storeNL and in the fix_peri_neigh. One of the only symptoms that I can see is that even at t=0 (when there is no migration), these pack/unpack functions are called using compute_nlcf and fix_storeNL, which I believe is wrong because they are not called at t=0 when using the PERI package. We were able to notice this difference by tracking (output to file) at what time step each pack/unpack function was called for both our code and the PERI package.

My questions are: Is it right that pack/unpack are called when there is no migration? When are these two functions called or which class/function triggers their call?

Thanks for your help.

Spencer.