From be5aa46df82b1aaa78250f3228465c6f7260c17a Mon Sep 17 00:00:00 2001
From: Trung Nguyen <ndactrung@gmail.com>
Date: Fri, 3 Sep 2021 17:32:41 -0500
Subject: [PATCH] Re-arranged the binsize call from the GPU lib in Atom so that
 the box bounds and bininv[xyz] are computed on the CPU side intact

---
 examples/amoeba/in.ubiquitin |  2 +-
 src/atom.cpp                 | 44 ++++++++++++------------------------
 2 files changed, 16 insertions(+), 30 deletions(-)

diff --git a/examples/amoeba/in.ubiquitin b/examples/amoeba/in.ubiquitin
index f017d8f122..2491493c45 100644
--- a/examples/amoeba/in.ubiquitin
+++ b/examples/amoeba/in.ubiquitin
@@ -4,7 +4,7 @@ units           real
 boundary        p p p
 
 atom_style amoeba
-atom_modify sort 1000 7.0
+#atom_modify sort 1000 7.0
 bond_style      class2
 angle_style     amoeba
 dihedral_style  none
diff --git a/src/atom.cpp b/src/atom.cpp
index 86e2b1151b..71cb2e9f31 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -2195,6 +2195,21 @@ void Atom::setup_sort_bins()
     return;
   }
 
+#ifdef LMP_GPU
+  if (userbinsize == 0.0) {   
+    int ifix = modify->find_fix("package_gpu");
+    if (ifix >= 0) {
+      const double subx = domain->subhi[0] - domain->sublo[0];
+      const double suby = domain->subhi[1] - domain->sublo[1];
+      const double subz = domain->subhi[2] - domain->sublo[2];
+
+      FixGPU *fix = static_cast<FixGPU *>(modify->fix[ifix]);
+      binsize = fix->binsize(subx, suby, subz, atom->nlocal,
+                             0.5 * neighbor->cutneighmax);
+    }
+  }
+#endif
+
   double bininv = 1.0/binsize;
 
   // nbin xyz = local bins
@@ -2271,35 +2286,6 @@ void Atom::setup_sort_bins()
   }
 #endif
 
-#ifdef LMP_GPU
-  if (userbinsize == 0.0) {
-    int ifix = modify->find_fix("package_gpu");
-    if (ifix >= 0) {
-      const double subx = domain->subhi[0] - domain->sublo[0];
-      const double suby = domain->subhi[1] - domain->sublo[1];
-      const double subz = domain->subhi[2] - domain->sublo[2];
-
-      FixGPU *fix = static_cast<FixGPU *>(modify->fix[ifix]);
-      binsize = fix->binsize(subx, suby, subz, atom->nlocal,
-                             neighbor->cutneighmax);
-      bininv = 1.0 / binsize;
-
-      nbinx = static_cast<int> (ceil(subx * bininv));
-      nbiny = static_cast<int> (ceil(suby * bininv));
-      nbinz = static_cast<int> (ceil(subz * bininv));
-      if (domain->dimension == 2) nbinz = 1;
-
-      if (nbinx == 0) nbinx = 1;
-      if (nbiny == 0) nbiny = 1;
-      if (nbinz == 0) nbinz = 1;
-
-      bininvx = bininv;
-      bininvy = bininv;
-      bininvz = bininv;
-    }
-  }
-#endif
-
   if (1.0*nbinx*nbiny*nbinz > INT_MAX)
     error->one(FLERR,"Too many atom sorting bins");