(#291) Improved FPE in g2g.

* Added some NaN checks in g2g.
* Changed tolerance for Ec/Fockc calculation.

Author: fedepedron

.gitignore

@@ -71,6 +71,6 @@ test/integration_test/iteration/main_gpu
 ###########################################################
 test/tests_engine/*.pyc
 ###########################################################
-# Profiling
+# Others
 ###########################################################
-test/profiling/*
+.vscode

g2g/init.cpp

@@ -319,18 +319,12 @@ void g2g_iteration(bool compute_energy, double* fort_energy_ptr,
                    double* fort_forces_ptr) {
   Timers timers;
 #ifdef _DEBUG 
-  feenableexcept(FE_INVALID);
-  feenableexcept(FE_DIVBYZERO);
-  //feenableexcept(FE_UNDERFLOW);
-  feenableexcept(FE_OVERFLOW);
+  feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW);
 #endif
   partition.solve(timers, compute_rmm, lda, compute_forces, compute_energy,
                   fort_energy_ptr, fort_forces_ptr, fortran_vars.OPEN);
 #ifdef _DEBUG 
-  fedisableexcept(FE_INVALID);
-  fedisableexcept(FE_DIVBYZERO);
-  //fedisableexcept(FE_UNDERFLOW);
-  fedisableexcept(FE_OVERFLOW);
+  fedisableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW);
 #endif
 }
 

g2g/pointxc/pbeOS_corr.h

@@ -234,8 +234,15 @@ __host__ __device__ void pbeOS_corr(scalar_type rho, scalar_type rs,
   scalar_type PREF = hz - GZ * T2 * hT / G;
   scalar_type FACT5 = GZ * ((scalar_type)2.0f * hT + t * hTT) / G;
   COMM = COMM - PREF * zet - uu * hTT - vv * hT - ww * (hZT - FACT5);
-//  if ( COMM != COMM ) { printf("NaN in COMM2 - hRS %E hRST %E T2 %E hT %E hTT %E \n",
-// hRS, hRST, T2, hT, hTT); };
+
+  #ifdef _DEBUG
+  if ( COMM != COMM ) { 
+      printf("NaN in COMM2 - hRS %E hRST %E T2 %E hT %E hTT %E \n", hRS, hRST, T2, hT, hTT);
+      printf("hTT: G3 %E t %E B %E Q8 %E Q9 %E FACT3 %E \n", G3, t, B, Q8, Q9, FACT3);
+      printf("hRST: rsthrd %E T2 %E hBT %E BEC %E ECRS %E \n", rsthrd, T2, hBT, BEC, ECRS);
+      printf("Q8: Q4 %E Q5 %E T2 %E \n", Q4, Q5, T2);
+  };
+  #endif
   dvc_a = COMM + PREF;
   dvc_b = COMM - PREF;
   // PBE POTENTIAL DONE !!!!!

g2g/pointxc/pbeOS_main.h

@@ -56,6 +56,7 @@ __host__ __device__ void pbeOS_main(
    // The limit 1e-18 is the one set in PBE's original paper, which
    // works in double precision. For single precision, we take 1e-12
    // so that rho^3 is still within precision.
+   // For Ec calculation, the limit is a bit lower (1e-10).
 #if FULL_DOUBLE
    const scalar_type MINIMUM_DENSITY_VALUE = (scalar_type) 1e-18;
 #else
@@ -72,7 +73,7 @@ __host__ __device__ void pbeOS_main(
    vxpbe_a = (scalar_type)0.0f;
    vxpbe_b = (scalar_type)0.0f;
 
- 
+
   scalar_type twodens = (scalar_type)2.0f * dens_a;
   if (twodens > MINIMUM_DENSITY_VALUE) {
     twodens2 = twodens * twodens;
@@ -117,7 +118,7 @@ __host__ __device__ void pbeOS_main(
   // Construct total density and contribution to Ex
   scalar_type rho = dens_a + dens_b;
   expbe = (expbe_a * (dens_a / rho) + expbe_b * (dens_b/rho));
-  if (rho < MINIMUM_DENSITY_VALUE) { return; };
+  if (rho < (MINIMUM_DENSITY_VALUE * (scalar_type) 1E5)) { return; };
 
   /*-------------------------------------------------------------//
   // PBE correlation
@@ -187,14 +188,16 @@ __host__ __device__ void pbeOS_main(
   vcpbe_a = vc_a + dvc_a;
   vcpbe_b = vc_b + dvc_b;
 
-// if (expbe_a != expbe_a) { printf("NaN in expbe_a \n");};
-// if (expbe_b != expbe_b) { printf("NaN in expbe_b \n");};
-// if (ec != ec) { printf("NaN in ec \n");};
-// if (h != h) { printf("NaN in h \n");};
-// if (vc_a != vc_a) { printf("NaN in vc_a \n");};
-// if (dvc_a != dvc_a) { printf("NaN in dvc_a \n");};
-// if (vc_b != vc_b) { printf("NaN in vc_b \n");};
-// if (dvc_b!= dvc_b) { printf("NaN in dvc_b \n");};
+#ifdef _DEBUG
+if (expbe_a != expbe_a) { printf("NaN in expbe_a \n");};
+if (expbe_b != expbe_b) { printf("NaN in expbe_b \n");};
+if (ec != ec)           { printf("NaN in ec \n");};
+if (h != h)             { printf("NaN in h \n");};
+if (vc_a != vc_a)       { printf("NaN in vc_a \n");};
+if (dvc_a != dvc_a)     { printf("NaN in dvc_a \n");};
+if (vc_b != vc_b)       { printf("NaN in vc_b \n");};
+if (dvc_b!= dvc_b)      { printf("NaN in dvc_b \n");};
+#endif
 }  // pbeOS_main
 
 #undef EASYPBE_PI