Merge pull request #4123 from akohlmey/rigid-image-fix

Address minimum image issue when restarting fix rigid/small

src/RIGID/fix_rigid_small.cpp

@@ -2174,7 +2174,7 @@ void FixRigidSmall::setup_bodies_static()
     xgc = body[ibody].xgc;
     double delta[3];
     MathExtra::sub3(xgc,xcm,delta);
-    domain->minimum_image(delta);
+    domain->minimum_image_big(delta);
     MathExtra::transpose_matvec(ex,ey,ez,delta,body[ibody].xgc_body);
     MathExtra::add3(xcm,delta,xgc);
   }

src/domain.cpp

@@ -973,6 +973,9 @@ void Domain::subbox_too_small_check(double thresh)
    changed "if" to "while" to enable distance to
      far-away ghost atom returned by atom->map() to be wrapped back into box
      could be problem for looking up atom IDs when cutoff > boxsize
+   should be used for most cases where the difference in the image count
+     is small (usually 0 or 1)
+   use minimum_image_big() when a large difference between image counts is expected
 ------------------------------------------------------------------------- */
 
 static constexpr double MAXIMGCOUNT = 16;
@@ -1045,6 +1048,72 @@ void Domain::minimum_image(double &dx, double &dy, double &dz) const
   }
 }
 
+/* ----------------------------------------------------------------------
+   minimum image convention in periodic dimensions
+   use 1/2 of box size as test
+   for triclinic, also add/subtract tilt factors in other dims as needed
+   allow multiple box lengths to enable distance to
+     far-away ghost atom returned by atom->map() to be wrapped back into box
+     could be problem for looking up atom IDs when cutoff > boxsize
+   this should be used when there is a large image count difference possible
+     this applies for example to fix rigid/small
+------------------------------------------------------------------------- */
+
+void Domain::minimum_image_big(double &dx, double &dy, double &dz) const
+{
+  if (triclinic == 0) {
+    if (xperiodic) {
+      double dfactor = dx/xprd + 0.5;
+      if (dx < 0) dfactor -= 1.0;
+      if (dfactor > MAXSMALLINT)
+        error->one(FLERR, "Atoms have moved too far apart ({}) for minimum image\n", dx);
+      dx -= xprd * static_cast<int>(dfactor);
+    }
+    if (yperiodic) {
+      double dfactor = dy/yprd + 0.5;
+      if (dy < 0) dfactor -= 1.0;
+      if (dfactor > MAXSMALLINT)
+        error->one(FLERR, "Atoms have moved too far apart ({}) for minimum image\n", dy);
+      dy -= yprd * static_cast<int>(dfactor);
+    }
+    if (zperiodic) {
+      double dfactor = dz/zprd + 0.5;
+      if (dz < 0) dfactor -= 1.0;
+      if (dfactor > MAXSMALLINT)
+        error->one(FLERR, "Atoms have moved too far apart ({}) for minimum image\n", dz);
+      dz -= zprd * static_cast<int>(dfactor);
+    }
+
+  } else {
+    if (zperiodic) {
+      double dfactor = dz/zprd + 0.5;
+      if (dz < 0) dfactor -= 1.0;
+      if (dfactor > MAXSMALLINT)
+        error->one(FLERR, "Atoms have moved too far apart ({}) for minimum image\n", dz);
+      int factor = static_cast<int>(dfactor);
+      dz -= zprd * factor;
+      dy -= yz * factor;
+      dx -= xz * factor;
+    }
+    if (yperiodic) {
+      double dfactor = dy/yprd + 0.5;
+      if (dy < 0) dfactor -= 1.0;
+      if (dfactor > MAXSMALLINT)
+        error->one(FLERR, "Atoms have moved too far apart ({}) for minimum image\n", dy);
+      int factor = static_cast<int>(dfactor);
+      dy -= yprd * factor;
+      dx -= xy * factor;
+    }
+    if (xperiodic) {
+      double dfactor = dx/xprd + 0.5;
+      if (dx < 0) dfactor -= 1.0;
+      if (dfactor > MAXSMALLINT)
+        error->one(FLERR, "Atoms have moved too far apart ({}) for minimum image\n", dx);
+      dx -= xprd * static_cast<int>(dfactor);
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    return local index of atom J or any of its images that is closest to atom I
    if J is not a valid index like -1, just return it

src/domain.h

@@ -119,6 +119,8 @@ class Domain : protected Pointers {
   void subbox_too_small_check(double);
   void minimum_image(double &, double &, double &) const;
   void minimum_image(double *delta) const { minimum_image(delta[0], delta[1], delta[2]); }
+  void minimum_image_big(double &, double &, double &) const;
+  void minimum_image_big(double *delta) const { minimum_image_big(delta[0], delta[1], delta[2]); }
   int closest_image(int, int);
   int closest_image(const double *const, int);
   void closest_image(const double *const, const double *const, double *const);