bugfix for 2 recenty reported neighbor issues, also a ReaxFF fix species issue

doc/src/pair_yukawa_colloid.txt

@@ -35,7 +35,7 @@ cutoff.
 In contrast to "pair_style yukawa"_pair_yukawa.html, this functional
 form arises from the Coulombic interaction between two colloid
 particles, screened due to the presence of an electrolyte, see the
-book by "Safran"_#Safran for a derivation in the context of DVLO
+book by "Safran"_#Safran for a derivation in the context of DLVO
 theory.  "Pair_style yukawa"_pair_yukawa.html is a screened Coulombic
 potential between two point-charges and uses no such approximation.
 

src/QEQ/fix_qeq.cpp

@@ -286,9 +286,6 @@ void FixQEq::setup_pre_force(int vflag)
   if (force->newton_pair == 0)
     error->all(FLERR,"QEQ with 'newton pair off' not supported");
 
-  // should not be needed
-  // neighbor->build_one(list);
-
   deallocate_storage();
   allocate_storage();
 

src/USER-DPD/nbin_ssa.cpp

@@ -76,6 +76,7 @@ void NBinSSA::bin_atoms()
       if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
       if (mask[i] & bitmask) {
         ibin = coord2bin(x[i]);
+        atom2bin[i] = ibin;
         bins_ssa[i] = gbinhead_ssa[ibin];
         gbinhead_ssa[ibin] = i;
       }
@@ -84,12 +85,14 @@ void NBinSSA::bin_atoms()
     for (i = nall-1; i >= nlocal; i--) {
       if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
       ibin = coord2bin(x[i]);
+      atom2bin[i] = ibin;
       bins_ssa[i] = gbinhead_ssa[ibin];
       gbinhead_ssa[ibin] = i;
     }
   }
   for (i = nlocal-1; i >= 0; i--) {
     ibin = coord2bin(x[i]);
+    atom2bin[i] = ibin;
     bins_ssa[i] = binhead_ssa[ibin];
     binhead_ssa[ibin] = i;
   }

src/USER-DPD/npair_half_bin_newton_ssa.cpp

@@ -141,7 +141,7 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
       }
     }
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
 
     // loop over all local atoms in other bins in "half" stencil
 

src/USER-OMP/npair_full_bin_ghost_omp.cpp

