diff --git a/EZFIO.tar.gz b/EZFIO.tar.gz
index 8efc885..7075609 100644
Binary files a/EZFIO.tar.gz and b/EZFIO.tar.gz differ
diff --git a/configure b/configure
index 49ce85f..ba415a7 100755
--- a/configure
+++ b/configure
@@ -3248,7 +3248,7 @@ echo "export EPLF_HAS_MPI=$EPLF_HAS_MPI" >> $HOME/.eplfrc
 echo "export EPLF_MPIRUN=$MPIRUN" >> $HOME/.eplfrc
 echo "************************************"
 echo "To finish the installation:"
-echo "1) Add the following line to your $HOME/.basrhc file:"
+echo "1) Add the following line to your $HOME/.bashrc file:"
 echo ". $HOME/.eplfrc"
 echo "2) Execute"
 echo ". $HOME/.eplfrc"
diff --git a/configure.ac b/configure.ac
index bd8e382..8d28a8d 100644
--- a/configure.ac
+++ b/configure.ac
@@ -211,7 +211,7 @@ echo "export EPLF_HAS_MPI=$EPLF_HAS_MPI" >> $HOME/.eplfrc
 echo "export EPLF_MPIRUN=$MPIRUN" >> $HOME/.eplfrc
 echo "************************************"
 echo "To finish the installation:"
-echo "1) Add the following line to your $HOME/.basrhc file:"
+echo "1) Add the following line to your $HOME/.bashrc file:"
 echo ". $HOME/.eplfrc"
 echo "2) Execute"
 echo ". $HOME/.eplfrc"
diff --git a/src/density.irp.f b/src/density.irp.f
index 45ed3e0..394d738 100644
--- a/src/density.irp.f
+++ b/src/density.irp.f
@@ -37,40 +37,10 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ]
   density_alpha_value_p += mo_value_p(i)**2
  enddo
 
-! TODO vectorization
- integer :: k,j,l, ik, il
- real    :: buffer
- real    :: phase
- integer :: exc(4)
- PROVIDE det
- PROVIDE elec_alpha_num
- do k=1,det_num
-  do l=1,det_num
-
-   exc(1) = abs(det_exc(k,l,1))
-   exc(2) = abs(det_exc(k,l,2))
-   exc(3) = exc(1)+exc(2)
-   exc(4) = exc(1)*exc(2)
-   if (exc(4) /= 0) then
-    exc(4) = exc(4)/abs(exc(4))
-   else
-    exc(4) = 1
-   endif
-
-   phase = dble(exc(4))
-   if (exc(3) == 0) then
-     buffer = 0.
-     do i=1,elec_alpha_num-mo_closed_num
-      buffer += mo_value_p(det(i,k,1))*mo_value_p(det(i,l,1))
-     enddo
-     density_alpha_value_p += phase*det_coef(k)*det_coef(l)*buffer
-   else if ( (exc(3) == 1).and.(exc(1) == 1) ) then
-     call get_single_excitation(k,l,ik,il,1)
-     buffer = mo_value_p(ik)*mo_value_p(il)
-     density_alpha_value_p += phase*det_coef(k)*det_coef(l)*buffer
-   endif
-
-  enddo
+ do m=1,one_e_density_num
+   i = one_e_density_indice(1,m)
+   j = one_e_density_indice(2,m)
+   density_alpha_value_p += mo_value_prod_p(i,j) * one_e_density_value(1,m)
  enddo
 
 END_PROVIDER
@@ -87,39 +57,10 @@ BEGIN_PROVIDER [ real, density_beta_value_p ]
   density_beta_value_p += mo_value_p(i)**2
  enddo
 
-! TODO vectorization
- integer :: k,j,l, ik, il
- real    :: buffer
- real    :: phase
- integer :: exc(4)
- PROVIDE det
- PROVIDE elec_beta_num
-
- do k=1,det_num
-  do l=1,det_num
-   exc(1) = abs(det_exc(k,l,1))
-   exc(2) = abs(det_exc(k,l,2))
-   exc(3) = exc(1)+exc(2)
-   exc(4) = exc(1)*exc(2)
-   if (exc(4) /= 0) then
-    exc(4) = exc(4)/abs(exc(4))
-   else
-    exc(4) = 1
-   endif
-
-   phase = dble(exc(4))
-   if (exc(3) == 0) then
-     buffer = 0.
-     do i=1,elec_beta_num-mo_closed_num
-      buffer += mo_value_p(det(i,k,2))*mo_value_p(det(i,l,2))
-     enddo
-     density_beta_value_p += phase*det_coef(k)*det_coef(l)*buffer
-   else if ( (exc(3) == 1).and.(exc(2) == 1) ) then
-     call get_single_excitation(k,l,ik,il,2)
-     buffer = mo_value_p(ik)*mo_value_p(il)
-     density_beta_value_p += phase*det_coef(k)*det_coef(l)*buffer
-   endif
-  enddo
+ do m=1,one_e_density_num
+   i = one_e_density_indice(1,m)
+   j = one_e_density_indice(2,m)
+   density_beta_value_p += mo_value_prod_p(i,j) * one_e_density_value(2,m)
  enddo
 
 END_PROVIDER
diff --git a/src/det.irp.f b/src/det.irp.f
index 6634f3e..07d45cb 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -46,21 +46,22 @@ END_PROVIDER
 
 integer function det_exc(k,l,p)
  implicit none
-! Degree of excitation between two determinants. Indices are alpha, beta
+! Degree of excitation+1 between two determinants. Indices are alpha, beta
 ! The sign is the phase factor
 
  integer :: k,l,p
+ integer :: i, j, jmax
 
- integer :: i, j
+ jmax = elec_num_2(p)-mo_closed_num
  det_exc = 0
- do i=1,elec_num_2(p)-mo_closed_num
 
+ do i=1,jmax
   logical :: found
   found = .False.
-  do j=1,elec_num_2(p)-mo_closed_num
+  do j=1,jmax
    if (det(j,l,p) == det(i,k,p)) then
      found = .True.
-!    exit
+     exit
    endif
   enddo
   if (.not.found) then
@@ -73,7 +74,10 @@ integer function det_exc(k,l,p)
 
  integer :: nperm
  nperm = 0
- integer :: buffer(0:mo_num-mo_closed_num)
+ integer, allocatable, save :: buffer(:)
+ if (.not. allocated(buffer)) then
+   allocate (buffer(0:mo_num-mo_closed_num))
+ endif
  do i=1,elec_num_2(p)-mo_closed_num
   buffer(i) = det(i,k,p)
  enddo
@@ -90,9 +94,10 @@ integer function det_exc(k,l,p)
    buffer(0) = buffer(i)
    buffer(i) = det(m,l,p)
    buffer(m) = buffer(0)
-   nperm += 1
+   nperm += m-i
   endif
  enddo
+ det_exc += 1
  det_exc *= (1-2*mod( nperm, 2 ))
 
 end
@@ -268,8 +273,8 @@ BEGIN_PROVIDER [ integer, two_e_density_num_max ]
  integer :: det_exc
  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) = abs(det_exc(k,l,1))-1
+   exc(2) = abs(det_exc(k,l,2))-1
    exc(3) = exc(1)+exc(2)
 
    do p=1,2
@@ -319,14 +324,10 @@ END_PROVIDER
    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(1) = abs(exc(1))-1
+   exc(2) = abs(exc(2))-1
    exc(3) = exc(1)+exc(2)
-   if (exc(4) /= 0) then
-    exc(4) = exc(4)/abs(exc(4))
-   else
-    exc(4) = 1
-   endif
+   exc(4) = exc(4)/abs(exc(4))
    phase = dble(exc(4))
 
    det_kl = phase*det_coef(k)*det_coef(l)
@@ -510,3 +511,47 @@ END_SHELL
  
 END_PROVIDER
 
+BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ]
+ implicit none
+ BEGIN_DOC  
+! One electron spin density matrix in MO space
+ END_DOC
+ integer :: i,j,k,l,p, il, jl
+ do p=1,2
+  do i=1,mo_active_num
+   do j=1,mo_active_num
+    one_e_density_mo(j,i,p) = 0.
+   enddo
+  enddo
+ enddo
+ real :: ckl, phase
+ integer :: exc(4), 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))-1
+   exc(2) = abs(exc(2))-1
+   exc(3) = exc(1)+exc(2)
+   exc(4) = exc(4)/abs(exc(4))
+   phase = dble(exc(4))
+   ckl = det_coef(k)*det_coef(l)*phase
+   do p=1,2
+    if (exc(3) == 0) then
+      do i=1,elec_num_2(p)-mo_closed_num
+       il = det(i,k,p) - mo_closed_num
+       one_e_density_mo(il,il,p) += ckl
+      enddo
+    else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
+      call get_single_excitation(k,l,il,jl,p)
+      jl -= mo_closed_num
+      il -= mo_closed_num
+      one_e_density_mo(il,jl,p) += ckl
+      one_e_density_mo(jl,il,p) += ckl
+    endif
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index 83c7894..242faa7 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -350,10 +350,3 @@ double precision function mo_eplf_integral(i,j)
  
 end function
 
-
-
-
-
-
-
-
diff --git a/version b/version
index e2e9306..d35cc4e 100644
--- a/version
+++ b/version
@@ -1 +1 @@
-VERSION=1.0.4
+VERSION=1.0.5