diff --git a/external/qp2-dependencies b/external/qp2-dependencies
index ce14f57b..f40bde09 160000
--- a/external/qp2-dependencies
+++ b/external/qp2-dependencies
@@ -1 +1 @@
-Subproject commit ce14f57b50511825a9fedb096749200779d3f4d4
+Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8
diff --git a/src/basis_correction/TODO b/src/basis_correction/TODO
index e28d593a..64a6ddeb 100644
--- a/src/basis_correction/TODO
+++ b/src/basis_correction/TODO
@@ -1 +1,4 @@
 change all correlation functionals with the pbe_on_top general
+factor 2 in two-rdm involved in: 
+  on-top, mu(r), pbe-on-top, sc_basis_corr and so on 
+  casscf : state_av_act_2_rdm_spin_trace_mo
diff --git a/src/basis_correction/pbe_on_top.irp.f b/src/basis_correction/pbe_on_top.irp.f
index cc9cdec9..a25fd61b 100644
--- a/src/basis_correction/pbe_on_top.irp.f
+++ b/src/basis_correction/pbe_on_top.irp.f
@@ -41,13 +41,15 @@
    
    if(mu_of_r_potential == "cas_ful")then
     ! You take the on-top of the CAS wave function which is computed with mu(r) 
-    on_top = on_top_cas_mu_r(ipoint,istate)
+    ! factor 2 because convention  N(N-1)/ 2 in provider on_top_cas_mu_r
+    on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate)
    else
     ! You take the on-top of the CAS wave function computed separately
+    ! No factor 2 because convention  N(N-1) in provider total_cas_on_top_density
     on_top = total_cas_on_top_density(ipoint,istate)
    endif
-!  We take the extrapolated on-top pair density * 2 because of normalization
-   on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
+!  We take the extrapolated on-top pair density 
+   on_top_extrap = mu_correction_of_on_top(mu,on_top)
 
    call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
  
@@ -103,13 +105,15 @@
    
    if(mu_of_r_potential == "cas_ful")then
     ! You take the on-top of the CAS wave function which is computed with mu(r) 
-    on_top = on_top_cas_mu_r(ipoint,istate)
+    ! factor 2 because convention  N(N-1)/ 2 in provider on_top_cas_mu_r
+    on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate)
    else
     ! You take the on-top of the CAS wave function computed separately
+    ! No factor 2 because convention  N(N-1) in provider total_cas_on_top_density
     on_top = total_cas_on_top_density(ipoint,istate)
    endif
-!  We take the extrapolated on-top pair density * 2 because of normalization
-   on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
+!  We take the extrapolated on-top pair density 
+   on_top_extrap = mu_correction_of_on_top(mu,on_top)
 
    call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
  
@@ -165,13 +169,15 @@
    
    if(mu_of_r_potential == "cas_ful")then
     ! You take the on-top of the CAS wave function which is computed with mu(r) 
-    on_top = on_top_cas_mu_r(ipoint,istate)
+    ! factor 2 because convention  N(N-1)/ 2 in provider on_top_cas_mu_r
+    on_top = 1.d0 * on_top_cas_mu_r(ipoint,istate)
    else
     ! You take the on-top of the CAS wave function computed separately
+    ! No factor 2 because convention  N(N-1) in provider total_cas_on_top_density
     on_top = total_cas_on_top_density(ipoint,istate)
    endif
-!  We DO NOT take the extrapolated on-top pair density, but there is  * 2 because of normalization
-   on_top_extrap = 2.d0 * on_top
+!  We DO NOT take the extrapolated on-top pair density
+   on_top_extrap =  on_top
 
    call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
  
diff --git a/src/basis_correction/print_routine.irp.f b/src/basis_correction/print_routine.irp.f
index 67c5c6c2..c2558d22 100644
--- a/src/basis_correction/print_routine.irp.f
+++ b/src/basis_correction/print_routine.irp.f
@@ -18,7 +18,7 @@ subroutine print_basis_correction
  print*, ''
  print*, 'For more details look at Journal of Chemical Physics 149, 194301 1-15 (2018)        '
  print*, '                         Journal of Physical Chemistry Letters 10, 2931-2937 (2019) '
- print*, '                         ???REF SC?'
+ print*, '                         Journal of Chemical Physics 152, 174104 (2020)             '
  print*, '****************************************'
  print*, '****************************************'
  print*, 'mu_of_r_potential = ',mu_of_r_potential
@@ -56,14 +56,14 @@ subroutine print_basis_correction
    print*,''
    print*,'********************************************'
    print*,'********************************************'
-   print*,'+) PBE-on-top Ecmd functional : (??????? REF-SCF ??????????)' 
+   print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) ' 
    print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization'
    do istate = 1, N_states
     write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-OT        , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate)
    enddo
    print*,''
    print*,'********************************************'