@@ -97,7 +97,7 @@ void NPairFullBinGhostOmp::build(NeighList *list)
     // no molecular test when i = ghost atom
 
     if (i < nlocal) {
-      ibin = coord2bin(x[i]);
+      ibin = atom2bin[i];
       for (k = 0; k < nstencil; k++) {
         for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
           if (i == j) continue;

src/USER-OMP/npair_full_bin_omp.cpp

@@ -90,7 +90,7 @@ void NPairFullBinOmp::build(NeighList *list)
     // loop over all atoms in surrounding bins in stencil including self
     // skip i = j
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
 
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

src/USER-OMP/npair_full_multi_omp.cpp

@@ -94,7 +94,7 @@ void NPairFullMultiOmp::build(NeighList *list)
     // skip if i,j neighbor cutoff is less than bin distance
     // skip i = j
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     s = stencil_multi[itype];
     distsq = distsq_multi[itype];
     cutsq = cutneighsq[itype];

src/USER-OMP/npair_half_bin_newtoff_ghost_omp.cpp

@@ -103,7 +103,7 @@ void NPairHalfBinNewtoffGhostOmp::build(NeighList *list)
     // no molecular test when i = ghost atom
 
     if (i < nlocal) {
-      ibin = coord2bin(x[i]);
+      ibin = atom2bin[i];
 
       for (k = 0; k < nstencil; k++) {
         for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

src/USER-OMP/npair_half_bin_newtoff_omp.cpp

@@ -94,7 +94,7 @@ void NPairHalfBinNewtoffOmp::build(NeighList *list)
     // stores own/own pairs only once
     // stores own/ghost pairs on both procs
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
 
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

src/USER-OMP/npair_half_bin_newton_omp.cpp

@@ -130,7 +130,7 @@ void NPairHalfBinNewtonOmp::build(NeighList *list)
 
     // loop over all atoms in other bins in stencil, store every pair
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         jtype = type[j];

src/USER-OMP/npair_half_bin_newton_tri_omp.cpp

@@ -94,7 +94,7 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list)
     //         (equal zyx and j <= i)
     // latter excludes self-self interaction but allows superposed atoms
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         if (x[j][2] < ztmp) continue;

src/USER-OMP/npair_half_multi_newtoff_omp.cpp

@@ -97,7 +97,7 @@ void NPairHalfMultiNewtoffOmp::build(NeighList *list)
     // stores own/own pairs only once
     // stores own/ghost pairs on both procs
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     s = stencil_multi[itype];
     distsq = distsq_multi[itype];
     cutsq = cutneighsq[itype];

src/USER-OMP/npair_half_multi_newton_omp.cpp

@@ -131,7 +131,7 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list)
     // 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]);
+    ibin = atom2bin[i];
     s = stencil_multi[itype];
     distsq = distsq_multi[itype];
     cutsq = cutneighsq[itype];

src/USER-OMP/npair_half_multi_newton_tri_omp.cpp

@@ -99,7 +99,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list)
     //         (equal zyx and j <= i)
     // latter excludes self-self interaction but allows superposed atoms
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     s = stencil_multi[itype];
     distsq = distsq_multi[itype];
     cutsq = cutneighsq[itype];

src/USER-OMP/npair_half_respa_bin_newtoff_omp.cpp

@@ -117,7 +117,7 @@ void NPairHalfRespaBinNewtoffOmp::build(NeighList *list)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     if (moltemplate) {
       imol = molindex[i];
       iatom = molatom[i];

src/USER-OMP/npair_half_respa_bin_newton_omp.cpp

@@ -176,7 +176,7 @@ void NPairHalfRespaBinNewtonOmp::build(NeighList *list)
 
     // loop over all atoms in other bins in stencil, store every pair
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         jtype = type[j];

src/USER-OMP/npair_half_respa_bin_newton_tri_omp.cpp

@@ -128,7 +128,7 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list)
     //         (equal zyx and j <= i)
     // latter excludes self-self interaction but allows superposed atoms
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         if (x[j][2] < ztmp) continue;

src/USER-OMP/npair_half_size_bin_newtoff_omp.cpp

@@ -113,7 +113,7 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
 
     // loop over all atoms in surrounding bins in stencil including self
     // only store pair if i < j

src/USER-OMP/npair_half_size_bin_newton_omp.cpp

@@ -168,7 +168,7 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list)
 
     // loop over all atoms in other bins in stencil, store every pair
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         if (exclude && exclusion(i,j,type[i],type[j],mask,molecule)) continue;

src/USER-OMP/npair_half_size_bin_newton_tri_omp.cpp

@@ -84,7 +84,7 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
     //         (equal zyx and j <= i)
     // latter excludes self-self interaction but allows superposed atoms
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         if (x[j][2] < ztmp) continue;

src/USER-REAXC/fix_qeq_reax.cpp

@@ -416,9 +416,6 @@ void FixQEqReax::init_taper()
 
 void FixQEqReax::setup_pre_force(int vflag)
 {
-  // should not be needed
-  // neighbor->build_one(list);
-
   deallocate_storage();
   allocate_storage();
 

src/USER-REAXC/fix_reaxc_species.cpp

@@ -500,7 +500,7 @@ void FixReaxCSpecies::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp)
 
 AtomCoord FixReaxCSpecies::chAnchor(AtomCoord in1, AtomCoord in2)
 {
-  if (in1.x < in2.x)
+  if (in1.x <= in2.x)
     return in1;
   return in2;
 }

src/nbin.cpp

@@ -29,6 +29,7 @@ NBin::NBin(LAMMPS *lmp) : Pointers(lmp)
   maxbin = maxatom = 0;
   binhead = NULL;
   bins = NULL;
+  atom2bin = NULL;
 
   // geometry settings
 
@@ -42,6 +43,7 @@ NBin::~NBin()
 {
   memory->destroy(binhead);
   memory->destroy(bins);
+  memory->destroy(atom2bin);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -87,12 +89,15 @@ void NBin::bin_atoms_setup(int nall)
     memory->create(binhead,maxbin,"neigh:binhead");
   }
 
