From ef2f16a240aed17caff0ec6d5b52a83f50a63e7d Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Mon, 1 Aug 2016 11:52:07 -0400
Subject: [PATCH] more thorough checks on valid numbers when computing
 box/virial related properties

---
 src/KSPACE/ewald.cpp      |  9 +++++++++
 src/KSPACE/ewald_disp.cpp | 12 ++++++++++++
 src/fix_nh.cpp            |  6 ++++++
 3 files changed, 27 insertions(+)

diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp
index 39a98901cb..9d57e22f18 100644
--- a/src/KSPACE/ewald.cpp
+++ b/src/KSPACE/ewald.cpp
@@ -232,18 +232,27 @@ void Ewald::setup()
     kzmax = 1;
 
     err = rms(kxmax,xprd,natoms,q2);
+    if (!ISFINITE(err))
+      error->all(FLERR,"Non-numeric box dimensions - simulation unstable");
+
     while (err > accuracy) {
       kxmax++;
       err = rms(kxmax,xprd,natoms,q2);
     }
 
     err = rms(kymax,yprd,natoms,q2);
+    if (!ISFINITE(err))
+      error->all(FLERR,"Non-numeric box dimensions - simulation unstable");
+
     while (err > accuracy) {
       kymax++;
       err = rms(kymax,yprd,natoms,q2);
     }
 
     err = rms(kzmax,zprd_slab,natoms,q2);
+    if (!ISFINITE(err))
+      error->all(FLERR,"Non-numeric box dimensions - simulation unstable");
+
     while (err > accuracy) {
       kzmax++;
       err = rms(kzmax,zprd_slab,natoms,q2);
diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp
index 49b7cde12e..25b19c42d2 100644
--- a/src/KSPACE/ewald_disp.cpp
+++ b/src/KSPACE/ewald_disp.cpp
@@ -247,21 +247,33 @@ void EwaldDisp::setup()
   int kxmax = 1;
   int kymax = 1;
   int kzmax = 1;
+
   err = rms(kxmax,domain->h[0],natoms,q2,b2,M2);
+  if (!ISFINITE(err))
+    error->all(FLERR,"Non-numeric box dimensions - simulation unstable");
+
   while (err > accuracy) {
     kxmax++;
     err = rms(kxmax,domain->h[0],natoms,q2,b2,M2);
   }
+
   err = rms(kymax,domain->h[1],natoms,q2,b2,M2);
+  if (!ISFINITE(err))
+    error->all(FLERR,"Non-numeric box dimensions - simulation unstable");
   while (err > accuracy) {
     kymax++;
     err = rms(kymax,domain->h[1],natoms,q2,b2,M2);
   }
+
   err = rms(kzmax,domain->h[2]*slab_volfactor,natoms,q2,b2,M2);
+  if (!ISFINITE(err))
+    error->all(FLERR,"Non-numeric box dimensions - simulation unstable");
+
   while (err > accuracy) {
     kzmax++;
     err = rms(kzmax,domain->h[2]*slab_volfactor,natoms,q2,b2,M2);
   }
+
   nbox = MAX(kxmax,kymax);
   nbox = MAX(nbox,kzmax);
   double gsqxmx = unit[0]*unit[0]*kxmax*kxmax;
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index 725cb60d26..3ce7c4e302 100644
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@@ -1013,12 +1013,18 @@ void FixNH::couple()
     p_current[2] = tensor[2];
   }
 
+  if (!ISFINITE(p_current[0]) || !ISFINITE(p_current[1]) || !ISFINITE(p_current[2]))
+    error->all(FLERR,"Non-numeric pressure - simulation unstable");
+
   // switch order from xy-xz-yz to Voigt
 
   if (pstyle == TRICLINIC) {
     p_current[3] = tensor[5];
     p_current[4] = tensor[4];
     p_current[5] = tensor[3];
+
+    if (!ISFINITE(p_current[3]) || !ISFINITE(p_current[4]) || !ISFINITE(p_current[5]))
+      error->all(FLERR,"Non-numeric pressure - simulation unstable");
   }
 }