-   print*,'+) PBE-on-top no spin polarization Ecmd functional : (??????? REF-SCF ??????????)' 
+   print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)' 
    print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION'
    do istate = 1, N_states
     write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD SU-PBE-OT     , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)
diff --git a/src/mu_of_r/f_psi_i_a_v_utils.irp.f b/src/mu_of_r/f_psi_i_a_v_utils.irp.f
index 427da199..0d08e193 100644
--- a/src/mu_of_r/f_psi_i_a_v_utils.irp.f
+++ b/src/mu_of_r/f_psi_i_a_v_utils.irp.f
@@ -194,7 +194,7 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
   do b = 1, n_act_orb  ! 2
    do c = 1, n_act_orb  ! 1
     do d = 1, n_act_orb  ! 2
-     rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate)
+     rho = mos_array_act_r1(c) * mos_array_act_r2(d) * 0.5d0 * act_2_rdm_ab_mo(d,c,b,a,istate)
      rho_tilde(b,a) += rho 
      two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) 
     enddo
diff --git a/src/mu_of_r/f_psi_utils.irp.f b/src/mu_of_r/f_psi_utils.irp.f
index bdd76f18..95f10f36 100644
--- a/src/mu_of_r/f_psi_utils.irp.f
+++ b/src/mu_of_r/f_psi_utils.irp.f
@@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, full_occ_2_rdm_cntrctd_trans, (n_points_final_
      do j = 1, n_core_inact_act_orb
       do i = 1, n_core_inact_act_orb
                                              !                                 1 2 1 2 
-       full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint)
+       full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += 0.5d0 * full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint)
       enddo
      enddo
     enddo
diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f
index 5c41acdc..f9c3b3b3 100644
--- a/src/mu_of_r/mu_of_r_conditions.irp.f
+++ b/src/mu_of_r/mu_of_r_conditions.irp.f
@@ -110,7 +110,7 @@
  do istate = 1, N_states
   do ipoint = 1, n_points_final_grid
    f_psi  = f_psi_cas_ab(ipoint,istate) 
-   on_top = on_top_cas_mu_r(ipoint,istate)
+   on_top = on_top_cas_mu_r(ipoint,istate) 
    if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then
      w_psi   = 1.d+10
    else 
diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f
index 41b28aea..c550e991 100644
--- a/src/two_body_rdm/act_2_rdm.irp.f
+++ b/src/two_body_rdm/act_2_rdm.irp.f
@@ -4,21 +4,18 @@
  BEGIN_DOC
 !                             12 12
 !                 1 2 1 2 == <ij|kl>
-! act_2_rdm_ab_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons 
+! act_2_rdm_ab_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of alpha/beta+beta/alpha electrons 
 ! 
-! <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>
+!   <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>
+!
+! + <Psi_{istate}| a^{\dagger}_{i \beta} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \beta} |Psi_{istate}>
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act}
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is alpha, electron 2 is beta
-!
-!  act_2_rdm_ab_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta
-!                      
-!                      Therefore you don't necessary have symmetry between electron 1 and 2 
  END_DOC 
  integer :: ispin
  double precision :: wall_1, wall_2
@@ -44,6 +41,7 @@
  endif
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_ab_mo',wall_2 - wall_1
+ act_2_rdm_ab_mo *= 2.d0
  END_PROVIDER 
 
 
@@ -56,7 +54,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
  END_DOC 
@@ -84,6 +82,7 @@
 
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_aa_mo',wall_2 - wall_1
+ act_2_rdm_aa_mo *= 2.d0
  END_PROVIDER 
  
 
@@ -96,7 +95,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
  END_DOC 
@@ -124,6 +123,7 @@
 
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_bb_mo',wall_2 - wall_1
+ act_2_rdm_bb_mo *= 2.d0
  END_PROVIDER 
 
  BEGIN_PROVIDER [double precision, act_2_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
@@ -135,7 +135,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec}^{act} * (N_{elec}^{act} - 1)/2 
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec}^{act} * (N_{elec}^{act} - 1)
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
  END_DOC 
@@ -161,6 +161,7 @@
   call ezfio_set_two_body_rdm_io_two_body_rdm_spin_trace("Read")
  endif
 
+ act_2_rdm_spin_trace_mo *= 2.d0
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_spin_trace_mo',wall_2 - wall_1
  END_PROVIDER 
diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f
index de3d97b9..30e2685a 100644
--- a/src/two_body_rdm/example.irp.f
+++ b/src/two_body_rdm/example.irp.f
@@ -117,10 +117,10 @@ subroutine routine_active_only
        endif
  
 
-       wee_ab(istate)    += vijkl * rdmab
-       wee_aa(istate)    += vijkl * rdmaa
-       wee_bb(istate)    += vijkl * rdmbb
-       wee_tot(istate)    += vijkl * rdmtot
+       wee_ab(istate)    += 0.5d0 * vijkl * rdmab
+       wee_aa(istate)    += 0.5d0 * vijkl * rdmaa
+       wee_bb(istate)    += 0.5d0 * vijkl * rdmbb
+       wee_tot(istate)   += 0.5d0 * vijkl * rdmtot
 
       enddo
      enddo
@@ -144,13 +144,13 @@ subroutine routine_active_only
    print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate)
    print*,'--------------------------'
    print*,'accu_aa       = ',accu_aa
-   print*,'N_a (N_a-1)/2 = ', elec_alpha_num*(elec_alpha_num-1)*0.5
+   print*,'N_a (N_a-1)   = ', elec_alpha_num*(elec_alpha_num-1)
    print*,'accu_bb       = ',accu_bb
-   print*,'N_b (N_b-1)/2 = ', elec_beta_num*(elec_beta_num-1)*0.5
+   print*,'2 N_b (N_b-1) = ', elec_beta_num*(elec_beta_num-1)*2
    print*,'accu_ab       = ',accu_ab
    print*,'N_a N_b       = ', elec_beta_num*elec_alpha_num
    print*,'accu_tot      = ',accu_tot
-   print*,'Ne(Ne-1)/2    = ',(elec_num-1)*elec_num * 0.5
+   print*,'Ne(Ne-1)      = ',(elec_num-1)*elec_num 
   enddo
  wee_aa_st_av  = 0.d0
  wee_bb_st_av  = 0.d0
@@ -192,10 +192,10 @@ subroutine routine_active_only
       print*,spin_trace,rdm_tot_st_av,dabs(spin_trace - rdm_tot_st_av)
      endif
 
-     wee_aa_st_av += vijkl * rdm_aa_st_av
-     wee_bb_st_av += vijkl * rdm_bb_st_av
-     wee_ab_st_av += vijkl * rdm_ab_st_av
-     wee_tot_st_av += vijkl * rdm_tot_st_av
+     wee_aa_st_av  += 0.5d0 * vijkl * rdm_aa_st_av
+     wee_bb_st_av  += 0.5d0 * vijkl * rdm_bb_st_av
+     wee_ab_st_av  += 0.5d0 * vijkl * rdm_ab_st_av
+     wee_tot_st_av += 0.5d0 * vijkl * rdm_tot_st_av
     enddo
    enddo
   enddo
@@ -279,10 +279,10 @@ subroutine routine_full_mos
        rdmbb  =  full_occ_2_rdm_bb_mo(l,k,j,i,istate)
        rdmtot =  full_occ_2_rdm_spin_trace_mo(l,k,j,i,istate)
 
-       wee_ab(istate) += vijkl * rdmab
-       wee_aa(istate) += vijkl * rdmaa
-       wee_bb(istate) += vijkl * rdmbb
-       wee_tot(istate)+= vijkl * rdmtot
+       wee_ab(istate) += 0.5d0 * vijkl * rdmab
+       wee_aa(istate) += 0.5d0 * vijkl * rdmaa
+       wee_bb(istate) += 0.5d0 * vijkl * rdmbb
+       wee_tot(istate)+= 0.5d0 * vijkl * rdmtot
       enddo
      enddo
      aa_norm(istate) += full_occ_2_rdm_aa_mo(j,i,j,i,istate)
@@ -310,18 +310,19 @@ subroutine routine_full_mos
    print*,'Normalization of two-rdms '
    print*,''
    print*,'aa_norm(istate)         = ',aa_norm(istate)
-   print*,'N_alpha(N_alpha-1)/2    = ',elec_num_tab(1) * (elec_num_tab(1) - 1)/2
+   print*,'N_alpha(N_alpha-1)      = ',elec_num_tab(1) * (elec_num_tab(1) - 1)
    print*,''
    print*,'bb_norm(istate)         = ',bb_norm(istate)
-   print*,'N_alpha(N_alpha-1)/2    = ',elec_num_tab(2) * (elec_num_tab(2) - 1)/2
+   print*,'N_alpha(N_alpha-1)      = ',elec_num_tab(2) * (elec_num_tab(2) - 1)
    print*,''
    print*,'ab_norm(istate)         = ',ab_norm(istate)
-   print*,'N_alpha * N_beta        = ',elec_num_tab(1) * elec_num_tab(2)
+   print*,'N_alpha * N_beta *2     = ',elec_num_tab(1) * elec_num_tab(2) * 2
    print*,''
    print*,'tot_norm(istate)        = ',tot_norm(istate)
-   print*,'N(N-1)/2                = ',elec_num*(elec_num - 1)/2
+   print*,'N(N-1)                  = ',elec_num*(elec_num - 1)
   enddo
 
+! return 
  wee_aa_st_av  = 0.d0
  wee_bb_st_av  = 0.d0
  wee_ab_st_av  = 0.d0
@@ -341,10 +342,10 @@ subroutine routine_full_mos
      rdm_ab_st_av  = state_av_full_occ_2_rdm_ab_mo(l,k,j,i)
      rdm_tot_st_av = state_av_full_occ_2_rdm_spin_trace_mo(l,k,j,i)
 
-     wee_aa_st_av  += vijkl * rdm_aa_st_av
-     wee_bb_st_av  += vijkl * rdm_bb_st_av
-     wee_ab_st_av  += vijkl * rdm_ab_st_av
-     wee_tot_st_av += vijkl * rdm_tot_st_av
+     wee_aa_st_av  += 0.5d0 * vijkl * rdm_aa_st_av
+     wee_bb_st_av  += 0.5d0 * vijkl * rdm_bb_st_av
+     wee_ab_st_av  += 0.5d0 * vijkl * rdm_ab_st_av
+     wee_tot_st_av += 0.5d0 * vijkl * rdm_tot_st_av
     enddo
    enddo
   enddo
@@ -354,17 +355,12 @@ subroutine routine_full_mos
  print*,''
  print*,'STATE AVERAGE ENERGY '
    print*,'wee_aa_st_av            = ',wee_aa_st_av
-   print*,'wee_aa_st_av_2          = ',wee_aa_st_av_2
    print*,'wee_bb_st_av            = ',wee_bb_st_av
-   print*,'wee_bb_st_av_2          = ',wee_bb_st_av_2
    print*,'wee_ab_st_av            = ',wee_ab_st_av
-   print*,'wee_ab_st_av_2          = ',wee_ab_st_av_2
    print*,'Sum of components       = ',wee_aa_st_av   + wee_bb_st_av   + wee_ab_st_av
-   print*,'Sum of components_2     = ',wee_aa_st_av_2 + wee_bb_st_av_2 + wee_ab_st_av_2
    print*,''
    print*,'Full energy '
    print*,'wee_tot_st_av           = ',wee_tot_st_av
-   print*,'wee_tot_st_av_2         = ',wee_tot_st_av_2
    print*,'wee_tot_st_av_3         = ',wee_tot_st_av_3
 
 end
diff --git a/src/two_body_rdm/full_orb_2_rdm.irp.f b/src/two_body_rdm/full_orb_2_rdm.irp.f
index fba88172..78a28acf 100644
--- a/src/two_body_rdm/full_orb_2_rdm.irp.f
+++ b/src/two_body_rdm/full_orb_2_rdm.irp.f
@@ -4,22 +4,18 @@
  full_occ_2_rdm_ab_mo = 0.d0
  integer :: i,j,k,l,iorb,jorb,korb,lorb,istate
  BEGIN_DOC
-! full_occ_2_rdm_ab_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons 
+! full_occ_2_rdm_ab_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of alpha/beta + beta/alpha electrons 
 !
-! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
+!   <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
+!
+! + <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \beta} |Psi>
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta}
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} * 2
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA 
-!
-!  full_occ_2_rdm_ab_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta
-!                      
-!                      Therefore you don't necessary have symmetry between electron 1 and 2 
-!
 !  !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO ARE SET TO ZERO 
  END_DOC 
  full_occ_2_rdm_ab_mo = 0.d0
@@ -50,7 +46,7 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                     alph beta alph beta
-      full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = 2.d0 * one_e_dm_mo_beta(jorb,iorb,istate)
      enddo
     enddo
    enddo
@@ -64,7 +60,7 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                     alph beta alph beta
-      full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = 2.d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
      enddo
     enddo
    enddo
@@ -76,7 +72,7 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                     alph beta alph beta
-      full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0
+      full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 2.D0
      enddo
     enddo
 
@@ -93,7 +89,7 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                     alph beta alph beta
-       full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = 2.d0 * one_e_dm_mo_beta(jorb,iorb,istate)
       enddo
      enddo
     enddo
@@ -107,7 +103,7 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                     alph beta alph beta
-       full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = 2.d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
       enddo
      enddo
     enddo
@@ -119,7 +115,7 @@
      do k = 1, n_core_orb
       korb = list_core(k)
       !                                     alph beta alph beta
-      full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0
+      full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 2.D0
      enddo
     enddo
    endif
@@ -139,7 +135,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
@@ -172,11 +168,11 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                       1     2   1    2    : DIRECT TERM 
-      full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) +=  0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
-      full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) +=  0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) +=  1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) +=  1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
       !                                       1     2   1    2    : EXCHANGE TERM 
-      full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
-      full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
      enddo
     enddo
    enddo
@@ -186,8 +182,8 @@
     jorb = list_inact(j)
     do k = 1, n_inact_orb
      korb = list_inact(k)
-     full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) +=  0.5d0 
-     full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -=  0.5d0 
+     full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) +=  1.0d0 
+     full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -=  1.0d0 
     enddo
    enddo
 
@@ -203,11 +199,11 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                       1     2   1    2    : DIRECT TERM 
-       full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) +=  0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
-       full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) +=  0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) +=  1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) +=  1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
        !                                       1     2   1    2    : EXCHANGE TERM 
-       full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
-       full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
       enddo
      enddo
     enddo
@@ -217,8 +213,8 @@
      jorb = list_core(j)
      do k = 1, n_core_orb
       korb = list_core(k)
-      full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) +=  0.5d0 
-      full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -=  0.5d0 
+      full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) +=  1.0d0 
+      full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -=  1.0d0 
      enddo
     enddo
    endif
@@ -237,7 +233,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
@@ -270,11 +266,11 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                       1     2   1    2    : DIRECT TERM 
-      full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) +=  0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
-      full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) +=  0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) +=  1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) +=  1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
       !                                       1     2   1    2    : EXCHANGE TERM 
-      full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
-      full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
      enddo
     enddo
    enddo
@@ -284,8 +280,8 @@
     jorb = list_inact(j)
     do k = 1, n_inact_orb
      korb = list_inact(k)
-     full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) +=  0.5d0 
-     full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -=  0.5d0 
+     full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) +=  1.0d0 
+     full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -=  1.0d0 
     enddo
    enddo
 
@@ -301,11 +297,11 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                       1     2   1    2    : DIRECT TERM 
-       full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) +=  0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
-       full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) +=  0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) +=  1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) +=  1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
        !                                       1     2   1    2    : EXCHANGE TERM 
-       full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
-       full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
       enddo
      enddo
     enddo
@@ -315,8 +311,8 @@
      jorb = list_core(j)
      do k = 1, n_core_orb
       korb = list_core(k)
-      full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) +=  0.5d0 
-      full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -=  0.5d0 
+      full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) +=  1.0d0 
+      full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -=  1.0d0 
      enddo
     enddo
    endif
@@ -335,14 +331,14 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
 !
 !  !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
 !  !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero 
 ! The two-electron energy of each state can be computed as:
 !
-!   \sum_{i,j,k,l = 1, n_core_inact_act_orb} full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) * < ii jj | kk ll > 
+!   \sum_{i,j,k,l = 1, n_core_inact_act_orb} full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) * 1/2 * < ii jj | kk ll > 
 !
 !   with ii = list_core_inact_act(i), jj = list_core_inact_act(j), kk = list_core_inact_act(k), ll = list_core_inact_act(l)
  END_DOC 
@@ -377,11 +373,11 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                       1     2   1    2    : DIRECT TERM 
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) +=  0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
-      full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) +=  0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) +=  1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) +=  1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
       !                                       1     2   1    2    : EXCHANGE TERM 
-      full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
      enddo
     enddo
    enddo
@@ -390,8 +386,8 @@
     jorb = list_inact(j)
     do k = 1, n_inact_orb
      korb = list_inact(k)
-     full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) +=  0.5d0 
-     full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -=  0.5d0 
+     full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) +=  1.0d0 
+     full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -=  1.0d0 
     enddo
    enddo
    if (.not.no_core_density)then
@@ -403,11 +399,11 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                       1     2   1    2    : DIRECT TERM 
-       full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) +=  0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
-       full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) +=  0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) +=  1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) +=  1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
        !                                       1     2   1    2    : EXCHANGE TERM 
-       full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
-       full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
       enddo
      enddo
     enddo
@@ -416,8 +412,8 @@
      jorb = list_core(j)
      do k = 1, n_core_orb
       korb = list_core(k)
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) +=  0.5d0 
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -=  0.5d0 
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) +=  1.0d0 
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -=  1.0d0 
      enddo
     enddo
    endif
@@ -433,11 +429,11 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                       1     2   1    2    : DIRECT TERM 
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) +=  0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
-      full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) +=  0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) +=  1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) +=  1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
       !                                       1     2   1    2    : EXCHANGE TERM 
-      full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
      enddo
     enddo
    enddo
@@ -446,8 +442,8 @@
     jorb = list_inact(j)
     do k = 1, n_inact_orb
      korb = list_inact(k)
-     full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) +=  0.5d0 
-     full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -=  0.5d0 
+     full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) +=  1.0d0 
+     full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -=  1.0d0 
     enddo
    enddo
    if (.not.no_core_density)then
@@ -459,11 +455,11 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                       1     2   1    2    : DIRECT TERM 
-       full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) +=  0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
-       full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) +=  0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) +=  1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) +=  1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
        !                                       1     2   1    2    : EXCHANGE TERM 
-       full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
-       full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
       enddo
      enddo
     enddo
@@ -472,8 +468,8 @@
      jorb = list_core(j)
      do k = 1, n_core_orb
       korb = list_core(k)
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) +=  0.5d0 
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -=  0.5d0 
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) +=  1.0d0 
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -=  1.0d0 
      enddo
     enddo
    endif
