From 72d7ba85adaada362532dcee610dcfac894be7e5 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Thu, 4 Apr 2024 20:34:37 -0400
Subject: [PATCH] call alternate minimum image code from @stanmoore1 from fix
 rigid/small

---
 src/RIGID/fix_rigid_small.cpp |  2 +-
 src/domain.cpp                | 69 +++++++++++++++++++++++++++++++++++
 src/domain.h                  |  2 +
 3 files changed, 72 insertions(+), 1 deletion(-)

diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 9e185a4de2..96bd6019d7 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -2160,7 +2160,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);
   }
diff --git a/src/domain.cpp b/src/domain.cpp
index 7513d384e2..c825016b27 100644
--- a/src/domain.cpp
+++ b/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
diff --git a/src/domain.h b/src/domain.h
index ab054f1b50..ac4dd12e61 100644
--- a/src/domain.h
+++ b/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);