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

src/COLLOID/pair_colloid.cpp

@@ -348,7 +348,6 @@ double PairColloid::init_one(int i, int j)
   sigma6[j][i] = sigma6[i][j];
   diameter[j][i] = diameter[i][j];
   cut[j][i] = cut[i][j];
-  cutsq[j][i] = cutsq[i][j] = cut[i][j] * cut[i][j];
 
   double epsilon = a12[i][j]/144.0;
   lj1[j][i] = lj1[i][j] = 48.0 * epsilon * sigma6[i][j] * sigma6[i][j];
@@ -359,7 +358,8 @@ double PairColloid::init_one(int i, int j)
   offset[j][i] = offset[i][j] = 0.0;
   if (offset_flag) {
     double tmp;
-    offset[j][i] = offset[i][j] = single(0,0,i,j,cutsq[i][j],0.0,1.0,tmp);
+    offset[j][i] = offset[i][j] = 
+      single(0,0,i,j,cut[i][j]*cut[i][j],0.0,1.0,tmp);
   }
 
   return cut[i][j];

src/MANYBODY/pair_airebo.cpp

@@ -48,14 +48,17 @@ PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp)
 {
   single_enable = 0;
   one_coeff = 1;
+  ghostneigh = 1;
 
   maxlocal = 0;
   REBO_numneigh = NULL;
   REBO_firstneigh = NULL;
-  closestdistsq = NULL;
   maxpage = 0;
   pages = NULL;
   nC = nH = NULL;
+
+  time1 = 0.0;
+  time2 = 0.0;
 }
 
 /* ----------------------------------------------------------------------
@@ -68,13 +71,13 @@ PairAIREBO::~PairAIREBO()
   memory->sfree(REBO_firstneigh);
   for (int i = 0; i < maxpage; i++) memory->sfree(pages[i]);
   memory->sfree(pages);
-  memory->sfree(closestdistsq);
   memory->sfree(nC);
   memory->sfree(nH);
 
   if (allocated) {
     memory->destroy_2d_int_array(setflag);
     memory->destroy_2d_double_array(cutsq);
+    memory->destroy_2d_double_array(cutghost);
 
     memory->destroy_2d_double_array(cutljsq);
     memory->destroy_2d_double_array(lj1);
@@ -92,12 +95,18 @@ void PairAIREBO::compute(int eflag, int vflag)
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = vflag_atom = 0;
 
+  double tstart = MPI_Wtime();
   REBO_neigh();
+  double tmid = MPI_Wtime();
   FREBO(eflag,vflag);
   if (ljflag) FLJ(eflag,vflag);
   if (torflag) TORSION(eflag,vflag);
   
   if (vflag_fdotr) virial_compute();
+  double tstop = MPI_Wtime();
+
+  time1 += tmid-tstart;
+  time2 += tstop-tmid;
 }
 
 /* ----------------------------------------------------------------------
@@ -115,6 +124,7 @@ void PairAIREBO::allocate()
       setflag[i][j] = 0;
 
   cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
+  cutghost = memory->create_2d_double_array(n+1,n+1,"pair:cutghost");
 
   // only sized by C,H = 2 types
 
@@ -210,11 +220,12 @@ void PairAIREBO::init_style()
   if (force->newton_pair == 0)
     error->all("Pair style AIREBO requires newton pair on");
 
-  // need a full neighbor list
+  // need a full neighbor list, including neighbors of ghosts
 
   int irequest = neighbor->request(this);
   neighbor->requests[irequest]->half = 0;
   neighbor->requests[irequest]->full = 1;
+  neighbor->requests[irequest]->ghost = 1;
 
   // local REBO neighbor list memory
 
@@ -261,6 +272,7 @@ double PairAIREBO::init_one(int i, int j)
   //   rcLJmax + 2*rcmax, since I-J < rcLJmax and J-L,L-N = REBO distances
   //   long interaction = I-J with I = owned and J = ghost
   //   cutlj*sigma, since I-J < LJ cutoff
+  // cutghost = REBO cutoff used in REBO_neigh() for neighbors of ghosts
 
   double cutmax = cut3rebo;
   if (ljflag) {
@@ -268,14 +280,14 @@ double PairAIREBO::init_one(int i, int j)
     cutmax = MAX(cutmax,cutlj*sigma[0][0]);
   }
 
-  cutsq[i][j] = cutmax * cutmax;
+  cutghost[i][j] = rcmax[ii][jj];
   cutljsq[ii][jj] = cutlj*sigma[ii][jj] * cutlj*sigma[ii][jj];
   lj1[ii][jj] = 48.0 * epsilon[ii][jj] * pow(sigma[ii][jj],12.0);
   lj2[ii][jj] = 24.0 * epsilon[ii][jj] * pow(sigma[ii][jj],6.0);
   lj3[ii][jj] = 4.0 * epsilon[ii][jj] * pow(sigma[ii][jj],12.0);
   lj4[ii][jj] = 4.0 * epsilon[ii][jj] * pow(sigma[ii][jj],6.0);
 
-  cutsq[j][i] = cutsq[i][j];
+  cutghost[j][i] = cutghost[i][j];
   cutljsq[jj][ii] = cutljsq[ii][jj];
   lj1[jj][ii] = lj1[ii][jj];
   lj2[jj][ii] = lj2[ii][jj];
@@ -292,7 +304,7 @@ double PairAIREBO::init_one(int i, int j)
 
 void PairAIREBO::REBO_neigh()
 {
-  int i,j,ii,jj,m,n,inum,jnum,itype,jtype;
+  int i,j,ii,jj,m,n,allnum,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS;
   int *ilist,*jlist,*numneigh,**firstneigh;
   int *neighptr;
@@ -306,41 +318,28 @@ void PairAIREBO::REBO_neigh()
     maxlocal = atom->nmax;
     memory->sfree(REBO_numneigh);
     memory->sfree(REBO_firstneigh);
-    memory->sfree(closestdistsq);
     memory->sfree(nC);
     memory->sfree(nH);
     REBO_numneigh = (int *)
       memory->smalloc(maxlocal*sizeof(int),"AIREBO:numneigh");
     REBO_firstneigh = (int **)
       memory->smalloc(maxlocal*sizeof(int *),"AIREBO:firstneigh");
-    closestdistsq = (double *) 
-      memory->smalloc(maxlocal*sizeof(double),"AIREBO:closestdistsq");
     nC = (double *) memory->smalloc(maxlocal*sizeof(double),"AIREBO:nC");
     nH = (double *) memory->smalloc(maxlocal*sizeof(double),"AIREBO:nH");
   }
   
-  inum = list->inum;
+  allnum = list->inum + list->gnum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  // initialize ghost atom references to -1 and closest distance = cutljrebosq
+  // store all REBO neighs of owned and ghost atoms
-
-  for (i = nlocal; i < nall; i++) {
-    REBO_numneigh[i] = -1;
-    closestdistsq[i] = cutljrebosq;
-  }
-
-  // store all REBO neighs of owned atoms
   // scan full neighbor list of I
-  // if J is ghost and within LJ cutoff:
-  //   flag it via REBO_numneigh so its REBO neighbors will be stored below
-  //   REBO requires neighbors of neighbors of i,j in each i,j LJ interaction
 
   npage = 0;
   int npnt = 0;
 
-  for (ii = 0; ii < inum; ii++) {
+  for (ii = 0; ii < allnum; ii++) {
     i = ilist[ii];
 
     if (pgsize - npnt < oneatom) {
@@ -367,82 +366,6 @@ void PairAIREBO::REBO_neigh()
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
 
-      if (rsq < rcmaxsq[itype][jtype]) {
-	neighptr[n++] = j;
-	if (jtype == 0)
-	  nC[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS);
-	else
-	  nH[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS);
-      }
-      if (j >= nlocal && rsq < closestdistsq[j]) {
-	REBO_numneigh[j] = i;
-	closestdistsq[j] = rsq;
-      }
-    }
-
-    REBO_firstneigh[i] = neighptr;
-    REBO_numneigh[i] = n;
-    npnt += n;
-    if (npnt >= pgsize)
-      error->one("Neighbor list overflow, boost neigh_modify one or page");
-  }
-
-  // store REBO neighs of ghost atoms I that have been flagged in REBO_numneigh
-  // find by scanning full neighbor list of owned atom M = closest neigh of I
-
-  for (i = nlocal; i < nall; i++) {
-
-    if (pgsize - npnt < oneatom) {
-      npnt = 0;
-      npage++;
-      if (npage == maxpage) add_pages(npage);
-    }
-    neighptr = &pages[npage][npnt];
-    n = 0;
-
-    xtmp = x[i][0];
-    ytmp = x[i][1];
-    ztmp = x[i][2];
-    itype = map[type[i]];
-    nC[i] = nH[i] = 0.0;
-    m = REBO_numneigh[i];
-    if (m < 0) {
-      REBO_firstneigh[i] = neighptr;
-      REBO_numneigh[i] = 0;
-      continue;
-    }
-
-    jtype = map[type[m]];
-    delx = xtmp - x[m][0];
-    dely = ytmp - x[m][1];
-    delz = ztmp - x[m][2];
-    rsq = delx*delx + dely*dely + delz*delz;
-
-    // add M as neigh of I if close enough
-
-    if (rsq < rcmaxsq[itype][jtype]) {
-      neighptr[n++] = m;
-      if (jtype == 0)
-	nC[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS);
-      else
-	nH[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS);
-    }
-
-    jlist = firstneigh[m];
-    jnum = numneigh[m];
-
-    // add M's neighbors as neighs of I if close enough
-    // skip I itself
-
-    for (jj = 0; jj < jnum; jj++) {
-      j = jlist[jj];
-      if (j == i) continue;
-      jtype = map[type[j]];
-      delx = xtmp - x[j][0];
-      dely = ytmp - x[j][1];
-      delz = ztmp - x[j][2];
-      rsq = delx*delx + dely*dely + delz*delz;
-
       if (rsq < rcmaxsq[itype][jtype]) {
 	neighptr[n++] = j;
 	if (jtype == 0)

src/MANYBODY/pair_airebo.h

@@ -36,6 +36,8 @@ class PairAIREBO : public Pair {
   double memory_usage();
 
  private:
+  double time1,time2;
+
   int me;
   int ljflag,torflag;              // 0/1 if LJ,torsion terms included
   int maxlocal;                    // size of numneigh, firstneigh arrays

src/PERI/pair_peri_lps.cpp

@@ -409,11 +409,6 @@ double PairPeriLPS::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all("All pair coeffs are not set");
 
-  cutsq[i][j] = cut[i][j] * cut[i][j];
-  cutsq[j][i] = cutsq[i][j];
-
-  // set other j,i parameters
-
   bulkmodulus[j][i] = bulkmodulus[i][j];		
   shearmodulus[j][i] = shearmodulus[i][j];		
   s00[j][i] = s00[i][j];		

src/PERI/pair_peri_pmb.cpp

@@ -346,11 +346,6 @@ double PairPeriPMB::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all("All pair coeffs are not set");
 
-  cutsq[i][j] = cut[i][j] * cut[i][j];
-  cutsq[j][i] = cutsq[i][j];
-
-  // set other j,i parameters
-
   kspring[j][i] = kspring[i][j];
   alpha[j][i] = alpha[i][j];
   s00[j][i] = s00[i][j];

src/USER-REAXC/fix_reax_c.cpp

@@ -152,5 +152,5 @@ void FixReaxC::unpack_comm(int n, int first, double *buf)
   m = 0;
   last = first + n;
   for (i = first; i < last; i++)
-    num_bonds[i] = buf[m++];
+    num_bonds[i] = static_cast<int> (buf[m++]);
 }

src/USER-REAXC/reaxc_types.h

@@ -30,7 +30,6 @@
 #include "string.h"
 #include "sys/time.h"
 #include "time.h"
-#include "zlib.h"
 
 /************* SOME DEFS - crucial for reax_types.h *********/
 

