diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp
index 4215584e4e..61ab0cd0d2 100644
--- a/src/compute_pair_local.cpp
+++ b/src/compute_pair_local.cpp
@@ -138,12 +138,14 @@ void ComputePairLocal::compute_local()
 int ComputePairLocal::compute_pairs(int flag)
 {
   int i,j,m,n,ii,jj,inum,jnum,itype,jtype;
+  tagint itag,jtag;
   double xtmp,ytmp,ztmp,delx,dely,delz;
   double rsq,eng,fpair,factor_coul,factor_lj;
   int *ilist,*jlist,*numneigh,**firstneigh;
   double *ptr;
 
   double **x = atom->x;
+  tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -162,6 +164,9 @@ int ComputePairLocal::compute_pairs(int flag)
 
   // loop over neighbors of my atoms
   // skip if I or J are not in group
+  // for newton = 0 and J = ghost atom,
+  //   need to insure I,J pair is only output by one proc
+  //   use same itag,jtag lobic as in Neighbor::neigh_half_nsq()
   // for flag = 0, just count pair interactions within force cutoff
   // for flag = 1, calculate requested output fields
 
@@ -176,6 +181,7 @@ int ComputePairLocal::compute_pairs(int flag)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
+    itag = tag[i];
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
@@ -187,7 +193,23 @@ int ComputePairLocal::compute_pairs(int flag)
       j &= NEIGHMASK;
 
       if (!(mask[j] & groupbit)) continue;
-      if (newton_pair == 0 && j >= nlocal) continue;
+
+      // itag = jtag is possible for long cutoffs that include images of self
+
+      if (newton_pair == 0 && j >= nlocal) {
+        jtag = tag[j];
+        if (itag > jtag) {
+          if ((itag+jtag) % 2 == 0) continue;
+        } else if (itag < jtag) {
+          if ((itag+jtag) % 2 == 1) continue;
+        } else {
+          if (x[j][2] < ztmp) continue;
+          if (x[j][2] == ztmp) {
+            if (x[j][1] < ytmp) continue;
+            if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+          }
+        }
+      }
 
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp
index 3589762a45..17ce63c775 100644
--- a/src/neigh_half_nsq.cpp
+++ b/src/neigh_half_nsq.cpp
@@ -251,8 +251,8 @@ void Neighbor::half_nsq_no_newton_ghost(NeighList *list)
 
 void Neighbor::half_nsq_newton(NeighList *list)
 {
-  int i,j,n,itype,jtype,itag,jtag,which,bitmask,imol,iatom,moltemplate;
-  tagint tagprev;
+  int i,j,n,itype,jtype,which,bitmask,imol,iatom,moltemplate;
+  tagint itag,jtag,tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr;