diff --git a/plugins/Properties/EZFIO.cfg b/plugins/Properties/EZFIO.cfg
index 02f42ad8..40ccd8b9 100644
--- a/plugins/Properties/EZFIO.cfg
+++ b/plugins/Properties/EZFIO.cfg
@@ -2,4 +2,10 @@
 type: double precision
 doc: z point on which the integrated delta rho is calculated
 interface: ezfio,provider,ocaml
-default: 3.9
\ No newline at end of file
+default: 3.9
+
+[threshld_two_bod_dm]
+type: double precision
+doc: threshold for the values of the alpha/beta two body dm evaluation
+interface: ezfio,provider,ocaml
+default: 0.000001
diff --git a/plugins/Properties/eval_ab_two_body_dm.irp.f b/plugins/Properties/eval_ab_two_body_dm.irp.f
new file mode 100644
index 00000000..595b0087
--- /dev/null
+++ b/plugins/Properties/eval_ab_two_body_dm.irp.f
@@ -0,0 +1,80 @@
+program two_bod_ab_dm
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call routine
+
+end
+
+subroutine routine
+ implicit none
+ integer :: i,j,k,l
+ double precision :: dx1,dx2
+ double precision :: interval_1,interval_2
+ integer :: nx1,nx2
+ double precision :: r1(3),r2(3)
+ double precision :: xmin_1,xmax_1,xmin_2,xmax_2
+ double precision :: compute_extra_diag_two_body_dm_ab,two_bod_extra_diag
+ double precision :: compute_diag_two_body_dm_ab,two_bod_diag
+
+ double precision,allocatable :: mos_array_r1(:),mos_array_r2(:)
+ double precision :: test_diag, test_extra_diag
+ double precision, allocatable :: x_array(:),y_array(:),z_diag_array(:),z_extra_diag_array(:),z_total_array(:)
+ allocate(mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num))
+ 
+
+
+ ! He triplet S
+ ! r1 = x
+ ! r2 = z
+ nx1 = 100
+ nx2 = 100
+ allocate(x_array(nx1),y_array(nx2),z_diag_array(nx1),z_extra_diag_array(nx1),z_total_array(nx1))
+ xmin_1 = -2.d0
+ xmax_1 =  2.d0
+ xmin_2 = -2.d0
+ xmax_2 =  2.d0
+ interval_1 = xmax_1 - xmin_1
+ interval_2 = xmax_2 - xmin_2
+ dx1 = interval_1/dble(nx1)
+ dx2 = interval_2/dble(nx2)
+ r1 = 0.d0
+ r2 = 0.d0
+ 
+ double precision :: x_tmp,y_tmp
+ x_tmp = xmin_1
+ do i = 1, nx1
+  x_array(i) = x_tmp 
+  write(i_unit_x_two_body_dm_ab,'(10000(F16.10,X))')x_array(i)
+  x_tmp += dx1
+ enddo
+ x_tmp = xmin_2
+ do i = 1, nx1
+  y_array(i) = x_tmp 
+  write(i_unit_y_two_body_dm_ab,'(10000(F16.10,X))')x_array(i)
+  x_tmp += dx2
+ enddo
+
+
+ ! initialization 
+ r1(1) = xmin_1
+ do i = 1, nx1
+  r2(3) = xmin_2
+  do j = 1, nx2
+   two_bod_extra_diag = compute_extra_diag_two_body_dm_ab(r1,r2)
+   two_bod_diag = compute_diag_two_body_dm_ab(r1,r2)
+   z_diag_array(j) = two_bod_diag 
+   z_extra_diag_array(j) = two_bod_extra_diag 
+   z_total_array(j) = two_bod_extra_diag + two_bod_diag
+!  write(i_unit_two_body_dm_ab,'(100(F16.10,X))')r1(1),r2(3),two_bod_diag,two_bod_extra_diag,two_bod_diag+two_bod_extra_diag
+   r2(3) += dx2
+  enddo
+  write(i_unit_z_two_body_diag_dm_ab,'(10000(F16.10,X))')z_diag_array(:)
+  write(i_unit_z_two_body_extra_diag_dm_ab,'(10000(F16.10,X))')z_extra_diag_array(:)
+  write(i_unit_z_two_body_total_dm_ab,'(10000(F16.10,X))')z_total_array(:)
+   r1(1) += dx1
+ enddo
+
+ deallocate(mos_array_r1,mos_array_r2)
+
+end
diff --git a/plugins/Properties/iunit_two_bod.irp.f b/plugins/Properties/iunit_two_bod.irp.f
new file mode 100644
index 00000000..e14d9893
--- /dev/null
+++ b/plugins/Properties/iunit_two_bod.irp.f
@@ -0,0 +1,45 @@
+BEGIN_PROVIDER [integer, i_unit_x_two_body_dm_ab]
+ implicit none
+ integer :: getUnitAndOpen
+ character*(128) :: file_name
+ file_name = trim(trim(ezfio_filename)//'/properties/two_body_dm_x')
+ i_unit_x_two_body_dm_ab = getUnitAndOpen(file_name,'w')
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, i_unit_y_two_body_dm_ab]
+ implicit none
+ integer :: getUnitAndOpen
+ character*(128) :: file_name
+ file_name = trim(trim(ezfio_filename)//'/properties/two_body_dm_y')
+ i_unit_y_two_body_dm_ab = getUnitAndOpen(file_name,'w')
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, i_unit_z_two_body_extra_diag_dm_ab]
+ implicit none
+ integer :: getUnitAndOpen
+ character*(128) :: file_name
+ file_name = trim(trim(ezfio_filename)//'/properties/two_body_dm_extra_diag')
+ i_unit_z_two_body_extra_diag_dm_ab = getUnitAndOpen(file_name,'w')
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, i_unit_z_two_body_diag_dm_ab]
+ implicit none
+ integer :: getUnitAndOpen
+ character*(128) :: file_name
+ file_name = trim(trim(ezfio_filename)//'/properties/two_body_dm_diag')
+ i_unit_z_two_body_diag_dm_ab = getUnitAndOpen(file_name,'w')
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, i_unit_z_two_body_total_dm_ab]
+ implicit none
+ integer :: getUnitAndOpen
+ character*(128) :: file_name
+ file_name = trim(trim(ezfio_filename)//'/properties/two_body_dm_total')
+ i_unit_z_two_body_total_dm_ab = getUnitAndOpen(file_name,'w')
+
+END_PROVIDER 
+
diff --git a/plugins/Properties/test_two_body_dm.irp.f b/plugins/Properties/test_two_body_dm.irp.f
new file mode 100644
index 00000000..6fc02abf
--- /dev/null
+++ b/plugins/Properties/test_two_body_dm.irp.f
@@ -0,0 +1,67 @@
+program test_two_bod
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call routine
+end
+subroutine routine
+ implicit none
+ integer :: i,j,k,l
+ double precision :: accu,get_two_body_dm_ab_map_element,get_mo_bielec_integral_schwartz
+ accu = 0.d0
+
+ ! Diag part of the two body dm
+ do i = 1, n_act_orb
+  do j = 1, n_act_orb
+   accu += two_body_dm_ab_diag(i,j) * mo_bielec_integral_jj(i,j)
+  enddo
+ enddo
+ print*,'BI ELECTRONIC   =  ',accu
+
+ double precision :: accu_extra_diag
+ accu_extra_diag = 0.d0
+ do l = 1, n_act_orb  ! p2 
+  do k = 1, n_act_orb  ! h2 
+   do j = 1, n_act_orb  ! p1 
+    do i = 1,n_act_orb   ! h1 
+     accu_extra_diag += two_body_dm_ab_big_array(i,j,k,l) * get_mo_bielec_integral_schwartz(i,k,j,l,mo_integrals_map)
+    enddo
+   enddo
+  enddo
+ enddo
+ print*,'extra_diag = ',accu_extra_diag
+ double precision :: average_mono
+ call get_average(mo_mono_elec_integral,one_body_dm_mo,average_mono)
+ print*,'BI ELECTRONIC   =  ',accu+accu_extra_diag
+ print*,'MONO ELECTRONIC =  ',average_mono
+ print*,'Total elec      =  ',accu+average_mono + accu_extra_diag
+ print*,'Total           =  ',accu+average_mono+nuclear_repulsion +accu_extra_diag
+ double precision :: e_0,hij
+ call u0_H_u_0(e_0,psi_coef,n_det,psi_det,N_int)
+ print*,'<Psi| H |Psi>   =  ',e_0 + nuclear_repulsion
+ integer :: degree
+ integer                        :: exc(0:2,2,2)
+ integer                        :: h1,h2,p1,p2,s1,s2
+ double precision               :: phase
+ integer :: n_elements
+ n_elements = 0
+ accu = 0.d0
+ do i = 1, N_det
+  do j = i+1, N_det
+   call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+   if(degree.gt.2)cycle
+!  if(degree.ne.1)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)
+   if(s1.eq.s2)cycle
+   n_elements += 1
+   call i_H_j(psi_det(1,1,i),psi_det(1,1,j),N_int,hij)
+   accu += 2.d0 * hij * psi_coef(i,1) * psi_coef(j,1)
+  enddo
+ enddo
+ print*,'n_elements  = ',n_elements
+ print*,'<Psi| extra diag   ',accu
+ print*,'dm                 ',accu_extra_diag
+ 
+
+end
diff --git a/src/Determinants/two_body_dm_map.irp.f b/src/Determinants/two_body_dm_map.irp.f
index 38467040..c89b1125 100644
--- a/src/Determinants/two_body_dm_map.irp.f
+++ b/src/Determinants/two_body_dm_map.irp.f
@@ -209,7 +209,7 @@ BEGIN_PROVIDER [double precision, two_body_dm_ab_diag, (mo_tot_num, mo_tot_num)]
    contrib = psi_coef(i,1)**2
    do j = 1, elec_beta_num  
     k = occ(j,2)
