diff --git a/src/Util.irp.f b/src/Util.irp.f
index eae7655..3907f3c 100644
--- a/src/Util.irp.f
+++ b/src/Util.irp.f
@@ -15,6 +15,12 @@
 !
 !end function
 
+
+double precision function binom(k,n)
+implicit double precision(a-h,o-z)
+ binom=fact(n)/(fact(k)*fact(n-k))
+end
+
 double precision function fact2(n)
   implicit none
   integer :: n
@@ -80,7 +86,7 @@ double precision function rintgauss(n)
 
   integer :: n
 
-  rintgauss = sqrt(pi)
+  rintgauss = dsqrt_pi
   if ( n == 0 ) then
     return
   else if ( n == 1 ) then
diff --git a/src/constants.F b/src/constants.F
index 7d671f3..0b2e29d 100644
--- a/src/constants.F
+++ b/src/constants.F
@@ -1 +1,2 @@
    double precision, parameter :: pi=3.14159265359d0
+   double precision, parameter :: dsqrt_pi=1.77245385091d0
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index eb63902..fa6ad29 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -5,9 +5,12 @@ BEGIN_PROVIDER  [ real, eplf_gamma ]
  END_DOC
  include 'constants.F'
  real            :: eps, N
- N = 0.1
- eps = -real(dlog(tiny(1.d0)))
- eplf_gamma = (4./3.*pi*density_value_p/N)**(2./3.) * eps
+! N = 0.1
+! eps = -real(dlog(tiny(1.d0)))
+! eplf_gamma = (4./3.*pi*density_value_p/N)**(2./3.) * eps
+ N=0.001
+ eps = 50.
+ eplf_gamma = eps**2 *0.5d0*(4./9.*pi*density_value_p * N**(-1./3.))**(2./3.)
 !eplf_gamma = 1.e10
 !eplf_gamma = 1.e5
 END_PROVIDER
@@ -339,155 +342,67 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
   double precision , intent(in) :: gmma  ! eplf_gamma 
   real, intent(in)    :: xa,xb,xr ! Centers
   integer, intent(in) :: i,j   ! Powers of xa and xb
-  integer :: ii, jj, kk, ll
-  real :: xp1,xp
-  real :: p1,p
-  double precision :: S(0:i+1,0:j+1)
-  double precision :: inv_p(2), di(max(i,j)), dj(j), c
+  integer :: ii, jj, kk
  
   ASSERT (a>0.)
   ASSERT (b>0.)
   ASSERT (i>=0)
   ASSERT (j>=0)
  
-  ! Gaussian product
-  real :: t(2), xab(2), ab(2)
-  inv_p(1) = 1.d0/(a+b)
-  p1 = a+b
-  ab(1) = a*b
-  inv_p(2) = 1.d0/(p1+gmma)
-  t(1) = (a*xa+b*xb)
-  xab(1) = xa-xb
-  xp1 = t(1)*inv_p(1)
-  p = p1+gmma
-  ab(2) = p1*gmma
-  t(2) = (p1*xp1+gmma*xr)
-  xab(2) = xp1-xr
-  xp = t(2)*inv_p(2)
-  c = dble(ab(1)*inv_p(1)*xab(1)**2 + &
-              ab(2)*inv_p(2)*xab(2)**2)
+  double precision :: alpha, beta, delta, c
+  alpha = a*xa*xa + b*xb*xb + gmma*xr*xr
+  delta = a*xa    + b*xb    + gmma*xr
+  beta  = a       + b       + gmma
+  c = alpha-delta**2/beta
   
-! double precision, save :: c_accu(2)
-! c_accu(1) += c
-! c_accu(2) += 1.d0
-! print *, c_accu(1)/c_accu(2)
-
   if ( c > 32.d0 ) then   ! Cut-off on exp(-32)
    ao_eplf_integral_primitive_oneD = 0.d0
    return
   endif
 
+  double precision :: u,v,w
   c = exp(-c) 
