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

src/OPT/pair_eam_opt.cpp

@@ -16,7 +16,7 @@
      James Fischer, High Performance Technologies, Inc.
      Charles Cornwell, High Performance Technologies, Inc.
      David Richie, Stone Ridge Technology
-     Vincent Natol, Stone Ridge Technology
+     Vincent Natoli, Stone Ridge Technology
 ------------------------------------------------------------------------- */
 
 #include "pair_eam_opt.h"

src/OPT/pair_lj_charmm_coul_long_opt.cpp

@@ -15,7 +15,7 @@
    Contributing authors:  
      James Fischer, High Performance Technologies, Inc.
      David Richie, Stone Ridge Technology
-     Vincent Natol, Stone Ridge Technology
+     Vincent Natoli, Stone Ridge Technology
 ------------------------------------------------------------------------- */
 
 #include "pair_lj_charmm_coul_long_opt.h"

src/OPT/pair_lj_cut_opt.cpp

@@ -15,7 +15,7 @@
    Contributing authors:  
      James Fischer, High Performance Technologies, Inc.
      David Richie, Stone Ridge Technology
-     Vincent Natol, Stone Ridge Technology
+     Vincent Natoli, Stone Ridge Technology
 ------------------------------------------------------------------------- */
 
 #include "pair_lj_cut_opt.h"

src/OPT/pair_morse_opt.cpp

@@ -15,7 +15,7 @@
    Contributing authors:  
      James Fischer, High Performance Technologies, Inc.
      David Richie, Stone Ridge Technology
-     Vincent Natol, Stone Ridge Technology
+     Vincent Natoli, Stone Ridge Technology
 ------------------------------------------------------------------------- */
 
 #include "pair_morse_opt.h"

src/comm.cpp

@@ -39,6 +39,8 @@ using namespace LAMMPS_NS;
 #define MAX(a,b) ((a) > (b) ? (a) : (b))
 #define BIG 1.0e20
 
+enum{SINGLE,MULTI};
+
 /* ----------------------------------------------------------------------
    setup MPI and allocate buffer space 
 ------------------------------------------------------------------------- */
@@ -50,6 +52,9 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp)
 
   user_procgrid[0] = user_procgrid[1] = user_procgrid[2] = 0;
 
+  style = SINGLE;
+  multilo = multihi = NULL;
+
   // initialize comm buffers & exchange memory
 
   maxsend = BUFMIN;
@@ -78,6 +83,11 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp)
 Comm::~Comm()
 {
   free_swap();
+  if (style == MULTI) {
+    free_multi();
+    memory->destroy_2d_double_array(cutghostmulti);
+  }
+
   memory->sfree(maxsendlist);
   if (sendlist) for (int i = 0; i < maxswap; i++) memory->sfree(sendlist[i]);
   memory->sfree(sendlist);
@@ -118,6 +128,18 @@ void Comm::init()
   }
 
   if (force->newton == 0) maxreverse = 0;
+
+  // memory for multi-style communication
+
+  if (style == MULTI && multilo == NULL) {
+    allocate_multi(maxswap);
+    cutghostmulti = 
+      memory->create_2d_double_array(atom->ntypes+1,3,"comm:cutghostmulti");
+  }
+  if (style == SINGLE && multilo) {
+    free_multi();
+    memory->destroy_2d_double_array(cutghostmulti);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -165,42 +187,62 @@ void Comm::set_procs()
 }
 
 /* ----------------------------------------------------------------------
-   setup spatial-decomposition communication patterns 
+   setup spatial-decomposition communication patterns
-   function of neighbor cutoff and current box size 
+   function of neighbor cutoff(s) and current box size
+   single style sets slab boundaries (slablo,slabhi) based on max cutoff
+   multi style sets type-dependent slab boundaries (multilo,multihi)
 ------------------------------------------------------------------------- */
 
 void Comm::setup()
 {
-  // cutghost = distance at which ghost atoms need to be acquired
+  // cutghost = max distance at which ghost atoms need to be acquired
   // for orthogonal:
   //   cutghost is in box coords = neigh->cutghost in all 3 dims
   // for triclinic:
   //   neigh->cutghost = distance between tilted planes in box coords
   //   cutghost is in lamda coords = distance between those planes
+  // for multi:
+  //   cutghostmulti = same as cutghost, only for each atom type
 
+  int i;
+  int ntypes = atom->ntypes;
   double *prd,*prd_border,*sublo,*subhi;
-
+  
   if (triclinic == 0) {
     prd = prd_border = domain->prd;
     sublo = domain->sublo;
     subhi = domain->subhi;
     cutghost[0] = cutghost[1] = cutghost[2] = neighbor->cutghost;
+    if (style == MULTI) {
+      double *cuttype = neighbor->cuttype;
+      for (i = 1; i <= ntypes; i++)
+	cutghostmulti[i][0] = cutghostmulti[i][1] = cutghostmulti[i][2] = 
+	  cuttype[i];
+    }
   } else {
     prd = domain->prd;
     prd_border = domain->prd_lamda;
     sublo = domain->sublo_lamda;
     subhi = domain->subhi_lamda;
     double *h_inv = domain->h_inv;
-    double length;
+    double length0,length1,length2;
-    length = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]);
+    length0 = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]);
-    cutghost[0] = neighbor->cutghost * length;
+    cutghost[0] = neighbor->cutghost * length0;
-    length = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]);
+    length1 = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]);
-    cutghost[1] = neighbor->cutghost * length;
+    cutghost[1] = neighbor->cutghost * length1;
-    length = h_inv[2];
+    length2 = h_inv[2];
-    cutghost[2] = neighbor->cutghost * length;
+    cutghost[2] = neighbor->cutghost * length2;
+    if (style == MULTI) {
+      double *cuttype = neighbor->cuttype;
+      for (i = 1; i <= ntypes; i++) {
+	cutghostmulti[i][0] = cuttype[i] * length0;
+	cutghostmulti[i][1] = cuttype[i] * length1;
+	cutghostmulti[i][2] = cuttype[i] * length2;
+      }
+    }
   }
 
-  // need = # of procs I need atoms from in each dim
+  // need = # of procs I need atoms from in each dim based on max cutoff
   // for 2d, don't communicate in z
 
   need[0] = static_cast<int> (cutghost[0] * procgrid[0] / prd_border[0]) + 1;
@@ -219,41 +261,57 @@ void Comm::setup()
   // allocate comm memory
 
   nswap = 2 * (need[0]+need[1]+need[2]);
-  if (nswap > maxswap) grow_swap();
+  if (nswap > maxswap) grow_swap(nswap);
 
   // setup parameters for each exchange:
   // sendproc = proc to send to at each swap
   // recvproc = proc to recv from at each swap
-  // slablo/slabhi = boundaries for slab of atoms to send at each swap
+  // for style SINGLE:
+  //   slablo/slabhi = boundaries for slab of atoms to send at each swap
   //   use -BIG/midpt/BIG to insure all atoms included even if round-off occurs
   //   if round-off, atoms recvd across PBC can be < or > than subbox boundary
   //   note that borders() only loops over subset of atoms during each swap
-  // set slablo > slabhi for swaps across non-periodic boundaries
+  //   set slablo > slabhi for swaps across non-periodic boundaries
-  //   this insures no atoms are swapped
+  //     this insures no atoms are swapped
-  //   only for procs owning sub-box at non-periodic end of global box
+  //     only for procs owning sub-box at non-periodic end of global box
-  // pbc_flag: 0 = nothing across a boundary, 1 = somthing across a boundary
+  // for style MULTI:
+  //   multilo/multihi is same as slablo/slabhi, only for each atom type
+  // pbc_flag: 0 = nothing across a boundary, 1 = something across a boundary
   // pbc = -1/0/1 for PBC factor in each of 3/6 orthog/triclinic dirs
   // for triclinic, slablo/hi and pbc_border will be used in lamda (0-1) coords
   // 1st part of if statement is sending to the west/south/down
   // 2nd part of if statement is sending to the east/north/up
 
+  int dim,ineed;
+
   int iswap = 0;