-  // bins = per-atom vector
+  // bins and atom2bin = per-atom vectors
+  // for both local and ghost atoms
 
   if (nall > maxatom) {
     maxatom = nall;
     memory->destroy(bins);
     memory->create(bins,maxatom,"neigh:bins");
+    memory->destroy(atom2bin);
+    memory->create(atom2bin,maxatom,"neigh:atom2bin");
   }
 }
 
@@ -148,6 +153,6 @@ bigint NBin::memory_usage()
 {
   bigint bytes = 0;
   bytes += maxbin*sizeof(int);
-  bytes += maxatom*sizeof(int);
+  bytes += 2*maxatom*sizeof(int);
   return bytes;
 }

src/nbin.h

@@ -33,6 +33,7 @@ class NBin : protected Pointers {
 
   int *binhead;                // index of first atom in each bin
   int *bins;                   // index of next atom in same bin
+  int *atom2bin;               // bin assignment for each atom (local+ghost)
 
   double cutoff_custom;        // cutoff set by requestor
 

src/nbin_standard.cpp

@@ -211,12 +211,14 @@ void NBinStandard::bin_atoms()
     for (i = nall-1; i >= nlocal; i--) {
       if (mask[i] & bitmask) {
         ibin = coord2bin(x[i]);
+        atom2bin[i] = ibin;
         bins[i] = binhead[ibin];
         binhead[ibin] = i;
       }
     }
     for (i = atom->nfirst-1; i >= 0; i--) {
       ibin = coord2bin(x[i]);
+      atom2bin[i] = ibin;
       bins[i] = binhead[ibin];
       binhead[ibin] = i;
     }
@@ -224,6 +226,7 @@ void NBinStandard::bin_atoms()
   } else {
     for (i = nall-1; i >= 0; i--) {
       ibin = coord2bin(x[i]);
+      atom2bin[i] = ibin;
       bins[i] = binhead[ibin];
       binhead[ibin] = i;
     }

src/neigh_request.cpp

@@ -50,6 +50,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   // default is no Intel-specific neighbor list build
   // default is no Kokkos neighbor list build
   // default is no Shardlow Splitting Algorithm (SSA) neighbor list build
+  // default is no list-specific cutoff
   // default is no storage of auxiliary floating point values
 
   occasional = 0;

src/neighbor.cpp

@@ -640,6 +640,24 @@ int Neighbor::init_pair()
   delete [] neigh_stencil;
   delete [] neigh_pair;
 
+  // error check on requests
+  // do not allow occasional, ghost, bin list
+  //   b/c it still uses variant of coord2bin() in NPair() method
+  //     instead of atom2bin, this could cause error b/c stoms have
+  //     moved out of proc domain by time occasional list is built
+  //   solution would be to use a different NBin variant
+  //     that used Npair::coord2bin(x,ix,iy,iz) (then delete it from NPair)
+  //     and stored the ix,iy,iz values for all atoms (including ghosts)
+  //     at time of binning when neighbor lists are rebuilt,
+  //     similar to what vanilla Nbin::coord2atom() does now in atom2bin
+
+  if (style == BIN) {
+    for (i = 0; i < nrequest; i++)
+      if (requests[i]->occasional && requests[i]->ghost)
+        error->all(FLERR,"Cannot request an occasional binned neighbor list "
+                   "with ghost info");
+  }
+
   // morph requests in various ways
   // purpose is to avoid duplicate or inefficient builds
   // may add new requests if a needed request to derive from does not exist
@@ -1669,7 +1687,6 @@ int Neighbor::choose_stencil(NeighRequest *rq)
   else if (rq->newton == 1) newtflag = 1;
   else if (rq->newton == 2) newtflag = 0;
 
-
   //printf("STENCIL RQ FLAGS: hff %d %d n %d g %d s %d newtflag %d\n",
   //       rq->half,rq->full,rq->newton,rq->ghost,rq->ssa,
   //       newtflag);
@@ -2087,7 +2104,7 @@ void Neighbor::build(int topoflag)
   }
 
   // bin atoms for all NBin instances
-  // not just NBin associated with perpetual lists
+  // not just NBin associated with perpetual lists, also occasional lists
   // b/c cannot wait to bin occasional lists in build_one() call
   // if bin then, atoms may have moved outside of proc domain & bin extent,
   //   leading to errors or even a crash
