From d8f8a3a36ab2fafdb63447fb2ffb841175fcc4b6 Mon Sep 17 00:00:00 2001
From: Vsevak <thevsevak@gmail.com>
Date: Sat, 28 May 2022 00:24:24 +0300
Subject: [PATCH] Handle inconsistent J molecules in tip4p/gpu

---
 lib/gpu/geryon/hip_device.h            |  2 +-
 lib/gpu/lal_lj_tip4p_long.cu           | 42 ++++++++++++++++++++------
 src/GPU/pair_lj_cut_tip4p_long_gpu.cpp |  7 +++--
 3 files changed, 38 insertions(+), 13 deletions(-)

diff --git a/lib/gpu/geryon/hip_device.h b/lib/gpu/geryon/hip_device.h
index fadeec8711..f809323ee7 100644
--- a/lib/gpu/geryon/hip_device.h
+++ b/lib/gpu/geryon/hip_device.h
@@ -394,7 +394,7 @@ UCL_Device::~UCL_Device() {
   clear();
 }
 
-int UCL_Device::set_platform(const int) {
+int UCL_Device::set_platform(const int pid) {
   clear();
   #ifdef UCL_DEBUG
   assert(pid<num_platforms());
diff --git a/lib/gpu/lal_lj_tip4p_long.cu b/lib/gpu/lal_lj_tip4p_long.cu
index 8b52f09d68..21fb1a91c8 100644
--- a/lib/gpu/lal_lj_tip4p_long.cu
+++ b/lib/gpu/lal_lj_tip4p_long.cu
@@ -217,26 +217,43 @@ __kernel void k_lj_tip4p_reneigh(const __global numtyp4 *restrict x_,
         // set iH1,iH2 to closest image to O
         iH1 = closest_image(i, iH1, sametag, x_);
         iH2 = closest_image(i, iH2, sametag, x_);
+
+        //printf("%d %f %f %f\n", (iH1 < 0 || iH2 < 0) ? 0:1, ix.x, ix.y, ix.z );
+        /*
+        if (iH1 < 0) {
+                printf("i=%d\ttag[i]=%d\tmap[tag[i]-1]=%d\tmap[tag[i]-2]=%d\n", i, tag[i], map[tag[i]+1], map[tag[i]+2]);
+        }
+        if (iH2 < 0) {
+                printf("i=%d\ttag[i]=%d\tmap[tag[i]-1]=%d\tmap[tag[i]-2]=%d\n", i, tag[i], map[tag[i]+1], map[tag[i]+2]);
+        }*/
+
         hneigh[i*4  ] = iH1;
         hneigh[i*4+1] = iH2;
         hneigh[i*4+2] = -1;
       }
     }
-    if (itype == typeH) {
+    if (i<inum && itype == typeH) {
       if (hneigh[i*4+2] != -1) {
         int iI, iH;
         iI = atom_mapping(map,tag[i] - 1);
-        numtyp4 iIx; fetch4(iIx,iI,pos_tex); //x_[iI];
+        iO  = closest_image(i,iI,sametag, x_);
+        //printf("%d %f %f %f\n", (iI < 0) ? 2:3, ix.x, ix.y, ix.z );
+        /*
+        // printf("iI = %d iO closest = %d\n",iI, iO);
+        if (iI < 0) {
+        printf("i=%d\ttag[i]=%d\tmap[tag[i]-1]=%d\tmap[tag[i]-2]=%d\n", i, tag[i], map[tag[i]-1],map[tag[i]-2]);
+        }*/
+        numtyp4 iIx; fetch4(iIx,iO,pos_tex); //x_[iI];
         if ((int)iIx.w == typeH) {
           iO = atom_mapping(map,tag[i] - 2);
           iO  = closest_image(i, iO, sametag, x_);
-          iH1 = closest_image(i, iI, sametag, x_);
-          iH2 = i;
+          //iH1 = closest_image(i, iI, sametag, x_);
+          //iH2 = i;
         } else { //if ((int)iIx.w == typeO)
-          iH = atom_mapping(map, tag[i] + 1);
-          iO  = closest_image(i,iI,sametag, x_);
-          iH1 = i;
-          iH2 = closest_image(i,iH,sametag, x_);
+          //iH = atom_mapping(map, tag[i] + 1);
+          //iO  = closest_image(i,iI,sametag, x_);
+          //iH1 = i;
+          //iH2 = closest_image(i,iH,sametag, x_);
         }
         hneigh[i*4+0] = iO;
         hneigh[i*4+1] += -1;
@@ -268,7 +285,14 @@ __kernel void k_lj_tip4p_newsite(const __global numtyp4 *restrict x_,
       iH2 = hneigh[i*4+1];
       iO  = i;
       numtyp qO; fetch(qO,iO,q_tex);
-      compute_newsite(iO,iH1,iH2, &m[iO], qO, alpha, x_);
+      if (iH1>=0 && iH2>=0) {
+      	compute_newsite(iO,iH1,iH2, &m[iO], qO, alpha, x_);
+      } else {
+        m[iO] = ix;
+        m[iO].w = qO;
+	hneigh[i*4] = iO;
+	hneigh[i*4+1] = iO;
+      }
     }
   }
 }
diff --git a/src/GPU/pair_lj_cut_tip4p_long_gpu.cpp b/src/GPU/pair_lj_cut_tip4p_long_gpu.cpp
index 5c2c391136..f6e528676d 100644
--- a/src/GPU/pair_lj_cut_tip4p_long_gpu.cpp
+++ b/src/GPU/pair_lj_cut_tip4p_long_gpu.cpp
@@ -186,8 +186,9 @@ void PairLJCutTIP4PLongGPU::init_style()
   alpha = qdist / (cos(0.5 * theta) * blen);
 
   cut_coulsq = cut_coul * cut_coul;
-  double cut_coulsqplus = (cut_coul + qdist + blen) * (cut_coul + qdist + blen);
-  if (maxcut < cut_coulsqplus) { cell_size = (cut_coul + qdist + blen) + neighbor->skin; }
+  double cut_coulplus = cut_coul + qdist + blen;
+  double cut_coulsqplus = cut_coulplus*cut_coulplus;
+  if (sqrt(maxcut) < cut_coulplus+blen) { cell_size = (cut_coulplus + blen) + neighbor->skin; }
   if (comm->cutghostuser < cell_size) {
     if (comm->me == 0)
       error->warning(FLERR,
@@ -205,7 +206,7 @@ void PairLJCutTIP4PLongGPU::init_style()
   GPU_EXTRA::check_flag(success, error, world);
   if (gpu_mode == GPU_FORCE) {
     auto req = neighbor->add_request(this, NeighConst::REQ_FULL);
-    req->set_cutoff(cut_coul + qdist + blen + neighbor->skin);
+    req->set_cutoff(cut_coulplus + neighbor->skin);
   }
 }