-    do l = 1, elec_beta_num  
+    do l = 1, elec_alpha_num
      m = occ(l,1)
      two_body_dm_ab_diag(k,m) += 0.5d0 * contrib
      two_body_dm_ab_diag(m,k) += 0.5d0 * contrib
@@ -276,43 +276,66 @@ subroutine insert_into_two_body_dm_big_array(big_array,dim1,dim2,dim3,dim4,contr
  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  
+ ! Two spin symmetry
+ big_array(h1,p1,h2,p2) +=  contrib  
+ big_array(h2,p2,h1,p1) +=  contrib  
+ ! Hermicity : hole-particle symmetry
+ big_array(p1,h1,p2,h2) +=  contrib  
+ big_array(p2,h2,p1,h1) +=  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
 
-double precision function compute_two_body_dm_ab(r1,r2)
+double precision function compute_extra_diag_two_body_dm_ab(r1,r2)
  implicit none
  double precision :: r1(3),r2(3)
  integer :: i,j,k,l
  double precision :: mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num)
  double precision :: contrib
- compute_two_body_dm_ab = 0.d0
- call give_all_mos_at_r(r1,mos_array_r1)
- call give_all_mos_at_r(r2,mos_array_r2)
+ compute_extra_diag_two_body_dm_ab = 0.d0
+!call give_all_mos_at_r(r1,mos_array_r1)
+!call give_all_mos_at_r(r2,mos_array_r2)
+ call give_all_act_mos_at_r(r1,mos_array_r1)
+ call give_all_act_mos_at_r(r2,mos_array_r2)
  do l = 1, n_act_orb  ! p2 
   contrib = mos_array_r2(l)
