From 3987b9794d6e17d00860d0543d12646d3e2d63f0 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Mon, 8 Jun 2020 12:52:08 +0200
Subject: [PATCH] removed a lot of floating point exceptions in DFT

---
 src/dft_utils_func/on_top_from_ueg.irp.f | 34 +++++++++++++++++++-----
 1 file changed, 27 insertions(+), 7 deletions(-)

diff --git a/src/dft_utils_func/on_top_from_ueg.irp.f b/src/dft_utils_func/on_top_from_ueg.irp.f
index 70560a7a..1965f37f 100644
--- a/src/dft_utils_func/on_top_from_ueg.irp.f
+++ b/src/dft_utils_func/on_top_from_ueg.irp.f
@@ -35,7 +35,11 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b)
  if (dabs(rho) > 1.d-12) then
   rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
   x = -d2*rs
-  g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*exp(x)
+  if(dabs(x).lt.50.d0)then
+   g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*dexp(x)
+  else
+   g0_UEG_mu_inf= 0.d0
+  endif
  else
   g0_UEG_mu_inf= 0.d0
  endif
@@ -67,7 +71,11 @@ double precision function g0_UEG_mu(mu,rho_a,rho_b)
  kf = (alpha*rs)**(-1d0)
  zeta = mu / kf
  x = -d2*rs*h_func(zeta)/ahd 
- g0_UEG_mu = (exp(x)/2d0) * (1d0- B*(h_func(zeta)/ahd)*rs + C*((h_func(zeta)**2d0)/(ahd**2d0))*(rs**2d0) + D*((h_func(zeta)**3d0)/(ahd**3d0))*(rs**3d0) + E*((h_func(zeta)**4d0)/(ahd**4d0))*(rs**4d0) )
+ if(dabs(x).lt.50.d0)then
+  g0_UEG_mu = (dexp(x)/2d0) * (1d0- B*(h_func(zeta)/ahd)*rs + C*((h_func(zeta)**2d0)/(ahd**2d0))*(rs**2d0) + D*((h_func(zeta)**3d0)/(ahd**3d0))*(rs**3d0) + E*((h_func(zeta)**4d0)/(ahd**4d0))*(rs**4d0) )
+ else
+  g0_UEG_mu = 0.d0
+ endif
  
 end
 
@@ -81,11 +89,11 @@ double precision function h_func(zeta)
  pi = 4d0 * datan(1d0)
  ahd = -0.36583d0
  alpha = (4d0/(9d0*pi))**(1d0/3d0)
- a1 = -(6d0*alpha/pi)*(1d0-log(2d0))
+ a1 = -(6d0*alpha/pi)*(1d0-dlog(2d0))
  b1 = 1.4919d0
  b3 = 1.91528d0
  a2 = ahd * b3
- b2 = (a1 - (b3*alpha/sqrt(pi)))/ahd
+ b2 = (a1 - (b3*alpha/dsqrt(pi)))/ahd
 
  h_func = (a1*zeta**2d0 + a2*zeta**3d0) / (1d0 + b1*zeta + b2*zeta**2d0 + b3*zeta**3d0)
 end
@@ -111,11 +119,23 @@ end
   D1 = -0.0127713d0
   E1 = 0.00185898d0
   B1 = 0.7317d0 - F1
-  rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0) 
+  if(dabs(rho).gt.1.d-20)then
+   rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0) 
+  else
+   rs = (3.d0 / (4.d0*pi*1.d-20))**(1.d0/3.d0) 
+  endif
    
   g0 = g0_UEG_mu_inf(rho_a, rho_b)
-  dg0drs = 0.5d0*((-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)-F1*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*exp(-F1*rs)
-  dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs
+  if(dabs(F1*rs).lt.50.d0)then
+   dg0drs = 0.5d0*((-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)-F1*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*dexp(-F1*rs)
+  else
+   dg0drs = 0.d0
+  endif
+  if(dabs(rho).gt.1.d-20)then
+   dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs
+  else
+   dg0drho = -((6.d0*dsqrt(pi)*1.d-40)**(-2.d0/3.d0))*dg0drs
+  endif
  
   end subroutine g0_dg0