@@ -489,14 +485,14 @@
       korb = list_inact(k)
       ! ALPHA INACTIVE - BETA ACTIVE
       !                                     alph beta alph beta
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
       !                                     beta alph beta alph
-      full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
       ! BETA INACTIVE - ALPHA ACTIVE
       !                                     beta alph beta alpha 
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
       !                                     alph beta alph beta 
-      full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+      full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
      enddo
     enddo
    enddo
@@ -506,8 +502,8 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                     alph beta alph beta
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0
-      full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0D0
+      full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 1.0D0
      enddo
     enddo
 
@@ -523,14 +519,14 @@
        korb = list_core(k)
        !! BETA ACTIVE - ALPHA CORE 
        !                                     alph beta alph beta
-       full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0D0 * one_e_dm_mo_beta(jorb,iorb,istate)
        !                                     beta alph beta alph 
-       full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0D0 * one_e_dm_mo_beta(jorb,iorb,istate)
        !! ALPHA ACTIVE - BETA CORE 
        !                                     alph beta alph beta
-       full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
        !                                     beta alph beta alph 
-       full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
+       full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
       enddo
      enddo
     enddo
@@ -540,8 +536,8 @@
      do k = 1, n_core_orb
       korb = list_core(k)
       !                                     alph beta alph beta
-      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0
-      full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0
+      full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0D0
+      full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 1.0D0
      enddo
     enddo
 
diff --git a/src/two_body_rdm/state_av_act_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f
index db793047..cd417a9d 100644
--- a/src/two_body_rdm/state_av_act_2rdm.irp.f
+++ b/src/two_body_rdm/state_av_act_2rdm.irp.f
@@ -2,21 +2,16 @@
  implicit none
  double precision, allocatable :: state_weights(:) 
  BEGIN_DOC
-! state_av_act_2_rdm_ab_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs
+! state_av_act_2_rdm_ab_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha/beta+beta/alpha electron pairs
 !
 !                                     = \sum_{istate} w(istate) * <Psi_{istate}| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi_{istate}>
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act}
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is alpha, electron 2 is beta
-!
-!  state_av_act_2_rdm_ab_mo(i,j,k,l) = i:alpha, j:beta, j:alpha, l:beta
-!                      
-!                      Therefore you don't necessary have symmetry between electron 1 and 2 
  END_DOC 
  allocate(state_weights(N_states))
  state_weights = state_average_weight
@@ -34,6 +29,8 @@
  call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  call wall_time(wall_2)
  print*,'Wall time to provide state_av_act_2_rdm_ab_mo',wall_2 - wall_1
+ ! factor 2 to have the correct normalization factor 
+ state_av_act_2_rdm_ab_mo *= 2.d0
 
  END_PROVIDER 
 
@@ -48,7 +45,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
  END_DOC 
@@ -63,6 +60,8 @@
  call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  call wall_time(wall_2)
  print*,'Wall time to provide state_av_act_2_rdm_aa_mo',wall_2 - wall_1
+ ! factor 2 to have the correct normalization factor 
+ state_av_act_2_rdm_aa_mo *= 2.d0
 
  END_PROVIDER 
 
@@ -76,7 +75,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
  END_DOC 
@@ -91,6 +90,8 @@
  call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  call wall_time(wall_2)
  print*,'Wall time to provide state_av_act_2_rdm_bb_mo',wall_2 - wall_1
+ ! factor 2 to have the correct normalization factor 
+ state_av_act_2_rdm_bb_mo *= 2.d0
 
  END_PROVIDER 
 
@@ -104,7 +105,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
  END_DOC
diff --git a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f
index b3a5fe65..2e44665d 100644
--- a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f
+++ b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f
@@ -4,22 +4,16 @@
  state_av_full_occ_2_rdm_ab_mo = 0.d0
  integer :: i,j,k,l,iorb,jorb,korb,lorb
  BEGIN_DOC
-! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) =  STATE AVERAGE physicist notation for 2RDM of alpha/beta electrons 
+! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) =  STATE AVERAGE physicist notation for 2RDM of alpha/beta + beta/alpha electrons 
 !
 !                                     = \sum_{istate} w(istate) * <Psi_{istate}| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi_{istate}>
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta}
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} * 2
 !
 !  !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 ! 
-!  !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA 
-!
-!  state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = i:alpha, j:beta, j:alpha, l:beta
-!                      
-!                      Therefore you don't necessary have symmetry between electron 1 and 2 
-!
 !  !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero 
  END_DOC 
  state_av_full_occ_2_rdm_ab_mo = 0.d0
@@ -47,7 +41,7 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                     alph beta alph beta
-      state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = 2.d0 * one_e_dm_mo_beta_average(jorb,iorb)
      enddo
     enddo
    enddo