@@ -2193,6 +2210,7 @@ void Neighbor::build_one(class NeighList *mylist, int preflag)
 
   // build the list
 
+  if (!mylist->copy) mylist->grow(atom->nlocal,atom->nlocal+atom->nghost);
   np->build_setup();
   np->build(mylist);
 }

src/npair.cpp

@@ -128,6 +128,7 @@ void NPair::copy_bin_info()
   bininvy = nb->bininvy;
   bininvz = nb->bininvz;
 
+  atom2bin = nb->atom2bin;
   bins = nb->bins;
   binhead = nb->binhead;
 }
@@ -198,53 +199,8 @@ int NPair::exclusion(int i, int j, int itype, int jtype,
 }
 
 /* ----------------------------------------------------------------------
-   convert atom coords into local bin #
+   same as coord2bin in Nbin, but also return ix,iy,iz offsets in each dim
-   for orthogonal, only ghost atoms will have coord >= bboxhi or coord < bboxlo
+   used by some of the ghost neighbor lists
-     take special care to insure ghosts are in correct bins even w/ roundoff
-     hi ghost atoms = nbin,nbin+1,etc
-     owned atoms = 0 to nbin-1
-     lo ghost atoms = -1,-2,etc
-     this is necessary so that both procs on either side of PBC
-       treat a pair of atoms straddling the PBC in a consistent way
-   for triclinic, doesn't matter since stencil & neigh list built differently
-------------------------------------------------------------------------- */
-
-int NPair::coord2bin(double *x)
-{
-  int ix,iy,iz;
-
-  if (!ISFINITE(x[0]) || !ISFINITE(x[1]) || !ISFINITE(x[2]))
-    error->one(FLERR,"Non-numeric positions - simulation unstable");
-
-  if (x[0] >= bboxhi[0])
-    ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx) + nbinx;
-  else if (x[0] >= bboxlo[0]) {
-    ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx);
-    ix = MIN(ix,nbinx-1);
-  } else
-    ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - 1;
-
-  if (x[1] >= bboxhi[1])
-    iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy) + nbiny;
-  else if (x[1] >= bboxlo[1]) {
-    iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy);
-    iy = MIN(iy,nbiny-1);
-  } else
-    iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - 1;
-
-  if (x[2] >= bboxhi[2])
-    iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz) + nbinz;
-  else if (x[2] >= bboxlo[2]) {
-    iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz);
-    iz = MIN(iz,nbinz-1);
-  } else
-    iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - 1;
-
-  return (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo);
-}
-
-/* ----------------------------------------------------------------------
-   same as coord2bin, but also return ix,iy,iz offsets in each dim
 ------------------------------------------------------------------------- */
 
 int NPair::coord2bin(double *x, int &ix, int &iy, int &iz)

src/npair.h

