diff --git a/Makefile b/Makefile
index a2e209d..2481d91 100644
--- a/Makefile
+++ b/Makefile
@@ -1,16 +1,17 @@
 # MPI-ifort
 IRPF90 = irpf90 -DMPI #-a -d 
 FC     = mpif90 -static-intel -static-libgcc -xT
-FCFLAGS= -O3 -g
+FCFLAGS= -O3
 
 # Gfortran
 #FC     = gfortran -ffree-line-length-none -static-libgcc -static
 #FCFLAGS= -O3 -ffast-math -L ~/QCIO/lib
 
 # Mono
-#IRPF90 = irpf90
-#FC     = ifort -static-intel -static-libgcc 
-#FCFLAGS= -O3 -axP 
+IRPF90 = irpf90
+FC     = ifort -static-intel -static-libgcc 
+FCFLAGS= -O3 -axP 
+FCFLAGS= -O3 -xT -g 
 
 SRC=
 OBJ=
diff --git a/Util.irp.f b/Util.irp.f
index 44261a9..dd06a5f 100644
--- a/Util.irp.f
+++ b/Util.irp.f
@@ -1,38 +1,16 @@
 
-double precision function Boys(x,n)
+recursive double precision function Boys(x,n) result(res)
  implicit none
  include 'constants.F'
 
  real, intent(in)     :: x
  integer, intent(in)  :: n
 
- integer :: k
-
- real, parameter    :: thr = 6.
- integer ,parameter :: Nmax = 20
-
- double precision :: fact,fact2
-
+ ASSERT (x > 0.)
  if (n == 0) then
-   if (x > thr) then
-     Boys = 0.5d0*sqrt(pi/x)
-   else
-     Boys = 1.d0/dble(2*n+1)
-     do k=1,Nmax
-       Boys = Boys + (-x)**k/dble(fact(k)*(2*n+2*k+1)) 
-     enddo
-     
-   endif
-
+   res = sqrt(pi/(4.*x))*erf(sqrt(x))
  else
-   if (x > thr) then
-     Boys = fact2(2*n-1)*0.5d0**(n+1)*sqrt(pi/x**(2*n+1))
-   else
-     Boys = 1.d0/dble(2*n+1)
-     do k=1,Nmax
-       Boys = Boys + (-x)**k/dble(fact(k)*(2*n+2*k+1)) 
-     enddo
-   endif
+   res = (dble(2*n-1) * Boys(x,(n-1)) - exp(-x) )/(2.*x)
  endif
 
 end function
@@ -96,3 +74,42 @@ double precision function fact(n)
 end function
 
 
+double precision function rintgauss(n)
+  implicit none
+  include 'constants.F'
+
+  integer :: n
+
+  rintgauss = sqrt(pi)
+  if ( n == 0 ) then
+    return
+  else if ( n == 1 ) then
+    rintgauss = 0.
+  else if ( mod(n,2) == 1) then
+    rintgauss = 0.
+  else
+    double precision :: fact2
+    rintgauss = rintgauss/(2.**(n/2))
+    rintgauss = rintgauss * fact2(n-1)
+  endif
+end function
+
+double precision function goverlap(gamA,gamB,nA)
+  implicit none
+
+  real    :: gamA, gamB
+  integer :: nA(3)
+
+  double precision :: gamtot
+  gamtot = gamA+gamB
+
+  goverlap=1.0
+
+  integer :: l
+  double precision :: rintgauss
+  do l=1,3
+    goverlap = goverlap * rintgauss(nA(l)+nA(l))/ (gamtot**(0.5+float(nA(l))))
+  enddo
+
+end function
+
diff --git a/ao.irp.f b/ao.irp.f
index 1181f6c..3c57c65 100644
--- a/ao.irp.f
+++ b/ao.irp.f
@@ -123,10 +123,10 @@ END_PROVIDER
 !$OMP END CRITICAL (qcio_critical)
 
  double precision :: norm, norm2
- double precision :: overlap
+ double precision :: goverlap
  do i=1,ao_num
   do j=1,ao_prim_num(i)
-   norm = overlap(ao_expo(j,i),ao_expo(j,i),ao_power(i,:))
+   norm = goverlap(ao_expo(j,i),ao_expo(j,i),ao_power(i,:))
    norm = sqrt(norm)
    ao_coef(j,i) = buffer(j,i)/norm
   enddo
@@ -134,43 +134,3 @@ END_PROVIDER
 
 END_PROVIDER
 