-  
-  ! Obara-Saika recursion
-  S(0,0) = 1.d0
- 
-  do ii=1,max(i,j)
-   di(ii) = 0.5d0*inv_p(2)*dble(ii)
-  enddo
- 
-  xab(1) = xp-xa
-  xab(2) = xp-xb
-  S(1,0) = xab(1) ! * S(0,0)
-  if (i>1) then
-   do ii=1,i-1
-    S(ii+1,0) = xab(1) * S(ii,0) + di(ii)*S(ii-1,0)
-   enddo
-  endif
- 
-  S(0,1) = xab(2) ! * S(0,0)
-  if (j>1) then
-   do jj=1,j-1
-    S(0,jj+1) = xab(2) * S(0,jj) + di(jj)*S(0,jj-1)
-   enddo
-  endif
- 
-  do jj=1,j
-   S(1,jj) = xab(1) * S(0,jj) + di(jj) * S(0,jj-1)
-   do ii=2,i
-    S(ii,jj) = xab(1) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
-   enddo
-  enddo
- 
-  ao_eplf_integral_primitive_oneD = dsqrt(pi*inv_p(2))*c*S(i,j)
- 
- end function
+  w = 1.d0/sqrt(beta)
+  u=delta/beta-xa
+  v=delta/beta-xb
+
+  double precision :: fact2, binom, f, ac
+  ac = 0.d0
+  integer :: istart(0:1)
+  istart = (/ 0, 1 /)
+  do ii=0,i
+   do jj=istart(mod(ii,2)),j,2
+     kk=(ii+jj)/2
+     f=binom(ii,i)*binom(jj,j)*fact2(ii+jj-1)
+     ac += u**(i-ii)*v**(j-jj)*w**(ii+jj)*f/dble(2**kk)
+   enddo
+  enddo
+  ao_eplf_integral_primitive_oneD = dsqrt_pi*w*c*ac
+ 
+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_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
-! implicit none
-! include 'constants.F'
-!!
-! real, intent(in)    :: a,b,gmma ! Exponents
-! double precision, intent(in)    :: gmma
-! real, intent(in)    :: xa,xb,xr ! Centers
-! integer, intent(in) :: i,j   ! Powers of xa and xb
-! integer :: ii, jj, kk, ll
-! real :: xp1,xp
-! real :: p1,p
-! double precision :: xpa, xpb
-! double precision :: inv_p(2),S00, c
-! double precision :: ObaraS
-!!
-! ASSERT (a>0.)
-! ASSERT (b>0.)
-! ASSERT (i>=0)
-! ASSERT (j>=0)
-!!
-! ! Inlined Gaussian products (same as call gaussian_product)
-! real :: t(2), xab(2), ab(2)
-! inv_p(1) = 1.d0/(a+b)
-! p1 = a+b
-! ab(1) = a*b
-! inv_p(2) = 1.d0/(p1+gmma)
-! t(1) = (a*xa+b*xb)
-! xab(1) = xa-xb
-! xp1 = t(1)*inv_p(1)
-! p = p1+gmma
-! ab(2) = p1*gmma
-! t(2) = (p1*xp1+gmma*xr)
-! xab(2) = xp1-xr
-! xp = t(2)*inv_p(2)
-! c = exp(- real(ab(1)*inv_p(1)*xab(1)**2 + &
-!             ab(2)*inv_p(2)*xab(2)**2) )
-!!
-! xpa = xp-xa
-! xpb = xp-xb
-! S00 = sqrt(real(pi*inv_p(2)))*c
-! ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p(2),S00)
-!!
-!end function
-!!
-!recursive double precision function ObaraS(i,j,xpa,xpb,inv_p,S00) result(res)
-! implicit none
-! integer, intent(in) :: i, j
-! double precision, intent(in) :: xpa, xpb, inv_p
-! double precision,intent(in) :: S00
-!!
-! if (i == 0) then
-!  if (j == 0) then
-!   res = S00
-!  else ! (j>0)
-!   res = xpb*ObaraS(0,j-1,xpa,xpb,inv_p,S00)
-!   if (j>1) then
-!    res += 0.5d0*dble(j-1)*inv_p*ObaraS(0,j-2,xpa,xpb,inv_p,S00)
-!   endif
-!  endif ! (i==0).and.(j>0)
-! else   ! (i>0)
-!  if (j==0) then
-!   res = xpa*ObaraS(i-1,0,xpa,xpb,inv_p,S00) 
-!   if (i>1) then
-!     res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,0,xpa,xpb,inv_p,S00) 
-!   endif
-!  else  ! (i>0).and.(j>0)
-!   res = xpa * ObaraS(i-1,j,xpa,xpb,inv_p,S00) 
-!   if (i>1) then
-!    res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00) 
-!   endif
-!   res += 0.5d0*dble(j)*inv_p*ObaraS(i-1,j-1,xpa,xpb,inv_p,S00) 
-!  endif  ! (i>0).and.(j>0)
-! endif  ! (i>0)
-!!
-!end function
 
 double precision function ao_eplf_integral(i,j,gmma,center)
  implicit none