Merge pull request #2151 from akohlmey/move-ubuf-to-lmptype

Move definition(s) of union ubuf to lmptype.h

doc/src/reset_ids.rst

@@ -8,7 +8,14 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   reset_ids
+   reset_ids keyword values ...
+
+   * zero or more keyword/value pairs may be appended
+   * keyword = *sort*
+
+.. parsed-literal::
+
+   *sort* value = *yes* or *no*
 
 Examples
 """"""""
@@ -16,6 +23,7 @@ Examples
 .. code-block:: LAMMPS
 
    reset_ids
+   reset_ids sort yes
 
 Description
 """""""""""
@@ -33,11 +41,32 @@ e.g. due to atoms moving outside a simulation box with fixed
 boundaries (see the "boundary command"), or due to evaporation (see
 the "fix evaporate" command).
 
-Note that the resetting of IDs is not really a compression, where gaps
+If the *sort* keyword is used with a setting of *yes*, then the
-in atom IDs are removed by decrementing atom IDs that are larger.
+assignment of new atom IDs will be the same no matter how many
-Instead the IDs for all atoms are erased, and new IDs are assigned so
+processors LAMMPS is running on.  This is done by first doing a
-that the atoms owned by an individual processor have consecutive IDs,
+spatial sort of all the atoms into bins and sorting them within each
-as the :doc:`create_atoms <create_atoms>` command explains.
+bin.  Because the set of bins is independent of the number of
+processors, this enables a consistent assignment of new IDs to each
+atom.
+
+This can be useful to do after using the "create_atoms" command and/or
+"replicate" command.  In general those commands do not guarantee
+assignment of the same atom ID to the same physical atom when LAMMPS
+is run on different numbers of processors.  Enforcing consistent IDs
+can be useful for debugging or comparing output from two different
+runs.
+
+Note that the spatial sort requires communication of atom IDs and
+coordinates between processors in an all-to-all manner.  This is done
+efficiently in LAMMPS, but it is more expensive than how atom IDs are
+reset without sorting.
+
+Note that whether sorting or not, the resetting of IDs is not a
+compression, where gaps in atom IDs are removed by decrementing atom
+IDs that are larger.  Instead the IDs for all atoms are erased, and
+new IDs are assigned so that the atoms owned by an individual
+processor have consecutive IDs, as the :doc:`create_atoms
+<create_atoms>` command explains.
 
 .. note::
 
@@ -59,4 +88,7 @@ Related commands
 
 :doc:`delete_atoms <delete_atoms>`
 
-**Default:** none
+Default
+"""""""
+
+By default, *sort* is no.

src/MANYBODY/pair_adp.cpp

@@ -554,7 +554,7 @@ void PairADP::read_file(char *filename)
       ValueTokenizer values = reader.next_values(1);
       file->nelements = values.next_int();
 
-      if (values.count() != file->nelements + 1)
+      if ((int)values.count() != file->nelements + 1)
         error->one(FLERR,"Incorrect element names in ADP potential file");
 
       file->elements = new char*[file->nelements];

src/MANYBODY/pair_airebo.cpp

@@ -3488,7 +3488,7 @@ void PairAIREBO::read_file(char *filename)
       // global parameters
       current_section = "global parameters";
 
-      for(int i = 0; i < params.size(); i++) {
+      for(int i = 0; i < (int)params.size(); i++) {
         *params[i] = reader.next_double();
       }
 

src/MANYBODY/pair_comb3.cpp

@@ -312,7 +312,7 @@ double PairComb3::init_one(int i, int j)
 
 void PairComb3::read_lib()
 {
-  int i,j,k,l,m;
+  int i,j,k,l;
   int ii,jj,kk,ll,mm,iii;
 
   // open library file on proc 0

src/MANYBODY/pair_eam.cpp

@@ -469,8 +469,6 @@ void PairEAM::read_file(char *filename)
     PotentialFileReader reader(lmp, filename, "EAM");
 
     try {
-      char * line = nullptr;
-
       reader.skip_line();
 
       ValueTokenizer values = reader.next_values(2);

src/MANYBODY/pair_eam_alloy.cpp

@@ -130,7 +130,7 @@ void PairEAMAlloy::read_file(char *filename)
       ValueTokenizer values = reader.next_values(1);
       file->nelements = values.next_int();
 
-      if (values.count() != file->nelements + 1)
+      if ((int)values.count() != file->nelements + 1)
         error->one(FLERR,"Incorrect element names in EAM potential file");
 
       file->elements = new char*[file->nelements];

src/MANYBODY/pair_eam_cd.cpp

@@ -518,7 +518,7 @@ void PairEAMCD::read_h_coeff(char *filename)
     int degree = values.next_int();
     nhcoeff = degree+1;
 
-    if (values.count() != nhcoeff + 1 || nhcoeff < 1)
+    if ((int)values.count() != nhcoeff + 1 || nhcoeff < 1)
       error->one(FLERR, "Failed to read h(x) function coefficients in EAM file.");
 
     hcoeff = new double[nhcoeff];

src/MANYBODY/pair_eam_fs.cpp

@@ -122,8 +122,6 @@ void PairEAMFS::read_file(char *filename)
     PotentialFileReader reader(lmp, filename, "EAM");
 
     try {
-      char * line = nullptr;
-
       reader.skip_line();
       reader.skip_line();
       reader.skip_line();
@@ -132,7 +130,7 @@ void PairEAMFS::read_file(char *filename)
       ValueTokenizer values = reader.next_values(1);
       file->nelements = values.next_int();
 
-      if (values.count() != file->nelements + 1)
+      if ((int)values.count() != file->nelements + 1)
         error->one(FLERR,"Incorrect element names in EAM potential file");
 
       file->elements = new char*[file->nelements];

src/RIGID/fix_rigid_small.cpp

@@ -72,6 +72,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   create_attribute = 1;
   dof_flag = 1;
   enforce2d_flag = 1;
+  stores_ids = 1;
 
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);

src/RIGID/fix_shake.cpp

@@ -64,6 +64,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   thermo_virial = 1;
   create_attribute = 1;
   dof_flag = 1;
+  stores_ids = 1;
 
   // error check
 

src/USER-OMP/pair_eam_alloy_omp.cpp

@@ -130,7 +130,7 @@ void PairEAMAlloyOMP::read_file(char *filename)
       ValueTokenizer values = reader.next_values(1);
       file->nelements = values.next_int();
 
-      if (values.count() != file->nelements + 1)
+      if ((int)values.count() != file->nelements + 1)
         error->one(FLERR,"Incorrect element names in EAM potential file");
 
       file->elements = new char*[file->nelements];

src/USER-OMP/pair_eam_fs_omp.cpp

@@ -130,7 +130,7 @@ void PairEAMFSOMP::read_file(char *filename)
       ValueTokenizer values = reader.next_values(1);
       file->nelements = values.next_int();
 
-      if (values.count() != file->nelements + 1)
+      if ((int)values.count() != file->nelements + 1)
         error->one(FLERR,"Incorrect element names in EAM potential file");
 
       file->elements = new char*[file->nelements];

src/atom_vec.h

@@ -209,25 +209,6 @@ class AtomVec : protected Pointers {
 
   bool *threads;
 
-  // union data struct for packing 32-bit and 64-bit ints into double bufs
-  // this avoids aliasing issues by having 2 pointers (double,int)
-  //   to same buf memory
-  // constructor for 32-bit int prevents compiler
-  //   from possibly calling the double constructor when passed an int
-  // copy to a double *buf:
-  //   buf[m++] = ubuf(foo).d, where foo is a 32-bit or 64-bit int
-  // copy from a double *buf:
-  //   foo = (int) ubuf(buf[m++]).i;, where (int) or (tagint) match foo
-  //   the cast prevents compiler warnings about possible truncation
-
-  union ubuf {
-    double d;
-    int64_t i;
-    ubuf(double arg) : d(arg) {}
-    ubuf(int64_t arg) : i(arg) {}
-    ubuf(int arg) : i(arg) {}
-  };
-
   // local methods
 
   void grow_nmax();

src/comm.cpp

@@ -256,7 +256,6 @@ void Comm::init_exchange()
   int nfix = modify->nfix;
   Fix **fix = modify->fix;
 
-  int onefix;
   maxexchange_fix = 0;
   for (int i = 0; i < nfix; i++)
     maxexchange_fix += fix[i]->maxexchange;

src/compute.h

@@ -161,17 +161,6 @@ class Compute : protected Pointers {
     return j >> SBBITS & 3;
   }
 
-  // union data struct for packing 32-bit and 64-bit ints into double bufs
-  // see atom_vec.h for documentation
-
-  union ubuf {
-    double d;
-    int64_t i;
-    ubuf(double arg) : d(arg) {}
-    ubuf(int64_t arg) : i(arg) {}
-    ubuf(int arg) : i(arg) {}
-  };
-
   // private methods
 
   void adjust_dof_fix();

src/fix.cpp

@@ -80,6 +80,7 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) :
   maxexchange = 0;
   maxexchange_dynamic = 0;
   pre_exchange_migrate = 0;
+  stores_ids = 0;
 
   scalar_flag = vector_flag = array_flag = 0;
   peratom_flag = local_flag = 0;

src/fix.h

@@ -64,6 +64,7 @@ class Fix : protected Pointers {
   int maxexchange;               // max # of per-atom values for Comm::exchange()
   int maxexchange_dynamic;       // 1 if fix sets maxexchange dynamically
   int pre_exchange_migrate;      // 1 if fix migrates atoms in pre_exchange()
+  int stores_ids;                 // 1 if fix stores atom IDs
 
   int scalar_flag;               // 0/1 if compute_scalar() function exists
   int vector_flag;               // 0/1 if compute_vector() function exists
@@ -241,17 +242,6 @@ class Fix : protected Pointers {
   void v_tally(int, int *, double, double *);
   void v_tally(int, double *);
   void v_tally(int, int, double);
-
-  // union data struct for packing 32-bit and 64-bit ints into double bufs
-  // see atom_vec.h for documentation
-
-  union ubuf {
-    double d;
-    int64_t i;
-    ubuf(double arg) : d(arg) {}
-    ubuf(int64_t arg) : i(arg) {}
-    ubuf(int arg) : i(arg) {}
-  };
 };
 
 namespace FixConst {

src/library.cpp

@@ -1667,7 +1667,7 @@ int lammps_style_name(void* ptr, char * category, int index, char * buffer, int
   Info info(lmp);
   auto styles = info.get_available_styles(category);
 
-  if (index < styles.size()) {
+  if (index < (int)styles.size()) {
     strncpy(buffer, styles[index].c_str(), max_size);
     return true;
   }

src/lmptype.h

@@ -171,6 +171,46 @@ typedef int bigint;
 
 #endif
 
+  /// Data structe for packing 32-bit and 64-bit integers
+  /// into double (communication) buffers
+  ///
+  /// Using this union avoids aliasing issues by having member types
+  /// (double, int) referencing the same buffer memory location.
+  ///
+  /// The explicit constructor for 32-bit integers prevents compilers
+  /// from (incorrectly) calling the double constructor when storing
+  ///  an int into a double buffer.
+  /*
+\verbatim embed:rst
+**Usage:**
+
+To copy an integer into a double buffer:
+
+.. code-block: c++
+
+   double buf[2];
+   int    foo =   1;
+   tagint bar = 2<<40;
+   buf[1] = ubuf(foo).d;
+   buf[2] = ubuf(bar).d;
+
+To copy from a double buffer back to an int:
+
+.. code-block: c++
+
+   foo = (int)    ubuf(buf[1]).i;
+   bar = (tagint) ubuf(buf[2]).i;
+
+The typecast prevents compiler warnings about possible truncations.
+\endverbatim
+  */
+  union ubuf {
+    double  d;
+    int64_t i;
+    ubuf(const double  &arg) : d(arg) {}
+    ubuf(const int64_t &arg) : i(arg) {}
+    ubuf(const int     &arg) : i(arg) {}
+  };
 }
 
 // preprocessor macros for compiler specific settings

src/pair.h

@@ -258,17 +258,6 @@ class Pair : protected Pointers {
                       double, double, double, double, double, double);
   void virial_fdotr_compute();
 
-  // union data struct for packing 32-bit and 64-bit ints into double bufs
-  // see atom_vec.h for documentation
-
-  union ubuf {
-    double d;
-    int64_t i;
-    ubuf(double arg) : d(arg) {}
-    ubuf(int64_t arg) : i(arg) {}
-    ubuf(int arg) : i(arg) {}
-  };
-
   inline int sbmask(int j) const {
     return j >> SBBITS & 3;
   }

src/reset_ids.cpp

@@ -13,10 +13,13 @@
 
 #include "reset_ids.h"
 #include <mpi.h>
+#include <cstring>
 #include "atom.h"
 #include "atom_vec.h"
 #include "domain.h"
 #include "comm.h"
+#include "modify.h"
+#include "fix.h"
 #include "special.h"
 #include "memory.h"
 #include "error.h"
@@ -24,25 +27,54 @@
 
 using namespace LAMMPS_NS;
 
+#if defined(LMP_QSORT)
+// allocate space for static class variable
+// prototype for non-class function
+ResetIDs::AtomRvous *ResetIDs::sortrvous;
+static int compare_coords(const void *, const void *);
+#else
+#include "mergesort.h"
+// prototype for non-class function
+static int compare_coords(const int, const int, void *);
+#endif
+
+#define PERBIN 10
+#define BIG 1.0e20
+
 /* ---------------------------------------------------------------------- */
 
 ResetIDs::ResetIDs(LAMMPS *lmp) : Pointers(lmp) {}
 
 /* ---------------------------------------------------------------------- */
 
-void ResetIDs::command(int narg, char ** /* arg */)
+void ResetIDs::command(int narg, char **arg)
 {
   if (domain->box_exist == 0)
     error->all(FLERR,"Reset_ids command before simulation box is defined");
-  if (narg != 0) error->all(FLERR,"Illegal reset_ids command");
   if (atom->tag_enable == 0)
     error->all(FLERR,"Cannot use reset_ids unless atoms have IDs");
 
-  // NOTE: check if any fixes exist which store atom IDs?
+  for (int i = 0; i < modify->nfix; i++)
-  // if so, this operation will mess up the fix
+    if (modify->fix[i]->stores_ids)
+      error->all(FLERR,"Cannot use reset_ids when a fix exists that stores atom IDs");
 
   if (comm->me == 0) utils::logmesg(lmp,"Resetting atom IDs ...\n");
 
+  // process args
+
+  int sortflag = 0;
+
+  int iarg = 0;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"sort") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal reset_ids command");
+      if (strcmp(arg[iarg+1],"yes") == 0) sortflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) sortflag = 0;
+      else error->all(FLERR,"Illegal reset_ids command");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal reset_ids command");
+  }
+
   // create an atom map if one doesn't exist already
 
   int mapflag = 0;
@@ -84,9 +116,12 @@ void ResetIDs::command(int narg, char ** /* arg */)
     tag[i] = 0;
   }
 
-  // assign new contiguous IDs to owned atoms via tag_extend()
+  // assign new contiguous IDs to owned atoms
+  // if sortflag = no: use fast tag_extend()
+  // if sortflag = yes: use slower full spatial sort plus rendezvous comm
 
-  atom->tag_extend();
+  if (sortflag == 0) atom->tag_extend();
+  else sort();
 
   // newIDs = copy of new IDs
   // restore old IDs, consistent with existing atom map
@@ -96,7 +131,7 @@ void ResetIDs::command(int narg, char ** /* arg */)
   memory->create(newIDs,nall,1,"reset_ids:newIDs");
 
   for (int i = 0; i < nlocal; i++) {
-    newIDs[i][0] = tag[i];
+    newIDs[i][0] = ubuf(tag[i]).d;
     tag[i] = oldIDs[i];
   }
 
@@ -228,7 +263,7 @@ void ResetIDs::command(int narg, char ** /* arg */)
 
   atom->map_clear();
   atom->nghost = 0;
-  for (int i = 0; i < nlocal; i++) tag[i] = static_cast<tagint> (newIDs[i][0]);
+  for (int i = 0; i < nlocal; i++) tag[i] = (tagint) ubuf(newIDs[i][0]).i;
   atom->map_init();
   atom->map_set();
 
@@ -251,3 +286,317 @@ void ResetIDs::command(int narg, char ** /* arg */)
   memory->destroy(oldIDs);
   memory->destroy(newIDs);
 }
+
+/* ----------------------------------------------------------------------
+   spatial sort of atoms followed by rendezvous comm to reset atom IDs
+------------------------------------------------------------------------- */
+
+void ResetIDs::sort()
+{
+  double mylo[3],myhi[3],bboxlo[3],bboxhi[3];
+
+  int me = comm->me;
+  int nprocs = comm->nprocs;
+  int dim = domain->dimension;
+
+  // bboxlo,bboxhi = bounding box on all atoms in system
+  // expanded by 0.01 percent
+  // bbox should work for orthogonal or triclinic system
+
+  double **x = atom->x;
+  int nlocal = atom->nlocal;
+
+  mylo[0] = mylo[1] = mylo[2] = BIG;
+  myhi[0] = myhi[1] = myhi[2] = -BIG;
+
+  for (int i = 0; i < nlocal; i++) {
+    mylo[0] = MIN(mylo[0],x[i][0]);
+    mylo[1] = MIN(mylo[1],x[i][1]);
+    mylo[2] = MIN(mylo[2],x[i][2]);
+    myhi[0] = MAX(myhi[0],x[i][0]);
+    myhi[1] = MAX(myhi[1],x[i][1]);
+    myhi[2] = MAX(myhi[2],x[i][2]);
+  }
+
+  if (dim == 2) mylo[2] = myhi[2] = 0.0;
+
+  MPI_Allreduce(mylo,bboxlo,3,MPI_DOUBLE,MPI_MIN,world);
+  MPI_Allreduce(myhi,bboxhi,3,MPI_DOUBLE,MPI_MAX,world);
+
+  bboxlo[0] -= 0.0001 * (bboxhi[0]-bboxlo[0]);
+  bboxlo[1] -= 0.0001 * (bboxhi[1]-bboxlo[1]);
+  bboxlo[2] -= 0.0001 * (bboxhi[2]-bboxlo[2]);
+  bboxhi[0] += 0.0001 * (bboxhi[0]-bboxlo[0]);
+  bboxhi[1] += 0.0001 * (bboxhi[1]-bboxlo[1]);
+  bboxhi[2] += 0.0001 * (bboxhi[2]-bboxlo[2]);
+
+  // nbin_estimate = estimate of total number of bins, each with PERBIN atoms
+  // binsize = edge length of a cubic bin
+  // nbin xyz = bin count in each dimension
+
+  bigint nbin_estimate = atom->natoms / PERBIN;
+  double vol;
+  if (dim == 2) vol = (bboxhi[0]-bboxlo[0]) * (bboxhi[1]-bboxlo[1]);
+  else vol = (bboxhi[0]-bboxlo[0]) * (bboxhi[1]-bboxlo[1]) * (bboxhi[2]-bboxlo[2]);
+  double binsize = pow(vol/nbin_estimate,1.0/dim);
+
+  int nbinx = static_cast<int> ((bboxhi[0]-bboxlo[0]) / binsize) + 1;
+  int nbiny = static_cast<int> ((bboxhi[1]-bboxlo[1]) / binsize) + 1;
+  int nbinz = static_cast<int> ((bboxhi[2]-bboxlo[2]) / binsize) + 1;
+
+  double invx = 1.0 / (bboxhi[0]-bboxlo[0]);
+  double invy = 1.0 / (bboxhi[1]-bboxlo[1]);
+  double invz;
+  if (dim == 2) invz = 0.0;
+  else invz = 1.0 / (bboxhi[2]-bboxlo[2]);
+
+  // nbins = actual total number of bins
+  // bins are split evenly across procs
+  // lo-numbered procs may have own one bin less than rest of procs
+  // nlo = # of bins/proc for lowest procs
+  // nhi = nlo+1 = # of bins/proc for highest procs
+  // nplo = # of procs in low group
+  // nbinlo = # of bins owned by low group procs
+  // binlo to binhi-1 = bin indices this proc will own in Rvous decomp
+  // bins are numbered from 0 to Nbins-1
+
+  bigint nbins = (bigint) nbinx*nbiny*nbinz;
+  int nlo = nbins / nprocs;
+  int nhi = nlo + 1;
+  int nplo = nprocs - (nbins % nprocs);
+  bigint nbinlo = nplo*nlo;
+
+  if (me < nplo) {
+    binlo = me * nlo;
+    binhi = (me+1) * nlo;
+  } else {
+    binlo = nbinlo + (me-nplo) * nhi;
+    binhi = nbinlo + (me+1-nplo) * nhi;
+  }
+
+  // fill atombuf with info on my atoms
+  // ibin = which bin the atom is in
+  // proclist = proc that owns ibin
+
+  int *proclist;
+  memory->create(proclist,nlocal,"special:proclist");
+  AtomRvous *atombuf = (AtomRvous *)
+    memory->smalloc((bigint) nlocal*sizeof(AtomRvous),"resetIDs:idbuf");
+
+  int ibinx,ibiny,ibinz,iproc;
+  bigint ibin;
+
+  for (int i = 0; i < nlocal; i++) {
+    ibinx = static_cast<int> ((x[i][0]-bboxlo[0])*invx * nbinx);
+    ibiny = static_cast<int> ((x[i][1]-bboxlo[1])*invy * nbiny);
+    ibinz = static_cast<int> ((x[i][2]-bboxlo[2])*invz * nbinz);
+    ibin = (bigint) ibinz*nbiny*nbinx + (bigint) ibiny*nbinx + ibinx;
+
+    if (ibin < nbinlo) iproc = ibin / nlo;
+    else iproc = nplo + (ibin-nbinlo) / nhi;
+    proclist[i] = iproc;
+
+    atombuf[i].ibin = ibin;
+    atombuf[i].proc = me;
+    atombuf[i].ilocal = i;
+    atombuf[i].x[0] = x[i][0];
+    atombuf[i].x[1] = x[i][1];
+    atombuf[i].x[2] = x[i][2];
+  }
+
+  // perform rendezvous operation, send atombuf to other procs
+
+  char *buf;
+  int nreturn = comm->rendezvous(1,nlocal,(char *) atombuf,sizeof(AtomRvous),
+				 0,proclist,
+				 sort_bins,0,buf,sizeof(IDRvous),(void *) this);
+  IDRvous *outbuf = (IDRvous *) buf;
+
+  memory->destroy(proclist);
+  memory->sfree(atombuf);
+
+  // set new ID for all owned atoms
+
+  int ilocal;
+
+  for (int i = 0; i < nreturn; i++) {
+    ilocal = outbuf[i].ilocal;
+    atom->tag[ilocal] = outbuf[i].newID;
+  }
+
+  memory->sfree(outbuf);
+}
+
+/* ----------------------------------------------------------------------
+   process data for atoms assigned to me in rendezvous decomposition
+   inbuf = list of N AtomRvous datums
+   outbuf = list of N IDRvous datums, sent back to sending proc
+------------------------------------------------------------------------- */
+
+int ResetIDs::sort_bins(int n, char *inbuf,
+			int &flag, int *&proclist, char *&outbuf,
+			void *ptr)
+{
+  int i,ibin,index;
+
+  ResetIDs *rptr = (ResetIDs *) ptr;
+  Memory *memory = rptr->memory;
+  Error *error = rptr->error;
+  MPI_Comm world = rptr->world;
+  bigint binlo = rptr->binlo;
+  bigint binhi = rptr->binhi;
+
+  // nbins = my subset of bins from binlo to binhi-1
+  // initialize linked lists of my Rvous atoms in each bin
+
+  int *next,*head,*last,*count;
+
+  int nbins = binhi - binlo;
+  memory->create(next,n,"resetIDs:next");
+  memory->create(head,nbins,"resetIDs:head");
+  memory->create(last,nbins,"resetIDs:last");
+  memory->create(count,nbins,"resetIDs:count");
+
+  for (i = 0; i < n; i++) next[i] = -1;
+
+  for (ibin = 0; ibin < nbins; ibin++) {
+    head[ibin] = last[ibin] = -1;
+    count[ibin] = 0;
+  }
+
+  AtomRvous *in = (AtomRvous *) inbuf;
+
+  for (i = 0; i < n; i++) {
+    if (in[i].ibin < binlo || in[i].ibin >= binhi) {
+      //printf("BAD me %d i %d %d ibin %d binlohi %d %d\n",
+      //	     rptr->comm->me,i,n,in[i].ibin,binlo,binhi);
+      error->one(FLERR,"Bad spatial bin assignment in reset_ids sort");
+    }
+    ibin = in[i].ibin - binlo;
+    if (head[ibin] < 0) head[ibin] = i;
+    if (last[ibin] >= 0) next[last[ibin]] = i;
+    last[ibin] = i;
+    count[ibin]++;
+  }
+
+  int maxcount = 0;
+  for (ibin = 0; ibin < nbins; ibin++)
+    maxcount = MAX(maxcount,count[ibin]);
+
+  int *order;
+  memory->create(order,maxcount,"resetIDs:order");
+
+  // sort atoms in each bin by their spatial position
+  // recreate linked list for each bin based on sorted order
+
+  for (ibin = 0; ibin < nbins; ibin++) {
+    index = head[ibin];
+    for (i = 0; i < count[ibin]; i++) {
+      order[i] = index;
+      index = next[index];
+    }
+
+
+#if defined(LMP_QSORT)
+    sortrvous = in;
+    qsort(order,count[ibin],sizeof(int),compare_coords);
+#else
+    merge_sort(order,count[ibin],(void *) in,compare_coords);
+#endif
+
+    head[ibin] = last[ibin] = -1;
+    for (i = 0; i < count[ibin]; i++) {
+      if (head[ibin] < 0) head[ibin] = order[i];
+      if (last[ibin] >= 0) next[last[ibin]] = order[i];
+      last[ibin] = order[i];
+    }
+  }
+
+  // MPI_Scan() to find out where my atoms begin in globally sorted list
+
+  tagint ntag = n;
+  tagint nprev;
+  MPI_Scan(&ntag,&nprev,1,MPI_LMP_TAGINT,MPI_SUM,world);
+  nprev -= n;
+
+  // loop over bins and sorted atoms in each bin, reset ids to be consecutive
+  // pass outbuf of IDRvous datums back to comm->rendezvous
+
+  int nout = n;
+  memory->create(proclist,nout,"resetIDs:proclist");
+  IDRvous *out = (IDRvous *)
+    memory->smalloc(nout*sizeof(IDRvous),"resetIDs:out");
+
+  tagint one = nprev + 1;
+  for (ibin = 0; ibin < nbins; ibin++) {
+    index = head[ibin];
+    for (i = 0; i < count[ibin]; i++) {
+      proclist[index] = in[index].proc;
+      out[index].newID = one++;
+      out[index].ilocal = in[index].ilocal;
+      index = next[index];
+    }
+  }
+
+  outbuf = (char *) out;
+
+  // clean up
+  // Comm::rendezvous will delete proclist and out (outbuf)
+
+  memory->destroy(next);
+  memory->destroy(head);
+  memory->destroy(last);
+  memory->destroy(count);
+  memory->destroy(order);
+
+  // flag = 2: new outbuf
+
+  flag = 2;
+  return nout;
+}
+
+#if defined(LMP_QSORT)
+
+/* ----------------------------------------------------------------------
+   comparison function invoked by qsort()
+   accesses static class member sortrvous, set before call to qsort()
+------------------------------------------------------------------------- */
+
+int compare_coords(const void *iptr, const void *jptr)
+{
+  int i = *((int *) iptr);
+  int j = *((int *) jptr);
+  ResetIDs::AtomRvous *rvous = ResetIDs::sortrvous;
+  double *xi = rvous[i].x;
+  double *xj = rvous[j].x;
+  if (xi[0] < xj[0]) return -1;
+  if (xi[0] > xj[0]) return 1;
+  if (xi[1] < xj[1]) return -1;
+  if (xi[1] > xj[1]) return 1;
+  if (xi[2] < xj[2]) return -1;
+  if (xi[2] > xj[2]) return 1;
+  return 0;
+}
+
+#else
+
+/* ----------------------------------------------------------------------
+   comparison function invoked by merge_sort()
+   void pointer contains sortrvous
+------------------------------------------------------------------------- */
+
+int compare_coords(const int i, const int j, void *ptr)
+{
+  ResetIDs::AtomRvous *rvous = (ResetIDs::AtomRvous *) ptr;
+  double *xi = rvous[i].x;
+  double *xj = rvous[j].x;
+  if (xi[0] < xj[0]) return -1;
+  if (xi[0] > xj[0]) return 1;
+  if (xi[1] < xj[1]) return -1;
+  if (xi[1] > xj[1]) return 1;
+  if (xi[2] < xj[2]) return -1;
+  if (xi[2] > xj[2]) return 1;
+  return 0;
+}
+
+#endif

src/reset_ids.h

@@ -26,10 +26,33 @@ namespace LAMMPS_NS {
 
 class ResetIDs : protected Pointers {
  public:
+  struct AtomRvous {
+    bigint ibin;
+    int proc,ilocal;
+    double x[3];
+  };
+
+  struct IDRvous {
+    tagint newID;
+    int ilocal;
+  };
+
+  #if defined(LMP_QSORT)
+  // static variable across all ResetID objects, for qsort callback
+  static AtomRvous *sortrvous;
+#endif
+
   ResetIDs(class LAMMPS *);
   void command(int, char **);
 
  private:
+  bigint binlo,binhi;
+
+  // callback functions for rendezvous communication
+
+  static int sort_bins(int, char *, int &, int *&, char *&, void *);
+
+  void sort();
 };
 
 }

src/table_file_reader.cpp

@@ -38,7 +38,7 @@ TableFileReader::~TableFileReader() {
 
 char * TableFileReader::find_section_start(const std::string & keyword) {
   char * line = nullptr;
-  while (line = reader->next_line()) {
+  while ((line = reader->next_line())) {
     ValueTokenizer values(line);
     std::string word = values.next_string();
 

src/utils.cpp

@@ -547,7 +547,7 @@ std::string utils::get_potential_date(const std::string & path, const std::strin
   reader.ignore_comments = false;
   char * line = nullptr;
 
-  while (line = reader.next_line()) {
+  while ((line = reader.next_line())) {
     ValueTokenizer values(line);
     while (values.has_next()) {
       std::string word = values.next_string();

unittest/utils/test_fmtlib.cpp

@@ -42,45 +42,63 @@ TEST(FmtLib, insert_neg_int) {
 }
 
 TEST(FmtLib, insert_bigint) {
-    if (sizeof(bigint) == 4) GTEST_SKIP();
+#if defined(LAMMPS_BIGBIG) || defined(LAMMPS_SMALLBIG)
     const bigint val = 9945234592L;
     auto text = fmt::format("word {}",val);
     ASSERT_THAT(text, Eq("word 9945234592"));
+#else
+    GTEST_SKIP();
+#endif
 }
 
 TEST(FmtLib, insert_neg_bigint) {
-    if (sizeof(bigint) == 4) GTEST_SKIP();
+#if defined(LAMMPS_BIGBIG) || defined(LAMMPS_SMALLBIG)
     const bigint val = -9945234592L;
     auto text = fmt::format("word {}",val);
     ASSERT_THAT(text, Eq("word -9945234592"));
+#else
+    GTEST_SKIP();
+#endif
 }
 
 TEST(FmtLib, insert_tagint) {
-    if (sizeof(tagint) == 4) GTEST_SKIP();
+#if defined(LAMMPS_BIGBIG)
     const tagint val = 9945234592L;
     auto text = fmt::format("word {}",val);
     ASSERT_THAT(text, Eq("word 9945234592"));
+#else
+    GTEST_SKIP();
+#endif
 }
 
 TEST(FmtLib, insert_neg_tagint) {
-    if (sizeof(tagint) == 4) GTEST_SKIP();
+#if defined(LAMMPS_BIGBIG)
     const tagint val = -9945234592L;
     auto text = fmt::format("word {}",val);
     ASSERT_THAT(text, Eq("word -9945234592"));
+#else
+    GTEST_SKIP();
+#endif
 }
 
 TEST(FmtLib, insert_imageint) {
-    if (sizeof(imageint) == 4) GTEST_SKIP();
+#if defined(LAMMPS_BIGBIG)
     const imageint val = 9945234592L;
     auto text = fmt::format("word {}",val);
     ASSERT_THAT(text, Eq("word 9945234592"));
+#else
+    GTEST_SKIP();
+#endif
 }
 
 TEST(FmtLib, insert_neg_imageint) {
-    if (sizeof(imageint) == 4) GTEST_SKIP();
+#if defined(LAMMPS_BIGBIG)
     const imageint val = -9945234592L;
     auto text = fmt::format("word {}",val);
     ASSERT_THAT(text, Eq("word -9945234592"));
+#else
+    GTEST_SKIP();
+#endif
 }
 
 TEST(FmtLib, insert_double) {