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

src/neigh_full.cpp

@@ -163,6 +163,7 @@ void Neighbor::full_nsq_ghost(NeighList *list)
 
     // loop over all atoms, owned and ghost
     // skip i = j
+    // no molecular test when i = ghost atom
 
     if (i < nlocal) {
       for (j = 0; j < nall; j++) {
@@ -195,8 +196,7 @@ void Neighbor::full_nsq_ghost(NeighList *list)
         delz = ztmp - x[j][2];
         rsq = delx*delx + dely*dely + delz*delz;
 
-        if (rsq <= cutneighghostsq[itype][jtype])
+        if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j;
-          neighptr[n++] = j;
       }
     }
 
@@ -371,6 +371,7 @@ void Neighbor::full_bin_ghost(NeighList *list)
     // loop over all atoms in surrounding bins in stencil including self
     // when i is a ghost atom, must check if stencil bin is out of bounds
     // skip i = j
+    // no molecular test when i = ghost atom
 
     if (i < nlocal) {
       ibin = coord2bin(x[i]);
@@ -418,8 +419,7 @@ void Neighbor::full_bin_ghost(NeighList *list)
           delz = ztmp - x[j][2];
           rsq = delx*delx + dely*dely + delz*delz;
 
-          if (rsq <= cutneighghostsq[itype][jtype])
+          if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j;
-            neighptr[n++] = j;
         }
       }
     }

src/neigh_half_bin.cpp

@@ -187,6 +187,7 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list)
     // stores own/own pairs only once
     // stores own/ghost pairs with owned atom only, on both procs
     // stores ghost/ghost pairs only once
+    // no molecular test when i = ghost atom
 
     if (i < nlocal) {
       ibin = coord2bin(x[i]);
@@ -235,8 +236,7 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list)
           delz = ztmp - x[j][2];
           rsq = delx*delx + dely*dely + delz*delz;
 
-          if (rsq <= cutneighghostsq[itype][jtype])
+          if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j;
-            neighptr[n++] = j;
         }
       }
     }

src/neigh_half_nsq.cpp

@@ -32,8 +32,6 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr;
 
-  // loop over each atom, storing neighbors
-
   int **special = atom->special;
   int **nspecial = atom->nspecial;
   int *tag = atom->tag;
@@ -59,6 +57,8 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
   int npage = 0;
   int npnt = 0;
 
