diff --git a/src/density.irp.f b/src/density.irp.f
index b8276a4..5cac985 100644
--- a/src/density.irp.f
+++ b/src/density.irp.f
@@ -38,15 +38,25 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ]
  enddo
 
 ! TODO vectorization
- integer :: k,j,l
+ integer :: k,j,l, ik, il
  real    :: buffer
+ PROVIDE det
+ PROVIDE elec_alpha_num
  do k=1,det_num
   do l=1,det_num
-   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 += det_coef(k)*det_coef(l)*buffer
+
+   if (det_exc(k,l,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 += det_coef(k)*det_coef(l)*buffer
+   else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,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 += det_coef(k)*det_coef(l)*buffer
+   endif
+
   enddo
  enddo
 
@@ -65,15 +75,23 @@ BEGIN_PROVIDER [ real, density_beta_value_p ]
  enddo
 
 ! TODO vectorization
- integer :: k,j,l
+ integer :: k,j,l, ik, il
  real    :: buffer
+ PROVIDE det
+ PROVIDE elec_beta_num
  do k=1,det_num
   do l=1,det_num
-   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 += det_coef(k)*det_coef(l)*buffer
+   if (det_exc(k,l,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 += det_coef(k)*det_coef(l)*buffer
+   else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,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 += det_coef(k)*det_coef(l)*buffer
+   endif
   enddo
  enddo
 
@@ -98,9 +116,12 @@ BEGIN_PROVIDER [ double precision, density_alpha_grad_p, (3) ]
    density_alpha_grad_p(3) += 2.*mo_grad_p(i,3)*mo_value_p(i)
  enddo
 
-! TODO vectorization
+! TODO Faux !! => Regles de Slater
  integer :: k,j,l
  real    :: buffer(3)
+ if (det_num > 1) then
+   call abrt(irp_here,'Need to implement Slater rules for gradient')
+ endif
  do k=1,det_num
   do l=1,det_num
    buffer(1) = 0.
@@ -142,6 +163,9 @@ BEGIN_PROVIDER [ double precision, density_beta_grad_p, (3) ]
  enddo
 
 ! TODO vectorization
+ if (det_num > 1) then
+   call abrt(irp_here,'Need to implement Slater rules for gradient')
+ endif
  integer :: k,j,l
  real    :: buffer(3)
  do k=1,det_num
@@ -189,6 +213,9 @@ BEGIN_PROVIDER [ double precision, density_alpha_lapl_p ]
 ! TODO vectorization
  integer :: k,j,l
  real    :: buffer
+ if (det_num > 1) then
+   call abrt(irp_here,'Need to implement Slater rules for gradient')
+ endif
  do k=1,det_num
   do l=1,det_num
    buffer = 0.
@@ -225,6 +252,9 @@ BEGIN_PROVIDER [ double precision, density_beta_lapl_p ]
  density_beta_lapl_p *= 2.
 
 ! TODO vectorization
+ if (det_num > 1) then
+   call abrt(irp_here,'Need to implement Slater rules for laplacian')
+ endif
  integer :: k,j,l
  real    :: buffer
  do k=1,det_num
diff --git a/src/det.irp.f b/src/det.irp.f
index ae61fba..422bfbb 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -33,4 +33,168 @@ BEGIN_PROVIDER  [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ]
 END_PROVIDER
 
 
+BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 3) ]
+ implicit none
+ BEGIN_DOC  
+! Degree of excitation between two determinants. The sign is the phase.
+ END_DOC
+
+ integer :: p
+ do p=1,2
+
+  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
+
+    ! 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
+
+    det_exc(l,k,p) = det_exc(k,l,p)
+   enddo
+  enddo
+
+ enddo
+
+ do l=1,det_num
+  do k=l+1,det_num
+   det_exc(k,l,3) = det_exc(k,l,1) + det_exc(k,l,2)
+  enddo
+ enddo
+
+ ! Phase
+ do p=1,2
+  do i=mo_closed_num,mo_num
+   integer :: det_pos(det_num)
+   do k=1,det_num
+    det_pos(k) = 0
+    do j=1,elec_num_2(p)-mo_closed_num
+     if (det(j,k,p) == i) then
+      det_pos(k) = j
+     endif
+    enddo
+   enddo
+   do k=1,det_num
+    do l=k+1,det_num
+     det_exc(k,l,3) *= -2*mod( (det_pos(k)+det_pos(l)), 2 )+1
+    enddo
+   enddo
+  enddo
+
+ enddo
+
+ do l=1,det_num
+  do k=l+1,det_num
+   det_exc(l,k,3) = det_exc(k,l,3)
+  enddo
+ enddo
+
+END_PROVIDER
+
+subroutine get_single_excitation(k,l,m,n,p)
+ implicit none
+ integer, intent(in)  :: k, l  ! determinant indices
+ integer, intent(out) :: m, n  ! m->n excitation
+ integer, intent(in)  :: p     ! spin
+ 
+ logical :: found
+ integer :: i,j
+ m=0
+ n=0
+ do j=1,elec_num_2(p)-mo_closed_num
+  found = .False.
+  do i=1,elec_num_2(p)-mo_closed_num
+   if (det(j,k,p) == det(i,l,p)) then
+    found = .True.
+    exit
+   endif
+  enddo
+  if (.not.found) then
+    m = det(j,k,p)
+    exit
+  endif
+ enddo
+
+ do i=1,elec_num_2(p)-mo_closed_num
+  found = .False.
+  do j=1,elec_num_2(p)-mo_closed_num
+   if (det(i,k,p) == det(i,l,p)) then
+    found = .True.
+    exit
+   endif
+  enddo
+  if (.not.found) then
+    n = det(i,l,p)
+    exit
+  endif
+ enddo
+
+end
+
+subroutine get_double_excitation(k,l,m,n,r,s,p)
+ implicit none
+ integer, intent(in)  :: k, l  ! determinant indices
+ integer, intent(out) :: m, n  ! m->n excitation
+ integer, intent(out) :: r, s  ! r->s excitation
+ integer, intent(in)  :: p     ! spin
+ 
+ logical :: found
+ integer :: i,j
+ m=0
+ n=0
+ r=0
+ s=0
+ do j=1,elec_num_2(p)-mo_closed_num
+  found = .False.
+  do i=1,elec_num_2(p)-mo_closed_num
+   if (det(j,k,p) == det(i,l,p)) then
+    found = .True.
+    exit
+   endif
+  enddo
+  if (.not.found) then
+   if (m == 0) then
+    m = det(j,k,p)
+   else
+    r = det(j,k,p)
+    exit
+   endif
+  endif
+ enddo
+
+ do i=1,elec_num_2(p)-mo_closed_num
+  found = .False.
+  do j=1,elec_num_2(p)-mo_closed_num
+   if (det(i,k,p) == det(i,l,p)) then
+    found = .True.
+    exit
+   endif
+  enddo
+  if (.not.found) then
+   if (n == 0) then
+    n = det(i,l,p)
+   else
+    s = det(i,l,p)
+    exit
+   endif
+  endif
+ enddo
+
+end
 
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index 9d03028..3538b89 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -11,6 +11,7 @@ BEGIN_PROVIDER  [ real, eplf_gamma ]
 !eplf_gamma = (4./(3.*N)*pi*density_value_p)**(2./3.) * eps
  eplf_gamma = density_value_p * eps
 !eplf_gamma = 1.e10
+!eplf_gamma = 1.e5
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
@@ -79,41 +80,116 @@ END_PROVIDER
 
  do j=1,mo_closed_num
   do i=1,mo_closed_num
-    eplf_up_up += 2.d0*mo_value_p(i)* (  &
+    eplf_up_up += 2.d0* mo_value_p(i)* (  &
       mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
       mo_value_p(j)*mo_eplf_integral_matrix(i,j) )
-    eplf_up_dn += 2.d0*mo_value_p(j)**2 * &
-        mo_eplf_integral_matrix(i,i)
+    eplf_up_dn += 2.d0* mo_value_p(i)*mo_value_p(i)* &
+        mo_eplf_integral_matrix(j,j)
   enddo
  enddo
 
  integer :: k,l,m,n,p
+ integer :: ik,il,jk,jl
  double precision :: ckl
+ double precision :: phase
+ integer :: exc
+
+ PROVIDE det
+ PROVIDE elec_num_2
+
  do k=1,det_num
   do l=1,det_num
 
    ckl = det_coef(k)*det_coef(l)
 
-   do p=1,2
-    do m=1,elec_num_2(p)-mo_closed_num
-     j = det(m,k,p)
-     do n=1,elec_num_2(p)-mo_closed_num
-       i = det(n,l,p)
-       eplf_up_up += ckl*mo_value_p(i)* (  &
-         mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
-         mo_value_p(j)*mo_eplf_integral_matrix(i,j) )
-     enddo
-    enddo
-   enddo
+   exc = det_exc(k,l,3)
 
-   do m=1,elec_beta_num-mo_closed_num
-    j = det(m,k,2)
-    do n=1,elec_alpha_num-mo_closed_num
-      i = det(n,k,1)
-      eplf_up_dn += ckl * ( mo_value_p(i)**2 * mo_eplf_integral_matrix(j,j) &
-                          + mo_value_p(j)**2 * mo_eplf_integral_matrix(i,i) )
-    enddo
-   enddo
+   if ( exc < 0 ) then
+    phase = -0.5d0
+    exc = -exc
+   else
+    phase = 0.5d0
+   endif
+
+   if ( exc == 0 ) then
+     ! Sum over all alpha-alpha and beta-beta interactions
+     do p=1,2
+      do m=1,elec_num_2(p)-mo_closed_num
+       jk = det(m,k,p)
+       jl = det(m,l,p)
+       do n=1,elec_num_2(p)-mo_closed_num
+         ik = det(n,k,p)
+         il = det(n,l,p)
+         eplf_up_up += phase*ckl*mo_value_p(ik)* (  &
+           mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
+           mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
+       enddo
+      enddo
+     enddo
+
+     ! Sum over all alpha-beta interactions
+     do m=1,elec_beta_num-mo_closed_num
+      jk = det(m,k,2)
+      jl = det(m,l,2)
+      do n=1,elec_alpha_num-mo_closed_num
+        ik = det(n,k,1)
+        il = det(n,l,1)
+        eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) &
+                            + mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) )
+      enddo
+     enddo
+
+   else if ( exc == 1 ) then
+    
+     do p=1,2
+      if ( det_exc(k,l,p) == 1 ) then
+      ! Sum over only the sigma-sigma interactions involving the excitation
+       call get_single_excitation(k,l,ik,il,p)
+       do m=1,elec_num_2(p)-mo_closed_num
+        jk = det(m,k,p)
+        jl = det(m,l,p)
+        eplf_up_up += phase*ckl*mo_value_p(ik)* (  &
+            mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
+            mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
+       enddo
+
+       ! Sum over only the sigma-(sigma_bar) interactions involving the excitation
+       integer :: p2
+       p2 = 1+mod(p,2)
+       do m=1,elec_num_2(p2)-mo_closed_num
+        jk = det(m,k,p2)
+        jl = det(m,l,p2)
+        eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) &
+                            + mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) )
+       enddo
+      endif
+     enddo
+
+   else if (exc == 2) then
+
+     if ( ( det_exc(k,l,1) == 2 ).or.( det_exc(k,l,2) == 2 ) ) then
+     
+      ! Consider only the double excitations of same-spin electrons
+      if ( det_exc(k,l,1) == 2 ) then
+       call get_double_excitation(k,l,ik,jk,il,jl,1)
+      else if ( det_exc(k,l,2) == 2 ) then
+       call get_double_excitation(k,l,ik,jk,il,jl,2)
+      endif
+
+      eplf_up_up += phase*ckl*mo_value_p(ik)* (  &
+           mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
+           mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
+
+     else if ( det_exc(k,l,1) == 1 ) then
+      ! Consider only the double excitations of opposite-spin electrons
+      call get_single_excitation(k,l,ik,jk,1)
+      call get_single_excitation(k,l,il,jl,2)
+
+      eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) &
+                  + mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) )
+     endif
+
+   endif
 
   enddo
  enddo