-  for (int dim = 0; dim < 3; dim++) {
+  for (dim = 0; dim < 3; dim++) {
-    for (int ineed = 0; ineed < 2*need[dim]; ineed++) {
+    for (ineed = 0; ineed < 2*need[dim]; ineed++) {
       pbc_flag[iswap] = 0;
       pbc[iswap][0] = pbc[iswap][1] = pbc[iswap][2] =
 	pbc[iswap][3] = pbc[iswap][4] = pbc[iswap][5] = 0;
-
+      
       if (ineed % 2 == 0) {
 	sendproc[iswap] = procneigh[dim][0];
 	recvproc[iswap] = procneigh[dim][1];
-	if (ineed < 2) slablo[iswap] = -BIG;
+	if (style == SINGLE) {
-	else slablo[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
+	  if (ineed < 2) slablo[iswap] = -BIG;
-	slabhi[iswap] = sublo[dim] + cutghost[dim];
+	  else slablo[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
+	  slabhi[iswap] = sublo[dim] + cutghost[dim];
+	} else {
+	  for (i = 1; i <= ntypes; i++) {
+	    if (ineed < 2) multilo[iswap][i] = -BIG;
+	    else multilo[iswap][i] = 0.5 * (sublo[dim] + subhi[dim]);
+	    multihi[iswap][i] = sublo[dim] + cutghostmulti[i][dim];
+	  }
+	}
 	if (myloc[dim] == 0) {
-	  if (periodicity[dim] == 0)
+	  if (periodicity[dim] == 0) {
-	    slabhi[iswap] = slablo[iswap] - 1.0;
+	    if (style == SINGLE) slabhi[iswap] = slablo[iswap] - 1.0;
-	  else {
+	    else 
+	      for (i = 1; i <= ntypes; i++)
+		multihi[iswap][i] = multilo[iswap][i] - 1.0;
+	  } else {
 	    pbc_flag[iswap] = 1;
 	    pbc[iswap][dim] = 1;
 	    if (triclinic) {
@@ -262,16 +320,28 @@ void Comm::setup()
 	    }
 	  }
 	}
+	
       } else {
 	sendproc[iswap] = procneigh[dim][1];
 	recvproc[iswap] = procneigh[dim][0];
-	slablo[iswap] = subhi[dim] - cutghost[dim];
+	if (style == SINGLE) {
-	if (ineed < 2) slabhi[iswap] = BIG;
+	  slablo[iswap] = subhi[dim] - cutghost[dim];
-	else slabhi[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
+	  if (ineed < 2) slabhi[iswap] = BIG;
+	  else slabhi[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
+	} else {
+	  for (i = 1; i <= ntypes; i++) {
+	    multilo[iswap][i] = subhi[dim] - cutghostmulti[i][dim];
+	    if (ineed < 2) multihi[iswap][i] = BIG;
+	    else multihi[iswap][i] = 0.5 * (sublo[dim] + subhi[dim]);
+	  }
+	}
 	if (myloc[dim] == procgrid[dim]-1) {
-	  if (periodicity[dim] == 0)
+	  if (periodicity[dim] == 0) {
-	    slabhi[iswap] = slablo[iswap] - 1.0;
+	    if (style == SINGLE) slabhi[iswap] = slablo[iswap] - 1.0;
-	  else {
+	    else
+	      for (i = 1; i <= ntypes; i++)
+		multihi[iswap][i] = multilo[iswap][i] - 1.0;
+	  } else {
 	    pbc_flag[iswap] = 1;
 	    pbc[iswap][dim] = -1;
 	    if (triclinic) {
@@ -281,7 +351,7 @@ void Comm::setup()
 	  }
 	}
       }
-
+      
       iswap++;
     }
   }
@@ -509,10 +579,11 @@ void Comm::exchange()
 
 void Comm::borders()
 {
-  int i,n,iswap,dim,ineed,maxneed,nsend,nrecv,nfirst,nlast,smax,rmax;
+  int i,n,itype,iswap,dim,ineed,maxneed,nsend,nrecv,nfirst,nlast,smax,rmax;
   double lo,hi;
+  int *type;
   double **x;
-  double *buf;
+  double *buf,*mlo,*mhi;
   MPI_Request request;
   MPI_Status status;
   AtomVec *avec = atom->avec;
@@ -538,20 +609,36 @@ void Comm::borders()
       //   for later swaps in a dim, only check newly arrived ghosts
       // store sent atom indices in list for use in future timesteps
 
-      lo = slablo[iswap];
-      hi = slabhi[iswap];
       x = atom->x;
+      if (style == SINGLE) {
+	lo = slablo[iswap];
+	hi = slabhi[iswap];
+      } else {
+	type = atom->type;
+	mlo = multilo[iswap];
+	mhi = multihi[iswap];
+      }
       if (ineed % 2 == 0) {
 	nfirst = nlast;
 	nlast = atom->nlocal + atom->nghost;
       }
 
       nsend = 0;
-      for (i = nfirst; i < nlast; i++)
+      if (style == SINGLE) {
-	if (x[i][dim] >= lo && x[i][dim] <= hi) {
+	for (i = nfirst; i < nlast; i++)
-	  if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
+	  if (x[i][dim] >= lo && x[i][dim] <= hi) {
-	  sendlist[iswap][nsend++] = i;
+	    if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
+	    sendlist[iswap][nsend++] = i;
+	  }
+      } else {
+	for (i = nfirst; i < nlast; i++) {
+	  itype = type[i];
+	  if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) {
+	    if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
+	    sendlist[iswap][nsend++] = i;
+	  }
 	}
+      }
       
       // pack up list of border atoms
 
@@ -817,6 +904,277 @@ void Comm::reverse_comm_compute(Compute *compute)
   }
 }
 
+/* ----------------------------------------------------------------------
+   communicate atoms to new owning procs via irregular communication
+   unlike exchange(), allows for atoms to move arbitrary distances
+   first setup irregular comm pattern, then invoke it
+------------------------------------------------------------------------- */
+
+void Comm::irregular()
+{
+  // clear global->local map since atoms move to new procs
+
+  if (map_style) atom->map_clear();
+
+  // subbox bounds for orthogonal or triclinic
+
+  double *sublo,*subhi;
+  if (triclinic == 0) {
+    sublo = domain->sublo;
+    subhi = domain->subhi;
+  } else {
+    sublo = domain->sublo_lamda;
+    subhi = domain->subhi_lamda;
+  }
+
+  // loop over atoms, flag any that are not in my sub-box
+  // fill buffer with atoms leaving my box, using < and >=
+  // when atom is deleted, fill it in with last atom
+
+  AtomVec *avec = atom->avec;
+  double **x = atom->x;
+  int nlocal = atom->nlocal;
+
+  int nsend = 0;
+  int nsendatom = 0;
+  int *sizes = new int[nlocal];
+  int *proclist = new int[nlocal];
+
+  int i = 0;
+  while (i < nlocal) {
+    if (x[i][0] < sublo[0] || x[i][0] >= subhi[0] ||
+	x[i][1] < sublo[1] || x[i][1] >= subhi[1] ||
+	x[i][2] < sublo[2] || x[i][2] >= subhi[2]) {
+      if (nsend > maxsend) grow_send(nsend,1);
+      sizes[nsendatom] = avec->pack_exchange(i,&buf_send[nsend]);
+      nsend += sizes[nsendatom];
+      proclist[nsendatom] = 0;
+      nsendatom++;
+      avec->copy(nlocal-1,i);
+      nlocal--;
+    } else i++;
+  }
+  atom->nlocal = nlocal;
+
+  // create irregular communication plan, perform comm, destroy plan
+  // returned nrecv = size of buffer needed for incoming atoms
+
+  int nrecv;
+  Plan *plan = irregular_create(nsendatom,sizes,proclist,&nrecv);
+  if (nrecv > maxrecv) grow_recv(nrecv);
+  irregular_perform(plan,buf_send,sizes,buf_recv);
+  irregular_destroy(plan);
+
+  delete [] sizes;
+  delete [] proclist;
+
+  // add received atoms to my list
+
+  int m = 0;
+  while (m < nrecv) m += avec->unpack_exchange(&buf_recv[m]);
+
+  // reset global->local map
+
+  if (map_style) atom->map_set();
+}
+
+/* ----------------------------------------------------------------------
+   create an irregular communication plan
+   n = # of atoms to send
+   sizes = # of doubles for each atom
+   proclist = proc to send each atom to
+   return nrecvsize = total # of doubles I will recv
+------------------------------------------------------------------------- */
+
+Comm::Plan *Comm::irregular_create(int n, int *sizes, int *proclist,
+				   int *nrecvsize)
+{
+  int i;
+
+  // allocate plan and work vectors
+
+  Plan *plan = (struct Plan *) memory->smalloc(sizeof(Plan),"comm:plan");
+  /*
+  int *list = new int[nprocs];
+  int *counts = new int[nprocs];
+
+  // nrecv = # of messages I receive
+
+  for (i = 0; i < nprocs; i++) {
+    list[i] = 0;
+    counts[i] = 1;
+  }
+  for (i = 0; i < n; i++) list[proclist[i]] = 1;
+
+  int nrecv;
+  MPI_Reduce_scatter(list,&nrecv,counts,MPI_INT,MPI_SUM,world);
+
+  // storage for recv info
+
+  int *procs_from = new int[nrecv];
+  int *lengths_from = new int[nrecv];
+  MPI_Request *request = new MPI_Request[nrecv];
+  MPI_Status *status = new MPI_Status[nrecv];
+
+  // nsend = # of messages I send
+
+  for (i = 0; i < nprocs; i++) list[i] = 0;
+  for (i = 0; i < n; i++) list[proclist[i]]++;
+
+  int nsend = 0;
+  for (i = 0; i < nprocs; i++)
+    if (list[i] > 0) nsend++;
+  if (nself) nsend--;
+
+  // storage for send info
+
+  int *procs_to = new int[nsend];
+  int *lengths_to = new int[nsend];
+  int *indices_to = new int[n];
+
+  // set send info in procs_to and lengths_to
+  // each proc begins with iproc > me, and continues until iproc = me
+  // store pointer to sends in list for later use in setting indices_to
+
+  int iproc = me;
+  int isend = 0;
+  for (i = 0; i < nprocs; i++) {
+    iproc++;
+    if (iproc == nprocs) iproc = 0;
+    if (list[iproc] > 0) {
+      procs_to[isend] = iproc;
+      lengths_to[isend] = list[iproc];
+      list[iproc] = isend;
+      isend++;
+    }
+  }
+
+  // tell receivers how many I'm sending
+  // sendmax = largest # of datums I send
+
+  int sendmax = 0;
+  for (i = 0; i < nsend; i++) {
+    MPI_Send(&lengths_to[i],1,MPI_INT,procs_to[i],0,world);
+    sendmax = MAX(sendmax,lengths_to[i]);
+  }
+
+  // receive sizes and sources of incoming data
+  // nrecvsize = total # of datums I recv
+
+  *nrecvsize = 0;
+  for (i = 0; i < nrecv; i++) {
+    MPI_Recv(&lengths_from[i],1,MPI_INT,MPI_ANY_SOURCE,0,world,status);
+    procs_from[i] = status->MPI_SOURCE;
+    *nrecvsize += lengths_from[i];
+  }
+
+  // barrier to insure all my MPI_ANY_SOURCE messages are received
+  // else some procs could proceed to comm_do and start sending to me
+
+  MPI_Barrier(world);
+
+  // setup indices_to
+  // counts = current offset into indices_to for each proc I send to
+
+  counts[0] = 0;
+  for (i = 1; i < nsend; i++) counts[i] = counts[i-1] + lengths_to[i-1];
+
+  for (i = 0; i < n; i++) {
+    isend = list[proclist[i]];
+    indices_to[counts[isend]++] = i;
+  }
+
+  // free work vectors
+
+  delete [] counts;
+  delete [] list;
+    
+  // initialize plan and return it
+
+  plan->nsend = nsend;
+  plan->nrecv = nrecv;
+  plan->sendmax = sendmax;
+
+  plan->procs_to = procs_to;
+  plan->lengths_to = lengths_to;
+  plan->indices_to = indices_to;
+  plan->procs_from = procs_from;
+  plan->lengths_from = lengths_from;
+  plan->request = request;
+  plan->status = status;
+  */
+  return plan;
+}
+
+/* ----------------------------------------------------------------------
+   perform irregular communication
+   sendbuf = list of atoms to send
+   sizes = # of doubles for each atom
+   recvbuf = received atoms
+------------------------------------------------------------------------- */
+
+void Comm::irregular_perform(Plan *plan, double *sendbuf, int *sizes,
+			     double *recvbuf)
+{
+  int i,m;
+
+  // post all receives
+
+  int recv_offset = 0;
+  for (int irecv = 0; irecv < plan->nrecv; irecv++) {
+    MPI_Irecv(&recvbuf[recv_offset],plan->length_from[irecv],MPI_DOUBLE,
+	      plan->proc_from[irecv],0,world,&plan->request[irecv]);
+    recv_offset += plan->length_from[irecv];
+  }
+
+  // allocate buf for largest send
+
+  double *buf = (double *) memory->smalloc(plan->sendmax*sizeof(double),
+					   "comm::irregular");
+
+  // send each message
+  // pack buf with list of datums (datum = one atom)
+  // m = index of datum in sendbuf
+
+  int send_offset;
+  int ndatum = 0;
+  for (int isend = 0; isend < plan->nsend; isend++) {
+    send_offset = 0;
+    for (i = 0; i < plan->datum_send[isend]; i++) {
+      m = plan->index_send[ndatum++];
+      memcpy(&buf[send_offset],&sendbuf[plan->offset_send[m]],
+	     sizes[m]*sizeof(double));
+      send_offset += sizes[m];
+    }
+    MPI_Send(buf,plan->length_send[isend],MPI_DOUBLE,
+	     plan->proc_send[isend],0,world);
+  }       
+
+  // free temporary send buffer
+
+  memory->sfree(buf);
+
+  // wait on all incoming messages
+
+  if (plan->nrecv) MPI_Waitall(plan->nrecv,plan->request,plan->status);
+}
+
+/* ----------------------------------------------------------------------
+   destroy an irregular communication plan
+------------------------------------------------------------------------- */
+
+void Comm::irregular_destroy(Plan *plan)
+{
+  delete [] plan->proc_send;
+  delete [] plan->length_send;
+  delete [] plan->datum_send;
+  delete [] plan->index_send;
+  delete [] plan->offset_send;
+  delete [] plan->proc_from;
+  delete [] plan->length_from;
+  memory->sfree(plan);
+}
+
 /* ----------------------------------------------------------------------
    assign nprocs to 3d xprd,yprd,zprd box so as to minimize surface area 
    area = surface area of each of 3 faces of simulation box
@@ -937,25 +1295,29 @@ void Comm::grow_list(int iswap, int n)
    realloc the buffers needed for swaps 
 ------------------------------------------------------------------------- */
 
-void Comm::grow_swap()
+void Comm::grow_swap(int n)
 {
   free_swap();
-  allocate_swap(nswap);
+  allocate_swap(n);
-
+  if (style == MULTI) {
-  sendlist = (int **) memory->srealloc(sendlist,nswap*sizeof(int *),
+    free_multi();
-				       "comm:sendlist");
+    allocate_multi(n);
-  maxsendlist = (int *) memory->srealloc(maxsendlist,nswap*sizeof(int),
-					 "comm:maxsendlist");
-  for (int i = maxswap; i < nswap; i++) {
-    maxsendlist[i] = BUFMIN;
-    sendlist[i] = (int *) memory->smalloc(BUFMIN*sizeof(int),
-					  "comm:sendlist[i]");
   }
-  maxswap = nswap;
+
+  sendlist = (int **)
+    memory->srealloc(sendlist,n*sizeof(int *),"comm:sendlist");
+  maxsendlist = (int *)
+    memory->srealloc(maxsendlist,n*sizeof(int),"comm:maxsendlist");
+  for (int i = maxswap; i < n; i++) {
+    maxsendlist[i] = BUFMIN;
+    sendlist[i] = (int *)
+      memory->smalloc(BUFMIN*sizeof(int),"comm:sendlist[i]");
+  }
+  maxswap = n;
 }
 
 /* ----------------------------------------------------------------------
-   initial allocation of swap info 
+   allocation of swap info 
 ------------------------------------------------------------------------- */
 
 void Comm::allocate_swap(int n)
@@ -974,6 +1336,17 @@ void Comm::allocate_swap(int n)
   pbc = (int **) memory->create_2d_int_array(n,6,"comm:pbc");
 }
 
+
+/* ----------------------------------------------------------------------
+   allocation of multi-type swap info
+------------------------------------------------------------------------- */
+
+void Comm::allocate_multi(int n)
+{
+  multilo = memory->create_2d_double_array(n,atom->ntypes+1,"comm:multilo");
+  multihi = memory->create_2d_double_array(n,atom->ntypes+1,"comm:multihi");
+}
+
 /* ----------------------------------------------------------------------
    free memory for swaps 
 ------------------------------------------------------------------------- */
@@ -994,6 +1367,29 @@ void Comm::free_swap()
   memory->destroy_2d_int_array(pbc);
 }
 
+/* ----------------------------------------------------------------------
+   free memory for multi-type swaps
+------------------------------------------------------------------------- */
+
+void Comm::free_multi()
+{
+  memory->destroy_2d_double_array(multilo);
+  memory->destroy_2d_double_array(multihi);
+}
+
+/* ----------------------------------------------------------------------
+   set communication style
+------------------------------------------------------------------------- */
+
+void Comm::set(int narg, char **arg)
+{
+  if (narg != 1) error->all("Illegal communicate command");
+
+  if (strcmp(arg[0],"single") == 0) style = SINGLE;
+  else if (strcmp(arg[0],"multi") == 0) style = MULTI;
+  else error->all("Illegal communicate command");
+}
+
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory 
 ------------------------------------------------------------------------- */

src/comm.h

@@ -21,6 +21,7 @@ namespace LAMMPS_NS {
 class Comm : protected Pointers {
  public:
   int me,nprocs;                    // proc info
+  int style;                        // single vs multi-type comm
   int procgrid[3];                  // assigned # of procs in each dim
   int user_procgrid[3];             // user request for procs in each dim
   int myloc[3];                     // which proc I am in each dim
@@ -49,6 +50,9 @@ class Comm : protected Pointers {
   void comm_compute(class Compute *);          // forward comm from a Compute
   void reverse_comm_compute(class Compute *);  // reverse comm from a Compute
 
+  void irregular();                 // irregular communication across all procs
+
+  void set(int, char **);           // set communication style
   int memory_usage();               // tally memory usage
 
  private:
@@ -60,6 +64,8 @@ class Comm : protected Pointers {
   int *size_reverse_send;           // # to send in each reverse comm
   int *size_reverse_recv;           // # to recv in each reverse comm
   double *slablo,*slabhi;           // bounds of slab to send at each swap
+  double **multilo,**multihi;       // bounds of slabs for multi-type swap
+  double **cutghostmulti;           // cutghost on a per-type basis
   int *pbc_flag;                    // general flag for sending atoms thru PBC
   int **pbc;                        // dimension flags for PBC adjustments
   int comm_x_only,comm_f_only;      // 1 if only exchange x,f in for/rev comm
@@ -74,6 +80,21 @@ class Comm : protected Pointers {
   int maxsend,maxrecv;              // current size of send/recv buffer
   int maxforward,maxreverse;        // max # of datums in forward/reverse comm
 
+  struct Plan {                // plan for irregular communication
+    int nsend;                 // # of messages to send
+    int nrecv;                 // # of messages to recv
+    int sendmax;               // # of doubles in largest send message
+    int *proc_send;            // procs to send to
+    int *length_send;          // # of doubles to send to each proc
+    int *datum_send;           // # of datums to send to each proc
+    int *index_send;           // list of datums to send to each proc
+    int *offset_send;          // where each datum starts in send buffer
+    int *proc_from;            // procs to recv from
+    int *length_from;          // # of doubles to recv from each proc
+    MPI_Request *request;      // MPI requests for posted recvs
+    MPI_Status *status;        // MPI statuses for WaitAll
+  };
+
   void procs2box();                 // map procs to 3d box
   void cross(double, double, double,
 	     double, double, double,
@@ -81,9 +102,15 @@ class Comm : protected Pointers {
   void grow_send(int,int);          // reallocate send buffer
   void grow_recv(int);              // free/allocate recv buffer
   void grow_list(int, int);         // reallocate one sendlist
-  void grow_swap();                 // grow swap arrays
+  void grow_swap(int);              // grow swap and multi arrays
   void allocate_swap(int);          // allocate swap arrays
+  void allocate_multi(int);         // allocate multi arrays
   void free_swap();                 // free swap arrays
+  void free_multi();                // free multi arrays
+
+  struct Plan *irregular_create(int, int *, int *, int *);
+  void irregular_perform(Plan *, double *, int *, double *);
+  void irregular_destroy(Plan *);
 };
 
 }

src/domain.cpp

@@ -116,12 +116,14 @@ void Domain::set_initial_box()
     if (yz != 0.0 && (!yperiodic || !zperiodic))
       error->all("Triclinic box must be periodic in skewed dimensions");
 
+    /*
     if (fabs(xy/(boxhi[1]-boxlo[1])) > 0.5)
       error->all("Triclinic box skew is too large");
     if (fabs(xz/(boxhi[2]-boxlo[2])) > 0.5)
       error->all("Triclinic box skew is too large");
     if (fabs(yz/(boxhi[2]-boxlo[2])) > 0.5)
       error->all("Triclinic box skew is too large");
+    */
   }
 
   // adjust box lo/hi for shrink-wrapped dims

src/fix_ave_spatial.cpp

@@ -12,7 +12,7 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing author: Pieter in't Veld (SNL)
+   Contributing author: Pieter in 't Veld (SNL)
 ------------------------------------------------------------------------- */
 
 #include "stdlib.h"

src/fix_ave_time.cpp

@@ -12,7 +12,7 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing author: Pieter in't Veld (SNL)
+   Contributing author: Pieter in 't Veld (SNL)
 ------------------------------------------------------------------------- */
 
 #include "stdlib.h"

src/input.cpp

@@ -406,6 +406,7 @@ int Input::execute_command()
   else if (!strcmp(command,"bond_coeff")) bond_coeff();
   else if (!strcmp(command,"bond_style")) bond_style();
   else if (!strcmp(command,"boundary")) boundary();
+  else if (!strcmp(command,"communicate")) communicate();
   else if (!strcmp(command,"compute")) compute();
   else if (!strcmp(command,"compute_modify")) compute_modify();
   else if (!strcmp(command,"dielectric")) dielectric();
@@ -701,6 +702,13 @@ void Input::boundary()
 
 /* ---------------------------------------------------------------------- */
 
+void Input::communicate()
+{
+  comm->set(narg,arg);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void Input::compute()
 {
   modify->add_compute(narg,arg);

src/input.h

@@ -67,6 +67,7 @@ class Input : protected Pointers {
   void bond_coeff();
   void bond_style();
   void boundary();
+  void communicate();
   void compute();
   void compute_modify();
   void dielectric();

src/lammps.cpp

@@ -329,7 +329,7 @@ void LAMMPS::destroy()
   delete output;
   delete modify;          // modify must come after output, force, update
                           //   since they delete fixes
-  delete atom;            // must come after modify
+  delete atom;            // must come after modify, neighbor
                           //   since fixes delete callbacks in atom
   delete timer;
 }

src/neigh_full.cpp

@@ -55,7 +55,8 @@ void Neighbor::full_nsq()
     ytmp = x[i][1];
     ztmp = x[i][2];
 
-    // loop over all atoms, owned and ghost, only skip i = j
+    // loop over all atoms, owned and ghost
+    // skip i = j
 
     for (j = 0; j < nall; j++) {
       if (i == j) continue;
@@ -83,7 +84,7 @@ void Neighbor::full_nsq()
 }
 
 /* ----------------------------------------------------------------------
-   binned search for all neighbors
+   binned neighbor list construction for all neighbors
    every neighbor pair appears in list of both atoms i and j
 ------------------------------------------------------------------------- */
 
@@ -125,10 +126,11 @@ void Neighbor::full_bin()
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
-    ibin = coord2bin(x[i]);
 
     // loop over all atoms in surrounding bins in stencil including self
-    // only skip i = j
+    // skip i = j
+
+    ibin = coord2bin(x[i]);
 
     for (k = 0; k < nstencil_full; k++) {
       for (j = binhead[ibin+stencil_full[k]]; j >= 0; j = bins[j]) {
@@ -157,3 +159,87 @@ void Neighbor::full_bin()
       error->one("Neighbor list overflow, boost neigh_modify one or page");
   }
 }
+
+/* ----------------------------------------------------------------------
+   binned neighbor list construction for all neighbors
+   multi-type stencil is itype dependent and is distance checked
+   every neighbor pair appears in list of both atoms i and j
+------------------------------------------------------------------------- */
+
+void Neighbor::full_bin_multi()
+{
+  int i,j,k,n,itype,jtype,ibin,which,ns;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *neighptr,*s;
+  double *cutsq,*distsq;
+
+  // bin local & ghost atoms
+
+  bin_atoms();
+
+  // loop over each atom, storing neighbors
+
+  double **x = atom->x;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int nlocal = atom->nlocal;
+  int nall = atom->nlocal + atom->nghost;
+  int molecular = atom->molecular;
+
+  int npage = 0;
+  int npnt = 0;
+
+  for (i = 0; i < nlocal; i++) {
+
+    if (pgsize - npnt < oneatom) {
+      npnt = 0;
+      npage++;
+      if (npage == maxpage_full) add_pages_full(npage);
+    }
+
+    neighptr = &pages_full[npage][npnt];
+    n = 0;
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+
+    // loop over all atoms in other bins in stencil, including self
+    // skip if i,j neighbor cutoff is less than bin distance
+    // skip i = j
+
+    ibin = coord2bin(x[i]);
+    s = stencil_full_multi[itype];
+    distsq = distsq_full_multi[itype];
+    cutsq = cutneighsq[itype];
+    ns = nstencil_full_multi[itype];
+    for (k = 0; k < ns; k++) {
+      for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
+	jtype = type[j];
+	if (cutsq[jtype] < distsq[k]) continue;
+	if (i == j) continue;
+	if (exclude && exclusion(i,j,type,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(i,j);
+	  else which = 0;
+	  if (which == 0) neighptr[n++] = j;
+	  else if (which > 0) neighptr[n++] = which*nall + j;
+	}
+      }
+    }
+
+    firstneigh_full[i] = neighptr;
+    numneigh_full[i] = n;
+    npnt += n;
+    if (npnt >= pgsize)
+      error->one("Neighbor list overflow, boost neigh_modify one or page");
+  }
+}

src/neigh_half.cpp

@@ -168,7 +168,7 @@ void Neighbor::half_nsq_newton()
 
 /* ----------------------------------------------------------------------
    binned neighbor list construction with partial Newton's 3rd law
-   each owned atom i checks own bin and surrounding bins in non-Newton stencil
+   each owned atom i checks own bin and other bins in stencil
    pair stored once if i,j are both owned and i < j
    pair stored by me if j is ghost (also stored by proc owning j)
 ------------------------------------------------------------------------- */
@@ -212,7 +212,7 @@ void Neighbor::half_bin_no_newton()
     ytmp = x[i][1];
     ztmp = x[i][2];
 
-    // loop over all atoms in surrounding bins in stencil including self
+    // loop over all atoms in other bins in stencil including self
     // only store pair if i < j
     // stores own/own pairs only once
     // stores own/ghost pairs on both procs
@@ -249,7 +249,7 @@ void Neighbor::half_bin_no_newton()
 
 /* ----------------------------------------------------------------------
    binned neighbor list construction with partial Newton's 3rd law
-   each owned atom i checks own bin and surrounding bins in non-Newton stencil
+   each owned atom i checks own bin and other bins in stencil
    multi-type stencil is itype dependent and is distance checked
    pair stored once if i,j are both owned and i < j
    pair stored by me if j is ghost (also stored by proc owning j)
@@ -257,10 +257,10 @@ void Neighbor::half_bin_no_newton()
 
 void Neighbor::half_bin_no_newton_multi()
 {
-  int i,j,k,n,itype,jtype,ibin,which,ms;
+  int i,j,k,n,itype,jtype,ibin,which,ns;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr,*s;
-  double *cut,*dist;
+  double *cutsq,*distsq;
 
   // bin local & ghost atoms
 
@@ -295,22 +295,22 @@ void Neighbor::half_bin_no_newton_multi()
     ytmp = x[i][1];
     ztmp = x[i][2];
 
-    // loop over all atoms in surrounding bins in stencil including self
+    // loop over all atoms in other bins in stencil including self
     // only store pair if i < j
     // skip if i,j neighbor cutoff is less than bin distance
     // stores own/own pairs only once
     // stores own/ghost pairs on both procs
 
     ibin = coord2bin(x[i]);
-    s = mstencil[itype];
+    s = stencil_multi[itype];
-    dist = mdist[itype];
+    distsq = distsq_multi[itype];
-    cut = cutneigh[itype];
+    cutsq = cutneighsq[itype];
-    ms = mstencils[itype];
+    ns = nstencil_multi[itype];
-    for (k = 0; k < ms; k++) {
+    for (k = 0; k < ns; k++) {
       for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
 	if (j <= i) continue;
 	jtype = type[j];
-	if (cut[jtype] < dist[k]) continue;
+	if (cutsq[jtype] < distsq[k]) continue;
 	if (exclude && exclusion(i,j,type,mask,molecule)) continue;
 
 	delx = xtmp - x[j][0];
@@ -446,10 +446,10 @@ void Neighbor::half_bin_newton()
 
 void Neighbor::half_bin_newton_multi()
 {
-  int i,j,k,n,itype,jtype,ibin,which,ms;
+  int i,j,k,n,itype,jtype,ibin,which,ns;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr,*s;
-  double *cut,*dist;
+  double *cutsq,*distsq;
 
   // bin local & ghost atoms
 
@@ -515,14 +515,14 @@ void Neighbor::half_bin_newton_multi()
     // skip if i,j neighbor cutoff is less than bin distance
 
     ibin = coord2bin(x[i]);
-    s = mstencil[itype];
+    s = stencil_multi[itype];
-    dist = mdist[itype];
+    distsq = distsq_multi[itype];
-    cut = cutneigh[itype];
+    cutsq = cutneighsq[itype];
-    ms = mstencils[itype];
+    ns = nstencil_multi[itype];
-    for (k = 0; k < ms; k++) {
+    for (k = 0; k < ns; k++) {
       for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
 	jtype = type[j];
-	if (cut[jtype] < dist[k]) continue;
+	if (cutsq[jtype] < distsq[k]) continue;
 	if (exclude && exclusion(i,j,type,mask,molecule)) continue;
 
 	delx = xtmp - x[j][0];
@@ -637,10 +637,10 @@ void Neighbor::half_bin_newton_tri()
 
 void Neighbor::half_bin_newton_multi_tri()
 {
-  int i,j,k,n,itype,jtype,ibin,which,ms;
+  int i,j,k,n,itype,jtype,ibin,which,ns;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr,*s;
-  double *cut,*dist;
+  double *cutsq,*distsq;
 
   // bin local & ghost atoms
 
@@ -681,14 +681,14 @@ void Neighbor::half_bin_newton_multi_tri()
     // pairs for atoms j below i are excluded
 
     ibin = coord2bin(x[i]);
-    s = mstencil[itype];
+    s = stencil_multi[itype];
-    dist = mdist[itype];
+    distsq = distsq_multi[itype];
-    cut = cutneigh[itype];
+    cutsq = cutneighsq[itype];
-    ms = mstencils[itype];
+    ns = nstencil_multi[itype];
-    for (k = 0; k < ms; k++) {
+    for (k = 0; k < ns; k++) {
       for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
 	jtype = type[j];
-	if (cut[jtype] < dist[k]) continue;
+	if (cutsq[jtype] < distsq[k]) continue;
 	if (x[j][2] < ztmp) continue;
 	if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
 	if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue;

src/neigh_stencil.cpp

@@ -0,0 +1,461 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "neighbor.h"
+#include "atom.h"
+#include "memory.h"
+
+using namespace LAMMPS_NS;
+
+enum{NSQ,BIN,MULTI};     // also in neighbor.cpp
+
+/* ----------------------------------------------------------------------
+   routines to create a stencil = list of bin offsets
+   stencil = bins whose closest corner to central bin is within cutoff
+   sx,sy,sz = bin bounds = furthest the stencil could possibly extend
+   3d creates xyz stencil, 2d creates xy stencil
+   for half neigh list with partial Newton:
+     stencil is all surrounding bins
+     stencil includes self
+   for half neigh list with full Newton:
+     stencil is bins to the "upper right" of central bin
+     stencil does not include self
+   for half neigh list with triclinic:
+     stencil is all bins in z-plane of self and above, but not below
+     in 2d is all bins in y-plane of self and above, but not below
+     stencil includes self
+   for full neigh list:
+     stencil is all surrounding bins including self
+     regardless of newton on/off or triclinic
+   for multi:
+     create one stencil for each atom type
+     stencil is same bin bounds as newton on/off, triclinic, half/full
+     cutoff is not cutneighmaxsq, but max cutoff for that atom type
+------------------------------------------------------------------------- */
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_allocate(int sx, int sy, int sz)
+{
+  int i;
+
+  int nmax = (2*sz+1) * (2*sy+1) * (2*sx+1);
+
+  if (half) {
+    if (style == BIN) {
+      if (nmax > maxstencil) {
+	maxstencil = nmax;
+	memory->sfree(stencil);
+	stencil = (int *) memory->smalloc(nmax*sizeof(int),"neigh:stencil");
+      }
+
+    } else {
+      int n = atom->ntypes;
+      if (nstencil_multi == NULL) {
+	maxstencil_multi = 0;
+	nstencil_multi = new int[n+1];
+	stencil_multi = new int*[n+1];
+	distsq_multi = new double*[n+1];
+	for (i = 1; i <= n; i++) {
+	  nstencil_multi[i] = 0;
+	  stencil_multi[i] = NULL;
+	  distsq_multi[i] = NULL;
+	}
+      }
+      if (nmax > maxstencil_multi) {
+	maxstencil_multi = nmax;
+	for (int i = 1; i <= n; i++) {
+	  memory->sfree(stencil_multi[i]);
+	  memory->sfree(distsq_multi[i]);
+	  stencil_multi[i] = (int *)
+	    memory->smalloc(nmax*sizeof(int),"neigh:stencil_multi");
+	  distsq_multi[i] = (double *) memory->smalloc(nmax*sizeof(double),
+						       "neigh:distsq_multi");
+	}
+      }
+    }
+  }
+
+  if (full) {
+    if (style == BIN) {
+      if (nmax > maxstencil_full) {
+	maxstencil_full = nmax;
+	memory->sfree(stencil_full);
+	stencil_full = (int *) memory->smalloc(nmax*sizeof(int),
+					       "neigh:stencil_full");
+      }
+
+    } else {
+      int n = atom->ntypes;
+      if (nstencil_full_multi == NULL) {
+	maxstencil_full_multi = 0;
+	nstencil_full_multi = new int[n+1];
+	stencil_full_multi = new int*[n+1];
+	distsq_full_multi = new double*[n+1];
+	for (i = 1; i <= n; i++) {
+	  nstencil_full_multi[i] = 0;
+	  stencil_full_multi[i] = NULL;
+	  distsq_full_multi[i] = NULL;
+	}
+      }
+      if (nmax > maxstencil_full_multi) {
+	maxstencil_full_multi = nmax;
+	for (int i = 1; i <= n; i++) {
+	  memory->sfree(stencil_full_multi[i]);
+	  memory->sfree(distsq_full_multi[i]);
+	  stencil_full_multi[i] = (int *)
+	    memory->smalloc(nmax*sizeof(int),"neigh:stencil_full_multi");
+	  distsq_full_multi[i] = 
+	    (double *) memory->smalloc(nmax*sizeof(double),
+				       "neigh:distsq_full_multi");
+	}
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_none(int sx, int sy, int sz) {}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_3d_no_newton(int sx, int sy, int sz)
+{
+  int i,j,k;
+  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)
+	  stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_3d_no_newton_multi(int sx, int sy, int sz)
+{
+  int i,j,k,n;
+  double rsq,typesq;
+  int *s;
+  double *distsq;
+
+  int ntypes = atom->ntypes;
+  for (int itype = 1; itype <= ntypes; itype++) {
+    typesq = cuttypesq[itype];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    n = 0;
+    for (k = -sz; k <= sz; k++)
+      for (j = -sy; j <= sy; j++)
+	for (i = -sx; i <= sx; i++) {
+	  rsq = bin_distance(i,j,k);
+	  if (rsq < typesq) {
+	    distsq[n] = rsq;
+	    s[n++] = k*mbiny*mbinx + j*mbinx + i;
+	  }
+	}
+    nstencil_multi[itype] = n;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_3d_newton(int sx, int sy, int sz)
+{
+  int i,j,k;
+  nstencil = 0;
+
+  for (k = 0; k <= sz; k++)
+    for (j = -sy; j <= sy; j++)
+      for (i = -sx; i <= sx; i++)
+	if (k > 0 || j > 0 || (j == 0 && i > 0))
+	  if (bin_distance(i,j,k) < cutneighmaxsq)
+	    stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_3d_newton_multi(int sx, int sy, int sz)
+{
+  int i,j,k,n;
+  double rsq,typesq;
+  int *s;
+  double *distsq;
+
+  int ntypes = atom->ntypes;
+  for (int itype = 1; itype <= ntypes; itype++) {
+    typesq = cuttypesq[itype];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    n = 0;
+    for (k = 0; k <= sz; k++)
+      for (j = -sy; j <= sy; j++)
+	for (i = -sx; i <= sx; i++)
+	  if (k > 0 || j > 0 || (j == 0 && i > 0)) {
+	    rsq = bin_distance(i,j,k);
+	    if (rsq < typesq) {
+	      distsq[n] = rsq;
+	      s[n++] = k*mbiny*mbinx + j*mbinx + i;
+	    }
+	  }
+    nstencil_multi[itype] = n;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_3d_newton_tri(int sx, int sy, int sz)
+{
+  int i,j,k;
+  nstencil = 0;
+
+  for (k = 0; k <= sz; k++)
+    for (j = -sy; j <= sy; j++)
+      for (i = -sx; i <= sx; i++)
+	if (bin_distance(i,j,k) < cutneighmaxsq)
+	  stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_3d_newton_multi_tri(int sx, int sy, int sz)
+{
+  int i,j,k,n;
+  double rsq,typesq;
+  int *s;
+  double *distsq;
+
+  int ntypes = atom->ntypes;
+  for (int itype = 1; itype <= ntypes; itype++) {
+    typesq = cuttypesq[itype];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    n = 0;
+    for (k = 0; k <= sz; k++)
+      for (j = -sy; j <= sy; j++)
+	for (i = -sx; i <= sx; i++) {
+	  rsq = bin_distance(i,j,k);
+	  if (rsq < typesq) {
+	    distsq[n] = rsq;
+	    s[n++] = k*mbiny*mbinx + j*mbinx + i;
+	  }
+	}
+    nstencil_multi[itype] = n;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_2d_no_newton(int sx, int sy, int sz)
+{
+  int i,j;
+  nstencil = 0;
+
+  for (j = -sy; j <= sy; j++)
+    for (i = -sx; i <= sx; i++)
+      if (bin_distance(i,j,0) < cutneighmaxsq)
+	stencil[nstencil++] = j*mbinx + i;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_2d_no_newton_multi(int sx, int sy, int sz)
+{
+  int i,j,n;
+  double rsq,typesq;
+  int *s;
+  double *distsq;
+
+  int ntypes = atom->ntypes;
+  for (int itype = 1; itype <= ntypes; itype++) {
+    typesq = cuttypesq[itype];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    n = 0;
+    for (j = -sy; j <= sy; j++)
+      for (i = -sx; i <= sx; i++) {
+	rsq = bin_distance(i,j,0);
+	if (rsq < typesq) {
+	  distsq[n] = rsq;
+	  s[n++] = j*mbinx + i;
+	}
+      }
+    nstencil_multi[itype] = n;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_2d_newton(int sx, int sy, int sz)
+{
+  int i,j;
+  nstencil = 0;
+
+  for (j = 0; j <= sy; j++)
+    for (i = -sx; i <= sx; i++)
+      if (j > 0 || (j == 0 && i > 0))
+	if (bin_distance(i,j,0) < cutneighmaxsq)
+	  stencil[nstencil++] = j*mbinx + i;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_2d_newton_multi(int sx, int sy, int sz)
+{
+  int i,j,n;
+  double rsq,typesq;
+  int *s;
+  double *distsq;
+
+  int ntypes = atom->ntypes;
+  for (int itype = 1; itype <= ntypes; itype++) {
+    typesq = cuttypesq[itype];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    n = 0;
+    for (j = 0; j <= sy; j++)
+      for (i = -sx; i <= sx; i++)
+	if (j > 0 || (j == 0 && i > 0)) {
+	  rsq = bin_distance(i,j,0);
+	  if (rsq < typesq) {
+	    distsq[n] = rsq;
+	    s[n++] = j*mbinx + i;
+	  }
+	}
+    nstencil_multi[itype] = n;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_2d_newton_tri(int sx, int sy, int sz)
+{
+  int i,j;
+  nstencil = 0;
+
+  for (j = 0; j <= sy; j++)
+    for (i = -sx; i <= sx; i++)
+      if (bin_distance(i,j,0) < cutneighmaxsq)
+	stencil[nstencil++] = j*mbinx + i;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_half_2d_newton_multi_tri(int sx, int sy, int sz)
+{
+  int i,j,n;
+  double rsq,typesq;
+  int *s;
+  double *distsq;
+
+  int ntypes = atom->ntypes;
+  for (int itype = 1; itype <= ntypes; itype++) {
+    typesq = cuttypesq[itype];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    n = 0;
+    for (j = 0; j <= sy; j++)
+      for (i = -sx; i <= sx; i++) {
+	rsq = bin_distance(i,j,0);
+	if (rsq < typesq) {
+	  distsq[n] = rsq;
+	  s[n++] = j*mbinx + i;
+	}
+      }
+    nstencil_multi[itype] = n;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_full_3d(int sx, int sy, int sz)
+{
+  int i,j,k;
+  nstencil_full = 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)
+	  stencil_full[nstencil_full++] = k*mbiny*mbinx + j*mbinx + i;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_full_3d_multi(int sx, int sy, int sz)
+{
+  int i,j,k,n;
+  double rsq,typesq;
+  int *s;
+  double *distsq;
+
+  int ntypes = atom->ntypes;
+  for (int itype = 1; itype <= ntypes; itype++) {
+    typesq = cuttypesq[itype];
+    s = stencil_full_multi[itype];
+    distsq = distsq_full_multi[itype];
+    n = 0;
+    for (k = -sz; k <= sz; k++)
+      for (j = -sy; j <= sy; j++)
+	for (i = -sx; i <= sx; i++) {
+	  rsq = bin_distance(i,j,k);
+	  if (rsq < typesq) {
+	    distsq[n] = rsq;
+	    s[n++] = k*mbiny*mbinx + j*mbinx + i;
+	  }
+	}
+    nstencil_full_multi[itype] = n;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_full_2d(int sx, int sy, int sz)
+{
+  int i,j;
+  nstencil_full = 0;
+
+  for (j = -sy; j <= sy; j++)
+    for (i = -sx; i <= sx; i++)
+      if (bin_distance(i,j,0) < cutneighmaxsq)
+	stencil_full[nstencil_full++] = j*mbinx + i;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Neighbor::stencil_full_2d_multi(int sx, int sy, int sz)
+{
+  int i,j,n;
+  double rsq,typesq;
+  int *s;
+  double *distsq;
+
+  int ntypes = atom->ntypes;
+  for (int itype = 1; itype <= ntypes; itype++) {
+    typesq = cuttypesq[itype];
+    s = stencil_full_multi[itype];
+    distsq = distsq_full_multi[itype];
+    n = 0;
+    for (j = -sy; j <= sy; j++)
+      for (i = -sx; i <= sx; i++) {
+	rsq = bin_distance(i,j,0);
+	if (rsq < typesq) {
+	  distsq[n] = rsq;
+	  s[n++] = j*mbinx + i;
+	}
+      }
+    nstencil_full_multi[itype] = n;
+  }
+}

src/neighbor.cpp

@@ -44,7 +44,7 @@ using namespace LAMMPS_NS;
 #define MIN(a,b) ((a) < (b) ? (a) : (b))
 #define MAX(a,b) ((a) > (b) ? (a) : (b))
 
-enum{NSQ,BIN,MULTI};
+enum{NSQ,BIN,MULTI};     // also in neigh_stencil.cpp
 
 /* ---------------------------------------------------------------------- */
 
@@ -62,8 +62,9 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
 
   maxlocal = 0;
     
-  cutneigh = NULL;
   cutneighsq = NULL;
+  cuttype = NULL;
+  cuttypesq = NULL;
   fixchecklist = NULL;
 
   // last neighbor info
@@ -90,27 +91,35 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
   maxbin = 0;
   bins = NULL;
 
-  maxstencil = 0;
+  nstencil = maxstencil = 0;
   stencil = NULL;
-  maxstencil_full = 0;
+  nstencil_full = maxstencil_full = 0;
   stencil_full = NULL;
 
-  //mstencil = NULL;
+  maxstencil_multi = 0;
-  
+  nstencil_multi = NULL;
+  stencil_multi = NULL;
+  distsq_multi = NULL;
+
+  maxstencil_full_multi = 0;
+  nstencil_full_multi = NULL;
+  stencil_full_multi = NULL;
+  distsq_full_multi = NULL;
+
   // half neighbor list info
 
   half = half_command = 0;
-  maxpage = 0;
   numneigh = NULL;
   firstneigh = NULL;
+  maxpage = 0;
   pages = NULL;
 
   // full neighbor list info
 
   full = 0;
-  maxpage_full = 0;
   numneigh_full = NULL;
   firstneigh_full = NULL;
+  maxpage_full = 0;
   pages_full = NULL;
 
   // shear history neighbor list info
@@ -118,19 +127,20 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
   fix_history = NULL;
   firsttouch = NULL;
   firstshear = NULL;
+  maxpage_history = 0;
   pages_touch = NULL;
   pages_shear = NULL;
 
   // multiple respa neighbor list info
 
   respa = 0;
-  maxpage_inner = 0;
-  maxpage_middle = 0;
   numneigh_inner = NULL;
   firstneigh_inner = NULL;
-  pages_inner = NULL;
   numneigh_middle = NULL;
   firstneigh_middle = NULL;
+  maxpage_inner = 0;
+  maxpage_middle = 0;
+  pages_inner = NULL;
   pages_middle = NULL;
 
   // bond list info
@@ -149,8 +159,10 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
 
 Neighbor::~Neighbor()
 {
-  memory->destroy_2d_double_array(cutneigh);
   memory->destroy_2d_double_array(cutneighsq);
+  delete [] cuttype;
+  delete [] cuttypesq;
+
   delete [] fixchecklist;
   memory->destroy_2d_double_array(xhold);
 
@@ -168,9 +180,29 @@ Neighbor::~Neighbor()
 
   memory->sfree(binhead);
   memory->sfree(bins);
+
   memory->sfree(stencil);
   memory->sfree(stencil_full);
 
+  if (nstencil_multi) {
+    for (int i = 1; i <= atom->ntypes; i++) {
+      memory->sfree(stencil_multi[i]);
+      memory->sfree(distsq_multi[i]);
+    }
+    delete [] nstencil_multi;
+    delete [] stencil_multi;
+    delete [] distsq_multi;
+  }
+  if (nstencil_full_multi) {
+    for (int i = 1; i <= atom->ntypes; i++) {
+      memory->sfree(stencil_full_multi[i]);
+      memory->sfree(distsq_full_multi[i]);
+    }
+    delete [] nstencil_full_multi;
+    delete [] stencil_full_multi;
+    delete [] distsq_full_multi;
+  }
+
   memory->destroy_2d_int_array(bondlist);
   memory->destroy_2d_int_array(anglelist);
   memory->destroy_2d_int_array(dihedrallist);
@@ -178,33 +210,22 @@ Neighbor::~Neighbor()
 
   memory->sfree(numneigh);
   memory->sfree(firstneigh);
-  for (int i = 0; i < maxpage; i++) memory->sfree(pages[i]);
+  clear_pages();
-  memory->sfree(pages);
 
   memory->sfree(numneigh_full);
   memory->sfree(firstneigh_full);
-  for (int i = 0; i < maxpage_full; i++) memory->sfree(pages_full[i]);
+  clear_pages_full();
-  memory->sfree(pages_full);
 
-  if (fix_history) {
+  memory->sfree(firsttouch);
-    memory->sfree(firsttouch);
+  memory->sfree(firstshear);
-    memory->sfree(firstshear);
+  clear_pages_history();
-    for (int i = 0; i < maxpage; i++) memory->sfree(pages_touch[i]);
-    for (int i = 0; i < maxpage; i++) memory->sfree(pages_shear[i]);
-    memory->sfree(pages_touch);
-    memory->sfree(pages_shear);
-  }
 
-  if (respa) {
+  memory->sfree(numneigh_inner);
-    memory->sfree(numneigh_inner);
+  memory->sfree(firstneigh_inner);
-    memory->sfree(firstneigh_inner);
+  memory->sfree(numneigh_middle);
-    for (int i = 0; i < maxpage_inner; i++) memory->sfree(pages_inner[i]);
+  memory->sfree(firstneigh_middle);
-    memory->sfree(pages_inner);
+  clear_pages_inner();
-    memory->sfree(numneigh_middle);
+  clear_pages_middle();
-    memory->sfree(firstneigh_middle);
-    for (int i = 0; i < maxpage_middle; i++) memory->sfree(pages_middle[i]);
-    memory->sfree(pages_middle);
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -235,31 +256,34 @@ void Neighbor::init()
   triggersq = 0.25*skin*skin;
 
   n = atom->ntypes;
-  if (cutneigh == NULL) {
+  if (cutneighsq == NULL) {
-    cutneigh = memory->create_2d_double_array(n+1,n+1,"neigh:cutneigh");
     cutneighsq = memory->create_2d_double_array(n+1,n+1,"neigh:cutneighsq");
+    cuttype = new double[n+1];
+    cuttypesq = new double[n+1];
   }
 
-  double cutoff,delta;
+  double cutoff,delta,cut;
-  for (i = 1; i <= n; i++)
+  cutneighmin = BIG;
-    for (j = i; j <= n; j++) {
+  cutneighmax = 0.0;
+
+  for (i = 1; i <= n; i++) {
+    cuttype[i] = cuttypesq[i] = 0.0;
+    for (j = 1; j <= n; j++) {
       if (force->pair) cutoff = sqrt(force->pair->cutsq[i][j]);
       else cutoff = 0.0;
       if (cutoff > 0.0) delta = skin;
       else delta = 0.0;
-      cutneigh[i][j] = cutneigh[j][i] = cutoff + delta;
+      cut = cutoff + delta;
-      cutneighsq[i][j] = cutneighsq[j][i] = (cutoff+delta) * (cutoff+delta);
-    }
 
-  cutneighmin = BIG;
+      cutneighsq[i][j] = cut*cut;
-  cutneighmax = 0.0;
+      cuttype[i] = MAX(cuttype[i],cut);
-  cutghost = skin;
+      cuttypesq[i] = MAX(cuttypesq[i],cut*cut);
-  for (i = 1; i <= n; i++)
+      cutneighmin = MIN(cutneighmin,cut);
-    for (j = i; j <= n; j++) {
+      cutneighmax = MAX(cutneighmax,cut);
-      cutneighmin = MIN(cutneighmin,cutneigh[i][j]);
-      cutneighmax = MAX(cutneighmax,cutneigh[i][j]);
-      cutghost = MAX(cutghost,cutneigh[i][j]);
     }
+  }
+  cutghost = MAX(cutneighmax,skin);
+  cutneighmaxsq = cutneighmax * cutneighmax;
 
   // check other classes that can induce reneighboring in decide()
 
@@ -313,6 +337,7 @@ void Neighbor::init()
   if (dist_check == 0) {
     memory->destroy_2d_double_array(xhold);
     maxhold = 0;
+    xhold = NULL;
   }
 
   if (style == NSQ) {
@@ -321,6 +346,8 @@ void Neighbor::init()
     memory->sfree(stencil);
     memory->sfree(stencil_full);
     maxbin = maxhead = maxstencil = maxstencil_full = 0;
+    binhead = NULL;
+    bins = NULL;
   }
 
   // 1st time allocation of xhold and bins
@@ -332,7 +359,7 @@ void Neighbor::init()
     }
   }
 
-  if (style == BIN) {
+  if (style != NSQ) {
     if (maxbin == 0) {
       maxbin = atom->nmax;
       bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins");
@@ -384,15 +411,11 @@ void Neighbor::init()
   }
 
   // ------------------------------------------------------------------
-  // half and full pairwise neighbor lists
+  // neighbor list flags and memory allocation/deallocation
 
   // determine whether to build half and full lists
   // query pair,fix,compute for their requirements
 
-  maxlocal = atom->nmax;
-  int half_previous = half;
-  int full_previous = full;
-
   half_once = full_once = 0;
   half_every = full_every = 0;
   if (force->pair) half_every = force->pair->neigh_half_every;
@@ -415,89 +438,19 @@ void Neighbor::init()
   if (full_every || full_once) full = 1;
   half = 1;
 
-  // setup/delete memory for half and full lists
+  // determine whether to build granular history lists
-
+  // fix_history = granular shear history fix
-  if (half == 0 && half_previous) {
-    memory->sfree(numneigh);
-    memory->sfree(firstneigh);
-    for (i = 0; i < maxpage; i++) memory->sfree(pages[i]);
-    memory->sfree(pages);
-    pages = NULL;
-    maxpage = 0;
-  } else if (half && half_previous == 0) {
-    numneigh =
-      (int *) memory->smalloc(maxlocal*sizeof(int),"neigh:numneigh");
-    firstneigh =
-      (int **) memory->smalloc(maxlocal*sizeof(int *),"neigh:firstneigh");
-    add_pages(0);
-  } else if (half && half_previous) {
-    memory->sfree(numneigh);
-    memory->sfree(firstneigh);
-    numneigh =
-      (int *) memory->smalloc(maxlocal*sizeof(int),"neigh:numneigh");
-    firstneigh =
-      (int **) memory->smalloc(maxlocal*sizeof(int *),"neigh:firstneigh");
-  }
-
-  if (full == 0 && full_previous) {
-    memory->sfree(numneigh_full);
-    memory->sfree(firstneigh_full);
-    for (i = 0; i < maxpage_full; i++) memory->sfree(pages_full[i]);
-    memory->sfree(pages_full);
-    pages_full = NULL;
-    maxpage_full = 0;
-  } else if (full && full_previous == 0) {
-    numneigh_full =
-      (int *) memory->smalloc(maxlocal*sizeof(int),"neigh:numneigh_full");
-    firstneigh_full =
-      (int **) memory->smalloc(maxlocal*sizeof(int *),"neigh:firstneigh_full");
-    add_pages_full(0);
-  } else if (full && full_previous) {
-    memory->sfree(numneigh_full);
-    memory->sfree(firstneigh_full);
-    numneigh_full =
-      (int *) memory->smalloc(maxlocal*sizeof(int),"neigh:numneigh_full");
-    firstneigh_full =
-      (int **) memory->smalloc(maxlocal*sizeof(int *),"neigh:firstneigh_full");
-  }
-
-  // setup/delete memory for shear history neighbor lists
-  // fix_history = granular shear history fix if it exists
   
-  FixShearHistory *fix_previous = fix_history;
   fix_history = NULL;
   if (force->pair_match("gran/history") || force->pair_match("gran/hertzian"))
     for (i = 0; i < modify->nfix; i++)
       if (strcmp(modify->fix[i]->style,"SHEAR_HISTORY") == 0) 
 	fix_history = (FixShearHistory *) modify->fix[i];
 
-  if (fix_history == NULL && fix_previous) {
+  // determine whether to build extra rRESPA lists
-    memory->sfree(firsttouch);
+  // respa = 1,2 if rRESPA requires inner,middle neighbor lists
-    memory->sfree(firstshear);
-    for (i = 0; i < maxpage; i++) memory->sfree(pages_touch[i]);
-    for (i = 0; i < maxpage; i++) memory->sfree(pages_shear[i]);
-    memory->sfree(pages_touch);
-    memory->sfree(pages_shear);
-    pages_touch = NULL;
-    pages_shear = NULL;
-  } else if (fix_history && fix_previous == NULL) {
-    firsttouch = (int **) memory->smalloc(maxlocal*sizeof(int *),"firsttouch");
-    firstshear = (double **)
-      memory->smalloc(maxlocal*sizeof(double *),"firstshear");
-    add_pages_history(0);
-  } else if (fix_history && fix_previous >= 0) {
-    memory->sfree(firsttouch);
-    memory->sfree(firstshear);
-    firsttouch = (int **) memory->smalloc(maxlocal*sizeof(int *),"firsttouch");
-    firstshear = (double **)
-      memory->smalloc(maxlocal*sizeof(double *),"firstshear");
-  }
-
-  // setup/delete memory for rRESPA neighbor lists
-  // respa = 1 if rRESPA requires extra neighbor lists
   // set neighbor cutoffs for multiple lists
 
-  int respa_previous = respa;
   respa = 0;
   if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) {
     if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
@@ -507,51 +460,6 @@ void Neighbor::init()
   if (respa && half_every == 0)
     error->all("Cannot use rRESPA with full neighbor lists");
 
-  if (respa == 0 && respa_previous) {
-    memory->sfree(numneigh_inner);
-    memory->sfree(firstneigh_inner);
-    for (i = 0; i < maxpage_inner; i++) memory->sfree(pages_inner[i]);
-    memory->sfree(pages_inner);
-    pages_inner = NULL;
-    maxpage_inner = 0;
-    if (respa_previous == 2) {
-      memory->sfree(numneigh_middle);
-      memory->sfree(firstneigh_middle);
-      for (i = 0; i < maxpage_middle; i++) memory->sfree(pages_middle[i]);
-      memory->sfree(pages_middle);
-      pages_middle = NULL;
-      maxpage_middle = 0;
-    }
-  } else if (respa && respa_previous == 0) {
-    numneigh_inner = (int *) 
-      memory->smalloc(maxlocal*sizeof(int),"neigh:numneigh_inner");
-    firstneigh_inner = (int **) 
-      memory->smalloc(maxlocal*sizeof(int *),"neigh:firstneigh_inner");
-    add_pages_inner(0);
-    if (respa == 2) {
-      numneigh_middle = (int *) 
-	memory->smalloc(maxlocal*sizeof(int),"neigh:numneigh_middle");
-      firstneigh_middle = (int **) 
-	memory->smalloc(maxlocal*sizeof(int *),"neigh:firstneigh_middle");
-      add_pages_middle(0);
-    }
-  } else if (respa && respa_previous) {
-    memory->sfree(numneigh_inner);
-    memory->sfree(firstneigh_inner);
-    numneigh_inner = (int *) 
-      memory->smalloc(maxlocal*sizeof(int),"neigh:numneigh_inner");
-    firstneigh_inner = (int **) 
-      memory->smalloc(maxlocal*sizeof(int *),"neigh:firstneigh_inner");
-    if (respa == 2) {
-      memory->sfree(numneigh_middle);
-      memory->sfree(firstneigh_middle);
-      numneigh_middle = (int *) 
-	memory->smalloc(maxlocal*sizeof(int),"neigh:numneigh_middle");
-      firstneigh_middle = (int **) 
-	memory->smalloc(maxlocal*sizeof(int *),"neigh:firstneigh_middle");
-    }
-  }
-
   if (respa) {
     double *cut_respa = ((Respa *) update->integrate)->cutoff;
     cut_inner_sq = (cut_respa[1] + skin) * (cut_respa[1] + skin);
@@ -559,8 +467,46 @@ void Neighbor::init()
     cut_middle_inside_sq = (cut_respa[0] - skin) * (cut_respa[0] - skin);
   }
 
-  // set ptrs to correct half, full, multi, triclinic build functions
+  // zero atom-length numneigh and firstneigh arrays
-  // cannot combine granular and rRESPA
+  // will be (re)allocated on first build()
+
+  maxlocal = 0;
+
+  memory->sfree(numneigh);
+  memory->sfree(firstneigh);
+  memory->sfree(numneigh_full);
+  memory->sfree(firstneigh_full);
+  memory->sfree(firsttouch);
+  memory->sfree(firstshear);
+  memory->sfree(numneigh_inner);
+  memory->sfree(firstneigh_inner);
+  memory->sfree(numneigh_middle);
+  memory->sfree(firstneigh_middle);
+
+  numneigh = numneigh_full = numneigh_inner = numneigh_middle = NULL;
+  firstneigh = firstneigh_full = NULL;
+  firsttouch = NULL;
+  firstshear = NULL;
+  firstneigh_inner = firstneigh_middle = NULL;
+
+  // clear old neighbor lists if no longer needed (whether exist or not)
+
+  if (half == 0) clear_pages();
+  if (full == 0) clear_pages_full();
+  if (fix_history == NULL) clear_pages_history();
+  if (respa == 0) clear_pages_inner();
+  if (respa == 0 || respa == 1) clear_pages_middle();
+
+  // setup new neighbor lists
+  // only if they don't exist so memory will persist from run to run
+
+  if (half && pages == NULL) add_pages(0);
+  if (full && pages_full == NULL) add_pages_full(0);
+  if (fix_history && pages_touch == NULL) add_pages_history(0);
+  if (respa >= 1 && pages_inner == NULL) add_pages_inner(0);
+  if (respa == 2 && pages_middle == NULL) add_pages_middle(0);
+
+  // set ptrs to half/full/multi/triclinic build & stencil functions
 
   if (half) {
     if (atom->check_style("granular")) {
@@ -569,24 +515,54 @@ void Neighbor::init()
 	  half_build = &Neighbor::granular_nsq_no_newton;
 	else half_build = &Neighbor::granular_nsq_newton;
       } else if (style == BIN) {
-	if (force->newton_pair == 0) 
+	if (force->newton_pair == 0) {
 	  half_build = &Neighbor::granular_bin_no_newton;
-	else if (triclinic)
+	  if (force->dimension == 3)
+	    half_stencil = &Neighbor::stencil_half_3d_no_newton;
+	  else
+	    half_stencil = &Neighbor::stencil_half_2d_no_newton;
+	} else if (triclinic) {
 	  half_build = &Neighbor::granular_bin_newton_tri;
-	else half_build = &Neighbor::granular_bin_newton;
+	  if (force->dimension == 3)
+	    half_stencil = &Neighbor::stencil_half_3d_newton_tri;
+	  else
+	    half_stencil = &Neighbor::stencil_half_2d_newton_tri;
+	} else {
+	  half_build = &Neighbor::granular_bin_newton;
+	  if (force->dimension == 3)
+	    half_stencil = &Neighbor::stencil_half_3d_newton;
+	  else
+	    half_stencil = &Neighbor::stencil_half_2d_newton;
+	}
       } else error->all("Neighbor multi not allowed with granular");
+
     } else if (respa) {
       if (style == NSQ) {
 	if (force->newton_pair == 0) 
 	  half_build = &Neighbor::respa_nsq_no_newton;
 	else half_build = &Neighbor::respa_nsq_newton;
       } else if (style == BIN) {
-	if (force->newton_pair == 0) 
+	if (force->newton_pair == 0) {
 	  half_build = &Neighbor::respa_bin_no_newton;
-	else if (triclinic)
+	  if (force->dimension == 3)
+	    half_stencil = &Neighbor::stencil_half_3d_no_newton;
+	  else
+	    half_stencil = &Neighbor::stencil_half_2d_no_newton;
+	} else if (triclinic) {
 	  half_build = &Neighbor::respa_bin_newton_tri;
-	else half_build = &Neighbor::respa_bin_newton;
+	  if (force->dimension == 3)
+	    half_stencil = &Neighbor::stencil_half_3d_newton_tri;
+	  else
+	    half_stencil = &Neighbor::stencil_half_2d_newton_tri;
+	} else {
+	  half_build = &Neighbor::respa_bin_newton;
+	  if (force->dimension == 3)
+	    half_stencil = &Neighbor::stencil_half_3d_newton;
+	  else
+	    half_stencil = &Neighbor::stencil_half_2d_newton;
+	}
       } else error->all("Neighbor multi not allowed with rRESPA");
+
     } else {
       if (style == NSQ) {
 	if (force->newton_pair == 0) {
@@ -598,31 +574,84 @@ void Neighbor::init()
 	}
       } else if (style == BIN) {
 	if (force->newton_pair == 0) {
-	  if (full_every) half_build = &Neighbor::half_full_no_newton;
+	  if (full_every) {
-	  else half_build = &Neighbor::half_bin_no_newton;
+	    half_build = &Neighbor::half_full_no_newton;
+	    half_stencil = &Neighbor::stencil_none;
+	  } else {
+	    half_build = &Neighbor::half_bin_no_newton;
+	    if (force->dimension == 3)
+	      half_stencil = &Neighbor::stencil_half_3d_no_newton;
+	    else
+	      half_stencil = &Neighbor::stencil_half_2d_no_newton;
+	  }
 	} else {
-	  if (full_every) half_build = &Neighbor::half_full_newton;
+	  if (full_every) {
-	  else if (triclinic) half_build = &Neighbor::half_bin_newton_tri;
+	    half_build = &Neighbor::half_full_newton;
-	  else half_build = &Neighbor::half_bin_newton;
+	    half_stencil = &Neighbor::stencil_none;
+	  } else if (triclinic) {
+	    half_build = &Neighbor::half_bin_newton_tri;
+	    if (force->dimension == 3)
+	      half_stencil = &Neighbor::stencil_half_3d_newton_tri;
+	    else
+	      half_stencil = &Neighbor::stencil_half_2d_newton_tri;
+	  } else {
+	    half_build = &Neighbor::half_bin_newton;
+	    if (force->dimension == 3)
+	      half_stencil = &Neighbor::stencil_half_3d_newton;
+	    else
+	      half_stencil = &Neighbor::stencil_half_2d_newton;
+	  }
 	}
       } else if (style == MULTI) {
 	if (force->newton_pair == 0) {
-	  if (full_every) half_build = &Neighbor::half_full_no_newton;
+	  if (full_every) {
-	  else half_build = &Neighbor::half_bin_no_newton_multi;
+	    half_build = &Neighbor::half_full_no_newton;
+	    half_stencil = &Neighbor::stencil_none;
+	  } else {
+	    half_build = &Neighbor::half_bin_no_newton_multi;
+	    if (force->dimension == 3)
+	      half_stencil = &Neighbor::stencil_half_3d_no_newton_multi;
+	    else
+	      half_stencil = &Neighbor::stencil_half_2d_no_newton_multi;
+	  }
 	} else {
-	  if (full_every) half_build = &Neighbor::half_full_newton;
+	  if (full_every) {
-	  else if (triclinic)
+	    half_build = &Neighbor::half_full_newton;
+	    half_stencil = &Neighbor::stencil_none;
+	  } else if (triclinic) {
 	    half_build = &Neighbor::half_bin_newton_multi_tri;
-	  else half_build = &Neighbor::half_bin_newton_multi;
+	    if (force->dimension == 3)
+	      half_stencil = &Neighbor::stencil_half_3d_newton_multi_tri;
+	    else
+	      half_stencil = &Neighbor::stencil_half_2d_newton_multi_tri;
+	  } else {
+	    half_build = &Neighbor::half_bin_newton_multi;
+	    if (force->dimension == 3)
+	      half_stencil = &Neighbor::stencil_half_3d_newton_multi;
+	    else
+	      half_stencil = &Neighbor::stencil_half_2d_newton_multi;
+	  }
 	}
       }
     }
+
   } else half_build = NULL;
 
   if (full) {
     if (style == NSQ) full_build = &Neighbor::full_nsq;
-    else if (style == BIN) full_build = &Neighbor::full_bin;
+    else if (style == BIN) {
-    else error->all("Neighbor multi not allowed with full neighbor lists");
+      full_build = &Neighbor::full_bin;
+      if (force->dimension == 3)
+	full_stencil = &Neighbor::stencil_full_3d;
+      else
+	full_stencil = &Neighbor::stencil_full_2d;
+    } else {
+      full_build = &Neighbor::full_bin_multi;
+      if (force->dimension == 3)
+	full_stencil = &Neighbor::stencil_full_3d_multi;
+      else
+	full_stencil = &Neighbor::stencil_full_2d_multi;
+    }
   } else full_build = NULL;
 
   // ------------------------------------------------------------------
@@ -663,7 +692,6 @@ void Neighbor::init()
   // bond_quartic sets bonds to 0
   // delete_bonds sets all interactions negative
 
-
   int bond_off = 0;
   int angle_off = 0;
   for (i = 0; i < modify->nfix; i++)
@@ -705,7 +733,7 @@ void Neighbor::init()
     }
   }
 
-  // set ptrs to correct intra-molecular build functions
+  // set ptrs to intra-molecular build functions
 
   if (bond_off) bond_build = &Neighbor::bond_partial;
   else bond_build = &Neighbor::bond_all;
@@ -1040,107 +1068,26 @@ void Neighbor::setup_bins()
     binhead = (int *) memory->smalloc(maxhead*sizeof(int),"neigh:binhead");
   }
 
-  // create stencil of bins whose closest corner to central bin
+  // create stencil of bins to search over in neighbor list construction
-  //   is within cutneighmax
+  // sx,sy,sz = max range of stencil extent
-  // stencil will be empty if cutneighmax = 0.0
+  // stencil is empty if cutneighmax = 0.0
-  // next(xyz) = how far the stencil could possibly extend
-  // for partial Newton (newton = 0)
-  //   stencil is all surrounding bins
-  //   stencil includes self
-  // for triclinic
-  //   stencil is all bins in z-plane of self and above, but not below
-  //   in 2d is all bins in y-plane of self and above, but not below
-  //   stencil includes self
-  // for full Newton (newton = 1)
-  //   stencil is bins to the "upper right" of central bin
-  //   stencil does not include self
-  // for full neighbor list (full = 1)
-  //   stencil is all surrounding bins including self, regardless of Newton
-  //   stored in stencil_full
-  // 3d creates xyz stencil, 2d is only xy
 
-  int nextx = static_cast<int> (cutneighmax*bininvx);
+  int sx = static_cast<int> (cutneighmax*bininvx);
-  if (nextx*binsizex < cutneighmax) nextx++;
+  if (sx*binsizex < cutneighmax) sx++;
-  int nexty = static_cast<int> (cutneighmax*bininvy);
+  int sy = static_cast<int> (cutneighmax*bininvy);
-  if (nexty*binsizey < cutneighmax) nexty++;
+  if (sy*binsizey < cutneighmax) sy++;
-  int nextz = static_cast<int> (cutneighmax*bininvz);
+  int sz = static_cast<int> (cutneighmax*bininvz);
-  if (nextz*binsizez < cutneighmax) nextz++;
+  if (sz*binsizez < cutneighmax) sz++;
 
-  int nmax = (2*nextz+1) * (2*nexty+1) * (2*nextx+1);
+  // allocate stencil memory and create stencil(s)
-  if (nmax > maxstencil) {
+  // check half/full instead of half_every/full_every so stencils will be
-    maxstencil = nmax;
+  //   allocated correctly whenever build_half() and build_full() are called
-    memory->sfree(stencil);
-    stencil = (int *) memory->smalloc(maxstencil*sizeof(int),"neigh:stencil");
-  }
 
-  int i,j,k;
+  stencil_allocate(sx,sy,sz);
-  nstencil = 0;
+  if (half) (this->*half_stencil)(sx,sy,sz);
-  double cutneighmaxsq = cutneighmax*cutneighmax;
+  if (full) (this->*full_stencil)(sx,sy,sz);
-
-  if (force->dimension == 3) {
-    if (force->newton_pair == 0) {
-      for (k = -nextz; k <= nextz; k++)
-	for (j = -nexty; j <= nexty; j++)
-	  for (i = -nextx; i <= nextx; i++)
-	    if (bin_distance(i,j,k) < cutneighmaxsq)
-	      stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
-    } else if (triclinic) {
-      for (k = 0; k <= nextz; k++)
-	for (j = -nexty; j <= nexty; j++)
-	  for (i = -nextx; i <= nextx; i++)
-	    if (bin_distance(i,j,k) < cutneighmaxsq)
-	      stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
-    } else {
-      for (k = 0; k <= nextz; k++)
-	for (j = -nexty; j <= nexty; j++)
-	  for (i = -nextx; i <= nextx; i++)
-	    if (k > 0 || j > 0 || (j == 0 && i > 0))
-	      if (bin_distance(i,j,k) < cutneighmaxsq)
-		stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
-    }
-  } else {
-    if (force->newton_pair == 0) {
-      for (j = -nexty; j <= nexty; j++)
-	for (i = -nextx; i <= nextx; i++)
-	  if (bin_distance(i,j,0) < cutneighmaxsq)
-	    stencil[nstencil++] = j*mbinx + i;
-    } else if (triclinic) {
-	for (j = 0; j <= nexty; j++)
-	  for (i = -nextx; i <= nextx; i++)
-	    if (bin_distance(i,j,0) < cutneighmaxsq)
-	      stencil[nstencil++] = j*mbinx + i;
-    } else {
-	for (j = 0; j <= nexty; j++)
-	  for (i = -nextx; i <= nextx; i++)
-	    if (j > 0 || (j == 0 && i > 0))
-	      if (bin_distance(i,j,0) < cutneighmaxsq)
-		stencil[nstencil++] = j*mbinx + i;
-    }
-  }
-
-  if (full) {
-    if (nmax > maxstencil_full) {
-      maxstencil_full = nmax;
-      memory->sfree(stencil_full);
-      stencil_full = (int *) memory->smalloc(maxstencil_full*sizeof(int),
-					     "neigh:stencil_full");
-    }
-    nstencil_full = 0;
-    if (force->dimension == 3) {
-      for (k = -nextz; k <= nextz; k++)
-	for (j = -nexty; j <= nexty; j++)
-	  for (i = -nextx; i <= nextx; i++)
-	    if (bin_distance(i,j,k) < cutneighmaxsq)
-	      stencil_full[nstencil_full++] = k*mbiny*mbinx + j*mbinx + i;
-    } else {
-      for (j = -nexty; j <= nexty; j++)
-	for (i = -nextx; i <= nextx; i++)
-	  if (bin_distance(i,j,0) < cutneighmaxsq)
-	    stencil_full[nstencil_full++] = j*mbinx + i;
-    }
-  }
 }
-      
+
 /* ----------------------------------------------------------------------
    compute closest distance between central bin (0,0,0) and bin (i,j,k)
 ------------------------------------------------------------------------- */
@@ -1149,26 +1096,17 @@ double Neighbor::bin_distance(int i, int j, int k)
 {
   double delx,dely,delz;
 
-  if (i > 0)
+  if (i > 0) delx = (i-1)*binsizex;
-    delx = (i-1)*binsizex;
+  else if (i == 0) delx = 0.0;
-  else if (i == 0)
+  else delx = (i+1)*binsizex;
-    delx = 0.0;
-  else
-    delx = (i+1)*binsizex;
 
-  if (j > 0)
+  if (j > 0) dely = (j-1)*binsizey;
-    dely = (j-1)*binsizey;
+  else if (j == 0) dely = 0.0;
-  else if (j == 0)
+  else dely = (j+1)*binsizey;
-    dely = 0.0;
-  else
-    dely = (j+1)*binsizey;
 
-  if (k > 0)
+  if (k > 0) delz = (k-1)*binsizez;
-    delz = (k-1)*binsizez;
+  else if (k == 0) delz = 0.0;
-  else if (k == 0)
+  else delz = (k+1)*binsizez;
-    delz = 0.0;
-  else
-    delz = (k+1)*binsizez;
  
   return (delx*delx + dely*dely + delz*delz);
 }
@@ -1317,8 +1255,8 @@ int Neighbor::memory_usage()
   if (fix_history) {
     bytes += maxlocal * sizeof(int *);
     bytes += maxlocal * sizeof(double *);
-    bytes += maxpage*pgsize * sizeof(int);
+    bytes += maxpage_history*pgsize * sizeof(int);
-    bytes += maxpage*pgsize*3 * sizeof(double);
+    bytes += maxpage_history*pgsize*3 * sizeof(double);
   }
 
   if (respa) {
@@ -1341,7 +1279,7 @@ int Neighbor::memory_usage()
 }
 
 /* ----------------------------------------------------------------------
-   add pages to half or full neighbor list, starting at npage
+   add pages to half/full/granular/rRESPA neighbor lists, starting at npage
 ------------------------------------------------------------------------- */
 
 void Neighbor::add_pages(int npage)
@@ -1363,18 +1301,16 @@ void Neighbor::add_pages_full(int npage)
       (int *) memory->smalloc(pgsize*sizeof(int),"neigh:pages_full[i]");
 }
 
-/* ----------------------------------------------------------------------
-   add pages to granular neighbor list, starting at npage
-------------------------------------------------------------------------- */
-
 void Neighbor::add_pages_history(int npage)
 {
+  maxpage_history += PGDELTA;
   pages_touch = (int **)
-    memory->srealloc(pages_touch,maxpage*sizeof(int *),"neigh:pages_touch");
+    memory->srealloc(pages_touch,maxpage_history*sizeof(int *),
+		     "neigh:pages_touch");
   pages_shear = (double **)
-    memory->srealloc(pages_shear,maxpage*sizeof(double *),
+    memory->srealloc(pages_shear,maxpage_history*sizeof(double *),
 		     "neigh:pages_shear");
-  for (int i = npage; i < maxpage; i++) {
+  for (int i = npage; i < maxpage_history; i++) {
     pages_touch[i] = (int *)
       memory->smalloc(pgsize*sizeof(int),"neigh:pages_touch[i]");
     pages_shear[i] = (double *)
@@ -1382,10 +1318,6 @@ void Neighbor::add_pages_history(int npage)
   }
 }
 
-/* ----------------------------------------------------------------------
-   add pages to rRESPA inner neighbor list, starting at npage_inner
-------------------------------------------------------------------------- */
-
 void Neighbor::add_pages_inner(int npage_inner)
 {
   maxpage_inner += PGDELTA;
@@ -1397,10 +1329,6 @@ void Neighbor::add_pages_inner(int npage_inner)
       (int *) memory->smalloc(pgsize*sizeof(int),"neigh:pages_inner[i]");
 }
 
-/* ----------------------------------------------------------------------
-   add pages to rRESPA middle neighbor list, starting at npage_middle
-------------------------------------------------------------------------- */
-
 void Neighbor::add_pages_middle(int npage_middle)
 {
   maxpage_middle += PGDELTA;
@@ -1412,6 +1340,53 @@ void Neighbor::add_pages_middle(int npage_middle)
       (int *) memory->smalloc(pgsize*sizeof(int),"neigh:pages_middle[i]");
 }
 
+/* ----------------------------------------------------------------------
+   clear half/full/granular/rRESPA neighbor lists
+------------------------------------------------------------------------- */
+
+void Neighbor::clear_pages()
+{
+  for (int i = 0; i < maxpage; i++) memory->sfree(pages[i]);
+  memory->sfree(pages);
+  pages = NULL;
+  maxpage = 0;
+}
+
+void Neighbor::clear_pages_full()
+{
+  for (int i = 0; i < maxpage_full; i++) memory->sfree(pages_full[i]);
+  memory->sfree(pages_full);
+  pages_full = NULL;
+  maxpage_full = 0;
+}
+
+void Neighbor::clear_pages_history()
+{
+  for (int i = 0; i < maxpage_history; i++) memory->sfree(pages_touch[i]);
+  for (int i = 0; i < maxpage_history; i++) memory->sfree(pages_shear[i]);
+  memory->sfree(pages_touch);
+  memory->sfree(pages_shear);
+  pages_touch = NULL;
+  pages_shear = NULL;
+  maxpage_history = 0;
+}
+
+void Neighbor::clear_pages_inner()
+{
+  for (int i = 0; i < maxpage_inner; i++) memory->sfree(pages_inner[i]);
+  memory->sfree(pages_inner);
+  pages_inner = NULL;
+  maxpage_inner = 0;
+}
+
+void Neighbor::clear_pages_middle()
+{
+  for (int i = 0; i < maxpage_middle; i++) memory->sfree(pages_middle[i]);
+  memory->sfree(pages_middle);
+  pages_middle = NULL;
+  maxpage_middle = 0;
+}
+
 /* ----------------------------------------------------------------------
    determine if atom j is in special list of atom i
    if it is not, return 0

src/neighbor.h

@@ -30,6 +30,7 @@ class Neighbor : protected Pointers {
 
   double skin;                     // skin distance
   double cutghost;                 // distance for acquiring ghost atoms
+  double *cuttype;                 // for each type, max neigh cut w/ others
 
   int ncalls;                      // # of times build has been called
   int ndanger;                     // # of dangerous builds
@@ -94,11 +95,13 @@ class Neighbor : protected Pointers {
   int fix_check;                   // # of fixes that induce reneigh
   int *fixchecklist;               // which fixes to check
 
-  double **cutneigh;               // neighbor cutoff for type pairs
+  double **cutneighsq;             // neighbor cutneigh sq for each type pair
-  double **cutneighsq;             // cutneigh squared
-  double cutneighmax;              // max neighbor cutoff for all type pairs
   double cutneighmin;              // min neighbor cutoff (for cutforce > 0)
-  double triggersq;                // trigger build when atom moves this dist
+  double cutneighmax;              // max neighbor cutoff for all type pairs
+  double cutneighmaxsq;            // cutneighmax squared
+  double *cuttypesq;               // cuttype squared
+
+  double triggersq;                // trigger = build when atom moves this dist
 
   double **xhold;                  // atom coords at last neighbor build
   int maxhold;                     // size of xhold array
@@ -122,15 +125,21 @@ class Neighbor : protected Pointers {
 
   int nstencil;                    // # of bins in half neighbor stencil
   int *stencil;                    // list of bin offsets
-  int maxstencil;                  // size of stencil
+  int maxstencil;                  // max size of stencil
 
   int nstencil_full;               // # of bins in full neighbor stencil
   int *stencil_full;               // list of bin offsets
-  int maxstencil_full;             // size of stencil
+  int maxstencil_full;             // max size of stencil
 
-  int *mstencils;                  // # bins in each type-based multi stencil
+  int *nstencil_multi;             // # bins in each type-based multi stencil
-  int **mstencil;                  // list of bin offsets in each stencil
+  int **stencil_multi;             // list of bin offsets in each stencil
-  double **mdist;                  // list of distances to bins in each stencil
+  double **distsq_multi;           // sq distances to bins in each stencil
+  int maxstencil_multi;            // max sizes of stencils
+
+  int *nstencil_full_multi;        // # bins in full type-based multi stencil
+  int **stencil_full_multi;        // list of bin offsets in each stencil
+  double **distsq_full_multi;      // sq distances to bins in each stencil
+  int maxstencil_full_multi;       // max sizes of stencils
 
   int **pages;                     // half neighbor list pages
   int maxpage;                     // # of half pages currently allocated
@@ -143,6 +152,7 @@ class Neighbor : protected Pointers {
                                          // else is ptr to fix shear history
   int **pages_touch;               // pages of touch flags
   double **pages_shear;            // pages of shear values
+  int maxpage_history;             // # of history pages currently allocated
 
                                    // rRESPA neighbor lists
   int respa;                       // 0 = single neighbor list
@@ -183,17 +193,20 @@ class Neighbor : protected Pointers {
 
   void half_nsq_no_newton();          // fns for half neighbor lists
   void half_nsq_newton();
+
   void half_bin_no_newton();
   void half_bin_no_newton_multi();
   void half_bin_newton();
   void half_bin_newton_multi();
   void half_bin_newton_tri();
   void half_bin_newton_multi_tri();
+
   void half_full_no_newton();
   void half_full_newton();
 
   void full_nsq();                    // fns for full neighbor lists
   void full_bin();
+  void full_bin_multi();
 
   void granular_nsq_no_newton();      // fns for granular neighbor lists
   void granular_nsq_newton();
@@ -223,17 +236,51 @@ class Neighbor : protected Pointers {
   void improper_all();                // improper list with all impropers
   void improper_partial();            // exclude certain impropers
 
-  void add_pages(int);                // add pages to half neigh list
+  void add_pages(int);                // add neigh list pages
-  void add_pages_full(int);           // add pages to full neigh list
+  void add_pages_full(int);
-  void add_pages_history(int);        // add pages to granular neigh list
+  void add_pages_history(int);
-  void add_pages_inner(int);          // add pages to respa inner list
+  void add_pages_inner(int);
-  void add_pages_middle(int);         // add pages to respa middle list
+  void add_pages_middle(int);
+
+  void clear_pages();                 // clear all neigh list pages
+  void clear_pages_full();
+  void clear_pages_history();
+  void clear_pages_inner();
+  void clear_pages_middle();
 
   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 find_special(int, int);           // look for special neighbor
   int exclusion(int, int, int *, int *, int *);  // test for pair exclusion
+
+  typedef void (Neighbor::*SFnPtr)(int, int, int);
+  SFnPtr half_stencil;                // ptr to half stencil functions
+  SFnPtr full_stencil;                // ptr to full stencil functions
+
+  void stencil_allocate(int, int, int);
+
+  void stencil_none(int, int, int);   // fns for stencil creation
+
+  void stencil_half_3d_no_newton(int, int, int);
+  void stencil_half_3d_no_newton_multi(int, int, int);
+  void stencil_half_3d_newton(int, int, int);
+  void stencil_half_3d_newton_multi(int, int, int);
+  void stencil_half_3d_newton_tri(int, int, int);
+  void stencil_half_3d_newton_multi_tri(int, int, int);
+
+  void stencil_half_2d_no_newton(int, int, int);
+  void stencil_half_2d_no_newton_multi(int, int, int);
+  void stencil_half_2d_newton(int, int, int);
+  void stencil_half_2d_newton_multi(int, int, int);
+  void stencil_half_2d_newton_tri(int, int, int);
+  void stencil_half_2d_newton_multi_tri(int, int, int);
+
+  void stencil_full_3d(int, int, int);
+  void stencil_full_3d_multi(int, int, int);
+
+  void stencil_full_2d(int, int, int);
+  void stencil_full_2d_multi(int, int, int);
 };
 
 }

src/region_prism.cpp

@@ -12,7 +12,7 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing author: Pieter in't Veld (SNL)
+   Contributing author: Pieter in 't Veld (SNL)
 ------------------------------------------------------------------------- */
 
 #include "stdlib.h"