@@ -61,7 +55,7 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                     alph beta alph beta
-      state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = 2.d0 * one_e_dm_mo_alpha_average(jorb,iorb)
      enddo
     enddo
    enddo
@@ -73,7 +67,7 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                     alph beta alph beta
-      state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 1.D0
+      state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0
      enddo
     enddo
 
@@ -90,7 +84,7 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                     alph beta alph beta
-       state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = 2.d0 * one_e_dm_mo_beta_average(jorb,iorb)
       enddo
      enddo
     enddo
@@ -104,7 +98,7 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                     alph beta alph beta
-       state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = 2.d0 * one_e_dm_mo_alpha_average(jorb,iorb)
       enddo
      enddo
     enddo
@@ -116,7 +110,7 @@
      do k = 1, n_core_orb
       korb = list_core(k)
       !                                     alph beta alph beta
-      state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 1.D0
+      state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0
      enddo
     enddo
    endif
@@ -135,7 +129,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)
 !
 !  !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
@@ -167,11 +161,11 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                       1     2   1    2    : DIRECT TERM 
-      state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) +=  0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
-      state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) +=  0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) +=  1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) +=  1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
       !                                       1     2   1    2    : EXCHANGE TERM 
-      state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
-      state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
      enddo
     enddo
    enddo
@@ -181,8 +175,8 @@
     jorb = list_inact(j)
     do k = 1, n_inact_orb
      korb = list_inact(k)
-     state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) +=  0.5d0 
-     state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -=  0.5d0 
+     state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) +=  1.0d0 
+     state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -=  1.0d0 
     enddo
    enddo
 
@@ -198,11 +192,11 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                       1     2   1    2    : DIRECT TERM 
-       state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) +=  0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
-       state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) +=  0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) +=  1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) +=  1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
        !                                       1     2   1    2    : EXCHANGE TERM 
-       state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
-       state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
       enddo
      enddo
     enddo
@@ -212,8 +206,8 @@
      jorb = list_core(j)
      do k = 1, n_core_orb
       korb = list_core(k)
-      state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) +=  0.5d0 
-      state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -=  0.5d0 
+      state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) +=  1.0d0 
+      state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -=  1.0d0 
      enddo
     enddo
    endif
@@ -231,7 +225,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)
 !
 ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
@@ -263,11 +257,11 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                       1     2   1    2    : DIRECT TERM 
-      state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) +=  0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
-      state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) +=  0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) +=  1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) +=  1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
       !                                       1     2   1    2    : EXCHANGE TERM 
-      state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
-      state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
      enddo
     enddo
    enddo
@@ -277,8 +271,8 @@
     jorb = list_inact(j)
     do k = 1, n_inact_orb
      korb = list_inact(k)
-     state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) +=  0.5d0 
-     state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -=  0.5d0 
+     state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) +=  1.0d0 
+     state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -=  1.0d0 
     enddo
    enddo
 
@@ -294,11 +288,11 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                       1     2   1    2    : DIRECT TERM 
-       state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) +=  0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
-       state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) +=  0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) +=  1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) +=  1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
        !                                       1     2   1    2    : EXCHANGE TERM 
-       state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
-       state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
       enddo
      enddo
     enddo
@@ -308,8 +302,8 @@
      jorb = list_core(j)
      do k = 1, n_core_orb
       korb = list_core(k)
-      state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) +=  0.5d0 
-      state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -=  0.5d0 
+      state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) +=  1.0d0 
+      state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -=  1.0d0 
      enddo
     enddo
    endif
@@ -328,7 +322,7 @@
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
-! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2
+! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
 !
 !  !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !
@@ -364,11 +358,11 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                       1     2   1    2    : DIRECT TERM 
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) +=  0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
-      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) +=  0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) +=  1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) +=  1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
       !                                       1     2   1    2    : EXCHANGE TERM 
-      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
      enddo
     enddo
    enddo
@@ -377,8 +371,8 @@
     jorb = list_inact(j)
     do k = 1, n_inact_orb
      korb = list_inact(k)
-     state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) +=  0.5d0 
-     state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  0.5d0 
+     state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) +=  1.0d0 
+     state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  1.0d0 
     enddo
    enddo
    if (.not.no_core_density)then
@@ -390,11 +384,11 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                       1     2   1    2    : DIRECT TERM 
-       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) +=  0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
-       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) +=  0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) +=  1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) +=  1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
        !                                       1     2   1    2    : EXCHANGE TERM 
-       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
-       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
       enddo
      enddo
     enddo
@@ -403,8 +397,8 @@
      jorb = list_core(j)
      do k = 1, n_core_orb
       korb = list_core(k)
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) +=  0.5d0 
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  0.5d0 
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) +=  1.0d0 
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  1.0d0 
      enddo
     enddo
    endif
