diff --git a/bin/to_ezfio.exe b/bin/to_ezfio.exe
index 2d1d745..222b347 100755
Binary files a/bin/to_ezfio.exe and b/bin/to_ezfio.exe differ
diff --git a/src/det.irp.f b/src/det.irp.f
index 967cc7f..0276d4f 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -32,80 +32,70 @@ BEGIN_PROVIDER  [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ]
 
 END_PROVIDER
 
-
-BEGIN_PROVIDER [ integer*1, det_exc, (det_num, det_num, 2) ]
+BEGIN_PROVIDER [ real, mo_occ, (mo_tot_num) ]
  implicit none
- BEGIN_DOC  
-! Degree of excitation between two determinants. Indices are alpha, beta
-! The sign is the phase factor
+ BEGIN_DOC
+! Occupation numbers of molecular orbitals
  END_DOC
 
- integer :: p
- do p=1,2
+ call get_mo_basis_mo_occ(mo_occ)
 
-  integer :: k, l
-  do l=1,det_num
-   det_exc(l,l,p) = 0
-   do k=l+1,det_num
-    det_exc(k,l,p) = 0
+END_PROVIDER
 
-    ! Excitation degree
-    integer :: i, j
-    do i=1,elec_num_2(p)-mo_closed_num
 
-     logical :: found
-     found = .False.
-     do j=1,elec_num_2(p)-mo_closed_num
-      if (det(j,l,p) == det(i,k,p)) then
-        found = .True.
-        exit
-      endif
-     enddo
-     if (.not.found) then
-       det_exc(k,l,p) += 1
-     endif
 
-    enddo
+integer function det_exc(k,l,p)
+ implicit none
+! Degree of excitation between two determinants. Indices are alpha, beta
+! The sign is the phase factor
 
-   enddo
+ integer :: k,l,p
+
+ integer :: i, j
+ det_exc = 0
+ do i=1,elec_num_2(p)-mo_closed_num
+
+  logical :: found
+  found = .False.
+  do j=1,elec_num_2(p)-mo_closed_num
+   if (det(j,l,p) == det(i,k,p)) then
+     found = .True.
+!    exit
+   endif
   enddo
+  if (.not.found) then
+    det_exc += 1
+  endif
 
  enddo
 
  ! Phase
 
- do l=1,det_num
-  do k=l+1,det_num
-   integer :: nperm
-   nperm = 0
-   do p=1,2
-    integer :: buffer(0:mo_num-mo_closed_num)
-    do i=1,elec_num_2(p)-mo_closed_num
-     buffer(i) = det(i,k,p)
-    enddo
-    do i=1,elec_num_2(p)-mo_closed_num
-     if (buffer(i) /= det(i,l,p)) then
-      integer :: m
-      m=elec_num_2(p)-mo_closed_num
-      do j=i+1,elec_num_2(p)-mo_closed_num
-       if (buffer(i) == det(j,l,p)) then  ! found
-        m=j
-        exit
-       endif 
-      enddo
-      buffer(0) = buffer(i)
-      buffer(i) = det(m,l,p)
-      buffer(m) = buffer(0)
-      nperm += 1
-     endif
-    enddo
-    det_exc(k,l,p) *= (1-2*mod( nperm, 2 ))
-    det_exc(l,k,p) = det_exc(k,l,p)
-   enddo
-  enddo
+ integer :: nperm
+ nperm = 0
+ integer :: buffer(0:mo_num-mo_closed_num)
+ do i=1,elec_num_2(p)-mo_closed_num
+  buffer(i) = det(i,k,p)
  enddo
+ do i=1,elec_num_2(p)-mo_closed_num
+  if (buffer(i) /= det(i,l,p)) then
+   integer :: m
+   m=elec_num_2(p)-mo_closed_num
+   do j=i+1,elec_num_2(p)-mo_closed_num
+    if (buffer(i) == det(j,l,p)) then  ! found
+     m=j
+     exit
+    endif 
+   enddo
+   buffer(0) = buffer(i)
+   buffer(i) = det(m,l,p)
+   buffer(m) = buffer(0)
+   nperm += 1
+  endif
+ enddo
+ det_exc *= (1-2*mod( nperm, 2 ))
 
-END_PROVIDER
+end
 
 subroutine get_single_excitation(k,l,m,n,p)
  implicit none
@@ -198,3 +188,70 @@ subroutine get_double_excitation(k,l,m,n,r,s,p)
 
 end
 
