Fix: bug of NAN in SDFT

Qianruipku · Aug 27, 2023 · e9975c8 · e9975c8
1 parent db71b1c
commit e9975c8
Showing 1 changed file with 24 additions and 15 deletions.
diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_func.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_func.cpp
@@ -68,32 +68,38 @@ REAL Sto_Func<REAL>:: nxfd(REAL rawe)
 template<typename REAL>
 REAL Sto_Func<REAL>:: fdlnfd(REAL e)
 {
-    REAL e_mu = (e - mu) / this->tem ;
-    if(e_mu > 36)
+    REAL e_mu = (e - mu) / this->tem;
+    if (e_mu > 36)
         return 0;
-    else if(e_mu < -36)
+    else if (e_mu < -36)
         return 0;
     else
     {
         REAL f = 1 / (1 + exp(e_mu));
-        return (f * log(f) + (1.0-f) * log(1.0-f));
+        if (f == 0 || f == 1)
+            return 0;
+        else
+            return (f * log(f) + (1.0 - f) * log(1.0 - f));
     }
 }
 
 template<typename REAL>
 REAL Sto_Func<REAL>:: nfdlnfd(REAL rawe)
 {
-    REAL Ebar = (Emin + Emax)/2;
-	REAL DeltaE = (Emax - Emin)/2;
-    REAL ne_mu = (rawe * DeltaE + Ebar - mu) / this->tem ;
-    if(ne_mu > 36)
+    REAL Ebar = (Emin + Emax) / 2;
+    REAL DeltaE = (Emax - Emin) / 2;
+    REAL ne_mu = (rawe * DeltaE + Ebar - mu) / this->tem;
+    if (ne_mu > 36)
         return 0;
-    else if(ne_mu < -36)
+    else if (ne_mu < -36)
         return 0;
     else
     {
         REAL f = 1 / (1 + exp(ne_mu));
-        return f * log(f) + (1-f) * log(1-f);
+        if (f == 0 || f == 1)
+            return 0;
+        else
+            return f * log(f) + (1 - f) * log(1 - f);
     }
 }
 
@@ -103,14 +109,17 @@ REAL Sto_Func<REAL>:: n_root_fdlnfd(REAL rawe)
     REAL Ebar = (Emin + Emax)/2;
 	REAL DeltaE = (Emax - Emin)/2;
     REAL ne_mu = (rawe * DeltaE + Ebar - mu) / this->tem ;
-    if(ne_mu > 72)
+    if(ne_mu > 36)
         return 0;
-    else if(ne_mu < -72)
+    else if(ne_mu < -36)
         return 0;
     else
     {
         REAL f = 1 / (1 + exp(ne_mu));
-        return sqrt(-f * log(f) - (1-f) * log(1-f));
+        if (f == 0 || f == 1)
+            return 0;
+        else
+            return sqrt(-f * log(f) - (1-f) * log(1-f));
     }
 }
 
@@ -170,7 +179,7 @@ REAL Sto_Func<REAL>::ngauss(REAL rawe)
 	REAL DeltaE = (Emax - Emin)/2;
     REAL e = rawe * DeltaE + Ebar;
     REAL a = pow((targ_e-e),2)/2.0/pow(sigma,2);
-    if(a > 32)
+    if(a > 72)
         return 0;
     else
         return  exp(-a) /sqrt(TWOPI) / sigma ;
@@ -183,7 +192,7 @@ REAL Sto_Func<REAL>::nroot_gauss(REAL rawe)
 	REAL DeltaE = (Emax - Emin)/2;
     REAL e = rawe * DeltaE + Ebar;
     REAL a = pow((targ_e-e),2)/4.0/pow(sigma,2);
-    if(a > 32)
+    if(a > 72)
         return 0;
     else
         return  exp(-a) /sqrt(sqrt(TWOPI) * sigma) ;