git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15329 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/neigh_list.cpp

@@ -57,6 +57,14 @@ NeighList::NeighList(LAMMPS *lmp) :
   listcopy = NULL;
   listskip = NULL;
 
+  // USER-DPD package
+  ndxAIR_ssa = NULL;
+  maxbin_ssa = 0;
+  bins_ssa = NULL;
+  maxhead_ssa = 0;
+  binhead_ssa = NULL;
+  gbinhead_ssa = NULL;
+
   maxstencil = ghostflag = 0;
   stencil = NULL;
   stencilxyz = NULL;
@@ -89,6 +97,12 @@ NeighList::~NeighList()
 
   if (maxstencil) memory->destroy(stencil);
   if (ghostflag) memory->destroy(stencilxyz);
+  if (ndxAIR_ssa) memory->sfree(ndxAIR_ssa);
+  if (maxbin_ssa) memory->destroy(bins_ssa);
+  if (maxhead_ssa) {
+    memory->destroy(binhead_ssa);
+    memory->destroy(gbinhead_ssa);
+  }
 
   if (maxstencil_multi) {
     for (int i = 1; i <= atom->ntypes; i++) {
@@ -148,6 +162,11 @@ void NeighList::grow(int nmax)
   if (dnum)
     firstdouble = (double **) memory->smalloc(maxatoms*sizeof(double *),
                                               "neighlist:firstdouble");
+  if (ssaflag) {
+    if (ndxAIR_ssa) memory->sfree(ndxAIR_ssa);
+    ndxAIR_ssa = (uint16_t (*)[8]) memory->smalloc(sizeof(uint16_t)*8*maxatoms,
+      "neighlist:ndxAIR_ssa");
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -229,6 +248,7 @@ void NeighList::print_attributes()
   printf("  %d = grow flag\n",growflag);
   printf("  %d = stencil flag\n",stencilflag);
   printf("  %d = ghost flag\n",ghostflag);
+  printf("  %d = ssa flag\n",ssaflag);
   printf("\n");
   printf("  %d = pair\n",rq->pair);
   printf("  %d = fix\n",rq->fix);
@@ -290,6 +310,12 @@ bigint NeighList::memory_usage()
 
   if (maxstencil) bytes += memory->usage(stencil,maxstencil);
   if (ghostflag) bytes += memory->usage(stencilxyz,maxstencil,3);
+  if (ndxAIR_ssa) bytes += sizeof(uint16_t) * 8 * maxatoms;
+  if (maxbin_ssa) bytes += memory->usage(bins_ssa,maxbin_ssa);
+  if (maxhead_ssa) {
+    bytes += memory->usage(binhead_ssa,maxhead_ssa);
+    bytes += memory->usage(gbinhead_ssa,maxhead_ssa);
+  }
 
   if (maxstencil_multi) {
     bytes += memory->usage(stencil_multi,atom->ntypes,maxstencil_multi);

src/neigh_list.h

@@ -65,6 +65,16 @@ class NeighList : protected Pointers {
   NeighList *listcopy;          // me = copy list, point to list I copy from
   NeighList *listskip;          // me = skip list, point to list I skip from
 
+  // USER-DPD package and Shardlow Splitting Algorithm (SSA) support
+
+  int ssaflag;               // 1 if the list has the ndxAIR_ssa array
+  uint16_t (*ndxAIR_ssa)[8]; // for each atom, last neighbor index of each AIR
+  int *bins_ssa;             // index of next atom in each bin
+  int maxbin_ssa;            // size of bins_ssa array
+  int *binhead_ssa;          // index of 1st local atom in each bin
+  int *gbinhead_ssa;         // index of 1st ghost atom in each bin
+  int maxhead_ssa;           // size of binhead_ssa and gbinhead_ssa arrays
+
   // stencils of bin indices for neighbor finding
 
   int maxstencil;                  // max size of stencil

src/neigh_shardlow.cpp

@@ -107,6 +107,20 @@ void Neighbor::stencil_half_bin_3d_ssa(NeighList *list,
   }
 }
 
+// space for static variable ssaAIRptr so it
+// can be used in qsort's compair function "cmp_ssaAIR()"
+static int *ssaAIRptr;
+
+static int cmp_ssaAIR(const void *iptr, const void *jptr)
+{
+  int i = *((int *) iptr);
+  int j = *((int *) jptr);
+  if (ssaAIRptr[i] < ssaAIRptr[j]) return -1;
+  if (ssaAIRptr[i] > ssaAIRptr[j]) return 1;
+  return 0;
+}
+
+
 /* ----------------------------------------------------------------------
    build half list from full list for use by Shardlow Spliting Algorithm
    pair stored once if i,j are both owned and i < j
@@ -138,6 +152,7 @@ void Neighbor::half_from_full_newton_ssa(NeighList *list)
   // loop over parent full list
 
   for (ii = 0; ii < inum_full; ii++) {
+    int AIRct[8] = { 0 };
     n = 0;
     neighptr = ipage->vget();
 
@@ -153,8 +168,10 @@ void Neighbor::half_from_full_newton_ssa(NeighList *list)
       j = joriginal & NEIGHMASK;
       if (j < nlocal) {
         if (i > j) continue;
+        ++(AIRct[0]);
       } else {
-        if (ssaAIR[j] <= 0) continue;
+        if (ssaAIR[j] < 2) continue; // skip ghost atoms not in AIR
+        ++(AIRct[ssaAIR[j] - 1]);
       }
       neighptr[n++] = joriginal;
     }
@@ -165,6 +182,16 @@ void Neighbor::half_from_full_newton_ssa(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+
+    // sort the locals+ghosts in the neighbor list by their ssaAIR number
+    ssaAIRptr = atom->ssaAIR;
+    qsort(&(neighptr[0]), n, sizeof(int), cmp_ssaAIR);
+
+    // Do a prefix sum on the counts to turn them into indexes.
+    list->ndxAIR_ssa[i][0] = AIRct[0];
+    for (int ndx = 1; ndx < 8; ++ndx) {
+      list->ndxAIR_ssa[i][ndx] = AIRct[ndx] + list->ndxAIR_ssa[i][ndx - 1];
+    }
   }
 
   list->inum = inum;
@@ -218,22 +245,22 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
     exclude ghost atoms that are not in the Active Interaction Regions (AIR)
 ------------------------------------------------------------------------- */
 
-    if (mbins > maxhead_ssa) {
+    if (mbins > list->maxhead_ssa) {
-      maxhead_ssa = mbins;
+      list->maxhead_ssa = mbins;
-      memory->destroy(gbinhead_ssa);
+      memory->destroy(list->gbinhead_ssa);
-      memory->destroy(binhead_ssa);
+      memory->destroy(list->binhead_ssa);
-      memory->create(binhead_ssa,maxhead_ssa,"binhead_ssa");
+      memory->create(list->binhead_ssa,list->maxhead_ssa,"binhead_ssa");
-      memory->create(gbinhead_ssa,maxhead_ssa,"gbinhead_ssa");
+      memory->create(list->gbinhead_ssa,list->maxhead_ssa,"gbinhead_ssa");
     }
     for (i = 0; i < mbins; i++) {
-      gbinhead_ssa[i] = -1;
+      list->gbinhead_ssa[i] = -1;
-      binhead_ssa[i] = -1;
+      list->binhead_ssa[i] = -1;
     }
 
-    if (maxbin > maxbin_ssa) {
+    if (maxbin > list->maxbin_ssa) {
-      maxbin_ssa = maxbin;
+      list->maxbin_ssa = maxbin;
-      memory->destroy(bins_ssa);
+      memory->destroy(list->bins_ssa);
-      memory->create(bins_ssa,maxbin_ssa,"bins_ssa");
+      memory->create(list->bins_ssa,list->maxbin_ssa,"bins_ssa");
     }
 
     // bin in reverse order so linked list will be in forward order
@@ -241,26 +268,26 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
     if (includegroup) {
       int bitmask = group->bitmask[includegroup];
       for (i = nall-1; i >= nlocal; i--) {
-        if (ssaAIR[i] <= 0) continue; // skip ghost atoms not in AIR
+        if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
         if (mask[i] & bitmask) {
           ibin = coord2bin(x[i]);
-          bins_ssa[i] = gbinhead_ssa[ibin];
+          list->bins_ssa[i] = list->gbinhead_ssa[ibin];
-          gbinhead_ssa[ibin] = i;
+          list->gbinhead_ssa[ibin] = i;
         }
       }
       nlocal = atom->nfirst; // This is important for the code that follows!
     } else {
       for (i = nall-1; i >= nlocal; i--) {
-        if (ssaAIR[i] <= 0) continue; // skip ghost atoms not in AIR
+        if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
         ibin = coord2bin(x[i]);
-        bins_ssa[i] = gbinhead_ssa[ibin];
+        list->bins_ssa[i] = list->gbinhead_ssa[ibin];
-        gbinhead_ssa[ibin] = i;
+        list->gbinhead_ssa[ibin] = i;
       }
     }
     for (i = nlocal-1; i >= 0; i--) {
       ibin = coord2bin(x[i]);
-      bins_ssa[i] = binhead_ssa[ibin];
+      list->bins_ssa[i] = list->binhead_ssa[ibin];
-      binhead_ssa[ibin] = i;
+      list->binhead_ssa[ibin] = i;
     }
   } /* else reuse previous binning. See Neighbor::build_one comment. */
 
@@ -269,6 +296,7 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
   // loop over owned atoms, storing half of the neighbors
 
   for (i = 0; i < nlocal; i++) {
+    int AIRct[8] = { 0 };
     n = 0;
     neighptr = ipage->vget();
 
@@ -285,7 +313,7 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
     // loop over rest of local atoms in i's bin
     // just store them, since j is beyond i in linked list
 
-    for (j = bins_ssa[i]; j >= 0; j = bins_ssa[j]) {
+    for (j = list->bins_ssa[i]; j >= 0; j = list->bins_ssa[j]) {
 
       jtype = type[j];
       if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
@@ -316,7 +344,7 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
 
     // loop over all local atoms in other bins in "half" stencil
     for (k = 0; k < nstencil; k++) {
-      for (j = binhead_ssa[ibin+stencil[k]]; j >= 0; j = bins_ssa[j]) {
+      for (j = list->binhead_ssa[ibin+stencil[k]]; j >= 0; j = list->bins_ssa[j]) {
 
         jtype = type[j];
         if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
@@ -343,13 +371,15 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
         }
       }
     }
+    AIRct[0] = n;
 
     // loop over AIR ghost atoms in all bins in "full" stencil
     // Note: the non-AIR ghost atoms have already been filtered out
     // That is a significant time savings because of the "full" stencil
-    for (k = 0; k < maxstencil; k++) {
+    // Note2: only non-pure locals can have ghosts as neighbors
+    if (ssaAIR[i] == 1) for (k = 0; k < maxstencil; k++) {
       if (stencil[k] > mbins) break; /* Check if ghost stencil bins are exhausted */
-      for (j = gbinhead_ssa[ibin+stencil[k]]; j >= 0; j = bins_ssa[j]) {
+      for (j = list->gbinhead_ssa[ibin+stencil[k]]; j >= 0; j = list->bins_ssa[j]) {
 
         jtype = type[j];
         if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
@@ -368,11 +398,20 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
                                    onemols[imol]->nspecial[iatom],
                                    tag[j]-tagprev);
             else which = 0;
-            if (which == 0) neighptr[n++] = j;
+            if (which == 0) {
-            else if (domain->minimum_image_check(delx,dely,delz))
               neighptr[n++] = j;
-            else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
+              ++(AIRct[ssaAIR[j] - 1]);
-          } else neighptr[n++] = j;
+            } else if (domain->minimum_image_check(delx,dely,delz)) {
+              neighptr[n++] = j;
+              ++(AIRct[ssaAIR[j] - 1]);
+            } else if (which > 0) {
+              neighptr[n++] = j ^ (which << SBBITS);
+              ++(AIRct[ssaAIR[j] - 1]);
+            }
+          } else {
+            neighptr[n++] = j;
+            ++(AIRct[ssaAIR[j] - 1]);
+          }
         }
       }
     }
@@ -383,6 +422,16 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+
+    // sort the ghosts in the neighbor list by their ssaAIR number
+    ssaAIRptr = atom->ssaAIR;
+    qsort(&(neighptr[AIRct[0]]), n - AIRct[0], sizeof(int), cmp_ssaAIR);
+
+    // Do a prefix sum on the counts to turn them into indexes.
+    list->ndxAIR_ssa[i][0] = AIRct[0];
+    for (int ndx = 1; ndx < 8; ++ndx) {
+      list->ndxAIR_ssa[i][ndx] = AIRct[ndx] + list->ndxAIR_ssa[i][ndx - 1];
+    }
   }
 
   list->inum = inum;

src/neighbor.cpp

@@ -104,14 +104,6 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
   maxbin = 0;
   bins = NULL;
 
-  // USER-DPD SSA AIR binning
-
-  maxbin_ssa = 0;
-  bins_ssa = NULL;
-  binhead_ssa = NULL;
-  gbinhead_ssa = NULL;
-  maxhead_ssa = 0;
-
   // pair exclusion list info
 
   includegroup = 0;
@@ -186,10 +178,6 @@ Neighbor::~Neighbor()
   memory->destroy(binhead);
   memory->destroy(bins);
 
-  memory->destroy(gbinhead_ssa);
-  memory->destroy(binhead_ssa);
-  memory->destroy(bins_ssa);
-
   memory->destroy(ex1_type);
   memory->destroy(ex2_type);
   memory->destroy(ex_type);
@@ -391,15 +379,6 @@ void Neighbor::init()
     maxbin = maxhead = 0;
     binhead = NULL;
     bins = NULL;
-
-    // for USER-DPD Shardlow Splitting Algorithm (SSA)
-    memory->destroy(bins_ssa);
-    memory->destroy(binhead_ssa);
-    memory->destroy(gbinhead_ssa);
-    maxbin_ssa = maxhead_ssa = 0;
-    bins_ssa = NULL;
-    binhead_ssa = NULL;
-    gbinhead_ssa = NULL;
   }
 
   // 1st time allocation of xhold and bins
@@ -724,6 +703,9 @@ void Neighbor::init()
         lists[i]->ghostflag = 0;
         if (requests[i]->ghost) lists[i]->ghostflag = 1;
         if (requests[i]->ghost && !requests[i]->occasional) anyghostlist = 1;
+
+        lists[i]->ssaflag = 0;
+        if (requests[i]->ssa) lists[i]->ssaflag = 1;
       } else init_list_flags1_kokkos(i);
     }
 
@@ -2254,9 +2236,6 @@ bigint Neighbor::memory_usage()
   if (style != NSQ) {
     bytes += memory->usage(bins,maxbin);
     bytes += memory->usage(binhead,maxhead);
-    bytes += memory->usage(bins_ssa,maxbin_ssa);
-    bytes += memory->usage(binhead_ssa,maxhead_ssa);
-    bytes += memory->usage(gbinhead_ssa,maxhead_ssa);
   }
 
   for (int i = 0; i < nrequest; i++)

src/neighbor.h

@@ -176,14 +176,6 @@ class Neighbor : protected Pointers {
 
   double *zeroes;              // vector of zeroes for shear history init
 
-  // USER-DPD package
-
-  int *bins_ssa;             // ptr to next atom in each bin used by SSA
-  int maxbin_ssa;            // size of bins array used by SSA
-  int *binhead_ssa;          // ptr to 1st atom in each bin used by SSA
-  int *gbinhead_ssa;         // ptr to 1st ghost atom in each bin used by SSA
-  int maxhead_ssa;           // size of binhead array used by SSA
-
   // methods
 
   void bin_atoms();                     // bin all atoms