+BEGIN_PROVIDER [ real, ci_mo, (mo_num,mo_num,3) ]
+ implicit none
+ BEGIN_DOC  
+! Spin Density matrix in the AO basis
+ END_DOC
+ integer :: i,j,k,l,m,ispin, ik,il
+ do ispin=1,3
+
+  do j=1,mo_num
+   do i=1,mo_num
+    ci_mo(i,j,ispin) = 0.
+   enddo
+  enddo
+
+  do l=1,det_num
+   do m=1,det_num
+    real :: factor
+    factor = 2.*det_coef(l)*det_coef(m) 
+    do il=1,mo_closed_num
+     do ik=1,mo_closed_num
+      ci_mo(ik,il,ispin) += factor
+     enddo
+    enddo
+   enddo
+  enddo
+
+ enddo
+
+
+ do l=1,det_num
+  do m=1,det_num
+   factor = det_coef(l)*det_coef(m) 
+   do ispin=1,2
+    do j=mo_closed_num+1,elec_num_2(ispin)
+     ik = det(j-mo_closed_num,l,ispin)
+     do il=1,mo_closed_num
+      ci_mo(ik,il,ispin) += factor
+      ci_mo(il,ik,ispin) += factor
+     enddo
+     do i=mo_closed_num+1,elec_num_2(ispin)
+      il = det(i-mo_closed_num,m,ispin)
+      ci_mo(ik,il,ispin) += factor
+      ci_mo(il,ik,ispin) += factor
+     enddo
+    enddo
+   enddo
+
+   ispin=3
+   do j=mo_closed_num+1,elec_num_2(1)
+    ik = det(j-mo_closed_num,l,1)
+    do il=1,mo_closed_num
+     ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m) 
+     ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m) 
+    enddo
+    do i=mo_closed_num+1,elec_num_2(2)
+     il = det(i-mo_closed_num,m,2)
+     ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m) 
+     ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m) 
+    enddo
+   enddo
+
+  enddo
+ enddo
+
+END_PROVIDER
+
+
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index c3df9d1..28c9bac 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -68,6 +68,7 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
  enddo
 END_PROVIDER
 
+
  BEGIN_PROVIDER [ double precision, eplf_up_up ]
 &BEGIN_PROVIDER [ double precision, eplf_up_dn ]
  implicit none
@@ -84,40 +85,151 @@ END_PROVIDER
 
  do i=1,mo_closed_num
   do j=1,mo_closed_num
-    eplf_up_up += mo_value_p(i)* (  &
-      mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
-      mo_value_p(j)*mo_eplf_integral_matrix(j,i) )
-    eplf_up_dn += mo_value_p(i)*mo_value_p(i)* &
-        mo_eplf_integral_matrix(j,j)
+    double precision :: temp
+    temp = mo_value_prod_p(i,i)*mo_eplf_integral_matrix(j,j)
+    eplf_up_up += temp - mo_value_prod_p(j,i)*mo_eplf_integral_matrix(j,i) 
+    eplf_up_dn += temp
   enddo
  enddo
  eplf_up_up *= 2.d0
  eplf_up_dn *= 2.d0 
 
+ integer :: k,l,m
+
+ do m=1,eplf_factor_num_max
+  i=eplf_factor_indice(1,m)
+  j=eplf_factor_indice(2,m)
+  k=eplf_factor_indice(3,m)
+  l=eplf_factor_indice(4,m)
+  temp = mo_value_prod_p(i,j)*mo_eplf_integral_matrix(k,l)
+  eplf_up_up += eplf_factor_value(1,m)*temp
+  eplf_up_dn += eplf_factor_value(2,m)*temp
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, eplf_factor_num_max ]
+ implicit none
+ BEGIN_DOC  
+! Number of factors containing the Slater rules 
+ END_DOC
+
+ eplf_factor_num_max = 0
+
+ integer :: k,l
+ integer :: exc(3), nact, nact2, p, p2
+ integer :: det_exc
+ do k=1,det_num
+  do l=k,det_num
+   exc(1) = det_exc(k,l,1)
+   exc(2) = det_exc(k,l,2)
+   exc(4) = exc(1)*exc(2)
+   exc(1) = abs(exc(1))
+   exc(2) = abs(exc(2))
+   exc(3) = exc(1)+exc(2)
+
+   do p=1,2
+     p2 = 1+mod(p,2)
+     nact  = elec_num_2(p) -mo_closed_num
+     nact2 = elec_num_2(p2)-mo_closed_num
+     if ( exc(3) == 0 ) then
+      eplf_factor_num_max += 2*nact*mo_num
+     else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
+      eplf_factor_num_max += 2*mo_num
+     else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
+      eplf_factor_num_max += 2
+     else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
+      eplf_factor_num_max += 1
+     endif
+   enddo
+
+  enddo
+ enddo
+ 
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, eplf_factor_indice, (4,eplf_factor_num_max) ]
+&BEGIN_PROVIDER [ real, eplf_factor_value, (2,eplf_factor_num_max) ]
+ implicit none
+ BEGIN_DOC  
+! Compact representation of eplf factors
+ END_DOC
+
+ integer :: i,j,k,l,m
+
+ m=1
+ do i=1,mo_num
+  do j=1,mo_num
+   do k=1,mo_num
+    do l=1,mo_num
+     if ( (eplf_factor(1,l,k,j,i) /= 0.).or. &
+          (eplf_factor(2,l,k,j,i) /= 0.) ) then
+       eplf_factor_indice(1,m) = l
+       eplf_factor_indice(2,m) = k
+       eplf_factor_indice(3,m) = j
+       eplf_factor_indice(4,m) = i
+       eplf_factor_value(1,m) = eplf_factor(1,l,k,j,i)
+       eplf_factor_value(2,m) = eplf_factor(2,l,k,j,i)
+       m += 1
+     endif
+    enddo
+   enddo
+  enddo
+ enddo
+ FREE eplf_factor
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ real, eplf_factor, (2,mo_num,mo_num,mo_num,mo_num) ]
+ implicit none
+ BEGIN_DOC  
+! Factors containing the Slater rules 
+ END_DOC
+
+ integer :: i, j
+
  integer :: k,l,m,n,p,p2
  integer :: ik,il,jk,jl