@@ -420,11 +414,11 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                       1     2   1    2    : DIRECT TERM 
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) +=  0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
-      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) +=  0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) +=  1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) +=  1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
       !                                       1     2   1    2    : EXCHANGE TERM 
-      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
      enddo
     enddo
    enddo
@@ -433,8 +427,8 @@
     jorb = list_inact(j)
     do k = 1, n_inact_orb
      korb = list_inact(k)
-     state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) +=  0.5d0 
-     state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  0.5d0 
+     state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) +=  1.0d0 
+     state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  1.0d0 
     enddo
    enddo
    if (.not.no_core_density)then
@@ -446,11 +440,11 @@
       do k = 1, n_core_orb
        korb = list_core(k)
        !                                       1     2   1    2    : DIRECT TERM 
-       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) +=  0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
-       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) +=  0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) +=  1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) +=  1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
        !                                       1     2   1    2    : EXCHANGE TERM 
-       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
-       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
       enddo
      enddo
     enddo
@@ -459,8 +453,8 @@
      jorb = list_core(j)
      do k = 1, n_core_orb
       korb = list_core(k)
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) +=  0.5d0 
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  0.5d0 
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) +=  1.0d0 
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  1.0d0 
      enddo
     enddo
    endif
@@ -476,14 +470,14 @@
       korb = list_inact(k)
       ! ALPHA INACTIVE - BETA ACTIVE
       !                                     alph beta alph beta
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
       !                                     beta alph beta alph
-      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
       ! BETA INACTIVE - ALPHA ACTIVE
       !                                     beta alph beta alpha 
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
       !                                     alph beta alph beta 
-      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
+      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
      enddo
     enddo
    enddo
@@ -493,8 +487,8 @@
      do k = 1, n_inact_orb
       korb = list_inact(k)
       !                                     alph beta alph beta
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0
-      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 0.5D0
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
+      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0d0
      enddo
     enddo
 
@@ -510,14 +504,14 @@
        korb = list_core(k)
        !! BETA ACTIVE - ALPHA CORE 
        !                                     alph beta alph beta
-       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0D0 * one_e_dm_mo_beta_average(jorb,iorb)
        !                                     beta alph beta alph 
-       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0D0 * one_e_dm_mo_beta_average(jorb,iorb)
        !! ALPHA ACTIVE - BETA CORE 
        !                                     alph beta alph beta
-       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0D0 * one_e_dm_mo_alpha_average(jorb,iorb)
        !                                     beta alph beta alph 
-       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb)
+       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0D0 * one_e_dm_mo_alpha_average(jorb,iorb)
       enddo
      enddo
     enddo
@@ -527,8 +521,8 @@
      do k = 1, n_core_orb
       korb = list_core(k)
       !                                     alph beta alph beta
-      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0
-      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 0.5D0
+      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0D0
+      state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0D0
      enddo
     enddo
 
diff --git a/src/two_body_rdm/test_2_rdm.irp.f b/src/two_body_rdm/test_2_rdm.irp.f
index 123261d8..4eb8f9f0 100644
--- a/src/two_body_rdm/test_2_rdm.irp.f
+++ b/src/two_body_rdm/test_2_rdm.irp.f
@@ -2,7 +2,7 @@ program test_2_rdm
  implicit none
  read_wf = .True.
  touch read_wf 
- call routine_active_only
+! call routine_active_only
  call routine_full_mos
 end
 
diff --git a/src/two_body_rdm/two_e_dm_mo.irp.f b/src/two_body_rdm/two_e_dm_mo.irp.f
index a4dea15f..7e35fc7b 100644
--- a/src/two_body_rdm/two_e_dm_mo.irp.f
+++ b/src/two_body_rdm/two_e_dm_mo.irp.f
@@ -29,7 +29,8 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)]
      enddo
     enddo
    enddo
-   two_e_dm_mo(:,:,:,:) = two_e_dm_mo(:,:,:,:) * 2.d0
+   two_e_dm_mo(:,:,:,:) = two_e_dm_mo(:,:,:,:) 
+!   * 2.d0
 
  END_PROVIDER
 
diff --git a/src/two_rdm_routines/davidson_like_2rdm.irp.f b/src/two_rdm_routines/davidson_like_2rdm.irp.f
index 2e5aa4d1..ad7a3b21 100644
--- a/src/two_rdm_routines/davidson_like_2rdm.irp.f
+++ b/src/two_rdm_routines/davidson_like_2rdm.irp.f
@@ -4,8 +4,8 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
    BEGIN_DOC
    ! if ispin == 1 :: alpha/alpha 2rdm 
    !          == 2 :: beta /beta  2rdm 
-   !          == 3 :: alpha/beta  2rdm 
-   !          == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba))
+   !          == 3 :: alpha/beta + beta/alpha  2rdm 
+   !          == 4 :: spin traced 2rdm :: aa + bb + ab + ba
    !
    ! Assumes that the determinants are in psi_det
    !