Merge pull request #695 from timattox/USER-DPD_small_updates

USER-DPD small updates

src/KOKKOS/npair_ssa_kokkos.cpp

@@ -70,6 +70,7 @@ void NPairSSAKokkos<DeviceType>::copy_neighbor_info()
   k_ex2_bit = neighborKK->k_ex2_bit;
   k_ex_mol_group = neighborKK->k_ex_mol_group;
   k_ex_mol_bit = neighborKK->k_ex_mol_bit;
+  k_ex_mol_intra = neighborKK->k_ex_mol_intra;
 }
 
 /* ----------------------------------------------------------------------
@@ -217,8 +218,12 @@ int NPairSSAKokkosExecute<DeviceType>::exclusion(const int &i,const int &j,
 
   if (nex_mol) {
     for (m = 0; m < nex_mol; m++)
-      if (mask(i) & ex_mol_bit(m) && mask(j) & ex_mol_bit(m) &&
-          molecule(i) == molecule(j)) return 1;
+      if (ex_mol_intra[m]) { // intra-chain: exclude i-j pair if on same molecule
+        if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] &&
+            molecule[i] == molecule[j]) return 1;
+      } else                 // exclude i-j pair if on different molecules
+        if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] &&
+            molecule[i] != molecule[j]) return 1;
   }
 
   return 0;
@@ -418,6 +423,7 @@ fprintf(stdout, "tota%03d total %3d could use %6d inums, expected %6d inums. inu
          nex_mol,
          k_ex_mol_group.view<DeviceType>(),
          k_ex_mol_bit.view<DeviceType>(),
+         k_ex_mol_intra.view<DeviceType>(),
          bboxhi,bboxlo,
          domain->xperiodic,domain->yperiodic,domain->zperiodic,
          domain->xprd_half,domain->yprd_half,domain->zprd_half);
@@ -432,6 +438,7 @@ fprintf(stdout, "tota%03d total %3d could use %6d inums, expected %6d inums. inu
   k_ex2_bit.sync<DeviceType>();
   k_ex_mol_group.sync<DeviceType>();
   k_ex_mol_bit.sync<DeviceType>();
+  k_ex_mol_intra.sync<DeviceType>();
   k_bincount.sync<DeviceType>();
   k_bins.sync<DeviceType>();
   k_gbincount.sync<DeviceType>();

src/KOKKOS/npair_ssa_kokkos.h

@@ -76,6 +76,7 @@ class NPairSSAKokkos : public NPair {
   DAT::tdual_int_1d k_ex1_bit,k_ex2_bit;
   DAT::tdual_int_1d k_ex_mol_group;
   DAT::tdual_int_1d k_ex_mol_bit;
+  DAT::tdual_int_1d k_ex_mol_intra;
 
   // data from NBinSSA class
 
@@ -123,6 +124,7 @@ class NPairSSAKokkosExecute
   const int nex_mol;
   const typename AT::t_int_1d_const ex_mol_group;
   const typename AT::t_int_1d_const ex_mol_bit;
+  const typename AT::t_int_1d_const ex_mol_intra;
 
   // data from NBinSSA class
 
@@ -233,6 +235,7 @@ class NPairSSAKokkosExecute
         const int & _nex_mol,
         const typename AT::t_int_1d_const & _ex_mol_group,
         const typename AT::t_int_1d_const & _ex_mol_bit,
+        const typename AT::t_int_1d_const & _ex_mol_intra,
         const X_FLOAT *_bboxhi, const X_FLOAT* _bboxlo,
         const int & _xperiodic, const int & _yperiodic, const int & _zperiodic,
         const int & _xprd_half, const int & _yprd_half, const int & _zprd_half):
@@ -266,6 +269,7 @@ class NPairSSAKokkosExecute
     ex1_group(_ex1_group),ex2_group(_ex2_group),
     ex1_bit(_ex1_bit),ex2_bit(_ex2_bit),nex_mol(_nex_mol),
     ex_mol_group(_ex_mol_group),ex_mol_bit(_ex_mol_bit),
+    ex_mol_intra(_ex_mol_intra),
     xperiodic(_xperiodic),yperiodic(_yperiodic),zperiodic(_zperiodic),
     xprd_half(_xprd_half),yprd_half(_yprd_half),zprd_half(_zprd_half) {
 

src/USER-DPD/nbin_ssa.h

@@ -23,6 +23,8 @@ NBinStyle(ssa,
 #define LMP_NBIN_SSA_H
 
 #include "nbin_standard.h"
+#include "math.h"
+#include "error.h"
 
 namespace LAMMPS_NS {
 
@@ -47,79 +49,14 @@ class NBinSSA : public NBinStandard {
 
   bigint memory_usage();
 
-  inline
-  int coord2bin(const double & x,const double & y,const double & z) const
-  {
-    int ix,iy,iz;
-
-    if (x >= bboxhi_[0])
-      ix = static_cast<int> ((x-bboxhi_[0])*bininvx) + nbinx;
-    else if (x >= bboxlo_[0]) {
-      ix = static_cast<int> ((x-bboxlo_[0])*bininvx);
-      ix = MIN(ix,nbinx-1);
-    } else
-      ix = static_cast<int> ((x-bboxlo_[0])*bininvx) - 1;
-
-    if (y >= bboxhi_[1])
-      iy = static_cast<int> ((y-bboxhi_[1])*bininvy) + nbiny;
-    else if (y >= bboxlo_[1]) {
-      iy = static_cast<int> ((y-bboxlo_[1])*bininvy);
-      iy = MIN(iy,nbiny-1);
-    } else
-      iy = static_cast<int> ((y-bboxlo_[1])*bininvy) - 1;
-
-    if (z >= bboxhi_[2])
-      iz = static_cast<int> ((z-bboxhi_[2])*bininvz) + nbinz;
-    else if (z >= bboxlo_[2]) {
-      iz = static_cast<int> ((z-bboxlo_[2])*bininvz);
-      iz = MIN(iz,nbinz-1);
-    } else
-      iz = static_cast<int> ((z-bboxlo_[2])*bininvz) - 1;
-
-    return (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo);
-  }
-
-  inline
-  int coord2bin(const double & x,const double & y,const double & z, int* i) const
-  {
-    int ix,iy,iz;
-
-    if (x >= bboxhi_[0])
-      ix = static_cast<int> ((x-bboxhi_[0])*bininvx) + nbinx;
-    else if (x >= bboxlo_[0]) {
-      ix = static_cast<int> ((x-bboxlo_[0])*bininvx);
-      ix = MIN(ix,nbinx-1);
-    } else
-      ix = static_cast<int> ((x-bboxlo_[0])*bininvx) - 1;
-
-    if (y >= bboxhi_[1])
-      iy = static_cast<int> ((y-bboxhi_[1])*bininvy) + nbiny;
-    else if (y >= bboxlo_[1]) {
-      iy = static_cast<int> ((y-bboxlo_[1])*bininvy);
-      iy = MIN(iy,nbiny-1);
-    } else
-      iy = static_cast<int> ((y-bboxlo_[1])*bininvy) - 1;
-
-    if (z >= bboxhi_[2])
-      iz = static_cast<int> ((z-bboxhi_[2])*bininvz) + nbinz;
-    else if (z >= bboxlo_[2]) {
-      iz = static_cast<int> ((z-bboxlo_[2])*bininvz);
-      iz = MIN(iz,nbinz-1);
-    } else
-      iz = static_cast<int> ((z-bboxlo_[2])*bininvz) - 1;
-
-    i[0] = ix - mbinxlo;
-    i[1] = iy - mbinylo;
-    i[2] = iz - mbinzlo;
-
-    return (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo);
-  }
-
   inline
   int coord2bin(const double & x,const double & y,const double & z, int &ixo, int &iyo, int &izo) const
   {
     int ix,iy,iz;
 
+    if (!ISFINITE(x) || !ISFINITE(y) || !ISFINITE(z))
+      error->one(FLERR,"Non-numeric positions - simulation unstable");
+
     if (x >= bboxhi_[0])
       ix = static_cast<int> ((x-bboxhi_[0])*bininvx) + nbinx;
     else if (x >= bboxlo_[0]) {