Updating logic for multi

src/OPENMP/npair_multi_omp.cpp

@@ -118,15 +118,16 @@ void NPairMultiOmp<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
       if (icollection == jcollection) jbin = ibin;
       else jbin = coord2bin(x[i], jcollection);
 
-      // loop over all atoms in surrounding bins in stencil including self
-      // skip i = j
-      // use full stencil for all collection combinations
-
       s = stencil_multi[icollection][jcollection];
       ns = nstencil_multi[icollection][jcollection];
 
       for (k = 0; k < ns; k++) {
         js = binhead_multi[jcollection][jbin + s[k]];
+
+        // own-bin for half stencil
+        if (HALF)
+          if (flag_half_multi[icollection][jcollection] && s[k] == 0) js = bins[i];
+
         for (j = js; j >= 0; j = bins[j]) {
           if (!HALF) {
             // Full neighbor list
@@ -144,40 +145,52 @@ void NPairMultiOmp<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
             // below = lower z or (equal z and lower y) or (equal zy and lower x)
             //         (equal zyx and j <= i)
             // latter excludes self-self interaction but allows superposed atoms
-            if (x[j][2] < ztmp) continue;
+
-            if (x[j][2] == ztmp) {
+            // if same size (same collection), use half stencil
-              if (x[j][1] < ytmp) continue;
+            if (flag_half_multi[icollection][jcollection]) {
-              if (x[j][1] == ytmp) {
+              if (x[j][2] < ztmp) continue;
-                if (x[j][0] < xtmp) continue;
+              if (x[j][2] == ztmp) {
-                if (x[j][0] == xtmp && j <= i) continue;
+                if (x[j][1] < ytmp) continue;
+                if (x[j][1] == ytmp) {
+                  if (x[j][0] < xtmp) continue;
+                  if (x[j][0] == xtmp && j <= i) continue;
+                }
               }
             }
           } else {
             // Half neighbor list, newton on, orthonormal
-            // store every pair for every bin in stencil,except for i's bin
+            // if same size: uses half stencil so includes a check of the central bin
+            if (flag_half_multi[icollection][jcollection]){
+              if (s[k] == 0) {
+                // if same collection,
+                //   if j is owned atom, store it, since j is beyond i in linked list
+                //   if j is ghost, only store if j coords are "above and to the right" of i
 
-            if (stencil[k] == 0) {
+                // if different collections,
-              // if j is owned atom, store it, since j is beyond i in linked list
+                //   if j is owned atom, store it if j > i
-              // if j is ghost, only store if j coords are "above and to the "right" of i
+                //   if j is ghost, only store if j coords are "above and to the right" of i
-              if (j >= nlocal) {
+
-                if (x[j][2] < ztmp) continue;
+                if ((icollection != jcollection) && (j < i)) continue;
-                if (x[j][2] == ztmp) {
+
-                  if (x[j][1] < ytmp) continue;
+                if (j >= nlocal) {
-                  if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+                  if (x[j][2] < ztmp) continue;
+                  if (x[j][2] == ztmp) {
+                    if (x[j][1] < ytmp) continue;
+                    if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+                  }
                 }
               }
             }
           }
 
           jtype = type[j];
-          if (exclude && exclusion(i,j,itype,jtype,mask,molecule))  continue;
+          if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
 
           delx = xtmp - x[j][0];
           dely = ytmp - x[j][1];
           delz = ztmp - x[j][2];
           rsq = delx*delx + dely*dely + delz*delz;
 
-
           if (SIZE) {
             radsum = radius[i] + radius[j];
             cut = radsum + skin;

src/npair.cpp

@@ -174,6 +174,8 @@ void NPair::copy_stencil_info()
 
   nstencil_multi = ns->nstencil_multi;
   stencil_multi = ns->stencil_multi;
+
+  flag_half_multi = ns->flag_half_multi;
 }
 
 /* ----------------------------------------------------------------------

src/npair.h

@@ -96,6 +96,7 @@ class NPair : protected Pointers {
   int *nstencil_multi_old;
   int **stencil_multi_old;
   double **distsq_multi_old;
+  bool **flag_half_multi;
 
   int **nstencil_multi;
   int ***stencil_multi;

src/npair_multi.cpp

@@ -112,6 +112,11 @@ void NPairMulti<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
 
       for (k = 0; k < ns; k++) {
         js = binhead_multi[jcollection][jbin + s[k]];
+
+        // own-bin for half stencil
+        if (HALF)
+          if (flag_half_multi[icollection][jcollection] && s[k] == 0) js = bins[i];
+
         for (j = js; j >= 0; j = bins[j]) {
           if (!HALF) {
             // Full neighbor list
@@ -129,26 +134,39 @@ void NPairMulti<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
             // below = lower z or (equal z and lower y) or (equal zy and lower x)
             //         (equal zyx and j <= i)
             // latter excludes self-self interaction but allows superposed atoms
-            if (x[j][2] < ztmp) continue;
+
-            if (x[j][2] == ztmp) {
+            // if same size (same collection), use half stencil
-              if (x[j][1] < ytmp) continue;
+            if (flag_half_multi[icollection][jcollection]) {
-              if (x[j][1] == ytmp) {
+              if (x[j][2] < ztmp) continue;
-                if (x[j][0] < xtmp) continue;
+              if (x[j][2] == ztmp) {
-                if (x[j][0] == xtmp && j <= i) continue;
+                if (x[j][1] < ytmp) continue;
+                if (x[j][1] == ytmp) {
+                  if (x[j][0] < xtmp) continue;
+                  if (x[j][0] == xtmp && j <= i) continue;
+                }
               }
             }
           } else {
             // Half neighbor list, newton on, orthonormal
-            // store every pair for every bin in stencil,except for i's bin
+            // if same size: uses half stencil so includes a check of the central bin
+            if (flag_half_multi[icollection][jcollection]) {
+              if (s[k] == 0) {
+                // if same collection,
+                //   if j is owned atom, store it, since j is beyond i in linked list
+                //   if j is ghost, only store if j coords are "above and to the right" of i
 
-            if (stencil[k] == 0) {
+                // if different collections,
-              // if j is owned atom, store it, since j is beyond i in linked list
+                //   if j is owned atom, store it if j > i
-              // if j is ghost, only store if j coords are "above and to the "right" of i
+                //   if j is ghost, only store if j coords are "above and to the right" of i
-              if (j >= nlocal) {
+
-                if (x[j][2] < ztmp) continue;
+                if ((icollection != jcollection) && (j < i)) continue;
-                if (x[j][2] == ztmp) {
+
-                  if (x[j][1] < ytmp) continue;
+                if (j >= nlocal) {
-                  if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+                  if (x[j][2] < ztmp) continue;
+                  if (x[j][2] == ztmp) {
+                    if (x[j][1] < ytmp) continue;
+                    if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+                  }
                 }
               }
             }

src/nstencil_multi.cpp

@@ -87,13 +87,13 @@ void NStencilMulti<HALF, DIM_3D, TRI>::create()
 
       bin_collection = bin_collection_multi[icollection][jcollection];
       cutsq = cutcollectionsq[icollection][jcollection];
-
       half_flag = flag_half_multi[icollection][jcollection];
 
       // Half and ortho stencils include central bin first
       // This preserves the historical order of the neighbor list
       // as the old npair classes used to separately parse the central bin first
-      if (half_flag && (!TRI)) stencil[nstencil++] = 0;
+      if (HALF && (!TRI))
+        if (half_flag) stencil_multi[icollection][jcollection][ns++] = 0;
 
       // For half stencils, only the upper plane is needed
       int sy_min = sy;
@@ -106,15 +106,16 @@ void NStencilMulti<HALF, DIM_3D, TRI>::create()
       for (k = -sz_min; k <= sz; k++) {
         for (j = -sy_min; j <= sy; j++) {
           for (i = -sx; i <= sx; i++) {
-
             // Now only include "upper right" bins for half and ortho stencils
-            if (HALF) {
+            if (HALF && (!TRI)) {
-              if (half_flag && (!DIM_3D) && (!TRI))
+              if (half_flag) {
-                if (! (j > 0 || (j == 0 && i > 0))) continue;
+                if (DIM_3D) {
-              if (half_flag && DIM_3D && (!TRI))
+                  if (! (k > 0 || j > 0 || (j == 0 && i > 0))) continue;
-                if (! (k > 0 || j > 0 || (j == 0 && i > 0))) continue;
+                } else {
+                  if (! (j > 0 || (j == 0 && i > 0))) continue;
+                }
+              }
             }
-
             if (bin_distance_multi(i,j,k,bin_collection) < cutsq)
               stencil_multi[icollection][jcollection][ns++] = k * mbiny * mbinx + j * mbinx + i;
           }