Catch single atom pair constraint failures 65/4665/1
authorJim Phillips <jim@ks.uiuc.edu>
Mon, 8 Oct 2018 19:43:20 +0000 (14:43 -0500)
committerJim Phillips <jim@ks.uiuc.edu>
Mon, 8 Oct 2018 19:43:20 +0000 (14:43 -0500)
The non-iterative solution for constraining a single pair of
bonded atoms contains a sqrt() that returns a NaN if the argument
is negative.  Catch this case as for the interative method to
prevent NaNs from propagating to velocities, temperature,
pressure, and positions; and eventually causing a segfault that
looks to the user like a bug rather than unstable dynamics.

Change-Id: Ifcd34cd8af18ef337b214ea3f5a07c3622f12cc9

src/HomePatch.C
src/Settle.C
src/Settle.h

index 152b4d3..87ffd60 100644 (file)
@@ -2356,9 +2356,9 @@ int HomePatch::rattle1(const BigReal timestep, Tensor *virial,
     if (icnt == 1) {
       rattlePair<1>(&rattleParam[posParam],
         refx, refy, refz,
-        posx, posy, posz);
+        posx, posy, posz,
+        consFailure);
       done = true;
-      consFailure = false;
     } else {
       rattleN(icnt, &rattleParam[posParam],
         refx, refy, refz,
index ae16e8c..6894d74 100644 (file)
@@ -545,7 +545,7 @@ void settle1_SIMD(const Vector *ref, Vector *pos,
 template <int veclen>
 void rattlePair(const RattleParam* rattleParam,
   const BigReal *refx, const BigReal *refy, const BigReal *refz,
-  BigReal *posx, BigReal *posy, BigReal *posz) {
+  BigReal *posx, BigReal *posy, BigReal *posz, bool& consFailure) {
 
   int a = rattleParam[0].ia;
   int b = rattleParam[0].ib;
@@ -566,7 +566,15 @@ void rattlePair(const RattleParam* rattleParam,
   BigReal rma = rattleParam[0].rma;
   BigReal rmb = rattleParam[0].rmb;
 
-  BigReal gab = (-rpab + sqrt(rpab*rpab + refsq*diffsq))/(refsq*(rma + rmb));
+  BigReal gab;
+  BigReal sqrtarg = rpab*rpab + refsq*diffsq;
+  if ( sqrtarg < 0. ) {
+    consFailure = true;
+    gab = 0.;
+  } else {
+    consFailure = false;
+    gab = (-rpab + sqrt(sqrtarg))/(refsq*(rma + rmb));
+  }
 
   BigReal dpx = rabx * gab;
   BigReal dpy = raby * gab;
@@ -633,7 +641,7 @@ void rattleN(const int icnt, const RattleParam* rattleParam,
 //
 template void rattlePair<1>(const RattleParam* rattleParam,
   const BigReal *refx, const BigReal *refy, const BigReal *refz,
-  BigReal *posx, BigReal *posy, BigReal *posz);
+  BigReal *posx, BigReal *posy, BigReal *posz, bool& consFailure);
 template void settle1_SIMD<2>(const Vector *ref, Vector *pos,
   BigReal mOrmT, BigReal mHrmT, BigReal ra,
   BigReal rb, BigReal rc, BigReal rra);
index 7fea736..715be6e 100644 (file)
@@ -63,7 +63,7 @@ struct RattleParam {
 template <int veclen>
 void rattlePair(const RattleParam* rattleParam,
   const BigReal *refx, const BigReal *refy, const BigReal *refz,
-  BigReal *posx, BigReal *posy, BigReal *posz);
+  BigReal *posx, BigReal *posy, BigReal *posz, bool& consFailure);
 
 void rattleN(const int icnt, const RattleParam* rattleParam,
   const BigReal *refx, const BigReal *refy, const BigReal *refz,