diff --git a/EZFIO.tar.gz b/EZFIO.tar.gz
index 35e6d6d..8efc885 100644
Binary files a/EZFIO.tar.gz and b/EZFIO.tar.gz differ
diff --git a/bin/to_ezfio.exe b/bin/to_ezfio.exe
index 60898b9..2d1d745 100755
Binary files a/bin/to_ezfio.exe and b/bin/to_ezfio.exe differ
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index fa6ad29..61b37b8 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -241,8 +241,8 @@ BEGIN_PROVIDER [ real, eplf_value_p ]
    aa = min(1.d0,aa)
    ab = min(1.d0,ab)
    aa = -dlog(aa)
-   ab = -dlog(ab)
    aa = dsqrt(aa)
+   ab = -dlog(ab)
    ab = dsqrt(ab)
    eplf_value_p = (aa-ab)/(aa+ab+eps)
  else
@@ -381,28 +381,6 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
  
 end function
 
-double precision function rint(i,xa,a,j,xb,b,xr,gmma)
- implicit double precision(a-h,o-z)
- include 'constants.F'
- beta=a+b+gmma
- w=1.d0/dsqrt(beta)
- alpha=a*xa**2+b*xb**2+gmma*xr**2
- delta=xa*a+xb*b+xr*gmma
- u=delta/beta-xa
- v=delta/beta-xb
- ac=0.d0
- do n=0,i
-  do m=0,j
-   if(mod(n+m,2).eq.0)then
-    k=(n+m)/2
-    f=binom(n,i)*binom(m,j)*fact2(n+m-1)
-    ac=ac+u**(i-n)*v**(j-m)*w**(n+m)*f/2.d0**k
-   endif
-  enddo
- enddo
- rint=ac*dsqrt(pi)*w*dexp(-alpha+delta**2/beta)
-end
-
 
 double precision function ao_eplf_integral(i,j,gmma,center)
  implicit none
diff --git a/src/overlap.irp.f b/src/overlap.irp.f
index ff62a8b..3027a31 100644
--- a/src/overlap.irp.f
+++ b/src/overlap.irp.f
@@ -141,7 +141,7 @@ double precision function primitive_overlap_oneD(a,xa,i,b,xb,j)
  real :: xp
  real :: p
  double precision :: S(0:i+1,0:j)
- double precision :: inv_p, di(max(i,j)), dj(j)
+ double precision :: inv_p, di(max(i,j)), dj(j), c
 
  ASSERT (a>0.)
  ASSERT (b>0.)
@@ -153,15 +153,24 @@ double precision function primitive_overlap_oneD(a,xa,i,b,xb,j)
  xp = (a*xa+b*xb)
  inv_p = 1.d0/p
  xp = xp*inv_p
- S(0,0) = dexp(-a*b*inv_p*(xa-xb)**2)* dsqrt(pi*inv_p)
+
+ c = a*b*inv_p*(xa-xb)**2
+!if ( c > 32.d0 ) then   ! Cut-off on exp(-32)
+! primitive_overlap_oneD = 0.d0
+! return
+!endif
+
+ c =  dexp(-c)
+
+ S(0,0) = 1.d0
 
  ! Obara-Saika recursion
 
  if (i>0) then
-   S(1,0) = (xp-xa) * S(0,0)
+   S(1,0) = xp-xa
  endif
  if (j>0) then
-   S(0,1) = (xp-xb) * S(0,0)
+   S(0,1) = xp-xb
  endif
 
  do ii=1,max(i,j)
@@ -187,7 +196,7 @@ double precision function primitive_overlap_oneD(a,xa,i,b,xb,j)
   enddo
  enddo
 
- primitive_overlap_oneD = S(i,j)
+ primitive_overlap_oneD = S(i,j)*c* dsqrt(pi*inv_p)
 
 end function
 
diff --git a/version b/version
index 4b3683c..e2e9306 100644
--- a/version
+++ b/version
@@ -1 +1 @@
-VERSION=1.0.3
+VERSION=1.0.4