diff --git a/src/dft_utils_func/routines_exc_sr_lda.irp.f b/src/dft_utils_func/routines_exc_sr_lda.irp.f
index d64d826d..ea1dcd69 100644
--- a/src/dft_utils_func/routines_exc_sr_lda.irp.f
+++ b/src/dft_utils_func/routines_exc_sr_lda.irp.f
@@ -138,6 +138,8 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b)
 !Density and kF
  rho_a_2=rho_a*2.D0
  akf = ckf*(rho_a_2**f13)
+ ! Avoid division by zero
+ if (akf == 0.d0) akf = 1.d-20
  a = mu/(z2*akf)
  a2 = a*a
  a3 = a2*a
@@ -169,6 +171,7 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b)
 !Density and kF
  rho_b_2= rho_b * 2.d0
  akf = ckf*(rho_b_2**f13)
+ if (akf == 0.d0) akf = 1.d-20
  a = mu/(z2*akf)
  a2 = a*a
  a3 = a2*a
diff --git a/src/utils/need.irp.f b/src/utils/need.irp.f
index d3bedc30..0e36351c 100644
--- a/src/utils/need.irp.f
+++ b/src/utils/need.irp.f
@@ -92,42 +92,47 @@
       end
 
       double precision function erf0(x)
-      implicit double precision (a-h,o-z)
-      if(x.lt.0.d0)then
-        erf0=-gammp(0.5d0,x**2)
+      implicit none
+      double precision, intent(in) :: x
+      double precision, external :: gammp
+      if(x < 0.d0)then
+         erf0=-gammp(0.5d0,x*x)
       else
-        erf0=gammp(0.5d0,x**2)
+         erf0=gammp(0.5d0,x*x)
       endif
       end
 
 
       double precision function gammp(a,x)
-      implicit double precision (a-h,o-z)
-      if(x.lt.0..or.a.le.0.)stop 'error in gammp'
-      if(x.lt.a+1.)then
+      implicit none
+      double precision, intent(in) :: a, x
+      double precision             :: gln, gammcf
+      if(x<0.d0.or.a<=0.d0) then
+         stop 'error in gammp'
+      endif
+      if(x < a+1.d0)then
         call gser(gammp,a,x,gln)
       else
         call gcf(gammcf,a,x,gln)
-        gammp=1.-gammcf
+        gammp=1.d0-gammcf
       endif
-      return
       end
 
 
       subroutine gser(gamser,a,x,gln)
       implicit double precision (a-h,o-z)
-      parameter (itmax=100,eps=3.e-7)
+      parameter (itmax=100,eps=3.d-7)
       gln=gammln(a)
-      if(x.le.0.)then
-        if(x.lt.0.) stop 'error in gser'
-        gamser=0.
+      if(x.le.0.d0)then
+        if(x.lt.0.d0) stop 'error in gser'
+        gamser=0.d0
         return
       endif
       ap=a
-      sum=1./a
+      sum=1.d0/a
       del=sum
       do 11 n=1,itmax
-        ap=ap+1.
+        ap=ap+1.d0
         del=del*x/ap
         sum=sum+del
         if(abs(del).lt.abs(sum)*eps)go to 1
@@ -139,14 +144,14 @@
 
       subroutine gcf(gammcf,a,x,gln)
       implicit double precision (a-h,o-z)
-      parameter (itmax=100,eps=3.e-7)
+      parameter (itmax=100,eps=3.d-7)
       gln=gammln(a)
-      gold=0.
-      a0=1.
+      gold=0.d0
+      a0=1.d0
       a1=x
-      b0=0.
-      b1=1.
-      fac=1.
+      b0=0.d0
+      b1=1.d0
+      fac=1.d0
       do 11 n=1,itmax
         an=float(n)
         ana=an-a
@@ -155,8 +160,8 @@
         anf=an*fac
         a1=x*a0+anf*a1
         b1=x*b0+anf*b1
-        if(a1.ne.0.)then
-          fac=1./a1
+        if(a1.ne.0.d0)then
+          fac=1.d0/a1
           g=b1*fac
           if(abs((g-gold)/g).lt.eps)go to 1
           gold=g
diff --git a/src/utils/prim_in_r.irp.f b/src/utils/prim_in_r.irp.f
index af68ff6d..e0d124da 100644
--- a/src/utils/prim_in_r.irp.f
+++ b/src/utils/prim_in_r.irp.f
@@ -24,8 +24,9 @@ double precision function primitive_value_explicit(power_prim,center_prim,alpha,
 end
 
 double precision function give_pol_in_r(r,pol,center, alpha,iorder, max_dim)
- double precision    :: r(3), center(3), alpha,pol(0:max_dim,3)
+ implicit none
  integer, intent(in) :: iorder(3), max_dim
+ double precision    :: r(3), center(3), alpha,pol(0:max_dim,3)
  integer :: i,m
  double precision :: gauss(3), x
  gauss  = 0.d0
@@ -33,7 +34,7 @@ double precision function give_pol_in_r(r,pol,center, alpha,iorder, max_dim)
  do m = 1, 3
   x = r(m) - center(m)
   do i = 0, iorder(m)
-   gauss(m) += pol(i,m) * dexp(-alpha *x**2 ) * x**i 
+   gauss(m) += pol(i,m) * dexp(-alpha *x*x ) * x**i 
   enddo
  enddo
  give_pol_in_r = gauss(1) * gauss(2) * gauss(3)