diff --git a/src/USER-DPD/nbin_ssa.h b/src/USER-DPD/nbin_ssa.h
index 2a0175081ea02edab1647cdfaf210acc647300e0..12a572b94caab73125eb01887b89849189c11826 100644
--- a/src/USER-DPD/nbin_ssa.h
+++ b/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]) {