diff --git a/Makefile b/Makefile
index ae6fb28..18b10b7 100644
--- a/Makefile
+++ b/Makefile
@@ -1,7 +1,6 @@
-IRPF90 = irpf90 # -a #-d
-FC     = ifort 
-FCFLAGS= -O2 #-xT -fast-transcendentals
-
+IRPF90 = irpf90  #-a -d
+FC     = ifort
+FCFLAGS= -O3 -xT
 
 SRC=
 OBJ=
diff --git a/debug.irp.f b/debug.irp.f
index 3a3a2f4..f2fa12a 100644
--- a/debug.irp.f
+++ b/debug.irp.f
@@ -1,9 +1,42 @@
 program debug
  implicit none
- double precision :: eplf_integral, ao_overlap
  PROVIDE ao_prim_num_max
  integer :: i,j
+ integer :: k
  read(*,*) i,j
- print *,  eplf_integral(i,j,eplf_gamma,point)
- print *,  ao_overlap(i,j)
+
+ print *,  ''
+ do k=1,nucl_num
+  print *,  nucl_coord(k,:)
+ enddo
+ print *,  ''
+ print *, 'AO ', i
+ print *, 'prim num:', ao_prim_num(i)
+ print *, 'powers  :', ao_power(i,:)
+ print *, 'center  :', ao_nucl(i)
+ print *, 'expo  /  coef'
+ do k=1,ao_prim_num(i)
+   print *, ao_expo(k,i), ao_coef(k,i)
+ enddo
+
+ print *,  ''
+ print *, 'AO ', j
+ print *, 'prim num:', ao_prim_num(j)
+ print *, 'powers  :', ao_power(j,:)
+ print *, 'center  :', ao_nucl(j)
+ print *, 'expo  /  coef'
+ do k=1,ao_prim_num(j)
+   print *, ao_expo(k,j), ao_coef(k,j)
+ enddo
+
+ double precision :: ao_overlap, ao_overlap_numeric
+ print *, ''
+ print *, 'Overlap integral   :', ao_overlap(i,j)
+ print *, 'Overlap integral N :', ao_overlap_numeric(i,j)
+
+ double precision :: eplf_integral, eplf_integral_numeric
+ print *,  ''
+ print *, 'EPLF gamma         : ', eplf_gamma
+ print *, 'EPLF integral      :', eplf_integral(i,j,eplf_gamma,point)
+ print *, 'EPLF integral N    :', eplf_integral_numeric(i,j,eplf_gamma,point)
 end
diff --git a/overlap.irp.f b/overlap.irp.f
index 9122e95..e6bcc9d 100644
--- a/overlap.irp.f
+++ b/overlap.irp.f
@@ -1,3 +1,112 @@
+
+BEGIN_PROVIDER [ double precision, ao_overlap_matrix, (ao_num,ao_num) ]
+ implicit none
+ BEGIN_DOC  
+! Overlap matrix between the Atomic Orbitals
+ END_DOC
+
+ integer :: i, j
+ double precision :: ao_overlap
+ do j=1,ao_num
+  do i=j,ao_num
+   ao_overlap_matrix(i,j) = ao_overlap(i,j)
+  enddo
+ enddo
+ do j=1,ao_num
+  do i=1,j-1
+   ao_overlap_matrix(i,j) = ao_overlap(j,i)
+  enddo
+ enddo
+END_PROVIDER
+
+
+double precision function primitive_overlap_oneD_numeric(a,xa,i,b,xb,j)
+ implicit none
+ include 'constants.F'
+
+ 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
+ real :: x, xmin, xmax, dx
+
+ ASSERT (a>0.)
+ ASSERT (b>0.)
+ ASSERT (i>=0)
+ ASSERT (j>=0)
+
+ xmin = min(xa,xb) - 10.
+ xmax = max(xa,xb) + 10.
+ dx = (xmax-xmin)/real(Npoints)
+
+ real :: dtemp
+ dtemp = 0.
+ x = xmin
+ integer :: k
+ do k=1,Npoints
+  dtemp = dtemp + &
+   (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2))
+  x = x+dx
+ enddo
+ primitive_overlap_oneD_numeric = dtemp*dx
+
+end function
+
+double precision function ao_overlap_numeric(i,j)
+ implicit none
+ integer, intent(in) :: i, j
+ integer :: p,q,k
+ double precision :: integral(ao_prim_num_max,ao_prim_num_max)
+
+ double precision :: primitive_overlap_oneD_numeric
+
+ ASSERT(i>0)
+ ASSERT(j>0)
+ ASSERT(i<=ao_num)
+ ASSERT(j<=ao_num)
+
+ do q=1,ao_prim_num(j)
+  do p=1,ao_prim_num(i)
+   integral(p,q) =                         &
+   primitive_overlap_oneD_numeric (        &
+    ao_expo(p,i),                          &
+    nucl_coord(ao_nucl(i),1),              &
+    ao_power(i,1),                         &
+    ao_expo(q,j),                          &
+    nucl_coord(ao_nucl(j),1),              &
+    ao_power(j,1) ) *                      &
+   primitive_overlap_oneD_numeric (        &
+    ao_expo(p,i),                          &
+    nucl_coord(ao_nucl(i),2),              &
+    ao_power(i,2),                         &
+    ao_expo(q,j),                          &
+    nucl_coord(ao_nucl(j),2),              &
+    ao_power(j,2) ) *                      &
+   primitive_overlap_oneD_numeric (        &
+    ao_expo(p,i),                          &
+    nucl_coord(ao_nucl(i),3),              &
+    ao_power(i,3),                         &
+    ao_expo(q,j),                          &
+    nucl_coord(ao_nucl(j),3),              &
+    ao_power(j,3) ) 
+  enddo
+ enddo
+
+ do q=1,ao_prim_num(j)
+  do p=1,ao_prim_num(i)
+   integral(p,q) = integral(p,q)*ao_coef(p,i)*ao_coef(q,j)
+  enddo
+ enddo
+
+ ao_overlap_numeric = 0.
+ do q=1,ao_prim_num(j)
+  do p=1,ao_prim_num(i)
+   ao_overlap_numeric = ao_overlap_numeric + integral(p,q)
+  enddo
+ enddo
+ 
+end function
+
 subroutine gaussian_product(a,xa,b,xb,k,p,xp)
  implicit none
 ! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
@@ -139,24 +248,3 @@ double precision function ao_overlap(i,j)
  
 end function
 
-
-BEGIN_PROVIDER [ double precision, ao_overlap_matrix, (ao_num,ao_num) ]
- implicit none
- BEGIN_DOC  
-! Overlap matrix between the Atomic Orbitals
- END_DOC
-
- integer :: i, j
- double precision :: ao_overlap
- do j=1,ao_num
-  do i=j,ao_num
-   ao_overlap_matrix(i,j) = ao_overlap(i,j)
-  enddo
- enddo
- do j=1,ao_num
-  do i=1,j-1
-   ao_overlap_matrix(i,j) = ao_overlap(j,i)
-  enddo
- enddo
-END_PROVIDER
-