diff --git a/src/Determinants/two_body_dm_map.irp.f b/src/Determinants/two_body_dm_map.irp.f
index eccb9741..f570e2bf 100644
--- a/src/Determinants/two_body_dm_map.irp.f
+++ b/src/Determinants/two_body_dm_map.irp.f
@@ -8,6 +8,7 @@ BEGIN_PROVIDER [ type(map_type), two_body_dm_ab_map ]
   END_DOC
   integer(key_kind)              :: key_max
   integer(map_size_kind)         :: sze
+  use map_module
   call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max)
   sze = key_max
   call map_init(two_body_dm_ab_map,sze)
@@ -129,51 +130,56 @@ subroutine add_values_to_two_body_dm_map(mask_ijkl)
     call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
     if(degree>2)cycle
     call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
-    contrib = psi_coef(i,1) * psi_coef(j,1) * phase
     call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
     if(degree==2)then  ! case of the DOUBLE EXCITATIONS  ************************************
 
      if(s1==s2)cycle  ! Only the alpha/beta two body density matrix
      ! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J
+     if(h1>p1)cycle
+     if(h2>p2)cycle
+!    if(s1.ne.1)cycle
      n_elements += 1
+     contrib = psi_coef(i,1) * psi_coef(j,1) * phase
      buffer_value(n_elements) = contrib
      !DEC$ FORCEINLINE
+!    call mo_bielec_integrals_index(h1,p1,h2,p2,buffer_i(n_elements))
      call mo_bielec_integrals_index(h1,h2,p1,p2,buffer_i(n_elements))
-     if (n_elements == size_buffer) then
-       call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
-           real(mo_integrals_threshold,integral_kind))
-       n_elements = 0
-     endif
+!    if (n_elements == size_buffer) then
+!      call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
+!          real(mo_integrals_threshold,integral_kind))
+!      n_elements = 0
+!    endif
 
     else ! case of the SINGLE EXCITATIONS  ***************************************************
+     cycle
 
-     if(s1==1)then  ! Mono alpha : 
-      do k = 1, elec_beta_num  
-       m = occ(k,2)
-       n_elements += 1
-       buffer_value(n_elements) = contrib
-       ! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
-       call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements))
-       if (n_elements == size_buffer) then
-         call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
-             real(mo_integrals_threshold,integral_kind))
-         n_elements = 0
-       endif
-      enddo
-     else  ! Mono Beta : 
-      do k = 1, elec_alpha_num
-       m = occ(k,1)
-       n_elements += 1
-       buffer_value(n_elements) = contrib
-       ! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
-       call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements))
-       if (n_elements == size_buffer) then
-         call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
-             real(mo_integrals_threshold,integral_kind))
-         n_elements = 0
-       endif
-      enddo
-     endif
+!    if(s1==1)then  ! Mono alpha : 
+!     do k = 1, elec_beta_num  
+!      m = occ(k,2)
+!      n_elements += 1
+!      buffer_value(n_elements) = contrib
+!      ! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
+!      call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements))
+!      if (n_elements == size_buffer) then
+!        call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
+!            real(mo_integrals_threshold,integral_kind))
+!        n_elements = 0
+!      endif
+!     enddo
+!    else  ! Mono Beta : 
+!     do k = 1, elec_alpha_num
+!      m = occ(k,1)
+!      n_elements += 1
+!      buffer_value(n_elements) = contrib
+!      ! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
+!      call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements))
+!      if (n_elements == size_buffer) then
+!        call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
+!            real(mo_integrals_threshold,integral_kind))
+!        n_elements = 0
+!      endif
+!     enddo
+!    endif
 
     endif
    enddo
@@ -181,6 +187,9 @@ subroutine add_values_to_two_body_dm_map(mask_ijkl)
   print*,'n_elements = ',n_elements
   call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
       real(mo_integrals_threshold,integral_kind))
+  call map_unique(two_body_dm_ab_map)
+
+  deallocate(buffer_i,buffer_value)
 
 end
 
@@ -192,8 +201,8 @@ BEGIN_PROVIDER [double precision, two_body_dm_ab_diag, (mo_tot_num, mo_tot_num)]
  double precision :: contrib
  BEGIN_DOC 
  ! two_body_dm_ab_diag(k,m) = <\Psi | n_(k\alpha) n_(m\beta) | \Psi>
-
  END_DOC
+
  two_body_dm_ab_diag = 0.d0
   do i = 1, N_det ! i == |I>
    call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
