more careful checks and print error messages to global root.

src/REPLICA/fix_alchemy.cpp

@@ -35,12 +35,12 @@ using namespace FixConst;
 
 FixAlchemy::FixAlchemy(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), commbuf(nullptr)
 {
-  if (narg != 4) error->all(FLERR, "Incorrect number of arguments for fix alchemy");
+  if (narg != 4) error->universe_all(FLERR, "Incorrect number of arguments for fix alchemy");
-  if (universe->nworlds != 2) error->all(FLERR, "Must use exactly two partitions");
+  if (universe->nworlds != 2) error->universe_all(FLERR, "Must use exactly two partitions");
   if (utils::strmatch(arg[3], "^v_"))
     id_lambda = arg[3] + 2;
   else
-    error->all(FLERR, "Must use variable as lambda argument to fix alchemy");
+    error->universe_all(FLERR, "Must use variable as lambda argument to fix alchemy");
 
   lambda = epot[0] = epot[1] = epot[2] = 0.0;
   progress = 0;
@@ -73,7 +73,9 @@ FixAlchemy::FixAlchemy(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
   int allfail = 0;
   MPI_Allreduce(&fail, &allfail, 1, MPI_INT, MPI_MAX, universe->uworld);
   if (allfail)
-    error->all(FLERR, "Number of atoms and domain decomposition must match for both partitions");
+    error->universe_all(FLERR,
+                        "Number of atoms and domain decomposition must be the same "
+                        "on all partitions");
 
   id_pe = std::string(id) + "_pe";
   pe = modify->add_compute(id_pe + " all pe");
@@ -130,17 +132,19 @@ void FixAlchemy::init()
   memory->create(commbuf, sizeof(double) * nmax, "alchemy:nmax");
 
   if (modify->get_fix_by_style("^balance").size() > 0)
-    error->all(FLERR, "Fix alchemy is not compatible with load balancing");
+    error->universe_all(FLERR, "Fix alchemy is not compatible with load balancing");
 
   if (modify->get_fix_by_style("^alchemy").size() > 1)
-    error->all(FLERR, "There may only one fix alchemy at a time");
+    error->universe_all(FLERR, "There may only one fix alchemy at a time");
+
+  if (utils::strmatch(update->integrate_style, "^respa"))
+    error->universe_all(FLERR, "Must not use run style respa with fix alchemy");
 
   ivar = input->variable->find(id_lambda.c_str());
   if (ivar < 0)
-    error->universe_one(FLERR, fmt::format("Variable {} for fix alchemy does not exist", id_lambda));
+    error->universe_one(FLERR, fmt::format("Fix alchemy variable {} does not exist", id_lambda));
   if (!input->variable->equalstyle(ivar))
-    error->universe_one(FLERR,
+    error->universe_one(FLERR, fmt::format("Fix alchemy variable {} is invalid style", id_lambda));
-                        fmt::format("Variable {} for fix alchemy is invalid style", id_lambda));
   lambda = input->variable->compute_equal(ivar);
 
   // synchronize box dimensions, determine if resync during run will be needed.
@@ -156,33 +160,60 @@ void FixAlchemy::init()
 
 void FixAlchemy::setup(int vflag)
 {
-  if (utils::strmatch(update->integrate_style, "^respa")) {
-    auto respa = dynamic_cast<Respa *>(update->integrate);
-    respa->copy_flevel_f(ilevel_respa);
-    post_force_respa(vflag, ilevel_respa, 0);
-    respa->copy_f_flevel(ilevel_respa);
-  } else {
-    post_force(vflag);
-  }
-
   if (universe->me == 0) {
-    double delta = update->ntimestep - update->beginstep;
+    progress = 0;
-    if ((delta != 0.0) && (update->beginstep != update->endstep))
+    auto msg = fmt::format("Starting alchemical run\n");
-          delta /= update->endstep - update->beginstep;
-    progress = static_cast<int>(delta*100.0);
-    auto msg = fmt::format("Starting alchemical transformation at {:>3d}%\n", progress);
     if (universe->uscreen) fmt::print(universe->uscreen, msg);
     if (universe->ulogfile) fmt::print(universe->ulogfile, msg);
   }
+
+  // recheck domain decomposition, atom ordering, and synchronize positions
+
+  post_integrate();
+
+  // mix initial forces
+
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixAlchemy::post_integrate()
 {
+  // re-check that we have the same domain decomposition on all ranks
+  const int nlocal = atom->nlocal;
+  int my_nlocal[2] = {0, 0};
+  int all_nlocal[2] = {0, 0};
+  my_nlocal[universe->iworld] = nlocal;
+  MPI_Allreduce(my_nlocal, all_nlocal, 2, MPI_INT, MPI_SUM, samerank);
+  int fail = (all_nlocal[0] == all_nlocal[1]) ? 0 : 1;
+  int allfail = 0;
+  MPI_Allreduce(&fail, &allfail, 1, MPI_INT, MPI_MAX, universe->uworld);
+  if (allfail)
+    error->universe_all(FLERR,
+                        "Number of atoms and domain decomposition must be the same on "
+                        "all partitions");
+
+  // check that we have the same atom order on all ranks
+  // re-use communication buffer for positions and forces
+
+  tagint *tagbuf = (tagint *) commbuf;
+  tagint *tag = atom->tag;
+  if (universe->iworld == 0) {
+    for (int i = 0; i < nlocal; ++i) tagbuf[i] = tag[i];
+  }
+  MPI_Bcast(tagbuf, nlocal, MPI_LMP_TAGINT, 0, samerank);
+  fail = allfail = 0;
+  if (universe->iworld > 0) {
+    for (int i = 0; i < nlocal; ++i)
+      if (tag[i] != tagbuf[i]) fail = 1;
+  }
+  MPI_Allreduce(&fail, &allfail, 1, MPI_INT, MPI_MAX, universe->uworld);
+  if (allfail) error->universe_all(FLERR, "Atoms must have the same order on all partitions");
+
   // synchronize atom positions
 
-  const int nall = atom->nlocal + atom->nghost;
+  const int nall = atom->nlocal;
   MPI_Bcast(&atom->x[0][0], 3 * nall, MPI_DOUBLE, 0, samerank);
 
   // synchronize box dimensions, if needed
@@ -228,10 +259,10 @@ void FixAlchemy::post_force(int /*vflag*/)
     double delta = update->ntimestep - update->beginstep;
     if ((delta != 0.0) && (update->beginstep != update->endstep))
       delta /= update->endstep - update->beginstep;
-    int status = static_cast<int>(delta*100.0);
+    int status = static_cast<int>(delta * 100.0);
     if ((status / 10) > (progress / 10)) {
       progress = status;
-      auto msg = fmt::format("  Alchemical transformation progress: {:>3d}%\n", progress);
+      auto msg = fmt::format("  Alchemical run progress: {:>3d}%\n", progress);
       if (universe->uscreen) fmt::print(universe->uscreen, msg);
       if (universe->ulogfile) fmt::print(universe->ulogfile, msg);
     }

src/variable.cpp

@@ -239,7 +239,7 @@ void Variable::set(int narg, char **arg)
     style[nvar] = WORLD;
     num[nvar] = narg - 2;
     if (num[nvar] != universe->nworlds)
-      error->all(FLERR,"World variable count doesn't match # of partitions");
+      error->universe_all(FLERR,"World variable count doesn't match # of partitions");
     which[nvar] = universe->iworld;
     pad[nvar] = 0;
     data[nvar] = new char*[num[nvar]];