-double precision function rintgauss(n)
-  implicit none
-
-  integer :: n
-  double precision :: pi
-  pi = acos(-1.)
-
-  rintgauss = sqrt(pi)
-  if ( n == 0 ) then
-    return
-  else if ( n == 1 ) then
-    rintgauss = 0.
-  else if ( mod(n,2) == 1) then
-    rintgauss = 0.
-  else
-    double precision :: fact2
-    rintgauss = rintgauss/(2.**(n/2))
-    rintgauss = rintgauss * fact2(n-1)
-  endif
-end function
-
-double precision function overlap(gamA,gamB,nA)
-  implicit none
-
-  real    :: gamA, gamB
-  integer :: nA(3)
-
-  double precision :: gamtot
-  gamtot = gamA+gamB
-
-  overlap=1.0
-
-  integer :: l
-  double precision :: rintgauss
-  do l=1,3
-    overlap = overlap * rintgauss(nA(l)+nA(l))/ (gamtot**(0.5+float(nA(l))))
-  enddo
-
-end function
-
diff --git a/eplf.irp.f b/eplf.irp.f
index 01981bf..e60bbfa 100644
--- a/eplf.irp.f
+++ b/eplf.irp.f
@@ -3,7 +3,7 @@ BEGIN_PROVIDER  [ real, eplf_gamma ]
  BEGIN_DOC  
 ! Value of the gaussian for the EPLF
  END_DOC
- eplf_gamma = 10.
+ eplf_gamma = 10000.
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
@@ -27,27 +27,33 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
 ! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF
  END_DOC
  integer :: i, j, k, l
+ double precision :: t
  PROVIDE ao_eplf_integral_matrix
  PROVIDE mo_coef
  do i=1,mo_num
   do j=i,mo_num
     mo_eplf_integral_matrix(j,i) = 0.
   enddo
- enddo
 
-  do i=1,mo_num
-   do k=1,ao_num
-    if (mo_coef(k,i) /= 0.) then
-     do j=i,mo_num
-       do l=1,ao_num
-        mo_eplf_integral_matrix(j,i) = mo_eplf_integral_matrix(j,i) + &
-          mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(l,k)
-       enddo
-     enddo
-    endif
+
+  do k=1,ao_num
+   if (abs(mo_coef(k,i)) /= 0.) then
+    do l=1,ao_num
+     t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k)
+     if (abs(ao_eplf_integral_matrix(l,k))>1.d-16) then
+      do j=i,mo_num
+       mo_eplf_integral_matrix(j,i) = mo_eplf_integral_matrix(j,i) + &
+         t*mo_coef_transp(j,l)
+      enddo
+     endif
+    enddo
+   endif
   enddo
 
-  do j=i,mo_num
+ enddo
+
+ do i=1,mo_num
+  do j=i+1,mo_num
     mo_eplf_integral_matrix(i,j) = mo_eplf_integral_matrix(j,i)
   enddo
  enddo
@@ -132,7 +138,7 @@ double precision function ao_eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,
  real, intent(in)    :: a,b,gmma ! Exponents
  real, intent(in)    :: xa,xb,xr ! Centers
  integer, intent(in) :: i,j   ! Powers of xa and xb
- integer,parameter :: Npoints=1000
+ integer,parameter :: Npoints=10000
  real :: x, xmin, xmax, dx
 
  ASSERT (a>0.)
@@ -165,12 +171,13 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center)
  integer :: p,q,k
  double precision :: integral
  double precision :: ao_eplf_integral_primitive_oneD_numeric
- real :: gmma, center(3)
+ real :: gmma, center(3), c
 
  
   ao_eplf_integral_numeric = 0.
   do q=1,ao_prim_num(j)
    do p=1,ao_prim_num(i)
+    c = ao_coef(p,i)*ao_coef(q,j)
     integral =                              &
     ao_eplf_integral_primitive_oneD_numeric(   &
      ao_expo(p,i),                          &
@@ -199,7 +206,7 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center)
      ao_power(j,3),                         &
      gmma,                                  &
      center(3))
-    ao_eplf_integral_numeric = ao_eplf_integral_numeric + integral*ao_coef(p,i)*ao_coef(q,j)
+    ao_eplf_integral_numeric = ao_eplf_integral_numeric + c*integral
    enddo
   enddo
   
diff --git a/overlap.irp.f b/overlap.irp.f
index e6bcc9d..a6b4001 100644
--- a/overlap.irp.f
+++ b/overlap.irp.f
@@ -27,7 +27,7 @@ double precision function primitive_overlap_oneD_numeric(a,xa,i,b,xb,j)
  real, intent(in)    :: a,b   ! Exponents
  real, intent(in)    :: xa,xb ! Centers
  integer, intent(in) :: i,j   ! Powers of xa and xb
- integer,parameter :: Npoints=1000
+ integer,parameter :: Npoints=10000
  real :: x, xmin, xmax, dx
 
  ASSERT (a>0.)
diff --git a/test_1d.irp.f b/test_1d.irp.f
index 660c188..7dfdfec 100644
--- a/test_1d.irp.f
+++ b/test_1d.irp.f
@@ -1,7 +1,7 @@
 program debug
  PROVIDE ao_prim_num_max
- read(*,*) eplf_gamma
- TOUCH eplf_gamma
+!read(*,*) eplf_gamma
+!TOUCH eplf_gamma
  call run()
 end