+  // loop over owned atoms, storing neighbors
+
   for (i = 0; i < nlocal; i++) {
 
     if (pgsize - npnt < oneatom) {
@@ -76,6 +76,7 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
     ztmp = x[i][2];
 
     // loop over remaining atoms, owned and ghost
+    // only store pair if i < j
 
     for (j = i+1; j < nall; j++) {
       if (includegroup && !(mask[j] & bitmask)) continue;
@@ -109,6 +110,119 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
   list->inum = inum;
 }
 
+/* ----------------------------------------------------------------------
+   N^2 / 2 search for neighbor pairs with partial Newton's 3rd law
+   include neighbors of ghost atoms, but no "special neighbors" for ghosts
+   pair stored once if i,j are both owned and i < j
+   pair stored by me if i owned and j ghost (also stored by proc owning j)
+   pair stored once if i,j are both ghost and i < j
+------------------------------------------------------------------------- */
+
+void Neighbor::half_nsq_no_newton_ghost(NeighList *list)
+{
+  int i,j,n,itype,jtype,which,bitmask;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *neighptr;
+
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
+  double **x = atom->x;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+  int molecular = atom->molecular;
+  if (includegroup) {
+    nlocal = atom->nfirst;
+    bitmask = group->bitmask[includegroup];
+  }
+
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+  int **pages = list->pages;
+
+  int inum = 0;
+  int npage = 0;
+  int npnt = 0;
+
+  // loop over owned & ghost atoms, storing neighbors
+
+  for (i = 0; i < nall; i++) {
+
+    if (pgsize - npnt < oneatom) {
+      npnt = 0;
+      npage++;
+      if (npage == list->maxpage) pages = list->add_pages();
+    }
+
+    neighptr = &pages[npage][npnt];
+    n = 0;
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+
+    // loop over remaining atoms, owned and ghost
+    // only store pair if i < j
+    // stores own/own pairs only once
+    // stores own/ghost pairs with owned atom only, on both procs
+    // stores ghost/ghost pairs only once
+    // no molecular test when i = ghost atom
+
+    if (i < nlocal) {
+      for (j = i+1; j < nall; j++) {
+        if (includegroup && !(mask[j] & bitmask)) continue;
+        jtype = type[j];
+        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 (rsq <= cutneighsq[itype][jtype]) {
+          if (molecular) {
+            which = find_special(special[i],nspecial[i],tag[j]);
+            if (which == 0) neighptr[n++] = j;
+            else if (domain->minimum_image_check(delx,dely,delz))
+              neighptr[n++] = j;
+            else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
+          } else neighptr[n++] = j;
+        }
+      }
+
+    } else {
+      for (j = i+1; j < nall; j++) {
+        if (includegroup && !(mask[j] & bitmask)) continue;
+        jtype = type[j];
+        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 (rsq <= cutneighsq[itype][jtype]) neighptr[n++] = j;
+      }
+    }
+    
+    ilist[inum++] = i;
+    firstneigh[i] = neighptr;
+    numneigh[i] = n;
+    npnt += n;
+    if (n > oneatom)
+      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+  }
+
+  list->inum = atom->nlocal;
+  list->gnum = inum - atom->nlocal;
+}
+
 /* ----------------------------------------------------------------------
    N^2 / 2 search for neighbor pairs with full Newton's 3rd law
    every pair stored exactly once by some processor

src/neigh_stencil.cpp

@@ -61,6 +61,28 @@ void Neighbor::stencil_half_bin_2d_no_newton(NeighList *list,
 
 /* ---------------------------------------------------------------------- */
 
+void Neighbor::stencil_half_ghost_bin_2d_no_newton(NeighList *list,
+                                                   int sx, int sy, int sz)
+{
+  int i,j,k;
+  int *stencil = list->stencil;
+  int **stencilxyz = list->stencilxyz;
+  int nstencil = 0;
+
+  for (j = -sy; j <= sy; j++)
+    for (i = -sx; i <= sx; i++)
+      if (bin_distance(i,j,0) < cutneighmaxsq) {
+        stencilxyz[nstencil][0] = i;
+        stencilxyz[nstencil][1] = j;
+        stencilxyz[nstencil][2] = 0;
+        stencil[nstencil++] = j*mbinx + i;
+      }
+
+  list->nstencil = nstencil;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void Neighbor::stencil_half_bin_3d_no_newton(NeighList *list,
                                              int sx, int sy, int sz)
 {

src/neighbor.cpp

@@ -826,24 +826,48 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
     } else if (rq->half) {
       if (style == NSQ) {
         if (rq->newton == 0) {
-          if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton;
+          if (newton_pair == 0) {
-          else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton;
+            if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton;
+            else if (includegroup)
+              error->all(FLERR,"Neighbor include group not allowed "
+                         "with ghost neighbors");
+            else pb = &Neighbor::half_nsq_no_newton_ghost;
+          } else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton;
         } else if (rq->newton == 1) {
           pb = &Neighbor::half_nsq_newton;
-        } else if (rq->newton == 2) pb = &Neighbor::half_nsq_no_newton;
+        } else if (rq->newton == 2) {
+          if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton;
+          else if (includegroup)
+            error->all(FLERR,"Neighbor include group not allowed "
+                       "with ghost neighbors");
+          else pb = &Neighbor::half_nsq_no_newton_ghost;
+        }
       } else if (style == BIN) {
         if (rq->newton == 0) {
-          if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton;
+          if (newton_pair == 0) {
-          else if (triclinic == 0) pb = &Neighbor::half_bin_newton;
+            if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton;
-          else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri;
+            else if (includegroup)
+              error->all(FLERR,"Neighbor include group not allowed "
+                         "with ghost neighbors");
+            else pb = &Neighbor::half_bin_no_newton_ghost;
+          } else if (triclinic == 0) {
+            pb = &Neighbor::half_bin_newton;
+          } else if (triclinic == 1) 
+            pb = &Neighbor::half_bin_newton_tri;
         } else if (rq->newton == 1) {
           if (triclinic == 0) pb = &Neighbor::half_bin_newton;
           else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri;
         } else if (rq->newton == 2) {
-          if (rq->ghost) pb = &Neighbor::half_bin_no_newton_ghost;
+          if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton;
-          else pb = &Neighbor::half_bin_no_newton;
+          else if (includegroup)
+            error->all(FLERR,"Neighbor include group not allowed "
+                       "with ghost neighbors");
+          else pb = &Neighbor::half_bin_no_newton_ghost;
         }
       } else if (style == MULTI) {
+        if (rq->ghost == 1)
+          error->all(FLERR,
+                     "Neighbor multi not yet enabled for ghost neighbors");
         if (rq->newton == 0) {
           if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton;
           else if (triclinic == 0) pb = &Neighbor::half_multi_newton;
@@ -860,17 +884,18 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
         else if (includegroup)
           error->all(FLERR,
                      "Neighbor include group not allowed with ghost neighbors");
-        else if (rq->ghost == 1) pb = &Neighbor::full_nsq_ghost;
+        else pb = &Neighbor::full_nsq_ghost;
       } else if (style == BIN) {
         if (rq->ghost == 0) pb = &Neighbor::full_bin;
         else if (includegroup)
           error->all(FLERR,
                      "Neighbor include group not allowed with ghost neighbors");
-        else if (rq->ghost == 1) pb = &Neighbor::full_bin_ghost;
+        else pb = &Neighbor::full_bin_ghost;
       } else if (style == MULTI) {
-        if (rq->ghost == 0) pb = &Neighbor::full_multi;
+        if (rq->ghost == 1)
-        else error->all(FLERR,
+          error->all(FLERR,
                      "Neighbor multi not yet enabled for ghost neighbors");
+        pb = &Neighbor::full_multi;
       }
 
     } else if (rq->gran) {
@@ -895,6 +920,9 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
       } else if (style == MULTI)
         error->all(FLERR,"Neighbor multi not yet enabled for rRESPA");
     }
+
+  // OMP versions of build methods
+
   } else {
 
     if (rq->copy) pb = &Neighbor::copy_from;
@@ -912,23 +940,52 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
     } else if (rq->half) {
       if (style == NSQ) {
         if (rq->newton == 0) {
-          if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton_omp;
+          if (newton_pair == 0) {
-          else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton_omp;
+            if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton_omp;
+            else if (includegroup)
+              error->all(FLERR,"Neighbor include group not allowed "
+                         "with ghost neighbors");
+            // NOTE: missing OMP version of this one
+            else pb = &Neighbor::half_nsq_no_newton_ghost;
+          } else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton_omp;
         } else if (rq->newton == 1) {
           pb = &Neighbor::half_nsq_newton_omp;
         } else if (rq->newton == 2) {
-          pb = &Neighbor::half_nsq_no_newton_omp;
+          if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton_omp;
+          else if (includegroup)
+            error->all(FLERR,"Neighbor include group not allowed "
+                       "with ghost neighbors");
+          // NOTE: missing OMP version of this one
+          else pb = &Neighbor::half_nsq_no_newton_ghost;
         }
       } else if (style == BIN) {
         if (rq->newton == 0) {
-          if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton_omp;
+          if (newton_pair == 0) {
-          else if (triclinic == 0) pb = &Neighbor::half_bin_newton_omp;
+            if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton_omp;
-          else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri_omp;
+            else if (includegroup)
+              error->all(FLERR,"Neighbor include group not allowed "
+                         "with ghost neighbors");
+            // NOTE: missing OMP version of this one
+            else pb = &Neighbor::half_bin_no_newton_ghost;
+          } else if (triclinic == 0) {
+            pb = &Neighbor::half_bin_newton_omp;
+          } else if (triclinic == 1) 
+            pb = &Neighbor::half_bin_newton_tri_omp;
         } else if (rq->newton == 1) {
           if (triclinic == 0) pb = &Neighbor::half_bin_newton_omp;
           else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri_omp;
-        } else if (rq->newton == 2) pb = &Neighbor::half_bin_no_newton_omp;
+        } else if (rq->newton == 2) {
+          if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton_omp;
+          else if (includegroup)
+            error->all(FLERR,"Neighbor include group not allowed "
+                       "with ghost neighbors");
+          // NOTE: missing OMP version of this one
+          else pb = &Neighbor::half_bin_no_newton_ghost;
+        }
       } else if (style == MULTI) {
+        if (rq->ghost == 1)
+          error->all(FLERR,
+                     "Neighbor multi not yet enabled for ghost neighbors");
         if (rq->newton == 0) {
           if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton_omp;
           else if (triclinic == 0) pb = &Neighbor::half_multi_newton_omp;
@@ -945,17 +1002,18 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
         else if (includegroup)
           error->all(FLERR,
                      "Neighbor include group not allowed with ghost neighbors");
-        else if (rq->ghost == 1) pb = &Neighbor::full_nsq_ghost_omp;
+        else pb = &Neighbor::full_nsq_ghost_omp;
       } else if (style == BIN) {
         if (rq->ghost == 0) pb = &Neighbor::full_bin_omp;
         else if (includegroup)
           error->all(FLERR,
                      "Neighbor include group not allowed with ghost neighbors");
-        else if (rq->ghost == 1) pb = &Neighbor::full_bin_ghost_omp;
+        else pb = &Neighbor::full_bin_ghost_omp;
       } else if (style == MULTI) {
-        if (rq->ghost == 0) pb = &Neighbor::full_multi_omp;
+        if (rq->ghost == 1)
-        else error->all(FLERR,
+          error->all(FLERR,
                      "Neighbor multi not yet enabled for ghost neighbors");
+        pb = &Neighbor::full_multi_omp;
       }
 
     } else if (rq->gran) {
@@ -1004,10 +1062,13 @@ void Neighbor::choose_stencil(int index, NeighRequest *rq)
     if (style == BIN) {
       if (rq->newton == 0) {
         if (newton_pair == 0) {
-          if (dimension == 2)
+          if (dimension == 2) {
-            sc = &Neighbor::stencil_half_bin_2d_no_newton;
+            if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_2d_no_newton;
-          else if (dimension == 3)
+            else sc = &Neighbor::stencil_half_bin_2d_no_newton;
-            sc = &Neighbor::stencil_half_bin_3d_no_newton;
+          } else if (dimension == 3) {
+            if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_3d_no_newton;
+            else sc = &Neighbor::stencil_half_bin_3d_no_newton;
+          }
         } else if (triclinic == 0) {
           if (dimension == 2)
             sc = &Neighbor::stencil_half_bin_2d_newton;
@@ -1033,7 +1094,8 @@ void Neighbor::choose_stencil(int index, NeighRequest *rq)
         }
       } else if (rq->newton == 2) {
         if (dimension == 2)
-          sc = &Neighbor::stencil_half_bin_2d_no_newton;
+          if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_2d_no_newton;
+          else sc = &Neighbor::stencil_half_bin_2d_no_newton;
         else if (dimension == 3) {
           if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_3d_no_newton;
           else sc = &Neighbor::stencil_half_bin_3d_no_newton;

src/neighbor.h

@@ -177,6 +177,7 @@ class Neighbor : protected Pointers {
   PairPtr *pair_build;
 
   void half_nsq_no_newton(class NeighList *);
+  void half_nsq_no_newton_ghost(class NeighList *);
   void half_nsq_newton(class NeighList *);
 
   void half_bin_no_newton(class NeighList *);
@@ -223,6 +224,7 @@ class Neighbor : protected Pointers {
   StencilPtr *stencil_create;
 
   void stencil_half_bin_2d_no_newton(class NeighList *, int, int, int);
+  void stencil_half_ghost_bin_2d_no_newton(class NeighList *, int, int, int);
   void stencil_half_bin_3d_no_newton(class NeighList *, int, int, int);
   void stencil_half_ghost_bin_3d_no_newton(class NeighList *, int, int, int);
   void stencil_half_bin_2d_newton(class NeighList *, int, int, int);