@@ -202,9 +211,83 @@ BEGIN_PROVIDER [double precision, two_body_dm_ab_diag, (mo_tot_num, mo_tot_num)]
     k = occ(j,2)
     do l = 1, elec_beta_num  
      m = occ(l,1)
-     two_body_dm_ab_diag(k,m) += contrib
+     two_body_dm_ab_diag(k,m) += 0.5d0 * contrib
+     two_body_dm_ab_diag(m,k) += 0.5d0 * contrib
     enddo
    enddo
   enddo
 END_PROVIDER 
 
+BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
+ implicit none
+ integer :: i,j,k,l,m
+ integer                        :: degree
+ PROVIDE mo_coef psi_coef psi_det
+ integer                        :: exc(0:2,2,2)
+ integer                        :: h1,h2,p1,p2,s1,s2
+ double precision               :: phase
+ double precision               :: contrib
+ integer                        :: occ(N_int*bit_kind_size,2)
+ integer                        :: n_occ_ab(2)
+ two_body_dm_ab_big_array = 0.d0
+ BEGIN_DOC
+! The alpha-beta energy can be computed thanks to 
+! sum_{h1,p1,h2,p2} two_body_dm_ab_big_array(h1,p1,h2,p2) * (h1p1|h2p2)
+ END_DOC
+
+ do i = 1, N_det ! i == |I>
+  call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
+  do j = i+1, N_det ! j == <J|
+   call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+   if(degree>2)cycle
+   call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
+   call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+   contrib = 0.5d0 * psi_coef(i,1) * psi_coef(j,1) * phase
+   if(degree==2)then  ! case of the DOUBLE EXCITATIONS  ************************************
+    if(s1==s2)cycle  ! Only the alpha/beta two body density matrix
+    ! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J
+    call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2)
+
+   else if(degree==1)then! case of the SINGLE EXCITATIONS  ***************************************************
+    if(s1==1)then  ! Mono alpha : 
+     do k = 1, elec_beta_num  
+      m = occ(k,2)
+      ! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
+      call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
+     enddo
+    else  ! Mono Beta : 
+     do k = 1, elec_alpha_num
+      m = occ(k,1)
+      ! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
+      call insert_into_two_body_dm_big_array(two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
+     enddo
+    endif
+
+   endif
+  enddo
+ enddo
+ print*,'Big array for density matrix provided !'
+
+
+
+END_PROVIDER 
+
+subroutine insert_into_two_body_dm_big_array(big_array,dim1,dim2,dim3,dim4,contrib,h1,p1,h2,p2)
+ implicit none
+ integer, intent(in) :: h1,p1,h2,p2
+ integer, intent(in) :: dim1,dim2,dim3,dim4
+ double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4)
+ double precision :: contrib
+ big_array(h1,p1,h2,p2) += 1.d0 * contrib  
+ big_array(p1,h1,h2,p2) += 1.d0 * contrib  
+ big_array(h1,p1,p2,h2) += 1.d0 * contrib  
+ big_array(p1,h1,p2,h2) += 1.d0 * contrib  
+
+!big_array(h2,p2,h1,p1) += 1.d0 * contrib  
+!big_array(p2,h2,h1,p1) += 1.d0 * contrib  
+!if(p2.ne.h2)then
+!big_array(h2,p2,p1,h1) += 1.d0 * contrib  
+!big_array(p2,h2,p1,h1) += 1.d0 * contrib  
+!endif
+ 
+end
diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f
index aa2feead..0f338877 100644
--- a/src/MO_Basis/utils.irp.f
+++ b/src/MO_Basis/utils.irp.f
@@ -268,3 +268,26 @@ subroutine mo_sort_by_observable(observable,label)
 end
 
 
+subroutine give_all_mos_at_r(r,mos_array)
+ implicit none
+ double precision, intent(in) :: r(3)
+ double precision, intent(out) :: mos_array(mo_tot_num)
+ call give_specific_mos_at_r(r,mos_array, mo_coef)
+end
+
+subroutine give_specific_mos_at_r(r,mos_array, mo_coef_specific)
+ implicit none
+ double precision, intent(in) :: r(3)
+ double precision, intent(in)  :: mo_coef_specific(ao_num_align, mo_tot_num)
+ double precision, intent(out) :: mos_array(mo_tot_num)
+ double precision :: aos_array(ao_num),accu
+ integer :: i,j
+ call give_all_aos_at_r(r,aos_array)
+ do i = 1, mo_tot_num
+  accu = 0.d0
+  do j = 1, ao_num
+   accu += mo_coef_specific(j,i) * aos_array(j) 
+  enddo
+  mos_array(i) = accu
+ enddo
+end