- double precision :: phase,dtemp(2)
+ real :: phase
  integer :: exc(4), nact, nact2
+ real :: det_kl
+ integer :: det_exc
+
+ do m=1,2
+  do i=1,mo_num
+   do j=1,mo_num
+    do k=1,mo_num
+     do l=1,mo_num
+      eplf_factor(m,l,k,j,i) = 0.
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
 
  PROVIDE det
- PROVIDE elec_num_2
- PROVIDE mo_value_prod_p
 
  do k=1,det_num
   do l=k,det_num
 
-   exc(1) = abs(det_exc(k,l,1))
-   exc(2) = abs(det_exc(k,l,2))
+   exc(1) = det_exc(k,l,1)
+   exc(2) = det_exc(k,l,2)
+   exc(4) = exc(1)*exc(2)
+   exc(1) = abs(exc(1))
+   exc(2) = abs(exc(2))
    exc(3) = exc(1)+exc(2)
-   exc(4) = det_exc(k,l,1)*det_exc(k,l,2)
    if (exc(4) /= 0) then
     exc(4) = exc(4)/abs(exc(4))
    else
     exc(4) = 1
    endif
+   phase = dble(exc(4))
+
+   det_kl = phase*det_coef(k)*det_coef(l)
+   if (k /= l) then
+    det_kl += det_kl
+   endif
 
-   dtemp(1) = 0.d0
-   dtemp(2) = 0.d0
    do p=1,2
      p2 = 1+mod(p,2)
      nact  = elec_num_2(p) -mo_closed_num
@@ -129,30 +241,27 @@ END_PROVIDER
 
        do i=1,mo_closed_num
         ! Closed-open shell interactions
-        dtemp(1) += (  &
-           mo_value_prod_p(jl,jk)*mo_eplf_integral_matrix(i,i) - &
-           mo_value_prod_p(i,jk)*mo_eplf_integral_matrix(i,jl) )
-        dtemp(2) += mo_value_prod_p(jl,jk)*mo_eplf_integral_matrix(i,i) 
+        eplf_factor(1,jk,jl,i,i) += det_kl
+        eplf_factor(2,jk,jl,i,i) += det_kl
+        eplf_factor(1,i,jl,jk,i) -= det_kl
 
         !- Open-closed shell interactions
-        dtemp(1) += (  &
-           mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk) - &
-           mo_value_prod_p(i,jl)*mo_eplf_integral_matrix(i,jk) )  
-        dtemp(2) += mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk) 
+        eplf_factor(1,i,i,jk,jl) += det_kl
+        eplf_factor(2,i,i,jk,jl) += det_kl
+        eplf_factor(1,jk,i,i,jl) -= det_kl
        enddo
 
        !- Open-open shell interactions
        do m=1,nact
          ik = det(m,k,p)
          il = det(m,l,p)
-         dtemp(1) += (  &
-           mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(jl,jk) - &
-           mo_value_prod_p(jl,ik)*mo_eplf_integral_matrix(il,jk) )
+         eplf_factor(1,ik,il,jk,jl) += det_kl
+         eplf_factor(1,jk,il,ik,jl) -= det_kl
        enddo
        do m=1,nact2
          ik = det(m,k,p2)
          il = det(m,l,p2)
-         dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(jl,jk) 
+         eplf_factor(2,ik,il,jk,jl) += det_kl
        enddo
 
       enddo
@@ -164,60 +273,46 @@ END_PROVIDER
 
       do i=1,mo_closed_num
        !- Open-closed shell interactions