src/atom.cpp

@@ -55,6 +55,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
 
   firstgroupname = NULL;
   sortfreq = 1000;
+  nextsort = 0;
   userbinsize = 0.0;
   maxbin = maxnext = 0;
   binhead = NULL;
@@ -1357,6 +1358,10 @@ void Atom::sort()
 {
   int i,m,n,ix,iy,iz,ibin,empty;
 
+  // set next timestep for sorting to take place
+
+  nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq;
+
   // re-setup sort bins if needed
 
   if (domain->box_change) setup_sort_bins();
@@ -1440,10 +1445,6 @@ void Atom::sort()
   //int flagall;
   //MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
   //if (flagall) error->all("Atom sort did not operate correctly");
-
-  // set next timestep for sorting to take place
-
-  nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq;
 }
 
 /* ----------------------------------------------------------------------

src/comm.cpp

@@ -77,7 +77,7 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp)
   allocate_swap(maxswap);
 
   sendlist = (int **) memory->smalloc(maxswap*sizeof(int *),"comm:sendlist");
-  maxsendlist = new int[maxswap];
+  memory->create(maxsendlist,maxswap,"comm:maxsendlist");
   for (int i = 0; i < maxswap; i++) {
     maxsendlist[i] = BUFMIN;
     memory->create(sendlist[i],BUFMIN,"comm:sendlist[i]");
@@ -101,7 +101,7 @@ Comm::~Comm()
 
   if (sendlist) for (int i = 0; i < maxswap; i++) memory->destroy(sendlist[i]);
   memory->sfree(sendlist);
-  delete [] maxsendlist;
+  memory->destroy(maxsendlist);
 
   memory->destroy(buf_send);
   memory->destroy(buf_recv);
@@ -1225,7 +1225,7 @@ void Comm::free_swap()
   memory->destroy(slabhi);
   memory->destroy(firstrecv);
   memory->destroy(pbc_flag);
-  memory->destroy_2d_int_array(pbc);
+  memory->destroy(pbc);
 }
 
 /* ----------------------------------------------------------------------

src/neigh_derive.cpp

@@ -161,6 +161,7 @@ void Neighbor::half_from_full_newton(NeighList *list)
    build skip list for subset of types from parent list
    iskip and ijskip flag which atom types and type pairs to skip
    this is for half and full lists
+   if ghostflag, also store neighbors of ghost atoms & set inum,gnum correctly
 ------------------------------------------------------------------------- */
 
 void Neighbor::skip_from(NeighList *list)
