diff --git a/EZFIO.tar.gz b/EZFIO.tar.gz
deleted file mode 100644
index 7075609..0000000
Binary files a/EZFIO.tar.gz and /dev/null differ
diff --git a/EZFIO.tar.gz b/EZFIO.tar.gz
new file mode 120000
index 0000000..16a9f1b
--- /dev/null
+++ b/EZFIO.tar.gz
@@ -0,0 +1 @@
+EZFIO.1.0.19.tar.gz
\ No newline at end of file
diff --git a/Makefile b/Makefile
index d99abaa..7e94492 100644
--- a/Makefile
+++ b/Makefile
@@ -4,6 +4,9 @@ include make.config
 
 all: bin/eplf
 
+EZFIO.tar.gz:
+	ln -s EZFIO.*.tar.gz EZFIO.tar.gz
+
 EZFIO/config/eplf.config: EZFIO.tar.gz
 	tar -zxf EZFIO.tar.gz
 	cd EZFIO ; IRPF90="`echo $(IRPF90) | cut -f 1`" ./configure
diff --git a/src/debug_eplf.irp.f b/src/debug_eplf.irp.f
deleted file mode 100644
index a02f679..0000000
--- a/src/debug_eplf.irp.f
+++ /dev/null
@@ -1,53 +0,0 @@
-program debug
- implicit none
- PROVIDE ao_prim_num_max
- integer :: i,j
- integer :: k
- print *,  ''
-!print *,  'Occupation numbers'
-!do k=1,mo_num
-!  print *, k, mo_occ(k)
-!enddo
- read(*,*) 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 :: ao_eplf_integral, ao_eplf_integral_numeric
- print *,  ''
- print *, 'EPLF gamma         : ', eplf_gamma
- print *, 'EPLF integral      :', ao_eplf_integral(i,j,eplf_gamma,point)
- print *, 'EPLF integral N    :', ao_eplf_integral_numeric(i,j,eplf_gamma,point)
- 
- print *, ''
- print *, 'EPLF grid Npoints  :', grid_x_num, grid_y_num, grid_z_num
- print *, 'EPLF grid step     :', grid_step(:)
- print *, 'EPLF grid origin   :', grid_origin(:)
-
-end
diff --git a/src/density.irp.f b/src/density.irp.f
index 4fc1d34..8b6964d 100644
--- a/src/density.irp.f
+++ b/src/density.irp.f
@@ -34,12 +34,12 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ]
  density_alpha_value_p = 0.
 
  do i=1,mo_closed_num
-  density_alpha_value_p += mo_value_p(i)**2
+  density_alpha_value_p += mo_value_p(i)*mo_value_p(i)
  enddo
 
  do j=1,elec_alpha_num-mo_closed_num
   do i=1,elec_alpha_num-mo_closed_num
-   density_alpha_value_p += mo_value_prod_p(i,j) * one_e_density_mo(i,j,1)
+   density_alpha_value_p += mo_value_prod_p(i+mo_closed_num,j+mo_closed_num) * one_e_density_mo(i,j,1)
   enddo
  enddo
 
@@ -54,12 +54,12 @@ BEGIN_PROVIDER [ real, density_beta_value_p ]
  density_beta_value_p = 0.
 
  do i=1,mo_closed_num
-  density_beta_value_p += mo_value_p(i)**2
+  density_beta_value_p += mo_value_p(i)*mo_value_p(i)
  enddo
 
  do j=1,elec_beta_num-mo_closed_num
   do i=1,elec_beta_num-mo_closed_num
-   density_beta_value_p += mo_value_prod_p(i,j) * one_e_density_mo(i,j,2)
+   density_beta_value_p += mo_value_prod_p(i+mo_closed_num,j+mo_closed_num) * one_e_density_mo(i,j,2)
   enddo
  enddo
 END_PROVIDER
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index ea244e3..9590890 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -292,7 +292,6 @@ double precision function ao_eplf_integral(i,j,gmma,center)
  integer, intent(in) :: i, j
  real, intent(in)    :: center(3)
  double precision, intent(in) :: gmma
-!DEC$ ATTRIBUTES FORCEINLINE
  integer :: p,q,k
  double precision :: integral
  double precision :: ao_eplf_integral_primitive_oneD
diff --git a/src/overlap.irp.f b/src/overlap.irp.f
index c51f37b..a7c4b59 100644
--- a/src/overlap.irp.f
+++ b/src/overlap.irp.f
@@ -16,6 +16,23 @@ BEGIN_PROVIDER [ double precision, ao_overlap_matrix, (ao_num,ao_num) ]
 END_PROVIDER
 
 
+BEGIN_PROVIDER [ double precision, mo_overlap_matrix, (mo_tot_num,mo_tot_num) ]
+ implicit none
+ BEGIN_DOC  
+! Overlap matrix between the Atomic Orbitals
+ END_DOC
+
+ integer :: i, j
+ double precision :: mo_overlap
+ do j=1,mo_tot_num
+  do i=1,j
+   mo_overlap_matrix(i,j) = mo_overlap(i,j)
+   mo_overlap_matrix(j,i) = mo_overlap_matrix(i,j)
+  enddo
+ enddo
+END_PROVIDER
+
+
 double precision function primitive_overlap_oneD_numeric(a,xa,i,b,xb,j)
  implicit none
  include 'constants.F'
@@ -255,3 +272,23 @@ double precision function ao_overlap(i,j)
  
 end function
 
+double precision function mo_overlap(i,j)
+ implicit none
+ integer, intent(in) :: i, j
+ double precision :: tmp
+ integer :: p,q
+
+ mo_overlap = 0.d0
+ do q=1,ao_num
+  tmp = 0.d0
+  do p=1,ao_num
+   tmp += mo_coef(p,i)*ao_overlap_matrix(p,q)
+  enddo
+  mo_overlap += tmp*mo_coef(q,j)
+ enddo
+ 
+end function
+
+
+
+