diff --git a/src/grid.irp.f b/src/grid.irp.f
index 6c33399..51bdbe9 100644
--- a/src/grid.irp.f
+++ b/src/grid.irp.f
@@ -100,7 +100,7 @@ BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ]
   do iy=1,grid_y_num 
    point(2) = grid_origin(2)+(iy-1)*grid_step(2)
    do ix=1,grid_x_num 
-    icount = icount-1
+    icount -= 1
     if (icount == mpi_rank) then
      point(1) = grid_origin(1)+(ix-1)*grid_step(1)
      TOUCH point
@@ -116,7 +116,7 @@ BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ]
  IRP_IF MPI
    integer :: dim, ierr
    do iz=1,grid_z_num 
-    real   :: buffer(grid_x_num*grid_y_num)
+    real   :: buffer(grid_x_num*grid_y_num+1)
     icount = 0
     do iy=1,grid_y_num 
      do ix=1,grid_x_num 
@@ -127,6 +127,10 @@ BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ]
     dim = grid_x_num * grid_y_num
     call MPI_REDUCE(buffer,grid_$X(1,1,iz),dim,mpi_real, &
        mpi_sum,0,MPI_COMM_WORLD,ierr)
+    if (ierr /= MPI_SUCCESS) then
+      call abrt(irp_here,'Unable to fetch buffer')
+    endif
+    call barrier
    enddo
  IRP_ENDIF