@@ -169,6 +170,7 @@ void Neighbor::skip_from(NeighList *list)
   int *neighptr,*jlist;
 
   int *type = atom->type;
+  int nlocal = atom->nlocal;
   int nall = atom->nlocal + atom->nghost;
 
   int *ilist = list->ilist;
@@ -178,7 +180,8 @@ void Neighbor::skip_from(NeighList *list)
   int *ilist_skip = list->listskip->ilist;
   int *numneigh_skip = list->listskip->numneigh;
   int **firstneigh_skip = list->listskip->firstneigh;
-  int inum_skip = list->listskip->inum;
+  int num_skip = list->listskip->inum;
+  if (list->ghostflag) num_skip += list->listskip->gnum;
 
   int *iskip = list->iskip;
   int **ijskip = list->ijskip;
@@ -191,7 +194,7 @@ void Neighbor::skip_from(NeighList *list)
   // skip I atom entirely if iskip is set for type[I]
   // skip I,J pair if ijskip is set for type[I],type[J]
 
-  for (ii = 0; ii < inum_skip; ii++) {
+  for (ii = 0; ii < num_skip; ii++) {
     i = ilist_skip[ii];
     itype = type[i];
     if (iskip[itype]) continue;
@@ -226,6 +229,14 @@ void Neighbor::skip_from(NeighList *list)
   }
 
   list->inum = inum;
+  if (list->ghostflag) {
+    int num = 0;
+    for (i = 0; i < inum; i++)
+      if (ilist[i] < nlocal) num++;
+      else break;
+    list->inum = num;
+    list->gnum = inum - num;
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -486,6 +497,7 @@ void Neighbor::copy_from(NeighList *list)
   NeighList *listcopy = list->listcopy;
 
   list->inum = listcopy->inum;
+  list->gnum = listcopy->gnum;
   list->ilist = listcopy->ilist;
   list->numneigh = listcopy->numneigh;
   list->firstneigh = listcopy->firstneigh;

src/neigh_full.cpp

@@ -55,6 +55,8 @@ void Neighbor::full_nsq(NeighList *list)
   int npage = 0;
   int npnt = 0;
 
+  // loop over owned atoms, storing neighbors
+
   for (i = 0; i < nlocal; i++) {
 
     if (pgsize - npnt < oneatom) {
@@ -101,6 +103,109 @@ void Neighbor::full_nsq(NeighList *list)
   }
 
   list->inum = inum;
+  list->gnum = 0;
+}
+
+/* ----------------------------------------------------------------------
+   N^2 search for all neighbors
+   include neighbors of ghost atoms
+   every neighbor pair appears in list of both atoms i and j
+------------------------------------------------------------------------- */
+
+void Neighbor::full_nsq_ghost(NeighList *list)
+{
+  int i,j,n,itype,jtype,which;
+  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;
+
+  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 all atoms, owned and ghost
+    // skip i = j
+
+    if (i < nlocal) {
+      for (j = 0; j < nall; j++) {
+	if (i == j) 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]);
+	  else which = 0;
+	  if (which == 0) neighptr[n++] = j;
+	  else if (which > 0) neighptr[n++] = which*nall + j;
+	}
+      }
+    } else {
+      for (j = 0; j < nall; j++) {
+	if (i == j) 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 <= cutneighghostsq[itype][jtype]) {
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
+	  else which = 0;
+	  if (which == 0) neighptr[n++] = j;
+	  else if (which > 0) neighptr[n++] = which*nall + j;
+	}
+      }
+    }
+
+    ilist[inum++] = i;
+    firstneigh[i] = neighptr;
+    numneigh[i] = n;
+    npnt += n;
+    if (n > oneatom || npnt >= pgsize)
+      error->one("Neighbor list overflow, boost neigh_modify one or page");
+  }
+
+  list->inum = atom->nlocal;
+  list->gnum = inum - atom->nlocal;
 }
 
 /* ----------------------------------------------------------------------
@@ -114,12 +219,10 @@ void Neighbor::full_bin(NeighList *list)
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr;
 
-  // bin local & ghost atoms
+  // bin owned & ghost atoms
 
   bin_atoms();
 
-  // loop over each atom, storing neighbors
-
   int **special = atom->special;
   int **nspecial = atom->nspecial;
   int *tag = atom->tag;
@@ -144,6 +247,8 @@ void Neighbor::full_bin(NeighList *list)
   int npage = 0;
   int npnt = 0;
 
+  // loop over owned atoms, storing neighbors
+
   for (i = 0; i < nlocal; i++) {
 
     if (pgsize - npnt < oneatom) {
@@ -195,6 +300,135 @@ void Neighbor::full_bin(NeighList *list)
   }
 
   list->inum = inum;
+  list->gnum = 0;
+}
+
+/* ----------------------------------------------------------------------
+   binned neighbor list construction for all neighbors
+   include neighbors of ghost atoms
+   every neighbor pair appears in list of both atoms i and j
+------------------------------------------------------------------------- */
+
+void Neighbor::full_bin_ghost(NeighList *list)
+{
+  int i,j,k,n,itype,jtype,ibin,which;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int xbin,ybin,zbin,xbin2,ybin2,zbin2;
+  int *neighptr;
+
+  // bin owned & ghost atoms
+
+  bin_atoms();
+
+  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;
+
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+  int **pages = list->pages;
+  int nstencil = list->nstencil;
+  int *stencil = list->stencil;
+  int **stencilxyz = list->stencilxyz;
+
+  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 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
+
+    if (i < nlocal) {
+      ibin = coord2bin(x[i]);
+      for (k = 0; k < nstencil; k++) {
+	for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
+	  if (i == j) 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]);
+	    else which = 0;
+	    if (which == 0) neighptr[n++] = j;
+	    else if (which > 0) neighptr[n++] = which*nall + j;
+	  }
+	}
+      }
+
+    } else {
+      ibin = coord2bin(x[i],xbin,ybin,zbin);
+      for (k = 0; k < nstencil; k++) {
+	xbin2 = xbin + stencilxyz[k][0];
+	ybin2 = ybin + stencilxyz[k][1];
+	zbin2 = zbin + stencilxyz[k][2];
+	if (xbin2 < 0 || xbin2 >= mbinx ||
+	    ybin2 < 0 || ybin2 >= mbiny ||
+	    zbin2 < 0 || zbin2 >= mbinz) continue;
+	for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
+	  if (i == j) 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 <= cutneighghostsq[itype][jtype]) {
+	    if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
+	    else which = 0;
+	    if (which == 0) neighptr[n++] = j;
+	    else if (which > 0) neighptr[n++] = which*nall + j;
+	  }
+	}
+      }
+    }
+
+    ilist[inum++] = i;
+    firstneigh[i] = neighptr;
+    numneigh[i] = n;
+    npnt += n;
+    if (n > oneatom || npnt >= pgsize)
+      error->one("Neighbor list overflow, boost neigh_modify one or page");
+  }
+
+  list->inum = atom->nlocal;
+  list->gnum = inum - atom->nlocal;
 }
 
 /* ----------------------------------------------------------------------
@@ -297,4 +531,5 @@ void Neighbor::full_multi(NeighList *list)
   }
 
   list->inum = inum;
+  list->gnum = 0;
 }

src/neigh_list.cpp

@@ -30,10 +30,10 @@ enum{NSQ,BIN,MULTI};     // also in neighbor.cpp
 
 NeighList::NeighList(LAMMPS *lmp, int size) : Pointers(lmp)
 {
-  maxlocal = 0;
+  maxatoms = 0;
   pgsize = size;
 
-  inum = 0;
+  inum = gnum = 0;
   ilist = NULL;
   numneigh = NULL;
   firstneigh = NULL;
@@ -59,6 +59,7 @@ NeighList::NeighList(LAMMPS *lmp, int size) : Pointers(lmp)
 
   maxstencil = 0;
   stencil = NULL;
+  stencilxyz = NULL;
 
   maxstencil_multi = 0;
   nstencil_multi = NULL;
@@ -88,6 +89,8 @@ NeighList::~NeighList()
   memory->destroy(ijskip);
 
   if (maxstencil) memory->destroy(stencil);
+  if (ghostflag) memory->destroy(stencilxyz);
+
   if (maxstencil_multi) {
     for (int i = 1; i <= atom->ntypes; i++) {
       memory->destroy(stencil_multi[i]);
@@ -106,22 +109,27 @@ NeighList::~NeighList()
 
 void NeighList::grow(int nmax)
 {
-  // skip if grow not needed
+  // skip if grow not needed by this list
+  // each list stores own maxatoms, b/c list->grow() called at different times
+  // if list does not store neighbors of ghosts, compare nmax to maxatoms
+  // else compare nlocal+nghost to maxatoms
+  // if reset list size, set it to nmax
 
-  if (nmax <= maxlocal) return;
+  if (!ghostflag && atom->nlocal <= maxatoms) return;
-  maxlocal = nmax;
+  if (ghostflag && atom->nlocal+atom->nghost <= maxatoms) return;
+  maxatoms = nmax;
 
   memory->destroy(ilist);
   memory->destroy(numneigh);
   memory->sfree(firstneigh);
-  memory->create(ilist,maxlocal,"neighlist:ilist");
+  memory->create(ilist,maxatoms,"neighlist:ilist");
-  memory->create(numneigh,maxlocal,"neighlist:numneigh");
+  memory->create(numneigh,maxatoms,"neighlist:numneigh");
   firstneigh = (int **)
-    memory->smalloc(maxlocal*sizeof(int *),"neighlist:firstneigh");
+    memory->smalloc(maxatoms*sizeof(int *),"neighlist:firstneigh");
 
   if (dnum) 
     firstdouble = (double **)
-      memory->smalloc(maxlocal*sizeof(double *),"neighlist:firstdouble");
+      memory->smalloc(maxatoms*sizeof(double *),"neighlist:firstdouble");
 }
 
 /* ----------------------------------------------------------------------
@@ -138,6 +146,10 @@ void NeighList::stencil_allocate(int smax, int style)
       maxstencil = smax;
       memory->destroy(stencil);
       memory->create(stencil,maxstencil,"neighlist:stencil");
+      if (ghostflag) {
+	memory->destroy(stencilxyz);
+	memory->create(stencilxyz,maxstencil,3,"neighlist:stencilxyz");
+      }
     }
 
   } else {
@@ -220,6 +232,7 @@ void NeighList::print_attributes()
   printf("  %d = build flag\n",buildflag);
   printf("  %d = grow flag\n",growflag);
   printf("  %d = stencil flag\n",stencilflag);
+  printf("  %d = ghost flag\n",ghostflag);
   printf("\n");
   printf("  %d = pair\n",rq->pair);
   printf("  %d = fix\n",rq->fix);
@@ -237,6 +250,7 @@ void NeighList::print_attributes()
   printf("\n");
   printf("  %d = occasional\n",rq->occasional);
   printf("  %d = dnum\n",rq->dnum);
+  printf("  %d = ghost\n",rq->ghost);
   printf("  %d = copy\n",rq->copy);
   printf("  %d = skip\n",rq->skip);
   printf("  %d = otherlist\n",rq->otherlist);
@@ -246,24 +260,26 @@ void NeighList::print_attributes()
 
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory
-   if growflag = 0, maxlocal & maxpage will also be 0
+   if growflag = 0, maxatoms & maxpage will also be 0
    if stencilflag = 0, maxstencil * maxstencil_multi will also be 0
 ------------------------------------------------------------------------- */
 
 bigint NeighList::memory_usage()
 {
   bigint bytes = 0;
-  bytes += memory->usage(ilist,maxlocal);
+  bytes += memory->usage(ilist,maxatoms);
-  bytes += memory->usage(numneigh,maxlocal);
+  bytes += memory->usage(numneigh,maxatoms);
-  bytes += maxlocal * sizeof(int *);
+  bytes += maxatoms * sizeof(int *);
   bytes += memory->usage(pages,maxpage,pgsize);
 
   if (dnum) {
-    bytes += maxlocal * sizeof(double *);
+    bytes += maxatoms * sizeof(double *);
     bytes += memory->usage(dpages,maxpage,dnum*pgsize);
   }
 
   if (maxstencil) bytes += memory->usage(stencil,maxstencil);
+  if (ghostflag) bytes += memory->usage(stencilxyz,maxstencil,3);
+
   if (maxstencil_multi) {
     bytes += memory->usage(stencil_multi,atom->ntypes,maxstencil_multi);
     bytes += memory->usage(distsq_multi,atom->ntypes,maxstencil_multi);

src/neigh_list.h

@@ -26,10 +26,12 @@ class NeighList : protected Pointers {
   int buildflag;                   // 1 if pair_build invoked every reneigh
   int growflag;                    // 1 if stores atom-based arrays & pages
   int stencilflag;                 // 1 if stores stencil arrays
+  int ghostflag;                   // 1 if it stores neighbors of ghosts
 
   // data structs to store neighbor pairs I,J and associated values
 
   int inum;                        // # of I atoms neighbors are stored for
+  int gnum;                        // # of ghost atoms neighbors are stored for
   int *ilist;                      // local indices of I atoms
   int *numneigh;                   // # of J neighbors for each I atom
   int **firstneigh;                // ptr to 1st J int value of each I atom
@@ -45,7 +47,7 @@ class NeighList : protected Pointers {
   // iskip,ijskip are just ptrs to corresponding request
 
   int *iskip;         // iskip[i] = 1 if atoms of type I are not in list
-  int **ijskip;       // ijskip[i][j] =1 if pairs of type I,J are not in list
+  int **ijskip;       // ijskip[i][j] = 1 if pairs of type I,J are not in list
 
   // settings and pointers for related neighbor lists and fixes
 
@@ -64,6 +66,7 @@ class NeighList : protected Pointers {
   int maxstencil;                  // max size of stencil
   int nstencil;                    // # of bins in stencil
   int *stencil;                    // list of bin offsets
+  int **stencilxyz;                // bin offsets in xyz dims
 
   int maxstencil_multi;            // max sizes of stencils
   int *nstencil_multi;             // # bins in each type-based multi stencil
@@ -80,7 +83,7 @@ class NeighList : protected Pointers {
   bigint memory_usage();
 
  private:
-  int maxlocal;                    // size of allocated atom arrays
+  int maxatoms;                    // size of allocated atom arrays
 };
 
 }

src/neigh_request.cpp

@@ -37,6 +37,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   gran = granhistory = 0;
   respainner = respamiddle = respaouter = 0;
   half_from_full = 0;
+  ghost = 0;
 
   // default is use newton_pair setting in force
 
@@ -50,6 +51,10 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
 
   dnum = 0;
 
+  // default is no neighbors of ghosts
+
+  ghost = 0;
+
   // default is no copy or skip
 
   copy = 0;
@@ -96,6 +101,7 @@ int NeighRequest::identical(NeighRequest *other)
   if (newton != other->newton) same = 0;
   if (occasional != other->occasional) same = 0;
   if (dnum != other->dnum) same = 0;
+  if (ghost != other->ghost) same = 0;
 
   if (copy != other->copy) same = 0;
   if (same_skip(other) == 0) same = 0;
@@ -122,6 +128,7 @@ int NeighRequest::same_kind(NeighRequest *other)
   if (respaouter != other->respaouter) same = 0;
   if (half_from_full != other->half_from_full) same = 0;
   if (newton != other->newton) same = 0;
+  if (ghost != other->ghost) same = 0;
 
   return same;
 }
@@ -169,4 +176,5 @@ void NeighRequest::copy_request(NeighRequest *other)
 
   newton = other->newton;
   dnum = other->dnum;
+  ghost = other->ghost;
 }

src/neigh_request.h

@@ -31,7 +31,7 @@ class NeighRequest : protected Pointers {
   int compute;
   int command;
 
-  // kind of list requested, one flag is 1, others are 0
+  // kind of list requested
   // set by requesting class
 
   int half;              // 1 if half neigh list
@@ -63,6 +63,10 @@ class NeighRequest : protected Pointers {
 
   int dnum;
 
+  // 1 if also need neighbors of ghosts
+
+  int ghost;
+
   // set by neighbor and pair_hybrid after all requests are made
   // these settings do not change kind value
 

src/neigh_stencil.cpp

@@ -366,6 +366,28 @@ void Neighbor::stencil_full_bin_2d(NeighList *list,
 
 /* ---------------------------------------------------------------------- */
 
+void Neighbor::stencil_full_ghost_bin_2d(NeighList *list,
+					 int sx, int sy, int sz)
+{
+  int i,j;
+  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_full_bin_3d(NeighList *list,
 				   int sx, int sy, int sz)
 {
@@ -382,6 +404,28 @@ void Neighbor::stencil_full_bin_3d(NeighList *list,
   list->nstencil = nstencil;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_full_ghost_bin_3d(NeighList *list,
+					 int sx, int sy, int sz)
+{
+  int i,j,k;
+  int *stencil = list->stencil;
+  int **stencilxyz = list->stencilxyz;
+  int nstencil = 0;
+
+  for (k = -sz; k <= sz; k++)
+    for (j = -sy; j <= sy; j++)
+      for (i = -sx; i <= sx; i++)
+	if (bin_distance(i,j,k) < cutneighmaxsq) {
+	  stencilxyz[nstencil][0] = i;
+	  stencilxyz[nstencil][1] = j;
+	  stencilxyz[nstencil][2] = k;
+	  stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
+	}
+
+  list->nstencil = nstencil;
+}
 
 /* ---------------------------------------------------------------------- */
 

src/neighbor.cpp

@@ -73,6 +73,7 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
   build_once = 0;
 
   cutneighsq = NULL;
+  cutneighghostsq = NULL;
   cuttype = NULL;
   cuttypesq = NULL;
   fixchecklist = NULL;
@@ -139,6 +140,7 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
 Neighbor::~Neighbor()
 {
   memory->destroy(cutneighsq);
+  memory->destroy(cutneighghostsq);
   delete [] cuttype;
   delete [] cuttypesq;
   delete [] fixchecklist;
@@ -226,6 +228,7 @@ void Neighbor::init()
   n = atom->ntypes;
   if (cutneighsq == NULL) {
     memory->create(cutneighsq,n+1,n+1,"neigh:cutneighsq");
+    memory->create(cutneighghostsq,n+1,n+1,"neigh:cutneighghostsq");
     cuttype = new double[n+1];
     cuttypesq = new double[n+1];
   }
@@ -248,6 +251,11 @@ void Neighbor::init()
       cuttypesq[i] = MAX(cuttypesq[i],cut*cut);
       cutneighmin = MIN(cutneighmin,cut);
       cutneighmax = MAX(cutneighmax,cut);
+
+      if (force->pair && force->pair->ghostneigh) {
+	cut = force->pair->cutghost[i][j] + skin;
+	cutneighghostsq[i][j] = cut*cut;
+      }
     }
   }
   cutneighmaxsq = cutneighmax * cutneighmax;
@@ -463,7 +471,7 @@ void Neighbor::init()
     // skip: point this list at request->otherlist, copy skip info from request
     // half_from_full: point this list at preceeding full list
     // granhistory: set preceeding list's listgranhistory to this list
-    //   also set precedding list's ptr to FixShearHistory
+    //   also set preceeding list's ptr to FixShearHistory
     // respaouter: point this list at preceeding 1/2 inner/middle lists
     // pair and half: if there is a full non-occasional non-skip list
     //   change this list to half_from_full and point at the full list
@@ -551,11 +559,14 @@ void Neighbor::init()
       else stencil_create[i] = NULL;
     }
 
-    // set each list's build/grow/stencil flags based on neigh request
+    // set each list's build/grow/stencil/ghost flags based on neigh request
     // buildflag = 1 if its pair_build() invoked every reneighbor
     // growflag = 1 if it stores atom-based arrays and pages
     // stencilflag = 1 if it stores stencil arrays
+    // ghostflag = 1 if it stores neighbors of ghosts
+    // anyghostlist = 1 if any non-occasional list stores neighbors of ghosts
 
+    anyghostlist = 0;
     for (i = 0; i < nlist; i++) {
       lists[i]->buildflag = 1;
       if (pair_build[i] == NULL) lists[i]->buildflag = 0;
@@ -567,6 +578,10 @@ void Neighbor::init()
       lists[i]->stencilflag = 1;
       if (style == NSQ) lists[i]->stencilflag = 0;
       if (stencil_create[i] == NULL) lists[i]->stencilflag = 0;
+
+      lists[i]->ghostflag = 0;
+      if (requests[i]->ghost) lists[i]->ghostflag = 1;
+      if (requests[i]->ghost && !requests[i]->occasional) anyghostlist = 1;
     }
 
 #ifdef NEIGH_LIST_DEBUG
@@ -773,8 +788,8 @@ int Neighbor::request(void *requestor)
    determine which pair_build function each neigh list needs
    based on settings of neigh request
    copy -> copy_from function
-   skip -> granular function if gran with granhistory
+   skip -> granular function if gran with granhistory,
-           respa function if respaouter
+           respa function if respaouter,
 	   skip_from function for everything else
    half_from_full, half, full, gran, respaouter ->
      choose by newton and rq->newton and tri settings
@@ -804,10 +819,8 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
   } else if (rq->half) {
     if (style == NSQ) {
       if (rq->newton == 0) {
-	if (newton_pair == 0)
+	if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton;
-	  pb = &Neighbor::half_nsq_no_newton;
+	else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton;
-	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) {
@@ -821,9 +834,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
       } 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) {
+      } else if (rq->newton == 2) pb = &Neighbor::half_bin_no_newton;
-	pb = &Neighbor::half_bin_no_newton;
-      }
     } else if (style == MULTI) {
       if (rq->newton == 0) {
 	if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton;
@@ -832,15 +843,20 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
       } else if (rq->newton == 1) {
 	if (triclinic == 0) pb = &Neighbor::half_multi_newton;
 	else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri;
-      } else if (rq->newton == 2) {
+      } else if (rq->newton == 2) pb = &Neighbor::half_multi_no_newton;
-	pb = &Neighbor::half_multi_no_newton;
-      } 
     }
 
   } else if (rq->full) {
-    if (style == NSQ) pb = &Neighbor::full_nsq;
+    if (style == NSQ) {
-    else if (style == BIN) pb = &Neighbor::full_bin;
+      if (rq->ghost == 0) pb = &Neighbor::full_nsq;
-    else if (style == MULTI) pb = &Neighbor::full_multi;
+      else if (rq->ghost == 1) pb = &Neighbor::full_nsq_ghost;
+    } else if (style == BIN) {
+      if (rq->ghost == 0) pb = &Neighbor::full_bin;
+      else if (rq->ghost == 1) pb = &Neighbor::full_bin_ghost;
+    } else if (style == MULTI) {
+      if (rq->ghost == 0) pb = &Neighbor::full_multi;
+      else error->all("Neighbor multi not yet enabled for ghost neighbors");
+    }
 
   } else if (rq->gran) {
     if (style == NSQ) {
@@ -865,6 +881,11 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
       error->all("Neighbor multi not yet enabled for rRESPA");
   }
 
+  // general error check
+
+  if (rq->ghost && !rq->full)
+    error->all("Neighbors of ghost atoms only allowed for full neighbor lists");
+
   pair_build[index] = pb;
 }
 
@@ -961,8 +982,14 @@ void Neighbor::choose_stencil(int index, NeighRequest *rq)
 
   } else if (rq->full) {
     if (style == BIN) {
-      if (dimension == 2) sc = &Neighbor::stencil_full_bin_2d;
+      if (dimension == 2) {
-      else if (dimension == 3) sc = &Neighbor::stencil_full_bin_3d;
+	if (rq->ghost) sc = &Neighbor::stencil_full_ghost_bin_2d;
+	else sc = &Neighbor::stencil_full_bin_2d;
+      }
+      else if (dimension == 3) {
+	if (rq->ghost) sc = &Neighbor::stencil_full_ghost_bin_3d;
+	else sc = &Neighbor::stencil_full_bin_3d;
+      }
     } else if (style == MULTI) {
       if (dimension == 2) sc = &Neighbor::stencil_full_multi_2d;
       else if (dimension == 3) sc = &Neighbor::stencil_full_multi_3d;
@@ -1119,14 +1146,16 @@ void Neighbor::build()
     }
   }
 
-  // if necessary, extend atom arrays in pairwise lists
+  // if any lists store neighbors of ghosts:
-  // only for lists with growflag set and which are used every reneighbor
+  //   invoke grow() on all in case nlocal+nghost is now too big
+  // else only invoke grow() if nlocal has exceeded previous list size
+  // only for lists with growflag set and which are perpetual
 
-  if (atom->nlocal > maxlocal) {
+  if (anyghostlist || atom->nlocal > maxlocal) {
     maxlocal = atom->nmax;
     for (i = 0; i < nglist; i++) lists[glist[i]]->grow(maxlocal);
   }
-
+  
   // extend atom bin list if necessary
 
   if (style != NSQ && atom->nmax > maxbin) {
@@ -1612,6 +1641,42 @@ int Neighbor::coord2bin(double *x)
   return (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo);
 }
 
+/* ----------------------------------------------------------------------
+   same as coord2bin, but also return ix,iy,iz offsets in each dim
+------------------------------------------------------------------------- */
+
+int Neighbor::coord2bin(double *x, int &ix, int &iy, int &iz)
+{
+  if (x[0] >= bboxhi[0])
+    ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx) + nbinx;
+  else if (x[0] >= bboxlo[0]) {
+    ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx);
+    ix = MIN(ix,nbinx-1);
+  } else
+    ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - 1;
+  
+  if (x[1] >= bboxhi[1])
+    iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy) + nbiny;
+  else if (x[1] >= bboxlo[1]) {
+    iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy);
+    iy = MIN(iy,nbiny-1);
+  } else
+    iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - 1;
+  
+  if (x[2] >= bboxhi[2])
+    iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz) + nbinz;
+  else if (x[2] >= bboxlo[2]) {
+    iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz);
+    iz = MIN(iz,nbinz-1);
+  } else
+    iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - 1;
+
+  ix -= mbinxlo;
+  iy -= mbinylo;
+  iz -= mbinzlo;
+  return iz*mbiny*mbinx + iy*mbinx + ix;
+}
+
 /* ----------------------------------------------------------------------
    test if atom pair i,j is excluded from neighbor list
    due to type, group, molecule settings from neigh_modify command

src/neighbor.h

@@ -86,6 +86,7 @@ class Neighbor : protected Pointers {
   int *fixchecklist;               // which fixes to check
 
   double **cutneighsq;             // neighbor cutneigh sq for each type pair
+  double **cutneighghostsq;        // neighbor cutnsq for each ghost type pair
   double cutneighmaxsq;            // cutneighmax squared
   double *cuttypesq;               // cuttype squared
 
@@ -130,6 +131,9 @@ class Neighbor : protected Pointers {
 
   int special_flag[4];             // flags for 1-2, 1-3, 1-4 neighbors
 
+  int anyghostlist;                // 1 if any non-occasional list
+                                   // stores neighbors of ghosts
+
   int exclude;                     // 0 if no type/group exclusions, 1 if yes
 
   int nex_type;                    // # of entries in type exclusion list
@@ -155,6 +159,7 @@ class Neighbor : protected Pointers {
   void bin_atoms();                     // bin all atoms
   double bin_distance(int, int, int);   // distance between binx
   int coord2bin(double *);              // mapping atom coord to a bin
+  int coord2bin(double *, int &, int &, int&); // ditto
 
   int exclusion(int, int, int, int, int *, int *);  // test for pair exclusion
   void choose_build(int, class NeighRequest *);
@@ -210,7 +215,9 @@ class Neighbor : protected Pointers {
   void half_multi_newton_tri(class NeighList *);
 
   void full_nsq(class NeighList *);
+  void full_nsq_ghost(class NeighList *);
   void full_bin(class NeighList *);
+  void full_bin_ghost(class NeighList *);
   void full_multi(class NeighList *);
 
   void half_from_full_no_newton(class NeighList *);
@@ -252,7 +259,9 @@ class Neighbor : protected Pointers {
   void stencil_half_multi_3d_newton_tri(class NeighList *, int, int, int);
 
   void stencil_full_bin_2d(class NeighList *, int, int, int);
+  void stencil_full_ghost_bin_2d(class NeighList *, int, int, int);
   void stencil_full_bin_3d(class NeighList *, int, int, int);
+  void stencil_full_ghost_bin_3d(class NeighList *, int, int, int);
   void stencil_full_multi_2d(class NeighList *, int, int, int);
   void stencil_full_multi_3d(class NeighList *, int, int, int);
 

src/pair.cpp

@@ -55,6 +55,8 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
   respa_enable = 0;
   one_coeff = 0;
   no_virial_compute = 0;
+  ghostneigh = 0;
+
   nextra = 0;
   pvector = NULL;
 

src/pair.h

@@ -29,7 +29,7 @@ class Pair : protected Pointers {
   double *eatom,**vatom;         // accumulated per-atom energy/virial
 
   double cutforce;               // max cutoff for all atom pairs
-  double **cutsq;                // max cutoff sq for each atom pair
+  double **cutsq;                // cutoff sq for each atom pair
   int **setflag;                 // 0/1 = whether each i,j has been set
 
   int comm_forward;              // size of forward communication (0 if none)
@@ -39,6 +39,8 @@ class Pair : protected Pointers {
   int respa_enable;              // 1 if inner/middle/outer rRESPA routines
   int one_coeff;                 // 1 if allows only one coeff * * call
   int no_virial_compute;         // 1 if does not invoke virial_compute()
+  int ghostneigh;                // 1 if pair style needs neighbors of ghosts
+  double **cutghost;             // cutoff for each ghost pair
 
   int tail_flag;                 // pair_modify flag for LJ tail correction
   double etail,ptail;            // energy/pressure tail corrections

src/pair_hybrid.cpp

@@ -51,6 +51,7 @@ PairHybrid::~PairHybrid()
   if (allocated) {
     memory->destroy_2d_int_array(setflag);
     memory->destroy_2d_double_array(cutsq);
+    memory->destroy_2d_double_array(cutghost);
     memory->destroy_2d_int_array(nmap);
     memory->destroy_3d_int_array(map);
   }
@@ -157,6 +158,7 @@ void PairHybrid::allocate()
       setflag[i][j] = 0;
 
   cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
+  cutghost = memory->create_2d_double_array(n+1,n+1,"pair:cutghost");
 
   nmap = memory->create_2d_int_array(n+1,n+1,"pair:nmap");
   map = memory->create_3d_int_array(n+1,n+1,nstyles,"pair:map");
@@ -187,6 +189,7 @@ void PairHybrid::settings(int narg, char **arg)
   if (allocated) {
     memory->destroy_2d_int_array(setflag);
     memory->destroy_2d_double_array(cutsq);
+    memory->destroy_2d_double_array(cutghost);
     memory->destroy_2d_int_array(nmap);
     memory->destroy_3d_int_array(map);
   }
@@ -262,6 +265,7 @@ void PairHybrid::settings(int narg, char **arg)
   // single_enable = 0 if any sub-style = 0
   // respa_enable = 1 if any sub-style is set
   // no_virial_compute = 1 if any sub-style is set
+  // ghostneigh = 1 if any sub-style is set
 
   for (m = 0; m < nstyles; m++)
     if (styles[m]->single_enable == 0) single_enable = 0;
@@ -269,6 +273,8 @@ void PairHybrid::settings(int narg, char **arg)
     if (styles[m]->respa_enable) respa_enable = 1;
   for (m = 0; m < nstyles; m++)
     if (styles[m]->no_virial_compute) no_virial_compute = 1;
+  for (m = 0; m < nstyles; m++)
+    if (styles[m]->ghostneigh) ghostneigh = 1;
 }
 
 /* ----------------------------------------------------------------------
@@ -455,12 +461,14 @@ double PairHybrid::init_one(int i, int j)
   }
 
   // call init/mixing for all sub-styles of I,J
-  // set cutsq in sub-style just as pair::init_one() does
+  // set cutsq in sub-style just as Pair::init() does via call to init_one()
+  // set cutghost for I,J and J,I just as sub-style does
   // sum tail corrections for I,J
   // return max cutoff of all sub-styles assigned to I,J
   // if no sub-styles assigned to I,J (pair_coeff none), cutmax = 0.0 returned
 
   double cutmax = 0.0;
+  cutghost[i][j] = cutghost[j][i] = 0.0;
   if (tail_flag) etail_ij = ptail_ij = 0.0;
 
   nmap[j][i] = nmap[i][j];
@@ -470,6 +478,9 @@ double PairHybrid::init_one(int i, int j)
     double cut = styles[map[i][j][k]]->init_one(i,j);
     styles[map[i][j][k]]->cutsq[i][j] = 
       styles[map[i][j][k]]->cutsq[j][i] = cut*cut;
+    if (styles[map[i][j][k]]->ghostneigh)
+      cutghost[i][j] = cutghost[j][i] = 
+	MAX(cutghost[i][j],styles[map[i][j][k]]->cutghost[i][j]);
     if (tail_flag) {
       etail_ij += styles[map[i][j][k]]->etail_ij;
       ptail_ij += styles[map[i][j][k]]->ptail_ij;
@@ -491,7 +502,7 @@ void PairHybrid::modify_requests()
 
   // loop over pair requests only
   // if list is skip list and not copy, look for non-skip list of same kind
-  // if one exists, point at that one
+  // if one exists, point at that one via otherlist
   // else make new non-skip request of same kind and point at that one
   //   don't bother to set ID for new request, since pair hybrid ignores list
   // only exception is half_from_full: