diff --git a/src/DPD/atom_vec_dpd.cpp b/src/DPD/atom_vec_dpd.cpp
index c29dfe24c6fee9b7bfbcf2d1cbed3e1eed6dc5c6..f292e7f01d99cc2670257c9889ad9826bd8168e6 100644
--- a/src/DPD/atom_vec_dpd.cpp
+++ b/src/DPD/atom_vec_dpd.cpp
@@ -14,6 +14,7 @@
 #include "stdlib.h"
 #include "atom_vec_dpd.h"
 #include "atom.h"
+#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -58,9 +59,10 @@ void AtomVecDPD::zero_ghost(int n, int first)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecDPD::pack_comm(int n, int *list, double *buf,
-			  int pbc_flag, double *pbc_dist)
+			  int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -74,11 +76,20 @@ int AtomVecDPD::pack_comm(int n, int *list, double *buf,
       buf[m++] = v[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = v[j][0];
       buf[m++] = v[j][1];
       buf[m++] = v[j][2];
@@ -128,9 +139,10 @@ int AtomVecDPD::unpack_comm_one(int i, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecDPD::pack_border(int n, int *list, double *buf,
-			    int pbc_flag, double *pbc_dist)
+			    int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -147,11 +159,20 @@ int AtomVecDPD::pack_border(int n, int *list, double *buf,
       buf[m++] = v[j][2];
    }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
diff --git a/src/DPD/atom_vec_dpd.h b/src/DPD/atom_vec_dpd.h
index 29c3a449bda5f6e1777c9e22f8e227674113aa56..b86a6e13e77efdb73eafa62378fe1fc75d019d0e 100644
--- a/src/DPD/atom_vec_dpd.h
+++ b/src/DPD/atom_vec_dpd.h
@@ -22,11 +22,11 @@ class AtomVecDPD : public AtomVecAtomic {
  public:
   AtomVecDPD(class LAMMPS *, int, char **);
   void zero_ghost(int, int);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   int pack_comm_one(int, double *);
   void unpack_comm(int, int, double *);
   int unpack_comm_one(int, double *);
-  int pack_border(int, int *, double *, int, double *);
+  int pack_border(int, int *, double *, int, int *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
diff --git a/src/GRANULAR/atom_vec_granular.cpp b/src/GRANULAR/atom_vec_granular.cpp
index 2e5a5a9c8182b7b45120f7fc51a1bd51de177742..ad352f3a03e992fdee72b018649a94dae9410a5e 100644
--- a/src/GRANULAR/atom_vec_granular.cpp
+++ b/src/GRANULAR/atom_vec_granular.cpp
@@ -15,8 +15,9 @@
 #include "stdlib.h"
 #include "atom_vec_granular.h"
 #include "atom.h"
-#include "force.h"
+#include "domain.h"
 #include "modify.h"
+#include "force.h"
 #include "fix.h"
 #include "memory.h"
 #include "error.h"
@@ -177,9 +178,10 @@ void AtomVecGranular::copy(int i, int j)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecGranular::pack_comm(int n, int *list, double *buf,
-			       int pbc_flag, double *pbc_dist)
+			       int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -196,11 +198,20 @@ int AtomVecGranular::pack_comm(int n, int *list, double *buf,
       buf[m++] = phiv[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = v[j][0];
       buf[m++] = v[j][1];
       buf[m++] = v[j][2];
@@ -319,9 +330,10 @@ int AtomVecGranular::unpack_reverse_one(int i, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecGranular::pack_border(int n, int *list, double *buf,
-				 int pbc_flag, double *pbc_dist)
+				 int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -343,11 +355,20 @@ int AtomVecGranular::pack_border(int n, int *list, double *buf,
       buf[m++] = phiv[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
diff --git a/src/GRANULAR/atom_vec_granular.h b/src/GRANULAR/atom_vec_granular.h
index 15d2887a3a67e21fb86e73b592a5db47ce30e23c..25c26a67558bff22ced3cdd433f46b6c2c1051b3 100644
--- a/src/GRANULAR/atom_vec_granular.h
+++ b/src/GRANULAR/atom_vec_granular.h
@@ -27,7 +27,7 @@ class AtomVecGranular : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   int pack_comm_one(int, double *);
   void unpack_comm(int, int, double *);
   int unpack_comm_one(int, double *);
@@ -35,7 +35,7 @@ class AtomVecGranular : public AtomVec {
   int pack_reverse_one(int, double *);
   void unpack_reverse(int, int *, double *);
   int unpack_reverse_one(int, double *);
-  int pack_border(int, int *, double *, int, double *);
+  int pack_border(int, int *, double *, int, int *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp
index 7245fb1df14564462d229440b5671746088eb794..317d55be74c402e12e5bf9cf66fb0b03bd0bb423 100644
--- a/src/MANYBODY/pair_eam.cpp
+++ b/src/MANYBODY/pair_eam.cpp
@@ -817,8 +817,7 @@ void PairEAM::single_embed(int i, int itype, double &phi)
 
 /* ---------------------------------------------------------------------- */
 
-int PairEAM::pack_comm(int n, int *list, double *buf,
-		       int pbc_flag, double *pack_dist)
+int PairEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
 {
   int i,j,m;
 
diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h
index 865b751b3e01f2bf11185eafa17378a03e0d739d..29cebe433a6201d5640e9e3242859cdbb17c5b9c 100644
--- a/src/MANYBODY/pair_eam.h
+++ b/src/MANYBODY/pair_eam.h
@@ -31,7 +31,7 @@ class PairEAM : public Pair {
 
   void single_embed(int, int, double &);
 
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   int pack_reverse_comm(int, int, double *);
   void unpack_reverse_comm(int, int *, double *);
diff --git a/src/MOLECULE/atom_vec_angle.cpp b/src/MOLECULE/atom_vec_angle.cpp
index 2d5b2a137773da22d9201131353c985ff5db1b52..ad2243860deb9ab22b2d0d3f98ccc9b8db74831f 100644
--- a/src/MOLECULE/atom_vec_angle.cpp
+++ b/src/MOLECULE/atom_vec_angle.cpp
@@ -14,6 +14,7 @@
 #include "stdlib.h"
 #include "atom_vec_angle.h"
 #include "atom.h"
+#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -208,9 +209,10 @@ void AtomVecAngle::copy(int i, int j)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecAngle::pack_comm(int n, int *list, double *buf,
-			    int pbc_flag, double *pbc_dist)
+			    int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -221,11 +223,20 @@ int AtomVecAngle::pack_comm(int n, int *list, double *buf,
       buf[m++] = x[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
     }
   }
   return m;
@@ -280,9 +291,10 @@ void AtomVecAngle::unpack_reverse(int n, int *list, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecAngle::pack_border(int n, int *list, double *buf,
-			      int pbc_flag, double *pbc_dist)
+			      int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -297,11 +309,20 @@ int AtomVecAngle::pack_border(int n, int *list, double *buf,
       buf[m++] = molecule[j];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
diff --git a/src/MOLECULE/atom_vec_angle.h b/src/MOLECULE/atom_vec_angle.h
index 0ca4f087f346f2ae66877437372b04c97b6d4136..a1a5aea9257946d728be4f9baf29ab1c8965614b 100644
--- a/src/MOLECULE/atom_vec_angle.h
+++ b/src/MOLECULE/atom_vec_angle.h
@@ -26,11 +26,11 @@ class AtomVecAngle : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int, double *);
+  int pack_border(int, int *, double *, int, int *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
diff --git a/src/MOLECULE/atom_vec_bond.cpp b/src/MOLECULE/atom_vec_bond.cpp
index 22fbd9f1e289756db601d4299a3d170d8b020b2f..d276a7cea5c78ab5cc82c2a54b821fa688bf3bd7 100644
--- a/src/MOLECULE/atom_vec_bond.cpp
+++ b/src/MOLECULE/atom_vec_bond.cpp
@@ -14,6 +14,7 @@
 #include "stdlib.h"
 #include "atom_vec_bond.h"
 #include "atom.h"
+#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -177,9 +178,10 @@ void AtomVecBond::copy(int i, int j)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecBond::pack_comm(int n, int *list, double *buf,
-			   int pbc_flag, double *pbc_dist)
+			   int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -190,11 +192,20 @@ int AtomVecBond::pack_comm(int n, int *list, double *buf,
       buf[m++] = x[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
     }
   }
   return m;
@@ -249,9 +260,10 @@ void AtomVecBond::unpack_reverse(int n, int *list, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecBond::pack_border(int n, int *list, double *buf,
-			     int pbc_flag, double *pbc_dist)
+			     int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -266,11 +278,20 @@ int AtomVecBond::pack_border(int n, int *list, double *buf,
       buf[m++] = molecule[j];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
diff --git a/src/MOLECULE/atom_vec_bond.h b/src/MOLECULE/atom_vec_bond.h
index 9da8f73fd726fa02297dcf8483dc019c858e5bec..a5c5c77cb2bfac145f551895ea8f6bc2e6cc1d20 100644
--- a/src/MOLECULE/atom_vec_bond.h
+++ b/src/MOLECULE/atom_vec_bond.h
@@ -26,11 +26,11 @@ class AtomVecBond : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int, double *);
+  int pack_border(int, int *, double *, int, int *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
diff --git a/src/MOLECULE/atom_vec_full.cpp b/src/MOLECULE/atom_vec_full.cpp
index 2fa18dbd5adaaec21810427a6fe1ed44de7ceed1..4ae64fb7a5f2aeacd0c8911849e91f74a19143fe 100644
--- a/src/MOLECULE/atom_vec_full.cpp
+++ b/src/MOLECULE/atom_vec_full.cpp
@@ -14,6 +14,7 @@
 #include "stdlib.h"
 #include "atom_vec_full.h"
 #include "atom.h"
+#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -286,9 +287,10 @@ void AtomVecFull::copy(int i, int j)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecFull::pack_comm(int n, int *list, double *buf,
-			   int pbc_flag, double *pbc_dist)
+			   int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -299,11 +301,20 @@ int AtomVecFull::pack_comm(int n, int *list, double *buf,
       buf[m++] = x[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
     }
   }
   return m;
@@ -358,9 +369,10 @@ void AtomVecFull::unpack_reverse(int n, int *list, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecFull::pack_border(int n, int *list, double *buf,
-			     int pbc_flag, double *pbc_dist)
+			     int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -376,11 +388,20 @@ int AtomVecFull::pack_border(int n, int *list, double *buf,
       buf[m++] = molecule[j];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
diff --git a/src/MOLECULE/atom_vec_full.h b/src/MOLECULE/atom_vec_full.h
index b4ed2707638cbab7b3bceecf5d8520e9d50bebd8..22a35091d2f0c6f61d08f90994bbd52892f7e13f 100644
--- a/src/MOLECULE/atom_vec_full.h
+++ b/src/MOLECULE/atom_vec_full.h
@@ -26,11 +26,11 @@ class AtomVecFull : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int, double *);
+  int pack_border(int, int *, double *, int, int *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
diff --git a/src/MOLECULE/atom_vec_molecular.cpp b/src/MOLECULE/atom_vec_molecular.cpp
index 5c05b9fd638896b1bb784fb8e9632fa705712ae9..8cb26b1033c54ee34089b74a6e19c10019e590d0 100644
--- a/src/MOLECULE/atom_vec_molecular.cpp
+++ b/src/MOLECULE/atom_vec_molecular.cpp
@@ -14,6 +14,7 @@
 #include "stdlib.h"
 #include "atom_vec_molecular.h"
 #include "atom.h"
+#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -280,9 +281,10 @@ void AtomVecMolecular::copy(int i, int j)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMolecular::pack_comm(int n, int *list, double *buf, 
-				int pbc_flag, double *pbc_dist)
+				int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -293,11 +295,20 @@ int AtomVecMolecular::pack_comm(int n, int *list, double *buf,
       buf[m++] = x[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
     }
   }
   return m;
@@ -352,9 +363,10 @@ void AtomVecMolecular::unpack_reverse(int n, int *list, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMolecular::pack_border(int n, int *list, double *buf,
-				  int pbc_flag, double *pbc_dist)
+				  int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -369,11 +381,20 @@ int AtomVecMolecular::pack_border(int n, int *list, double *buf,
       buf[m++] = molecule[j];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
diff --git a/src/MOLECULE/atom_vec_molecular.h b/src/MOLECULE/atom_vec_molecular.h
index 2b2e3a49fe84b17cd32c422ea456cec810a56e4f..c436afa45c0180c972f377242c74d54be70c4893 100644
--- a/src/MOLECULE/atom_vec_molecular.h
+++ b/src/MOLECULE/atom_vec_molecular.h
@@ -26,11 +26,11 @@ class AtomVecMolecular : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int, double *);
+  int pack_border(int, int *, double *, int, int *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
diff --git a/src/atom_vec.h b/src/atom_vec.h
index 8278fdea5af34e5133c3171b3bf9f079753f6c0b..59972d8a07586248caa62889ac6face0ac41d70d 100644
--- a/src/atom_vec.h
+++ b/src/atom_vec.h
@@ -44,7 +44,7 @@ class AtomVec : protected Pointers {
   virtual void zero_ghost(int,int) {}
   virtual void copy(int, int) = 0;
 
-  virtual int pack_comm(int, int *, double *, int, double *) = 0;
+  virtual int pack_comm(int, int *, double *, int, int *) = 0;
   virtual int pack_comm_one(int, double *) {return 0;}
   virtual void unpack_comm(int, int, double *) = 0;
   virtual int unpack_comm_one(int, double *) {return 0;}
@@ -54,7 +54,7 @@ class AtomVec : protected Pointers {
   virtual void unpack_reverse(int, int *, double *) = 0;
   virtual int unpack_reverse_one(int, double *) {return 0;}
 
-  virtual int pack_border(int, int *, double *, int, double *) = 0;
+  virtual int pack_border(int, int *, double *, int, int *) = 0;
   virtual int pack_border_one(int, double *) {return 0;}
   virtual void unpack_border(int, int, double *) = 0;
   virtual int unpack_border_one(int, double *) {return 0;}
diff --git a/src/atom_vec_atomic.cpp b/src/atom_vec_atomic.cpp
index 9921b30c704b09485d2b55924541de0ef59d5d8b..562dea42d3994a5d88fe144ff19cce87c5e2916d 100644
--- a/src/atom_vec_atomic.cpp
+++ b/src/atom_vec_atomic.cpp
@@ -14,6 +14,7 @@
 #include "stdlib.h"
 #include "atom_vec_atomic.h"
 #include "atom.h"
+#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -102,9 +103,10 @@ void AtomVecAtomic::copy(int i, int j)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecAtomic::pack_comm(int n, int *list, double *buf,
-			     int pbc_flag, double *pbc_dist)
+			     int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -115,11 +117,20 @@ int AtomVecAtomic::pack_comm(int n, int *list, double *buf,
       buf[m++] = x[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
     }
   }
   return m;
@@ -174,9 +185,10 @@ void AtomVecAtomic::unpack_reverse(int n, int *list, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecAtomic::pack_border(int n, int *list, double *buf,
-			       int pbc_flag, double *pbc_dist)
+			       int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -190,11 +202,20 @@ int AtomVecAtomic::pack_border(int n, int *list, double *buf,
       buf[m++] = mask[j];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
diff --git a/src/atom_vec_atomic.h b/src/atom_vec_atomic.h
index 0cff9b0d0edb2545c20fbfe508d70c48a35cd58f..3158e22077f2fc6a8d4a4b88736887f29c789db8 100644
--- a/src/atom_vec_atomic.h
+++ b/src/atom_vec_atomic.h
@@ -25,11 +25,11 @@ class AtomVecAtomic : public AtomVec {
   void grow(int);
   void reset_ptrs();
   void copy(int, int);
-  virtual int pack_comm(int, int *, double *, int, double *);
+  virtual int pack_comm(int, int *, double *, int, int *);
   virtual void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  virtual int pack_border(int, int *, double *, int, double *);
+  virtual int pack_border(int, int *, double *, int, int *);
   virtual void unpack_border(int, int, double *);
   int pack_exchange(int, double *);
   int unpack_exchange(double *);
diff --git a/src/atom_vec_charge.cpp b/src/atom_vec_charge.cpp
index 20c36f28654199f3bb2eda3859b34baf122dd47e..631697eb7c2e56b932500eb4566ef9386b0b759e 100644
--- a/src/atom_vec_charge.cpp
+++ b/src/atom_vec_charge.cpp
@@ -14,6 +14,7 @@
 #include "stdlib.h"
 #include "atom_vec_charge.h"
 #include "atom.h"
+#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -134,9 +135,10 @@ void AtomVecCharge::copy(int i, int j)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecCharge::pack_comm(int n, int *list, double *buf,
-			     int pbc_flag, double *pbc_dist)
+			     int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -147,11 +149,20 @@ int AtomVecCharge::pack_comm(int n, int *list, double *buf,
       buf[m++] = x[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
     }
   }
   return m;
@@ -206,9 +217,10 @@ void AtomVecCharge::unpack_reverse(int n, int *list, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecCharge::pack_border(int n, int *list, double *buf,
-			       int pbc_flag, double *pbc_dist)
+			       int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -223,11 +235,20 @@ int AtomVecCharge::pack_border(int n, int *list, double *buf,
       buf[m++] = q[j];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
diff --git a/src/atom_vec_charge.h b/src/atom_vec_charge.h
index 955d74008ea58c118699b0136426d5cd5465c594..5b6efb9de7fe5ac05ba0415fbd7cf7f4a8f67107 100644
--- a/src/atom_vec_charge.h
+++ b/src/atom_vec_charge.h
@@ -26,11 +26,11 @@ class AtomVecCharge : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int, double *);
+  int pack_border(int, int *, double *, int, int *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
diff --git a/src/atom_vec_hybrid.cpp b/src/atom_vec_hybrid.cpp
index f5d01daf541cb49dd9ddb12c94035a4671e09649..93497484cd8f8782a7fd5379e06757d4bb924edc 100644
--- a/src/atom_vec_hybrid.cpp
+++ b/src/atom_vec_hybrid.cpp
@@ -15,6 +15,7 @@
 #include "string.h"
 #include "atom_vec_hybrid.h"
 #include "atom.h"
+#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -153,9 +154,10 @@ void AtomVecHybrid::copy(int i, int j)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecHybrid::pack_comm(int n, int *list, double *buf,
-			     int pbc_flag, double *pbc_dist)
+			     int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -167,11 +169,20 @@ int AtomVecHybrid::pack_comm(int n, int *list, double *buf,
       m += styles[hybrid[j]]->pack_comm_one(j,&buf[m]);
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       m += styles[hybrid[j]]->pack_comm_one(j,&buf[m]);
     }
   }
@@ -230,9 +241,10 @@ void AtomVecHybrid::unpack_reverse(int n, int *list, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecHybrid::pack_border(int n, int *list, double *buf,
-			       int pbc_flag, double *pbc_dist)
+			       int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -248,11 +260,20 @@ int AtomVecHybrid::pack_border(int n, int *list, double *buf,
       m += styles[hybrid[j]]->pack_border_one(j,&buf[m]);
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_dist[0];
-      buf[m++] = x[j][1] + pbc_dist[1];
-      buf[m++] = x[j][2] + pbc_dist[2];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
diff --git a/src/atom_vec_hybrid.h b/src/atom_vec_hybrid.h
index bb402d6e83d97769163c8e2f766181fbe75d4a18..50e688830dae9ab82a323c73c6cf9ed5ba7a0962 100644
--- a/src/atom_vec_hybrid.h
+++ b/src/atom_vec_hybrid.h
@@ -29,11 +29,11 @@ class AtomVecHybrid : public AtomVec {
   void grow(int);
   void reset_ptrs();
   void copy(int, int);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int, double *);
+  int pack_border(int, int *, double *, int, int *);
   void unpack_border(int, int, double *);
   int pack_exchange(int, double *);
   int unpack_exchange(double *);
diff --git a/src/comm.cpp b/src/comm.cpp
index 42f267ec08cae5b99d171975d9356322a9671848..bea87a64b3e3e66b1bb9c33558c083e033c164b7 100644
--- a/src/comm.cpp
+++ b/src/comm.cpp
@@ -171,11 +171,12 @@ void Comm::set_procs()
 
 void Comm::setup()
 {
-  // cutcomm = distance at which ghost atoms need to be acquired
-  // for orthogonal, cutcomm = box coords = cutneigh in all 3 dims
-  // for triclinic, cutneigh = distance between tilted planes in box coords
-  //   cutcomm = lamda coords = distance between those same planes
-  //                            can be different for each dim
+  // cutghost = 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
 
   double *prd,*prd_border,*sublo,*subhi;
 
@@ -183,7 +184,7 @@ void Comm::setup()
     prd = prd_border = domain->prd;
     sublo = domain->sublo;
     subhi = domain->subhi;
-    cutcomm[0] = cutcomm[1] = cutcomm[2] = neighbor->cutneigh;
+    cutghost[0] = cutghost[1] = cutghost[2] = neighbor->cutghost;
   } else {
     prd = domain->prd;
     prd_border = domain->prd_lamda;
@@ -192,19 +193,19 @@ void Comm::setup()
     double *h_inv = domain->h_inv;
     double length;
     length = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]);
-    cutcomm[0] = neighbor->cutneigh * length;
+    cutghost[0] = neighbor->cutghost * length;
     length = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]);
-    cutcomm[1] = neighbor->cutneigh * length;
+    cutghost[1] = neighbor->cutghost * length;
     length = h_inv[2];
-    cutcomm[2] = neighbor->cutneigh * length;
+    cutghost[2] = neighbor->cutghost * length;
   }
 
   // need = # of procs I need atoms from in each dim
   // for 2d, don't communicate in z
 
-  need[0] = static_cast<int> (cutcomm[0] * procgrid[0] / prd_border[0]) + 1;
-  need[1] = static_cast<int> (cutcomm[1] * procgrid[1] / prd_border[1]) + 1;
-  need[2] = static_cast<int> (cutcomm[2] * procgrid[2] / prd_border[2]) + 1;
+  need[0] = static_cast<int> (cutghost[0] * procgrid[0] / prd_border[0]) + 1;
+  need[1] = static_cast<int> (cutghost[1] * procgrid[1] / prd_border[1]) + 1;
+  need[2] = static_cast<int> (cutghost[2] * procgrid[2] / prd_border[2]) + 1;
   if (force->dimension == 2) need[2] = 0;
 
   // if non-periodic, do not communicate further than procgrid-1 away
@@ -230,8 +231,8 @@ void Comm::setup()
   // set slablo > slabhi for swaps across non-periodic boundaries
   //   this insures no atoms are swapped
   //   only for procs owning sub-box at non-periodic end of global box
-  // pbc_flag: 0 = not across a boundary, 1 = yes across a boundary
-  // pbc_dist/border: factors to add to atom coords across PBC for comm/borders
+  // pbc_flag: 0 = nothing across a boundary, 1 = somthing 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
@@ -240,36 +241,31 @@ void Comm::setup()
   for (int dim = 0; dim < 3; dim++) {
     for (int ineed = 0; ineed < 2*need[dim]; ineed++) {
       pbc_flag[iswap] = 0;
-      pbc_dist[iswap][0] = pbc_dist[iswap][1] = pbc_dist[iswap][2] = 0.0;
-      pbc_dist_border[iswap][0] = pbc_dist_border[iswap][1] = 
-	pbc_dist_border[iswap][2] = 0.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;
 	else slablo[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
-	slabhi[iswap] = sublo[dim] + cutcomm[dim];
+	slabhi[iswap] = sublo[dim] + cutghost[dim];
 	if (myloc[dim] == 0) {
 	  if (periodicity[dim] == 0)
 	    slabhi[iswap] = slablo[iswap] - 1.0;
 	  else {
 	    pbc_flag[iswap] = 1;
-	    pbc_dist[iswap][dim] = prd[dim];
-	    pbc_dist_border[iswap][dim] = prd_border[dim];
+	    pbc[iswap][dim] = 1;
 	    if (triclinic) {
-	      if (dim == 1) pbc_dist[iswap][0] += domain->xy;
-	      else if (dim == 2) {
-		pbc_dist[iswap][0] += domain->xz;
-		pbc_dist[iswap][1] += domain->yz;
-	      }
+	      if (dim == 1) pbc[iswap][5] = 1;
+	      else if (dim == 2) pbc[iswap][4] = pbc[iswap][3] = 1;
 	    }
 	  }
 	}
       } else {
 	sendproc[iswap] = procneigh[dim][1];
 	recvproc[iswap] = procneigh[dim][0];
-	slablo[iswap] = subhi[dim] - cutcomm[dim];
+	slablo[iswap] = subhi[dim] - cutghost[dim];
 	if (ineed < 2) slabhi[iswap] = BIG;
 	else slabhi[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
 	if (myloc[dim] == procgrid[dim]-1) {
@@ -277,14 +273,10 @@ void Comm::setup()
 	    slabhi[iswap] = slablo[iswap] - 1.0;
 	  else {
 	    pbc_flag[iswap] = 1;
-	    pbc_dist[iswap][dim] = -prd[dim];
-	    pbc_dist_border[iswap][dim] = -prd_border[dim];
+	    pbc[iswap][dim] = -1;
 	    if (triclinic) {
-	      if (dim == 1) pbc_dist[iswap][0] -= domain->xy;
-	      else if (dim == 2) {
-		pbc_dist[iswap][0] -= domain->xz;
-		pbc_dist[iswap][1] -= domain->yz;
-	      }
+	      if (dim == 1) pbc[iswap][5] = -1;
+	      else if (dim == 2) pbc[iswap][4] = pbc[iswap][3] = -1;
 	    }
 	  }
 	}
@@ -321,14 +313,14 @@ void Comm::communicate()
 	MPI_Irecv(buf,size_comm_recv[iswap],MPI_DOUBLE,
 		  recvproc[iswap],0,world,&request);
 	n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
-			    buf_send,pbc_flag[iswap],pbc_dist[iswap]);
+			    buf_send,pbc_flag[iswap],pbc[iswap]);
 	MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
 	MPI_Wait(&request,&status);
       } else {
 	MPI_Irecv(buf_recv,size_comm_recv[iswap],MPI_DOUBLE,
 		  recvproc[iswap],0,world,&request);
 	n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
-			    buf_send,pbc_flag[iswap],pbc_dist[iswap]);
+			    buf_send,pbc_flag[iswap],pbc[iswap]);
 	MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
 	MPI_Wait(&request,&status);
 	avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_recv);
@@ -338,11 +330,10 @@ void Comm::communicate()
       if (comm_x_only) {
 	if (sendnum[iswap])
 	  n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
-			      x[firstrecv[iswap]],
-			      pbc_flag[iswap],pbc_dist[iswap]);
+			      x[firstrecv[iswap]],pbc_flag[iswap],pbc[iswap]);
       } else {
 	n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
-			    buf_send,pbc_flag[iswap],pbc_dist[iswap]);
+			    buf_send,pbc_flag[iswap],pbc[iswap]);
 	avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_send);
       }
     }
@@ -567,8 +558,8 @@ void Comm::borders()
       if (nsend*size_border > maxsend)
 	grow_send(nsend*size_border,0);
       n = avec->pack_border(nsend,sendlist[iswap],buf_send,
-			    pbc_flag[iswap],pbc_dist_border[iswap]);
-
+			    pbc_flag[iswap],pbc[iswap]);
+      
       // swap atoms with other proc
       // put incoming ghosts at end of my atom arrays
       // if swapping with self, simply copy, no messages
@@ -635,7 +626,7 @@ void Comm::comm_pair(Pair *pair)
     // pack buffer
 
     n = pair->pack_comm(sendnum[iswap],sendlist[iswap],
-			buf_send,pbc_flag[iswap],pbc_dist[iswap]);
+			buf_send,pbc_flag[iswap],pbc[iswap]);
 
     // exchange with another proc
     // if self, set recv buffer to send buffer
@@ -704,7 +695,7 @@ void Comm::comm_fix(Fix *fix)
     // pack buffer
 
     n = fix->pack_comm(sendnum[iswap],sendlist[iswap],
-		       buf_send,pbc_flag[iswap],pbc_dist[iswap]);
+		       buf_send,pbc_flag[iswap],pbc[iswap]);
 
     // exchange with another proc
     // if self, set recv buffer to send buffer
@@ -773,7 +764,7 @@ void Comm::comm_compute(Compute *compute)
     // pack buffer
 
     n = compute->pack_comm(sendnum[iswap],sendlist[iswap],
-			   buf_send,pbc_flag[iswap],pbc_dist[iswap]);
+			   buf_send,pbc_flag[iswap],pbc[iswap]);
 
     // exchange with another proc
     // if self, set recv buffer to send buffer
@@ -980,9 +971,7 @@ void Comm::allocate_swap(int n)
   slabhi = (double *) memory->smalloc(n*sizeof(double),"comm:slabhi");
   firstrecv = (int *) memory->smalloc(n*sizeof(int),"comm:firstrecv");
   pbc_flag = (int *) memory->smalloc(n*sizeof(int),"comm:pbc_flag");
-  pbc_dist = (double **) memory->create_2d_double_array(n,3,"comm:pbc_dist");
-  pbc_dist_border = (double **)
-    memory->create_2d_double_array(n,3,"comm:pbc_dist_border");
+  pbc = (int **) memory->create_2d_int_array(n,6,"comm:pbc");
 }
 
 /* ----------------------------------------------------------------------
@@ -1002,8 +991,7 @@ void Comm::free_swap()
   memory->sfree(slabhi);
   memory->sfree(firstrecv);
   memory->sfree(pbc_flag);
-  memory->destroy_2d_double_array(pbc_dist);
-  memory->destroy_2d_double_array(pbc_dist_border);
+  memory->destroy_2d_int_array(pbc);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/comm.h b/src/comm.h
index 3725ebb9750ab7877dda1572c294ae1e9816cf35..92ca6e92d6c5cd52aa8baffbf84ef89b7a81e791 100644
--- a/src/comm.h
+++ b/src/comm.h
@@ -29,7 +29,7 @@ class Comm : protected Pointers {
   int need[3];                      // procs I need atoms from in each dim
   int maxforward_fix,maxreverse_fix;   // comm sizes called from Fix,Pair
   int maxforward_pair,maxreverse_pair;
-  double cutcomm[3];                // cutoffs used for acquiring ghost atoms
+  double cutghost[3];               // cutoffs used for acquiring ghost atoms
   
   Comm(class LAMMPS *);
   ~Comm();
@@ -60,9 +60,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
-  int *pbc_flag;                    // flags for sending atoms thru PBC
-  double **pbc_dist;                // distance to adjust atoms coords for PBC
-  double **pbc_dist_border;         // PBC distance for borders()
+  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
   int map_style;                    // non-0 if global->local mapping is done
 
diff --git a/src/compute.h b/src/compute.h
index 5e6c3cca86e4a1dacccefcfa914c58159e9b14f4..bae6bc1ab6e80cd0f07f1d17ec7231084096a848 100644
--- a/src/compute.h
+++ b/src/compute.h
@@ -55,7 +55,7 @@ class Compute : protected Pointers {
   virtual void compute_vector() {}
   virtual void compute_peratom() {}
 
-  virtual int pack_comm(int, int *, double *, int, double *) {return 0;}
+  virtual int pack_comm(int, int *, double *, int, int *) {return 0;}
   virtual void unpack_comm(int, int, double *) {}
   virtual int pack_reverse_comm(int, int, double *) {return 0;}
   virtual void unpack_reverse_comm(int, int *, double *) {}
diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp
index 2eecabd960d6dd5a9d1fb3b3e3c7d657a7ca8b29..660f31e2afff4c0a421452ac1949a280b9483c35 100644
--- a/src/delete_atoms.cpp
+++ b/src/delete_atoms.cpp
@@ -214,8 +214,8 @@ void DeleteAtoms::delete_overlap(int narg, char **arg, int *list)
 
   // error check on cutoff
 
-  if (cut > neighbor->cutneigh) 
-    error->all("Delete_atoms cutoff > neighbor cutoff");
+  if (cut > neighbor->cutghost) 
+    error->all("Delete_atoms cutoff > ghost cutoff");
 
   // create an atom map if one doesn't exist already
 
diff --git a/src/fix.h b/src/fix.h
index 97b18bb5780bfa23cb73d38ea642563ac9782c47..73049a200fbecaad04832b162c53e6fe05be0c8f 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -82,7 +82,7 @@ class Fix : protected Pointers {
   virtual int min_dof() {return 0;}
   virtual void min_xinitial(double *) {}
 
-  virtual int pack_comm(int, int *, double *, int, double *) {return 0;}
+  virtual int pack_comm(int, int *, double *, int, int *) {return 0;}
   virtual void unpack_comm(int, int, double *) {}
   virtual int pack_reverse_comm(int, int, double *) {return 0;}
   virtual void unpack_reverse_comm(int, int *, double *) {}
diff --git a/src/fix_orient_fcc.cpp b/src/fix_orient_fcc.cpp
index d6ac434d2c44a9c1dbd0695ea8ee7833722306a6..73b9f3b8958bec574ec5a20c374c5c1ca02015db 100644
--- a/src/fix_orient_fcc.cpp
+++ b/src/fix_orient_fcc.cpp
@@ -420,7 +420,7 @@ double FixOrientFCC::thermo(int n)
 /* ---------------------------------------------------------------------- */
 
 int FixOrientFCC::pack_comm(int n, int *list, double *buf,
-			    int pbc_flag, double *pair_dist)
+			    int pbc_flag, int *pbc)
 {
   int i,j,k,id,num;
   int *tag = atom->tag;
diff --git a/src/fix_orient_fcc.h b/src/fix_orient_fcc.h
index cdfc05673e699592af14299745ffbc622ad4d733..9b34e41535fb9d703ba71853f6c1b49cc9a69dc0 100644
--- a/src/fix_orient_fcc.h
+++ b/src/fix_orient_fcc.h
@@ -45,7 +45,7 @@ class FixOrientFCC : public Fix {
   void post_force(int);
   void post_force_respa(int, int, int);
   double thermo(int);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   int memory_usage();
 
diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp
index 263d742d84cecf59f4c50c639bc6be7bcfed5af4..0cca830821968a0781304b3484e2d6bc55f902b8 100644
--- a/src/fix_shake.cpp
+++ b/src/fix_shake.cpp
@@ -2252,10 +2252,10 @@ void FixShake::post_force_respa(int vflag_in, int ilevel, int iloop)
 
 /* ---------------------------------------------------------------------- */
 
-int FixShake::pack_comm(int n, int *list, double *buf,
-			int pbc_flag, double *pbc_dist)
+int FixShake::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
 {
   int i,j,m;
+  double dx,dy,dz;
 
   m = 0;
   if (pbc_flag == 0) {
@@ -2266,11 +2266,20 @@ int FixShake::pack_comm(int n, int *list, double *buf,
       buf[m++] = xshake[j][2];
     }
   } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = xshake[j][0] + pbc_dist[0];
-      buf[m++] = xshake[j][1] + pbc_dist[1];
-      buf[m++] = xshake[j][2] + pbc_dist[2];
+      buf[m++] = xshake[j][0] + dx;
+      buf[m++] = xshake[j][1] + dy;
+      buf[m++] = xshake[j][2] + dz;
     }
   }
   return 3;
diff --git a/src/fix_shake.h b/src/fix_shake.h
index 09e78dffc22d724cd7a6b821573e5f57bde68512..045661ccfa70e3648dd1ef63a6269bbf068eb46e 100644
--- a/src/fix_shake.h
+++ b/src/fix_shake.h
@@ -34,7 +34,7 @@ class FixShake : public Fix {
   void copy_arrays(int, int);
   int pack_exchange(int, double *);
   int unpack_exchange(int, double *);
-  int pack_comm(int, int *, double *, int, double *);
+  int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
 
   int dof(int);
diff --git a/src/force.h b/src/force.h
index 14008b0a0d70f51bac5b2ebaaa7e5b54e9ab75c9..e3876d53eba8d3f1f77152888881682c87a6ef79 100644
--- a/src/force.h
+++ b/src/force.h
@@ -58,21 +58,21 @@ class Force : protected Pointers {
   void init();
 
   void create_pair(char *);
-  Pair *new_pair(char *);
-  Pair *pair_match(char *);
+  class Pair *new_pair(char *);
+  class Pair *pair_match(char *);
 
   void create_bond(char *);
-  Bond *new_bond(char *);
-  Bond *bond_match(char *);
+  class Bond *new_bond(char *);
+  class Bond *bond_match(char *); 
 
   void create_angle(char *);
-  Angle *new_angle(char *);
+  class Angle *new_angle(char *);
 
   void create_dihedral(char *);
-  Dihedral *new_dihedral(char *);
+  class Dihedral *new_dihedral(char *);
 
   void create_improper(char *);
-  Improper *new_improper(char *);
+  class Improper *new_improper(char *);
 
   void create_kspace(int, char **);
 
diff --git a/src/input.cpp b/src/input.cpp
index 51425bbdc06743235b0ad097d7f33683d9d00bf0..5964e8b1710fea4a3bc37d451a97ac900bbcec7f 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -883,12 +883,7 @@ void Input::neigh_modify()
 
 void Input::neighbor_command()
 {
-  if (narg != 2) error->all("Illegal neighbor command");
-
-  neighbor->skin = atof(arg[0]);
-  if (strcmp(arg[1],"nsq") == 0) neighbor->style = 0;
-  else if (strcmp(arg[1],"bin") == 0) neighbor->style = 1;
-  else error->all("Illegal neighbor command");
+  neighbor->set(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/neigh_half.cpp b/src/neigh_half.cpp
index 37ca5ec665148db60ee11ef375481e85d01182b8..df7a19577e6da1b77f9a6c8a598e45fc2ef5751b 100644
--- a/src/neigh_half.cpp
+++ b/src/neigh_half.cpp
@@ -211,13 +211,14 @@ void Neighbor::half_bin_no_newton()
     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 store pair if i < j
     // stores own/own pairs only once
     // stores own/ghost pairs on both procs
 
+    ibin = coord2bin(x[i]);
+
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
 	if (j <= i) continue;
@@ -246,6 +247,94 @@ 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
+   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)
+------------------------------------------------------------------------- */
+
+void Neighbor::half_bin_no_newton_multi()
+{
+  int i,j,k,n,itype,jtype,ibin,which,ms;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *neighptr,*s;
+  double *cut,*dist;
+
+  // 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) add_pages(npage);
+    }
+
+    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
+    // 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];
+    dist = mdist[itype];
+    cut = cutneigh[itype];
+    ms = mstencils[itype];
+    for (k = 0; k < ms; 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 (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[i] = neighptr;
+    numneigh[i] = n;
+    npnt += n;
+    if (npnt >= pgsize)
+      error->one("Neighbor list overflow, boost neigh_modify one or page");
+  }
+}
+
 /* ----------------------------------------------------------------------
    binned neighbor list construction with full Newton's 3rd law
    each owned atom i checks its own bin and other bins in Newton stencil
@@ -348,6 +437,116 @@ void Neighbor::half_bin_newton()
   }
 }
 
+/* ----------------------------------------------------------------------
+   binned neighbor list construction with full Newton's 3rd law
+   each owned atom i checks its own bin and other bins in Newton stencil
+   multi-type stencil is itype dependent and is distance checked
+   every pair stored exactly once by some processor
+------------------------------------------------------------------------- */
+
+void Neighbor::half_bin_newton_multi()
+{
+  int i,j,k,n,itype,jtype,ibin,which,ms;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *neighptr,*s;
+  double *cut,*dist;
+
+  // 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) add_pages(npage);
+    }
+
+    neighptr = &pages[npage][npnt];
+    n = 0;
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+
+    // loop over rest of atoms in i's bin, ghosts are at end of linked list
+    // if j is owned atom, store it, since j is beyond i in linked list
+    // if j is ghost, only store if j coords are "above and to the right" of i
+
+    for (j = bins[i]; j >= 0; j = bins[j]) {
+      if (j >= nlocal) {
+	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;
+      }
+
+      if (exclude && exclusion(i,j,type,mask,molecule)) continue;
+
+      jtype = 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 <= 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;
+      }
+    }
+
+    // loop over all atoms in other bins in stencil, store every pair
+    // skip if i,j neighbor cutoff is less than bin distance
+
+    ibin = coord2bin(x[i]);
+    s = mstencil[itype];
+    dist = mdist[itype];
+    cut = cutneigh[itype];
+    ms = mstencils[itype];
+    for (k = 0; k < ms; k++) {
+      for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
+	jtype = type[j];
+	if (cut[jtype] < dist[k]) 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[i] = neighptr;
+    numneigh[i] = n;
+    npnt += n;
+    if (npnt >= pgsize)
+      error->one("Neighbor list overflow, boost neigh_modify one or page");
+  }
+}
+
 /* ----------------------------------------------------------------------
    binned neighbor list construction with Newton's 3rd law for triclinic
    each owned atom i checks its own bin and other bins in triclinic stencil
@@ -429,6 +628,94 @@ void Neighbor::half_bin_newton_tri()
   }
 }
 
+/* ----------------------------------------------------------------------
+   binned neighbor list construction with Newton's 3rd law for triclinic
+   each owned atom i checks its own bin and other bins in triclinic stencil
+   multi-type stencil is itype dependent and is distance checked
+   every pair stored exactly once by some processor
+------------------------------------------------------------------------- */
+
+void Neighbor::half_bin_newton_multi_tri()
+{
+  int i,j,k,n,itype,jtype,ibin,which,ms;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *neighptr,*s;
+  double *cut,*dist;
+
+  // 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) add_pages(npage);
+    }
+
+    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 bins, including self, in stencil
+    // skip if i,j neighbor cutoff is less than bin distance
+    // bins below self are excluded from stencil
+    // pairs for atoms j below i are excluded
+
+    ibin = coord2bin(x[i]);
+    s = mstencil[itype];
+    dist = mdist[itype];
+    cut = cutneigh[itype];
+    ms = mstencils[itype];
+    for (k = 0; k < ms; k++) {
+      for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
+	jtype = type[j];
+	if (cut[jtype] < dist[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;
+	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[i] = neighptr;
+    numneigh[i] = n;
+    npnt += n;
+    if (npnt >= pgsize)
+      error->one("Neighbor list overflow, boost neigh_modify one or page");
+  }
+}
+
 /* ----------------------------------------------------------------------
    build half list from full list
    pair stored once if i,j are both owned and i < j
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index 41f6efd49457440a10a3e150c4101866e03f76b5..67e6739869e95784904d40f208998b4862758c3a 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -39,10 +39,13 @@ using namespace LAMMPS_NS;
 #define LB_FACTOR 1.5
 #define SMALL 1.0e-6
 #define EXDELTA 1
+#define BIG 1.0e20
 
 #define MIN(a,b) ((a) < (b) ? (a) : (b))
 #define MAX(a,b) ((a) > (b) ? (a) : (b))
 
+enum{NSQ,BIN,MULTI};
+
 /* ---------------------------------------------------------------------- */
 
 Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
@@ -50,7 +53,7 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
 
-  style = 1;
+  style = BIN;
   every = 1;
   delay = 10;
   dist_check = 1;
@@ -59,6 +62,7 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
 
   maxlocal = 0;
     
+  cutneigh = NULL;
   cutneighsq = NULL;
   fixchecklist = NULL;
 
@@ -91,6 +95,8 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
   maxstencil_full = 0;
   stencil_full = NULL;
 
+  //mstencil = NULL;
+  
   // half neighbor list info
 
   half = half_command = 0;
@@ -143,6 +149,7 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
 
 Neighbor::~Neighbor()
 {
+  memory->destroy_2d_double_array(cutneigh);
   memory->destroy_2d_double_array(cutneighsq);
   delete [] fixchecklist;
   memory->destroy_2d_double_array(xhold);
@@ -217,37 +224,43 @@ void Neighbor::init()
   // ------------------------------------------------------------------
   // settings
 
-  // set cutneigh and trigger distance for reneighboring
+  // set neighbor cutoffs (force cutoff + skin)
+  // trigger determines when atoms migrate and neighbor lists are rebuilt
+  // cutneigh and cutneighsq determine what pairs go into neighbor list
+  //   set to 0 if cutforce = 0
+  // cutneighmin/max used for neighbor bin sizes
+  // cutghost determines comm distance = max of cutneigh & skin
+  //   may need ghosts for bonds even if all cutneigh = 0 (pair = NULL)
 
-  if (force->pair) cutneigh = force->pair->cutforce + skin;
-  else cutneigh = skin;
   triggersq = 0.25*skin*skin;
-  if (cutneighsq == NULL)
-    cutneighsq = memory->create_2d_double_array(atom->ntypes+1,atom->ntypes+1,
-						"neigh:cutneighsq");
-
-  // set neighbor cutoffs with skin included
-  // if no pair defined, cutneigh is just skin
 
   n = atom->ntypes;
-
-  if (force->pair) {
-    double cutoff;
-    double **cutsq = force->pair->cutsq;
-    for (i = 1; i <= n; i++)
-      for (j = i; j <= n; j++) {
-	cutoff = sqrt(cutsq[i][j]);
-	cutneighsq[i][j] = (cutoff+skin) * (cutoff+skin);
-	cutneighsq[j][i] = cutneighsq[i][j];
-      }
-  } else {
-    for (i = 1; i <= n; i++)
-      for (j = i; j <= n; j++) {
-	cutneighsq[i][j] = skin*skin;
-	cutneighsq[j][i] = cutneighsq[i][j];
-      }
+  if (cutneigh == 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");
   }
 
+  double cutoff,delta;
+  for (i = 1; i <= n; i++)
+    for (j = i; 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;
+      cutneighsq[i][j] = cutneighsq[j][i] = (cutoff+delta) * (cutoff+delta);
+    }
+
+  cutneighmin = BIG;
+  cutneighmax = 0.0;
+  cutghost = skin;
+  for (i = 1; i <= n; i++)
+    for (j = i; j <= n; j++) {
+      cutneighmin = MIN(cutneighmin,cutneigh[i][j]);
+      cutneighmax = MAX(cutneighmax,cutneigh[i][j]);
+      cutghost = MAX(cutghost,cutneigh[i][j]);
+    }
+
   // check other classes that can induce reneighboring in decide()
 
   restart_check = 0;
@@ -302,7 +315,7 @@ void Neighbor::init()
     maxhold = 0;
   }
 
-  if (style == 0) {
+  if (style == NSQ) {
     memory->sfree(bins);
     memory->sfree(binhead);
     memory->sfree(stencil);
@@ -319,7 +332,7 @@ void Neighbor::init()
     }
   }
 
-  if (style == 1) {
+  if (style == BIN) {
     if (maxbin == 0) {
       maxbin = atom->nmax;
       bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins");
@@ -546,36 +559,36 @@ void Neighbor::init()
     cut_middle_inside_sq = (cut_respa[0] - skin) * (cut_respa[0] - skin);
   }
 
-  // set ptrs to correct half, full, triclinic build functions
+  // set ptrs to correct half, full, multi, triclinic build functions
   // cannot combine granular and rRESPA
 
   if (half) {
     if (atom->check_style("granular")) {
-      if (style == 0) {
+      if (style == NSQ) {
 	if (force->newton_pair == 0) 
 	  half_build = &Neighbor::granular_nsq_no_newton;
 	else half_build = &Neighbor::granular_nsq_newton;
-      } else if (style == 1) {
+      } else if (style == BIN) {
 	if (force->newton_pair == 0) 
 	  half_build = &Neighbor::granular_bin_no_newton;
 	else if (triclinic)
 	  half_build = &Neighbor::granular_bin_newton_tri;
 	else half_build = &Neighbor::granular_bin_newton;
-      }
+      } else error->all("Neighbor multi not allowed with granular");
     } else if (respa) {
-      if (style == 0) {
+      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 == 1) {
+      } else if (style == BIN) {
 	if (force->newton_pair == 0) 
 	  half_build = &Neighbor::respa_bin_no_newton;
 	else if (triclinic)
 	  half_build = &Neighbor::respa_bin_newton_tri;
 	else half_build = &Neighbor::respa_bin_newton;
-      }
+      } else error->all("Neighbor multi not allowed with rRESPA");
     } else {
-      if (style == 0) {
+      if (style == NSQ) {
 	if (force->newton_pair == 0) {
 	  if (full_every) half_build = &Neighbor::half_full_no_newton;
 	  else half_build = &Neighbor::half_nsq_no_newton;
@@ -583,23 +596,33 @@ void Neighbor::init()
 	  if (full_every) half_build = &Neighbor::half_full_newton;
 	  else half_build = &Neighbor::half_nsq_newton;
 	}
-      } else if (style == 1) {
+      } else if (style == BIN) {
 	if (force->newton_pair == 0) {
 	  if (full_every) half_build = &Neighbor::half_full_no_newton;
 	  else half_build = &Neighbor::half_bin_no_newton;
 	} else {
 	  if (full_every) half_build = &Neighbor::half_full_newton;
-	  else if (triclinic)
-	    half_build = &Neighbor::half_bin_newton_tri;
+	  else if (triclinic) half_build = &Neighbor::half_bin_newton_tri;
 	  else half_build = &Neighbor::half_bin_newton;
 	}
+      } else if (style == MULTI) {
+	if (force->newton_pair == 0) {
+	  if (full_every) half_build = &Neighbor::half_full_no_newton;
+	  else half_build = &Neighbor::half_bin_no_newton_multi;
+	} else {
+	  if (full_every) half_build = &Neighbor::half_full_newton;
+	  else if (triclinic)
+	    half_build = &Neighbor::half_bin_newton_multi_tri;
+	  else half_build = &Neighbor::half_bin_newton_multi;
+	}
       }
     }
   } else half_build = NULL;
 
   if (full) {
-    if (style == 0) full_build = &Neighbor::full_nsq;
-    else full_build = &Neighbor::full_bin;
+    if (style == NSQ) full_build = &Neighbor::full_nsq;
+    else if (style == BIN) full_build = &Neighbor::full_bin;
+    else error->all("Neighbor multi not allowed with full neighbor lists");
   } else full_build = NULL;
 
   // ------------------------------------------------------------------
@@ -828,7 +851,7 @@ void Neighbor::build()
 
   // extend bin list if necessary
 
-  if (style && atom->nmax > maxbin) {
+  if (style != NSQ && atom->nmax > maxbin) {
     maxbin = atom->nmax;
     memory->sfree(bins);
     bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins");
@@ -890,33 +913,34 @@ void Neighbor::setup_bins()
 {
   // bbox lo/hi = bounding box of entire domain
   // bbox = size of bbox of entire domain
-  // for triclinic, bbox bounds all 8 corners of tilted box
-  // bsubbox lo/hi = bounding box of my subdomain extended by cutneigh
-  // for triclinic, subdomain with cutneigh extension is first computed
-  //   in lamda coords, including tilt factor via cutcomm,
-  //   domain->bbox then computes bbox of corner pts transformed to box coords
+  // bsubbox lo/hi = bounding box of my subdomain extended by ghost atoms
+  // for triclinic:
+  //   bbox bounds all 8 corners of tilted box
+  //   subdomain is in lamda coords
+  //   include dimension-dependent extension via comm->cutghost
+  //   domain->bbox() converts lamda extent to box coords and computes bbox
 
   double bbox[3],bsubboxlo[3],bsubboxhi[3];
 
   if (triclinic == 0) {
     bboxlo = domain->boxlo;
     bboxhi = domain->boxhi;
-    bsubboxlo[0] = domain->sublo[0] - cutneigh;
-    bsubboxlo[1] = domain->sublo[1] - cutneigh;
-    bsubboxlo[2] = domain->sublo[2] - cutneigh;
-    bsubboxhi[0] = domain->subhi[0] + cutneigh;
-    bsubboxhi[1] = domain->subhi[1] + cutneigh;
-    bsubboxhi[2] = domain->subhi[2] + cutneigh;
+    bsubboxlo[0] = domain->sublo[0] - cutghost;
+    bsubboxlo[1] = domain->sublo[1] - cutghost;
+    bsubboxlo[2] = domain->sublo[2] - cutghost;
+    bsubboxhi[0] = domain->subhi[0] + cutghost;
+    bsubboxhi[1] = domain->subhi[1] + cutghost;
+    bsubboxhi[2] = domain->subhi[2] + cutghost;
   } else {
     bboxlo = domain->boxlo_bound;
     bboxhi = domain->boxhi_bound;
     double lo[3],hi[3];
-    lo[0] = domain->sublo_lamda[0] - comm->cutcomm[0];
-    lo[1] = domain->sublo_lamda[1] - comm->cutcomm[1];
-    lo[2] = domain->sublo_lamda[2] - comm->cutcomm[2];
-    hi[0] = domain->subhi_lamda[0] + comm->cutcomm[0];
-    hi[1] = domain->subhi_lamda[1] + comm->cutcomm[1];
-    hi[2] = domain->subhi_lamda[2] + comm->cutcomm[2];
+    lo[0] = domain->sublo_lamda[0] - comm->cutghost[0];
+    lo[1] = domain->sublo_lamda[1] - comm->cutghost[1];
+    lo[2] = domain->sublo_lamda[2] - comm->cutghost[2];
+    hi[0] = domain->subhi_lamda[0] + comm->cutghost[0];
+    hi[1] = domain->subhi_lamda[1] + comm->cutghost[1];
+    hi[2] = domain->subhi_lamda[2] + comm->cutghost[2];
     domain->bbox(lo,hi,bsubboxlo,bsubboxhi);
   }
 
@@ -924,22 +948,30 @@ void Neighbor::setup_bins()
   bbox[1] = bboxhi[1] - bboxlo[1];
   bbox[2] = bboxhi[2] - bboxlo[2];
 
-  // test for too many global bins in any dimension due to huge domain
+  // for BIN style, binsize = 1/2 of max neighbor cutoff
+  // for MULTI style, binsize = 1/2 of min neighbor cutoff
+  // special case of all cutoffs = 0.0, binsize = box size
 
-  double cutneighinv = 1.0/cutneigh;
+  double binsize;
+  if (style == BIN) binsize = 0.5*cutneighmax;
+  else binsize = 0.5*cutneighmin;
+  if (binsize == 0.0) binsize = bbox[0];
+
+  // test for too many global bins in any dimension due to huge domain
 
-  if (2.0*bbox[0]*cutneighinv > INT_MAX || 2.0*bbox[1]*cutneighinv > INT_MAX ||
-      2.0*bbox[2]*cutneighinv > INT_MAX)
+  double binsizeinv = 1.0/binsize;
+  if (bbox[0]*binsizeinv > INT_MAX || bbox[1]*binsizeinv > INT_MAX ||
+      bbox[2]*binsizeinv > INT_MAX)
     error->all("Domain too large for neighbor bins");
 
   // divide box into bins
   // optimal size is roughly 1/2 the cutoff
   // use one bin even if cutoff >> bbox
 
-  nbinx = static_cast<int> (2.0*bbox[0]*cutneighinv);
-  nbiny = static_cast<int> (2.0*bbox[1]*cutneighinv);
+  nbinx = static_cast<int> (bbox[0]*binsizeinv);
+  nbiny = static_cast<int> (bbox[1]*binsizeinv);
   if (force->dimension == 3)
-    nbinz = static_cast<int> (2.0*bbox[2]*cutneighinv);
+    nbinz = static_cast<int> (bbox[2]*binsizeinv);
   else nbinz = 1;
 
   if (nbinx == 0) nbinx = 1;
@@ -1009,7 +1041,8 @@ void Neighbor::setup_bins()
   }
 
   // create stencil of bins whose closest corner to central bin
-  //   is within neighbor cutoff
+  //   is within cutneighmax
+  // stencil will be empty if cutneighmax = 0.0
   // next(xyz) = how far the stencil could possibly extend
   // for partial Newton (newton = 0)
   //   stencil is all surrounding bins
@@ -1026,12 +1059,12 @@ void Neighbor::setup_bins()
   //   stored in stencil_full
   // 3d creates xyz stencil, 2d is only xy
 
-  int nextx = static_cast<int> (cutneigh*bininvx);
-  if (nextx*binsizex < cutneigh) nextx++;
-  int nexty = static_cast<int> (cutneigh*bininvy);
-  if (nexty*binsizey < cutneigh) nexty++;
-  int nextz = static_cast<int> (cutneigh*bininvz);
-  if (nextz*binsizez < cutneigh) nextz++;
+  int nextx = static_cast<int> (cutneighmax*bininvx);
+  if (nextx*binsizex < cutneighmax) nextx++;
+  int nexty = static_cast<int> (cutneighmax*bininvy);
+  if (nexty*binsizey < cutneighmax) nexty++;
+  int nextz = static_cast<int> (cutneighmax*bininvz);
+  if (nextz*binsizez < cutneighmax) nextz++;
 
   int nmax = (2*nextz+1) * (2*nexty+1) * (2*nextx+1);
   if (nmax > maxstencil) {
@@ -1042,45 +1075,45 @@ void Neighbor::setup_bins()
 
   int i,j,k;
   nstencil = 0;
-  double cutsq = cutneigh*cutneigh;
+  double cutneighmaxsq = cutneighmax*cutneighmax;
 
   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) < cutsq)
+	    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) < cutsq)
+	    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) < cutsq)
+	      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) < cutsq)
+	  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) < cutsq)
+	    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) < cutsq)
+	      if (bin_distance(i,j,0) < cutneighmaxsq)
 		stencil[nstencil++] = j*mbinx + i;
     }
   }
@@ -1097,12 +1130,12 @@ void Neighbor::setup_bins()
       for (k = -nextz; k <= nextz; k++)
 	for (j = -nexty; j <= nexty; j++)
 	  for (i = -nextx; i <= nextx; i++)
-	    if (bin_distance(i,j,k) < cutsq)
+	    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) < cutsq)
+	  if (bin_distance(i,j,0) < cutneighmaxsq)
 	    stencil_full[nstencil_full++] = j*mbinx + i;
     }
   }
@@ -1140,6 +1173,23 @@ double Neighbor::bin_distance(int i, int j, int k)
   return (delx*delx + dely*dely + delz*delz);
 }
 
+/* ----------------------------------------------------------------------
+   set neighbor style and skin distance
+------------------------------------------------------------------------- */
+
+void Neighbor::set(int narg, char **arg)
+{
+  if (narg != 2) error->all("Illegal neighbor command");
+
+  skin = atof(arg[0]);
+  if (skin < 0.0) error->all("Illegal neighbor command");
+
+  if (strcmp(arg[1],"nsq") == 0) style = NSQ;
+  else if (strcmp(arg[1],"bin") == 0) style = BIN;
+  else if (strcmp(arg[1],"multi") == 0) style = MULTI;
+  else error->all("Illegal neighbor command");
+}
+
 /* ----------------------------------------------------------------------
    modify parameters of the pair-wise neighbor build
 ------------------------------------------------------------------------- */
@@ -1245,7 +1295,7 @@ int Neighbor::memory_usage()
 
   bytes += maxhold*3 * sizeof(double);
 
-  if (style == 0) {
+  if (style == NSQ) {
     bytes += maxbin * sizeof(int);
     bytes += maxhead * sizeof(int);
     bytes += maxstencil * sizeof(int);
diff --git a/src/neighbor.h b/src/neighbor.h
index 73c9afcf3c821719164eb75c4ffea76d46febc49..99d55d9c3b61bca420d17091e4377ef36c4177b9 100644
--- a/src/neighbor.h
+++ b/src/neighbor.h
@@ -29,7 +29,7 @@ class Neighbor : protected Pointers {
   int oneatom;                     // max # of neighbors for one atom
 
   double skin;                     // skin distance
-  double cutneigh;                 // neighbor cutoff
+  double cutghost;                 // distance for acquiring ghost atoms
 
   int ncalls;                      // # of times build has been called
   int ndanger;                     // # of dangerous builds
@@ -79,6 +79,7 @@ class Neighbor : protected Pointers {
   void build();               // create all neighbor lists (half,full,bond)
   void build_half();          // one-time creation of half neighbor list
   void build_full();          // one-time creation of full neighbor list
+  void set(int, char **);     // set neighbor style and skin distance
   void modify_params(int, char**);  // modify parameters of neighbor build
   int memory_usage();         // tally memory usage
   
@@ -93,8 +94,11 @@ class Neighbor : protected Pointers {
   int fix_check;                   // # of fixes that induce reneigh
   int *fixchecklist;               // which fixes to check
 
-  double **cutneighsq;             // neighbor cutoff squared for type pairs
-  double triggersq;                // sq distance to trigger build on
+  double **cutneigh;               // neighbor cutoff for type pairs
+  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 **xhold;                  // atom coords at last neighbor build
   int maxhold;                     // size of xhold array
@@ -106,10 +110,6 @@ class Neighbor : protected Pointers {
   int *binhead;                    // ptr to 1st atom in each bin
   int maxhead;                     // size of binhead array
 
-  int nstencil;                    // # of bins in stencil
-  int *stencil;                    // stencil list of bin offsets
-  int maxstencil;                  // size of stencil
-
   int mbins;                       // # of local bins and offset
   int mbinx,mbiny,mbinz;
   int mbinxlo,mbinylo,mbinzlo;
@@ -120,16 +120,24 @@ class Neighbor : protected Pointers {
   int triclinic;                   // 0 if domain is orthog, 1 if triclinic
   double *bboxlo,*bboxhi;          // copy of full domain bounding box
 
+  int nstencil;                    // # of bins in half neighbor stencil
+  int *stencil;                    // list of bin offsets
+  int maxstencil;                  // 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 *mstencils;                  // # bins in each type-based multi stencil
+  int **mstencil;                  // list of bin offsets in each stencil
+  double **mdist;                  // list of distances to bins in each stencil
+
   int **pages;                     // half neighbor list pages
   int maxpage;                     // # of half pages currently allocated
 
   int **pages_full;                // full neighbor list pages
   int maxpage_full;                // # of full pages currently allocated
 
-  int nstencil_full;               // # of bins in full neighbor stencil
-  int *stencil_full;               // full neighbor stencil list of bin offsets
-  int maxstencil_full;             // size of full neighbor stencil
-
                                          // granular neighbor list
   class FixShearHistory *fix_history;    // NULL if history not needed
                                          // else is ptr to fix shear history
@@ -173,25 +181,31 @@ class Neighbor : protected Pointers {
   FnPtr half_build;                   // ptr to half pair list functions
   FnPtr full_build;                   // ptr to full pair list functions
 
-  void half_nsq_no_newton();          // via nsq w/ pair newton off
-  void half_nsq_newton();             // via nsq w/ pair newton on
-  void half_bin_no_newton();          // via binning w/ pair newton off
-  void half_bin_newton();             // via binning w/ pair newton on
-  void half_full_no_newton();         // via full list w/ pair newton off
-  void half_full_newton();            // via full list w/ pair newton on
+  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();                    // 2 fns for full neighbor lists
+  void full_nsq();                    // fns for full neighbor lists
   void full_bin();
 
-  void granular_nsq_no_newton();      // 4 fns for granular systems
+  void granular_nsq_no_newton();      // fns for granular neighbor lists
   void granular_nsq_newton();
   void granular_bin_no_newton();
   void granular_bin_newton();
+  void granular_bin_newton_tri();
 
-  void respa_nsq_no_newton();         // 4 fns for multiple respa lists
+  void respa_nsq_no_newton();         // fns for respa neighbor lists
   void respa_nsq_newton();
   void respa_bin_no_newton();
   void respa_bin_newton();
+  void respa_bin_newton_tri();
 
   FnPtr bond_build;                   // ptr to bond list functions
   void bond_all();                    // bond list with all bonds
@@ -220,12 +234,6 @@ class Neighbor : protected Pointers {
   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
-
-  // PJV
-
-  void granular_bin_newton_tri();
-  void respa_bin_newton_tri();
-  void half_bin_newton_tri();
 };
 
 }
diff --git a/src/pair.h b/src/pair.h
index 5479078f94fcaf063b8aac060384baa01087a4dd..00f10c03972e5b0ccf3605809ac33b90ef597dcf 100644
--- a/src/pair.h
+++ b/src/pair.h
@@ -79,7 +79,7 @@ class Pair : protected Pointers {
   virtual void write_restart_settings(FILE *) {}
   virtual void read_restart_settings(FILE *) {}
 
-  virtual int pack_comm(int, int *, double *, int, double *) {return 0;}
+  virtual int pack_comm(int, int *, double *, int, int *) {return 0;}
   virtual void unpack_comm(int, int, double *) {}
   virtual int pack_reverse_comm(int, int, double *) {return 0;}
   virtual void unpack_reverse_comm(int, int *, double *) {}
diff --git a/src/velocity.cpp b/src/velocity.cpp
index 7f84f3e239eed6d191eeff1e17be3ccf87e58eb6..e2bf4b9b3a78745986a71a79e57631b3efa410f9 100644
--- a/src/velocity.cpp
+++ b/src/velocity.cpp
@@ -664,8 +664,8 @@ void Velocity::options(int narg, char **arg)
 void Velocity::triple(double *x, double *vx, double *vy, double *vz,
 		      int seed, RanPark *random)
 {
-  // lamda[3] = fractional position in box
-  // for triclinic box, just convert to lamda coords
+  // for orthogonal box, lamda = fractional position in box
+  // for triclinic box, convert to lamda coords
 
   double lamda[3];