-  if(dabs(contrib).lt.1.d-6)cycle
+! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
   do k = 1, n_act_orb  ! h2 
-    contrib *= mos_array_r2(k)
-   if(dabs(contrib).lt.1.d-6)cycle
+!  contrib *= mos_array_r2(k)
+!  if(dabs(contrib*mos_array_r2(k)).lt.threshld_two_bod_dm)cycle
    do j = 1, n_act_orb  ! p1 
-    contrib *= mos_array_r1(j)
-    if(dabs(contrib).lt.1.d-6)cycle
+!   contrib *= mos_array_r1(j)
+!   if(dabs(contrib).lt.threshld_two_bod_dm)cycle
     do i = 1,n_act_orb   ! h1 
-     compute_two_body_dm_ab += two_body_dm_ab_big_array(i,j,k,l) * mos_array_r1(i) * contrib
+     double precision :: contrib_tmp
+     contrib_tmp = mos_array_r1(i) * mos_array_r1(j) * mos_array_r2(k) * mos_array_r2(l)
+     compute_extra_diag_two_body_dm_ab += two_body_dm_ab_big_array(i,j,k,l) * contrib_tmp
     enddo
    enddo
   enddo
  enddo
 
 end
+
+double precision function compute_diag_two_body_dm_ab(r1,r2)
+ implicit none
+ double precision :: r1(3),r2(3)
+ integer :: i,j,k,l
+ double precision :: mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num)
+ double precision :: contrib,contrib_tmp
+ compute_diag_two_body_dm_ab = 0.d0
+ call give_all_act_mos_at_r(r1,mos_array_r1)
+ call give_all_act_mos_at_r(r2,mos_array_r2)
+ do l = 1, n_act_orb  ! 
+  contrib = mos_array_r2(l)*mos_array_r2(l)
+! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
+  do k = 1, n_act_orb  ! 
+    contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k)
+!  if(dabs(contrib).lt.threshld_two_bod_dm)cycle
+     compute_diag_two_body_dm_ab += two_body_dm_ab_diag(k,l) * contrib_tmp
+  enddo
+ enddo
+
+end
+
+