@@ -77,7 +77,7 @@ class NPair : protected Pointers {
   int mbinx,mbiny,mbinz;
   int mbinxlo,mbinylo,mbinzlo;
   double bininvx,bininvy,bininvz;
-  int *bins;
+  int *atom2bin,*bins;
   int *binhead;
   
   // data from NStencil class

src/npair_full_bin.cpp

@@ -80,7 +80,7 @@ void NPairFullBin::build(NeighList *list)
     // loop over all atoms in surrounding bins in stencil including self
     // skip i = j
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
 
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

src/npair_full_bin_atomonly.cpp

@@ -64,7 +64,7 @@ void NPairFullBinAtomonly::build(NeighList *list)
     // loop over all atoms in surrounding bins in stencil including self
     // skip i = j
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
 
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

src/npair_full_bin_ghost.cpp

@@ -87,7 +87,7 @@ void NPairFullBinGhost::build(NeighList *list)
     // no molecular test when i = ghost atom
 
     if (i < nlocal) {
-      ibin = coord2bin(x[i]);
+      ibin = atom2bin[i];
       for (k = 0; k < nstencil; k++) {
         for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
           if (i == j) continue;

src/npair_full_multi.cpp

@@ -83,7 +83,7 @@ void NPairFullMulti::build(NeighList *list)
     // skip if i,j neighbor cutoff is less than bin distance
     // skip i = j
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     s = stencil_multi[itype];
     distsq = distsq_multi[itype];
     cutsq = cutneighsq[itype];

src/npair_half_bin_atomonly_newton.cpp

@@ -90,7 +90,8 @@ void NPairHalfBinAtomonlyNewton::build(NeighList *list)
 
     // loop over all atoms in other bins in stencil, store every pair
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
+
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         jtype = type[j];

src/npair_half_bin_newtoff.cpp

@@ -84,7 +84,7 @@ void NPairHalfBinNewtoff::build(NeighList *list)
     // stores own/own pairs only once
     // stores own/ghost pairs on both procs
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
 
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

src/npair_half_bin_newtoff_ghost.cpp

@@ -92,7 +92,7 @@ void NPairHalfBinNewtoffGhost::build(NeighList *list)
     // no molecular test when i = ghost atom
 
     if (i < nlocal) {
-      ibin = coord2bin(x[i]);
+      ibin = atom2bin[i];
 
       for (k = 0; k < nstencil; k++) {
         for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

src/npair_half_bin_newton.cpp

@@ -119,7 +119,7 @@ void NPairHalfBinNewton::build(NeighList *list)
 
     // loop over all atoms in other bins in stencil, store every pair
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         jtype = type[j];

src/npair_half_bin_newton_tri.cpp

@@ -84,7 +84,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list)
     //         (equal zyx and j <= i)
     // latter excludes self-self interaction but allows superposed atoms
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         if (x[j][2] < ztmp) continue;

src/npair_half_multi_newtoff.cpp

@@ -87,7 +87,7 @@ void NPairHalfMultiNewtoff::build(NeighList *list)
     // stores own/own pairs only once
     // stores own/ghost pairs on both procs
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     s = stencil_multi[itype];
     distsq = distsq_multi[itype];
     cutsq = cutneighsq[itype];

src/npair_half_multi_newton.cpp

@@ -121,7 +121,7 @@ void NPairHalfMultiNewton::build(NeighList *list)
     // 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]);
+    ibin = atom2bin[i];
     s = stencil_multi[itype];
     distsq = distsq_multi[itype];
     cutsq = cutneighsq[itype];

src/npair_half_multi_newton_tri.cpp

@@ -88,7 +88,7 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
     //         (equal zyx and j <= i)
     // latter excludes self-self interaction but allows superposed atoms
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     s = stencil_multi[itype];
     distsq = distsq_multi[itype];
     cutsq = cutneighsq[itype];

src/npair_half_respa_bin_newtoff.cpp

@@ -101,7 +101,7 @@ void NPairHalfRespaBinNewtoff::build(NeighList *list)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     if (moltemplate) {
       imol = molindex[i];
       iatom = molatom[i];

src/npair_half_respa_bin_newton.cpp

@@ -160,7 +160,7 @@ void NPairHalfRespaBinNewton::build(NeighList *list)
 
     // loop over all atoms in other bins in stencil, store every pair
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         jtype = type[j];

src/npair_half_respa_bin_newton_tri.cpp

@@ -113,7 +113,7 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list)
     //         (equal zyx and j <= i)
     // latter excludes self-self interaction but allows superposed atoms
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         if (x[j][2] < ztmp) continue;

src/npair_half_size_bin_newtoff.cpp

@@ -105,7 +105,7 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list)
     ytmp = x[i][1];
     ztmp = x[i][2];
     radi = radius[i];
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
 
     // loop over all atoms in surrounding bins in stencil including self
     // only store pair if i < j

src/npair_half_size_bin_newton.cpp

@@ -156,7 +156,7 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
 
     // loop over all atoms in other bins in stencil, store every pair
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         if (exclude && exclusion(i,j,type[i],type[j],mask,molecule)) continue;

src/npair_half_size_bin_newton_tri.cpp

@@ -112,7 +112,7 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
     //         (equal zyx and j <= i)
     // latter excludes self-self interaction but allows superposed atoms
 
-    ibin = coord2bin(x[i]);
+    ibin = atom2bin[i];
     for (k = 0; k < nstencil; k++) {
       for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
         if (x[j][2] < ztmp) continue;