Creating nested neighbor lists in LAMMPS

LAMMPS LAMMPS Development
1

Hi,

I am trying to create two neighbor lists to iterate over: intermolecular (atoms not in the same molecule) and intramolecular (atoms in the same molecule). These neighbor lists are defined at a timestep = 0. These two lists will accessed every couple of timesteps during the simulation. I was wondering if there is a “LAMMPS-appropriate” way to construct these lists. Also, I have tried to browse all the various compute functions in LAMMPS, but I can not find a good example. Currently, I am using list<list> for the lists of list and list as a singular list. Please see below a code snippet where I construct these neighbor lists:

if (firstflag) {
nchunk = n; //this is a chunk compute
allocate(); //Memory allocation due a chunk compute
size_array_rows = nchunk;
if (nl_need) {
neighbor->build_one(list); //Building a neighbor list

  std::list<int> sgl_lst; //list construction

  inum = list->inum; // # of I atoms neighbors are stored for
  ilist = list->ilist; // local indices of I atoms
  numneigh = list->numneigh; //number of neighbors
  firstneigh = list->firstneigh; // ptr to 1st J int value of each I atom

  int **bondlist = neighbor->bondlist; // list of bonded neighbors
  int nbondlist = neighbor->nbondlist; //number of bonded neighbors



  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];  //picking an atom
    if (mask[i] & groupbit) {
      xtmp = x[i][0]; //positions
      ytmp = x[i][1];
      ztmp = x[i][2];
      jlist_er = firstneigh[i];
      jnum_er = numneigh[i];

      int ncount_ra = 0;
      for (jj_ra = 0; jj_ra < nbondlist; jj_ra++) {
        j_ra = bondlist[i][1];
        j_ra &= NEIGHMASK;
        sgl_lst.push_back(j_ra); //appending atom index to list
        ncount_ra++; 
        }

      n_intra_neigh->push_back(ncount_ra);
      ra_list.push_back(sgl_lst);

      std::unordered_set<int> set_neigh(sgl_lst.begin(), sgl_lst.end()); //defining a set 

      sgl_lst.clear(); //clearing list to prepare for intermolecular list

      int ncount_er = 0;
      for (jj = 0; jj < jnum_er; jj++) {
        j_er = jlist_er[jj];
        if (set_neigh.find(j) != set_neigh.end()) { //if atom index is not in intra_list then it is in inter_list
          j &= NEIGHMASK;
    
          delx = xtmp - x[j][0]; 
          dely = ytmp - x[j][1];
          delz = ztmp - x[j][2];
          rsq = delx * delx + dely * dely + delz * delz;
          if (rsq < cutsq ) { //radial cutoff to define intermolecular neighbor list
            sgl_lst.push_back(j_er);
            ncount_er++;
            }
          }
        }
      n_inter_neigh->push_back(ncount_er);
      er_list.push_back(sgl_lst);
      }
    }
  }

Also, in the header file:
std::list<std::list> ra_list, er_list;
std::vector *n_inter_neigh, *n_intra_neigh;

2

The conceptual issue here is that these neighbour lists are not actually nested – that is, atoms are not guaranteed to be “ordinary” neighbours simply because they are in the same molecule. Two examples:

  • your simulation could contain an arbitrarily long linear polymer.
  • less meaningfully, you could set every single molid in a simulation to 1 and, without any “molecular” potentials, it would return the exact same results.

In these cases it’s often easier for you to say what you want to do. Often when people articulate what they actually want to do it turns out that either (1) the feature is already present in LAMMPS, or (2) it is either not sensible or too difficult to do correctly in an MD simulation.