-       dtemp(1) += (  &
-           mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(i,i) - &
-           mo_value_prod_p(i,ik)*mo_eplf_integral_matrix(i,il) )
-       dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(i,i)
+       eplf_factor(1,ik,il,i,i) += det_kl
+       eplf_factor(2,ik,il,i,i) += det_kl
+       eplf_factor(1,i,il,ik,i) -= det_kl
 
        !- Closed-open shell interactions
-       dtemp(1) += (  &
-           mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk) - &
-           mo_value_prod_p(i,jl)*mo_eplf_integral_matrix(i,jk) )  
-       dtemp(2) += mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk)
+       eplf_factor(1,i,i,jk,jl) += det_kl
+       eplf_factor(2,i,i,jk,jl) += det_kl
+       eplf_factor(1,jk,i,i,jl) -= det_kl
       enddo
 
       !- Open-open shell interactions
       do m=1,nact
        jk = det(m,k,p)
        jl = det(m,l,p)
-       dtemp(1) += (  &
-           mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(jl,jk) - &
-           mo_value_prod_p(jl,ik)*mo_eplf_integral_matrix(il,jk) )
+       eplf_factor(1,ik,il,jk,jl) += det_kl
+       eplf_factor(1,jk,il,ik,jl) -= det_kl
       enddo
       do m=1,nact2
        jk = det(m,k,p2)
        jl = det(m,l,p2)
-       dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(jl,jk) 
+       eplf_factor(2,ik,il,jk,jl) += det_kl
       enddo
 
      else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
     
      ! Consider only the double excitations of same-spin electrons
       call get_double_excitation(k,l,ik,il,jk,jl,p)
-    
-      dtemp(1) += (  &
-             mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(jl,jk) - &
-             mo_value_prod_p(jl,ik)*mo_eplf_integral_matrix(il,jk) )
+      eplf_factor(1,ik,il,jk,jl) += det_kl
+      eplf_factor(1,jk,il,ik,jl) -= det_kl
     
      else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
     
-     ! Consider only the double excitations of opposite-spin electrons
-     call get_single_excitation(k,l,ik,il,p)
-     call get_single_excitation(k,l,jk,jl,p2)
-    
-      dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(jl,jk)
+      ! Consider only the double excitations of opposite-spin electrons
+      call get_single_excitation(k,l,ik,il,p)
+      call get_single_excitation(k,l,jk,jl,p2)
+      eplf_factor(2,ik,il,jk,jl) += det_kl
     
      endif
    enddo
 
-   phase = dble(exc(4))
-   eplf_up_up += phase * det_coef(k)*det_coef(l) * dtemp(1) 
-   eplf_up_dn += phase * det_coef(k)*det_coef(l) * dtemp(2)
-   if (k /= l) then
-    eplf_up_up += phase * det_coef(k)*det_coef(l) * dtemp(1)
-    eplf_up_dn += phase * det_coef(k)*det_coef(l) * dtemp(2)
-   endif
-
   enddo
  enddo
  
diff --git a/src/mo.irp.f b/src/mo.irp.f
index a83c022..4b6fac9 100644
--- a/src/mo.irp.f
+++ b/src/mo.irp.f
@@ -43,16 +43,6 @@ BEGIN_PROVIDER [ integer, mo_num ]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ real, mo_occ, (mo_tot_num) ]
- implicit none
- BEGIN_DOC
-! Occupation numbers of molecular orbitals
- END_DOC
-
- call get_mo_basis_mo_occ(mo_occ)
-
-END_PROVIDER
-
 BEGIN_PROVIDER [ real, mo_coef, (ao_num,mo_num) ]
  implicit none
  BEGIN_DOC
diff --git a/src/overlap.irp.f b/src/overlap.irp.f
index 3027a31..c51f37b 100644
--- a/src/overlap.irp.f
+++ b/src/overlap.irp.f
@@ -155,10 +155,10 @@ double precision function primitive_overlap_oneD(a,xa,i,b,xb,j)
  xp = xp*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
+ if ( c > 32.d0 ) then   ! Cut-off on exp(-32)
+  primitive_overlap_oneD = 0.d0
+  return
+ endif
 
  c =  dexp(-c)
 
diff --git a/src/test_1d.irp.f b/src/test_1d.irp.f
index 79e3112..ef845f6 100644
--- a/src/test_1d.irp.f
+++ b/src/test_1d.irp.f
@@ -10,7 +10,7 @@ subroutine run
  point(1) = 0.
  point(2) = 0.
  integer :: i
- do i=-40,40
+ do i=-40,60
   point(3) = real(i)/20.
   TOUCH point
   print *,  point(3), eplf_value_p, eplf_up_up, eplf_up_dn