diff --git a/config/gfortran.cfg b/config/gfortran.cfg
index c0aa875f..60e32235 100644
--- a/config/gfortran.cfg
+++ b/config/gfortran.cfg
@@ -35,14 +35,14 @@ OPENMP  : 1          ; Append OpenMP flags
 #                           -ffast-math and the Fortran-specific
 #                           -fno-protect-parens and -fstack-arrays.
 [OPT]
-FCFLAGS : -Ofast
+FCFLAGS : 
 
 # Profiling flags
 #################
 #
 [PROFILE]
 FC      : -p -g
-FCFLAGS : -Ofast 
+FCFLAGS : 
 
 # Debugging flags
 #################
diff --git a/plugins/FCIdump/NEEDED_CHILDREN_MODULES b/plugins/FCIdump/NEEDED_CHILDREN_MODULES
index 34de8ddb..8d60d3c7 100644
--- a/plugins/FCIdump/NEEDED_CHILDREN_MODULES
+++ b/plugins/FCIdump/NEEDED_CHILDREN_MODULES
@@ -1 +1 @@
-Determinants Davidson
+Determinants Davidson core_integrals
diff --git a/plugins/FCIdump/fcidump.irp.f b/plugins/FCIdump/fcidump.irp.f
index f93c1128..8d334fc5 100644
--- a/plugins/FCIdump/fcidump.irp.f
+++ b/plugins/FCIdump/fcidump.irp.f
@@ -1,21 +1,25 @@
 program fcidump
   implicit none
+  character*(128) :: output
+  integer :: i_unit_output,getUnitAndOpen
+  output=trim(ezfio_filename)//'.FCIDUMP'
+  i_unit_output = getUnitAndOpen(output,'w')
 
   integer :: i,j,k,l
-  integer :: ii(8), jj(8), kk(8),ll(8)
+  integer :: i1,j1,k1,l1
+  integer :: i2,j2,k2,l2
   integer*8 :: m
   character*(2), allocatable :: A(:)
 
-  print *, '&FCI NORB=', mo_tot_num, ', NELEC=', elec_num, &
+  write(i_unit_output,*) '&FCI NORB=', n_act_orb, ', NELEC=', elec_num-n_core_orb*2, &
    ', MS2=', (elec_alpha_num-elec_beta_num), ','
-  allocate (A(mo_tot_num))
+  allocate (A(n_act_orb))
   A = '1,'
-  print *, 'ORBSYM=', (A(i), i=1,mo_tot_num) 
-  print *,'ISYM=0,'
-  print *,'/'
+  write(i_unit_output,*) 'ORBSYM=', (A(i), i=1,n_act_orb) 
+  write(i_unit_output,*) 'ISYM=0,'
+  write(i_unit_output,*) '/'
   deallocate(A)
   
-  integer*8 :: i8, k1
   integer(key_kind), allocatable :: keys(:)
   double precision, allocatable  :: values(:)
   integer(cache_map_size_kind)   :: n_elements, n_elements_max
@@ -23,14 +27,18 @@ program fcidump
 
   double precision :: get_mo_bielec_integral, integral
 
-  do l=1,mo_tot_num
-   do k=1,mo_tot_num
-    do j=l,mo_tot_num
-     do i=k,mo_tot_num
-       if (i>=j) then
-          integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
+  do l=1,n_act_orb
+   l1 = list_act(l)
+   do k=1,n_act_orb
+    k1 = list_act(k)
+    do j=l,n_act_orb
+     j1 = list_act(j)
+     do i=k,n_act_orb
+      i1 = list_act(i)
+       if (i1>=j1) then
+          integral = get_mo_bielec_integral(i1,j1,k1,l1,mo_integrals_map)
           if (dabs(integral) > mo_integrals_threshold) then 
-            print *, integral, i,k,j,l
+            write(i_unit_output,*) integral, i,k,j,l
           endif
        end if
      enddo
@@ -38,13 +46,15 @@ program fcidump
    enddo
   enddo
 
-  do j=1,mo_tot_num
-   do i=j,mo_tot_num
-      integral = mo_mono_elec_integral(i,j)
+  do j=1,n_act_orb
+   j1 = list_act(j)
+   do i=j,n_act_orb
+    i1 = list_act(i)
+      integral = mo_mono_elec_integral(i1,j1) + core_fock_operator(i1,j1)
       if (dabs(integral) > mo_integrals_threshold) then 
-        print *, integral, i,j,0,0
+        write(i_unit_output,*) integral, i,j,0,0
       endif
    enddo
   enddo
-  print *, 0.d0, 0, 0, 0, 0
+  write(i_unit_output,*) core_energy, 0, 0, 0, 0
 end
diff --git a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES
index 16fce081..25d61c69 100644
--- a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES
+++ b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES
@@ -1 +1 @@
-Perturbation Selectors_no_sorted Hartree_Fock Davidson CISD
+Perturbation Selectors_no_sorted SCF_density Davidson CISD
diff --git a/plugins/FOBOCI/all_singles.irp.f b/plugins/FOBOCI/all_singles.irp.f
index 65d81e07..7c321b72 100644
--- a/plugins/FOBOCI/all_singles.irp.f
+++ b/plugins/FOBOCI/all_singles.irp.f
@@ -48,6 +48,7 @@ subroutine all_single(e_pt2)
     print*,'-----------------------'
     print*,'i = ',i
     call H_apply_just_mono(pt2, norm_pert, H_pert_diag,  N_st)
+    call make_s2_eigenfunction_first_order
     call diagonalize_CI
     print*,'N_det = ',N_det
     print*,'E        = ',CI_energy(1)
diff --git a/plugins/FOBOCI/create_1h_or_1p.irp.f b/plugins/FOBOCI/create_1h_or_1p.irp.f
index 41ec7b6c..c5205903 100644
--- a/plugins/FOBOCI/create_1h_or_1p.irp.f
+++ b/plugins/FOBOCI/create_1h_or_1p.irp.f
@@ -29,21 +29,13 @@ subroutine create_restart_and_1h(i_hole)
     enddo
    enddo
   enddo
+ 
 
  integer :: N_det_old
  N_det_old = N_det
- N_det += n_new_det
- allocate (new_det(N_int,2,n_new_det))
- if (psi_det_size < N_det) then
-   psi_det_size = N_det
-   TOUCH psi_det_size
- endif
- do i = 1, N_det_old 
-  do k = 1, N_int
-   psi_det(k,1,i) = old_psi_det(k,1,i)
-   psi_det(k,2,i) = old_psi_det(k,2,i)
-  enddo
- enddo
+
+ logical, allocatable :: duplicate(:)
+ allocate (new_det(N_int,2,n_new_det),duplicate(n_new_det))
 
  n_new_det = 0
   do j = 1, n_act_orb
@@ -58,19 +50,56 @@ subroutine create_restart_and_1h(i_hole)
      if(i_ok .ne. 1)cycle
      n_new_det +=1
      do k = 1, N_int
-      psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1) 
-      psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2) 
+      new_det(k,1,n_new_det) = key_tmp(k,1) 
+      new_det(k,2,n_new_det) = key_tmp(k,2) 
      enddo
-     psi_coef(n_det_old+n_new_det,:) = 0.d0
     enddo
    enddo
   enddo
 
-  SOFT_TOUCH N_det psi_det psi_coef
-  logical :: found_duplicates
-  if(n_act_orb.gt.1)then
-  call remove_duplicates_in_psi_det(found_duplicates)
+ integer :: i_test
+ duplicate = .False. 
+ do i = 1, n_new_det
+  if(duplicate(i))cycle
+  do j = i+1, n_new_det
+   i_test = 0
+   do ispin =1 ,2
+    do k = 1, N_int
+     i_test += popcnt(xor(new_det(k,ispin,i),new_det(k,ispin,j)))
+    enddo
+   enddo
+   if(i_test.eq.0)then
+    duplicate(j) = .True.
+   endif
+  enddo
+ enddo
+
+ integer :: n_new_det_unique
+ n_new_det_unique = 0
+ print*, 'uniq det'
+ do i = 1, n_new_det
+  if(.not.duplicate(i))then
+   n_new_det_unique += 1
   endif
+ enddo
+ print*, n_new_det_unique
+ N_det += n_new_det_unique
+ if (psi_det_size < N_det) then
+   psi_det_size = N_det
+   TOUCH psi_det_size
+ endif
+ do i = 1, n_new_det_unique
+  do ispin = 1, 2
+   do k = 1, N_int
+    psi_det(k,ispin,N_det_old+i) = new_det(k,ispin,i)
+   enddo
+  enddo
+  psi_coef(N_det_old+i,:) = 0.d0
+ enddo
+
+
+  SOFT_TOUCH N_det psi_det psi_coef
+  deallocate (new_det,duplicate)
 end
 
 subroutine create_restart_and_1p(i_particle)
@@ -107,18 +136,8 @@ subroutine create_restart_and_1p(i_particle)
 
  integer :: N_det_old
  N_det_old = N_det
- N_det += n_new_det
- allocate (new_det(N_int,2,n_new_det))
- if (psi_det_size < N_det) then
-   psi_det_size = N_det
-   TOUCH psi_det_size
- endif
- do i = 1, N_det_old 
-  do k = 1, N_int
-   psi_det(k,1,i) = old_psi_det(k,1,i)
-   psi_det(k,2,i) = old_psi_det(k,2,i)
-  enddo
- enddo
+ logical, allocatable :: duplicate(:)
+ allocate (new_det(N_int,2,n_new_det),duplicate(n_new_det))
 
  n_new_det = 0
   do j = 1, n_act_orb
@@ -133,17 +152,59 @@ subroutine create_restart_and_1p(i_particle)
      if(i_ok .ne. 1)cycle
      n_new_det +=1
      do k = 1, N_int
-      psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1) 
-      psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2) 
+      new_det(k,1,n_new_det) = key_tmp(k,1) 
+      new_Det(k,2,n_new_det) = key_tmp(k,2) 
      enddo
-     psi_coef(n_det_old+n_new_det,:) = 0.d0
     enddo
    enddo
   enddo
 
+ integer :: i_test
+ duplicate = .False. 
+ do i = 1, n_new_det
+  if(duplicate(i))cycle
+  call debug_det(new_det(1,1,i),N_int)
+  do j = i+1, n_new_det
+   i_test = 0
+   call debug_det(new_det(1,1,j),N_int)
+   do ispin =1 ,2
+    do k = 1, N_int
+     i_test += popcnt(xor(new_det(k,ispin,i),new_det(k,ispin,j)))
+    enddo
+   enddo
+   if(i_test.eq.0)then
+    duplicate(j) = .True.
+   endif
+  enddo
+ enddo
+
+ integer :: n_new_det_unique
+ n_new_det_unique = 0
+ print*, 'uniq det'
+ do i = 1, n_new_det
+  if(.not.duplicate(i))then
+   n_new_det_unique += 1
+  endif
+ enddo
+ print*, n_new_det_unique
+
+ N_det += n_new_det_unique
+ if (psi_det_size < N_det) then
+   psi_det_size = N_det
+   TOUCH psi_det_size
+ endif
+ do i = 1, n_new_det_unique
+  do ispin = 1, 2
+   do k = 1, N_int
+    psi_det(k,ispin,N_det_old+i) = new_det(k,ispin,i)
+   enddo
+  enddo
+  psi_coef(N_det_old+i,:) = 0.d0
+ enddo
+
   SOFT_TOUCH N_det psi_det psi_coef
-  logical :: found_duplicates
-  call remove_duplicates_in_psi_det(found_duplicates)
+  deallocate (new_det,duplicate)
+
 end
 
 subroutine create_restart_1h_1p(i_hole,i_part)
diff --git a/plugins/FOBOCI/density_matrix.irp.f b/plugins/FOBOCI/density_matrix.irp.f
index aaf80c4f..42138c00 100644
--- a/plugins/FOBOCI/density_matrix.irp.f
+++ b/plugins/FOBOCI/density_matrix.irp.f
@@ -32,6 +32,11 @@
     psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1) * inv_coef_ref_generators_restart
     norm_generators_restart += psi_coef_generators_restart(i,1)**2
    enddo
+   double precision :: inv_norm
+   inv_norm = 1.d0/dsqrt(norm_generators_restart)
+   do i = 1, N_det_generators_restart
+    psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1)  * inv_norm
+   enddo
 
 
      one_body_dm_mo_alpha_generators_restart = 0.d0
diff --git a/plugins/FOBOCI/dress_simple.irp.f b/plugins/FOBOCI/dress_simple.irp.f
index dd1ed221..fabbd834 100644
--- a/plugins/FOBOCI/dress_simple.irp.f
+++ b/plugins/FOBOCI/dress_simple.irp.f
@@ -175,6 +175,10 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
   endif
   do j = 1, Ndet_generators
    call i_h_j(psi_det_generators_input(1,1,j),psi_det_generators_input(1,1,i),N_int,hij)  ! Fill the zeroth order H matrix
+   if(i==j)then
+    call debug_det(psi_det_generators_input(1,1,i),N_int)
+    print*, hij
+   endif
    dressed_H_matrix(i,j) = hij
   enddo
  enddo
@@ -234,6 +238,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
   do i = 1, N_states
    i_state(i) = i
    E_ref(i) = eigvalues(i)
+   print*, 'E_ref(i)',E_ref(i)
   enddo
  endif
  do i = 1,N_states
@@ -287,7 +292,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
       accu += eigvectors(j,i) * psi_coef_ref(j,k)
      enddo
      print*,'accu = ',accu
-     if(dabs(accu).ge.0.72d0)then
+     if(dabs(accu).ge.0.60d0)then
       i_good_state(0) +=1
       i_good_state(i_good_state(0)) = i
      endif
diff --git a/plugins/FOBOCI/fobo_scf.irp.f b/plugins/FOBOCI/fobo_scf.irp.f
index 8a709154..8be36b8a 100644
--- a/plugins/FOBOCI/fobo_scf.irp.f
+++ b/plugins/FOBOCI/fobo_scf.irp.f
@@ -15,8 +15,6 @@ end
 
 subroutine run_prepare
  implicit none
-! no_oa_or_av_opt = .False.
-! touch no_oa_or_av_opt
   call damping_SCF
   call diag_inactive_virt_and_update_mos
 end
@@ -28,7 +26,7 @@ subroutine routine_fobo_scf
  print*,''
  character*(64) :: label
  label = "Natural"
- do i = 1, 5
+ do i = 1, 1
   print*,'*******************************************************************************'
   print*,'*******************************************************************************'
   print*,'FOBO-SCF Iteration ',i
@@ -54,7 +52,7 @@ subroutine routine_fobo_scf
   endif
   call FOBOCI_lmct_mlct_old_thr(i)
   call save_osoci_natural_mos
-  call damping_SCF
+! call damping_SCF
   call diag_inactive_virt_and_update_mos
   call clear_mo_map
   call provide_properties
diff --git a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f
index 46ca9662..3d8dfb08 100644
--- a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f
+++ b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f
@@ -55,6 +55,10 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter)
    call create_restart_and_1h(i_hole_osoci)
    call set_generators_to_psi_det
    print*,'Passed set generators'
+   integer :: m
+   do m = 1, N_det_generators
+    call debug_det(psi_det_generators(1,1,m),N_int)
+   enddo
    call set_bitmask_particl_as_input(reunion_of_bitmask)
    call set_bitmask_hole_as_input(reunion_of_bitmask)
    double precision :: e_pt2
@@ -82,7 +86,7 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter)
      call set_bitmask_particl_as_input(reunion_of_bitmask)
      call set_bitmask_hole_as_input(reunion_of_bitmask)
      call all_single(e_pt2)
-     call make_s2_eigenfunction_first_order
+!    call make_s2_eigenfunction_first_order
      threshold_davidson = 1.d-6
      soft_touch threshold_davidson davidson_criterion
      call diagonalize_ci
@@ -541,7 +545,6 @@ subroutine FOBOCI_lmct_mlct_old_thr_restart(iter)
      call print_generators_bitmasks_holes
      ! Impose that only the active part can be reached 
      call set_bitmask_hole_as_input(unpaired_bitmask)
-!!!  call all_single_h_core
      call create_restart_and_1p(i_particl_osoci)
 !!!  ! Update the generators 
      call set_generators_to_psi_det
diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f
index 7d194a54..26ce3b12 100644
--- a/plugins/FOBOCI/routines_foboci.irp.f
+++ b/plugins/FOBOCI/routines_foboci.irp.f
@@ -13,6 +13,8 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
  integer :: n_good_hole
  logical,allocatable :: is_a_ref_det(:)
  allocate(index_one_hole(n_det),index_one_hole_one_p(n_det),index_two_hole_one_p(N_det),index_two_hole(N_det),index_one_p(N_det),is_a_ref_det(N_det))
+ double precision, allocatable :: local_norm(:)
+ allocate(local_norm(N_states))
  
  n_one_hole = 0
  n_one_hole_one_p = 0
@@ -30,7 +32,6 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
    do k = 1, N_states
     inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k)
    enddo
-!  cycle
   endif
   
   ! Find all the determinants present in the reference wave function
@@ -59,10 +60,8 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
    enddo
   endif
  enddo
-!do k = 1, N_det
-! call debug_det(psi_det(1,1,k),N_int)
-! print*,'k,coef = ',k,psi_coef(k,1)/psi_coef(index_ref_generators_restart,1)
-!enddo
+
+
  print*,''
  print*,'n_good_hole = ',n_good_hole
  do k = 1,N_states
@@ -72,27 +71,37 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
   enddo
   print*,''
  enddo
- norm = 0.d0
 
- ! Set the wave function to the intermediate normalization
+ ! Set the wave function to the intermediate normalization 
  do k = 1, N_states
   do i = 1, N_det
    psi_coef(i,k) = psi_coef(i,k) * inv_coef_ref_generators_restart(k)
   enddo
  enddo
+
+
+ norm = 0.d0
  do k = 1,N_states
   print*,'state ',k
   do i = 1, N_det
-!!  print*,'psi_coef(i_ref) = ',psi_coef(i,1)
    if (is_a_ref_det(i))then
     print*,'i,psi_coef_ref = ',psi_coef(i,k)
-    cycle
    endif
    norm(k) += psi_coef(i,k) * psi_coef(i,k)
   enddo
   print*,'norm = ',norm(k)
  enddo
+ do k =1, N_states
+  local_norm(k) = 1.d0 / dsqrt(norm(k))
+ enddo
+ do k = 1,N_states
+  do i = 1, N_det
+   psi_coef(i,k) = psi_coef(i,k) * local_norm(k)
+  enddo
+ enddo
+ 
  deallocate(index_one_hole,index_one_hole_one_p,index_two_hole_one_p,index_two_hole,index_one_p,is_a_ref_det)
+ deallocate(local_norm)
  soft_touch psi_coef
 end
 
@@ -117,6 +126,8 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl)
  integer :: i_count
  allocate(index_one_hole(n_det),index_one_hole_one_p(n_det),index_two_hole_one_p(N_det),index_two_hole(N_det),index_one_p(N_det),is_a_ref_det(N_det))
  allocate(index_one_hole_two_p(n_det))
+ double precision, allocatable :: local_norm(:)
+ allocate(local_norm(N_states))
  
  n_one_hole = 0
  n_one_hole_one_p = 0
@@ -185,20 +196,29 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl)
    psi_coef(i,k) = psi_coef(i,k) * inv_coef_ref_generators_restart(k)
   enddo
  enddo
- do k = 1, N_states
+
+ norm = 0.d0
+ do k = 1,N_states
   print*,'state ',k
   do i = 1, N_det
-!! print*,'i = ',i, psi_coef(i,1)
    if (is_a_ref_det(i))then
     print*,'i,psi_coef_ref = ',psi_coef(i,k)
-    cycle
    endif
    norm(k) += psi_coef(i,k) * psi_coef(i,k)
   enddo
-  print*,'norm = ',norm
+  print*,'norm = ',norm(k)
+ enddo
+ do k =1, N_states
+  local_norm(k) = 1.d0 / dsqrt(norm(k))
+ enddo
+ do k = 1,N_states
+  do i = 1, N_det
+   psi_coef(i,k) = psi_coef(i,k) * local_norm(k)
+  enddo
  enddo
  soft_touch psi_coef
  deallocate(index_one_hole,index_one_hole_one_p,index_two_hole_one_p,index_two_hole,index_one_p,is_a_ref_det)
+ deallocate(local_norm)
 end
 
 
@@ -210,12 +230,60 @@ subroutine update_density_matrix_osoci
  END_DOC
  integer :: i,j
  integer :: iorb,jorb
+ ! active <--> inactive block
  do i = 1, mo_tot_num
   do j = 1, mo_tot_num
-   one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) + (one_body_dm_mo_alpha_average(i,j) - one_body_dm_mo_alpha_generators_restart(i,j))
-   one_body_dm_mo_beta_osoci(i,j) = one_body_dm_mo_beta_osoci(i,j) + (one_body_dm_mo_beta_average(i,j) - one_body_dm_mo_beta_generators_restart(i,j))
+   one_body_dm_mo_alpha_osoci(i,j) += one_body_dm_mo_alpha_average(i,j) - one_body_dm_mo_alpha_generators_restart(i,j)
+   one_body_dm_mo_beta_osoci(i,j) += one_body_dm_mo_beta_average(i,j) - one_body_dm_mo_beta_generators_restart(i,j)
   enddo
  enddo
+!do i = 1, n_act_orb
+! iorb = list_act(i)
+! do j = 1, n_inact_orb
+!  jorb = list_inact(j)
+!  one_body_dm_mo_alpha_osoci(iorb,jorb)+= one_body_dm_mo_alpha_average(iorb,jorb)
+!  one_body_dm_mo_alpha_osoci(jorb,iorb)+= one_body_dm_mo_alpha_average(jorb,iorb)
+!  one_body_dm_mo_beta_osoci(iorb,jorb) += one_body_dm_mo_beta_average(iorb,jorb) 
+!  one_body_dm_mo_beta_osoci(jorb,iorb) += one_body_dm_mo_beta_average(jorb,iorb) 
+! enddo
+!enddo
+
+!! active <--> virt block
+!do i = 1, n_act_orb
+! iorb = list_act(i)
+! do j = 1, n_virt_orb
+!  jorb = list_virt(j)
+!  one_body_dm_mo_alpha_osoci(iorb,jorb)+= one_body_dm_mo_alpha_average(iorb,jorb)
+!  one_body_dm_mo_alpha_osoci(jorb,iorb)+= one_body_dm_mo_alpha_average(jorb,iorb)
+!  one_body_dm_mo_beta_osoci(iorb,jorb) += one_body_dm_mo_beta_average(iorb,jorb) 
+!  one_body_dm_mo_beta_osoci(jorb,iorb) += one_body_dm_mo_beta_average(jorb,iorb) 
+! enddo
+!enddo
+
+!! virt <--> virt block
+!do j = 1, n_virt_orb
+!  jorb = list_virt(j)
+!  one_body_dm_mo_alpha_osoci(jorb,jorb)+= one_body_dm_mo_alpha_average(jorb,jorb)
+!  one_body_dm_mo_beta_osoci(jorb,jorb) += one_body_dm_mo_beta_average(jorb,jorb) 
+!enddo
+
+!! inact <--> inact block
+!do j = 1, n_inact_orb
+!  jorb = list_inact(j)
+!  one_body_dm_mo_alpha_osoci(jorb,jorb) -= one_body_dm_mo_alpha_average(jorb,jorb)
+!  one_body_dm_mo_beta_osoci(jorb,jorb)  -= one_body_dm_mo_beta_average(jorb,jorb) 
+!enddo
+ double precision :: accu_alpha, accu_beta
+ accu_alpha = 0.d0
+ accu_beta = 0.d0
+ do i = 1, mo_tot_num
+  accu_alpha += one_body_dm_mo_alpha_osoci(i,i) 
+  accu_beta += one_body_dm_mo_beta_osoci(i,i) 
+! write(*,'(I3,X,100(F16.10,X))') i,one_body_dm_mo_alpha_osoci(i,i),one_body_dm_mo_beta_osoci(i,i),one_body_dm_mo_alpha_osoci(i,i)+one_body_dm_mo_beta_osoci(i,i)
+ enddo 
+ print*, 'accu_alpha/beta',accu_alpha,accu_beta
+ 
+ 
 
 
 end
@@ -263,6 +331,12 @@ subroutine initialize_density_matrix_osoci
  implicit none
  one_body_dm_mo_alpha_osoci = one_body_dm_mo_alpha_generators_restart
  one_body_dm_mo_beta_osoci  = one_body_dm_mo_beta_generators_restart
+ integer :: i
+ print*, '8*********************'
+  print*, 'initialize_density_matrix_osoci'
+ do i = 1, mo_tot_num
+  print*,one_body_dm_mo_alpha_osoci(i,i),one_body_dm_mo_alpha_generators_restart(i,i)
+ enddo
 end
 
 subroutine rescale_density_matrix_osoci(norm)
@@ -438,6 +512,10 @@ subroutine save_osoci_natural_mos
    endif
   enddo
  enddo
+ print*, 'test'
+ print*, 'test'
+ print*, 'test'
+ print*, 'test'
  do i = 1, mo_tot_num
   do j = i+1, mo_tot_num
    if(dabs(tmp(i,j)).le.threshold_fobo_dm)then
@@ -445,7 +523,9 @@ subroutine save_osoci_natural_mos
       tmp(j,i) = 0.d0
    endif
   enddo
+  print*, tmp(i,i)
  enddo
+ 
 
  label = "Natural"
  
diff --git a/plugins/MRPT/MRPT_Utils.main.irp.f b/plugins/MRPT/MRPT_Utils.main.irp.f
index af3713c5..ab7a0ccb 100644
--- a/plugins/MRPT/MRPT_Utils.main.irp.f
+++ b/plugins/MRPT/MRPT_Utils.main.irp.f
@@ -18,11 +18,11 @@ subroutine routine_3
  print *,  'N_states = ', N_states
  do i = 1, N_States
   print*,'State',i
-  write(*,'(A12,X,I3,A3,XX,F16.10)')    '  PT2       ', i,' = ', second_order_pt_new(i)
-  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E         ', i,' = ', psi_ref_average_value(i)
-  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E+PT2     ', i,' = ', psi_ref_average_value(i)+second_order_pt_new(i)
-  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E dressed ', i,' = ', CI_dressed_pt2_new_energy(i)
-  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  S^2       ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i)
+  write(*,'(A12,X,I3,A3,XX,F20.16)')    '  PT2       ', i,' = ', second_order_pt_new(i)
+  write(*,'(A12,X,I3,A3,XX,F22.16)')    '  E         ', i,' = ', psi_ref_average_value(i)
+  write(*,'(A12,X,I3,A3,XX,F22.16)')    '  E+PT2     ', i,' = ', psi_ref_average_value(i)+second_order_pt_new(i)
+  write(*,'(A12,X,I3,A3,XX,F22.16)')    '  E dressed ', i,' = ', CI_dressed_pt2_new_energy(i)
+  write(*,'(A12,X,I3,A3,XX,F20.16)')    '  S^2       ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i)
   print*,'coef before and after'
   do j = 1, N_det_ref
    print*,psi_ref_coef(j,i),CI_dressed_pt2_new_eigenvectors(j,i)
diff --git a/plugins/MRPT/print_1h2p.irp.f b/plugins/MRPT/print_1h2p.irp.f
index b9f6575b..f20f12b6 100644
--- a/plugins/MRPT/print_1h2p.irp.f
+++ b/plugins/MRPT/print_1h2p.irp.f
@@ -5,6 +5,12 @@ program print_1h2p
  call routine
 end
 
+subroutine routine
+ implicit none
+ provide one_anhil_one_creat_inact_virt
+
+end
+
 subroutine routine_2
  implicit none
  integer :: i,j,degree
@@ -27,7 +33,7 @@ subroutine routine_2
 
 end
 
-subroutine routine
+subroutine routine_3
  implicit none
  double precision,allocatable :: matrix_1h2p(:,:,:) 
  double precision :: accu(2)
diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f
index 563b4bdf..f7e48e4f 100644
--- a/plugins/MRPT_Utils/energies_cas.irp.f
+++ b/plugins/MRPT_Utils/energies_cas.irp.f
@@ -22,6 +22,40 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange, (N_states)]
 END_PROVIDER 
 
 
+BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange_bis, (N_states)]
+ implicit none
+ integer :: i,j
+ double precision :: energies(N_states)
+ integer(bit_kind), allocatable :: psi_in_ref(:,:,:)
+ allocate (psi_in_ref(N_int,2,n_det_ref))
+ integer(bit_kind), allocatable :: psi_in_active(:,:,:)
+ allocate (psi_in_active(N_int,2,n_det_ref))
+ double precision, allocatable :: psi_ref_coef_in(:, :)
+ allocate(psi_ref_coef_in(N_det_ref, N_states))
+
+ do i = 1, N_det_ref
+  do j = 1, N_int
+   psi_in_ref(j,1,i) =  psi_ref(j,1,i) 
+   psi_in_ref(j,2,i) =  psi_ref(j,2,i) 
+
+   psi_in_active(j,1,i) =  psi_active(j,1,i) 
+   psi_in_active(j,2,i) =  psi_active(j,2,i) 
+  enddo
+  do j = 1, N_states
+   psi_ref_coef_in(i,j) = psi_ref_coef(i,j)
+  enddo
+ enddo 
+ do i = 1, N_states
+  call u0_H_dyall_u0_no_exchange_bis(energies,psi_in_ref,psi_ref_coef_in,n_det_ref,n_det_ref,n_det_ref,N_states,i)
+  energy_cas_dyall_no_exchange_bis(i) = energies(i)
+  print*,  'energy_cas_dyall(i)_no_exchange_bis',  energy_cas_dyall_no_exchange_bis(i)
+ enddo
+ deallocate (psi_in_ref)
+ deallocate (psi_in_active)
+ deallocate(psi_ref_coef_in)
+END_PROVIDER 
+
+
 
 BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
  implicit none
@@ -604,6 +638,8 @@ END_PROVIDER
  double precision, allocatable :: psi_in_out_coef(:,:)
  use bitmasks
  allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
+ integer(bit_kind), allocatable :: psi_in_active(:,:,:)
+ allocate (psi_in_active(N_int,2,n_det_ref))
 
  integer :: iorb,jorb,i_ok
  integer :: state_target
@@ -614,6 +650,9 @@ END_PROVIDER
 
 
  double precision :: thresh_norm
+ integer :: other_spin(2)
+ other_spin(1) = 2
+ other_spin(2) = 1
   
  thresh_norm = 1.d-20
 
@@ -644,14 +683,14 @@ END_PROVIDER
        print*, 'pb, i_ok ne 0 !!!'
       endif
       call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij)
+      integer                        :: exc(0:2,2,2)
+      double precision :: phase
+      call get_mono_excitation(psi_in_out(1,1,i),psi_ref(1,1,i),exc,phase,N_int)
       do j = 1, n_states
-        double precision ::  coef,contrib
-        coef = psi_ref_coef(i,j) !* psi_ref_coef(i,j)
-        psi_in_out_coef(i,j) = coef * hij 
-        if(orb_i == 5 .and. orb_v == 20)then
-!       if(orb_i == 2 .and. orb_v == 6 )then
+        psi_in_out_coef(i,j) = psi_ref_coef(i,j)* hij * phase
+!       if(orb_i == 5 .and. orb_v == 20)then
+        if(orb_i == 2 .and. orb_v == 6 )then
          print*, i, ispin
-         print*, coef * hij,coef,hij
         endif
         norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) 
       enddo
@@ -663,22 +702,31 @@ END_PROVIDER
        one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 0.d0
       else
        norm_no_inv(j,ispin) = norm(j,ispin)
-       one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin)
+!      one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin)
        norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
-       if(orb_i == 5 .and. orb_v == 20)then
-!      if(orb_i == 2 .and. orb_v == 6 )then
+!      if(orb_i == 5 .and. orb_v == 20)then
+       if(orb_i == 2 .and. orb_v == 6 )then
         print*,ispin ,norm(j,ispin)
        endif
       endif
      enddo
+     integer :: iorb_annil,hole_particle,spin_exc,orb
+     double precision :: norm_out_bis(N_states)
      do i = 1, N_det_ref
        do j = 1, N_states
         psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin)
         norm_bis(j,ispin) += psi_in_out_coef(i,j) *  psi_in_out_coef(i,j)
        enddo
+     enddo
+     
+     do i = 1, N_det_ref
        do j = 1, N_int
+ !      psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1))
+ !      psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1))
         psi_in_out(j,1,i) = psi_active(j,1,i)
         psi_in_out(j,2,i) = psi_active(j,2,i)
+ !      psi_in_out(j,1,i) = psi_ref(j,1,i)
+ !      psi_in_out(j,2,i) = psi_ref(j,2,i)
        enddo
      enddo
      do state_target = 1, N_states
@@ -686,29 +734,35 @@ END_PROVIDER
       if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
        call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
        energies_alpha_beta(state_target, ispin) +=  energies(state_target) 
-       if(orb_i == 5 .and. orb_v == 20)then
-!      if(orb_i == 2 .and. orb_v == 6 )then
-        print*, ispin, energy_cas_dyall_no_exchange(1) , energies_alpha_beta(state_target, ispin)
-        print*, ispin, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target, ispin)
-       endif
       endif
      enddo
     enddo ! ispin 
    do state_target = 1, N_states
     if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
      one_anhil_one_creat_inact_virt(iorb,vorb,state_target) =  energy_cas_dyall_no_exchange(state_target) - &
-      ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) & 
+!       0.5d0 * (energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2))
+      ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) )  & 
      /( norm_bis(state_target,1) +  norm_bis(state_target,2) )
     else  
      one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0
     endif
     if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.30584271462d0) < 1.d-11)then
+!   if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.29269686324d0) < 1.d-11)then
+     print*, ''
      print*, orb_i,orb_v
-     print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,1) / norm_bis(state_target,1)
-     print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,2) / norm_bis(state_target,2)
+     print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,1)  !/ norm_bis(state_target,1)
+     print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,2)  !/ norm_bis(state_target,2)
      print*, fock_core_inactive_total_spin_trace(orb_i,1)
      print*, fock_virt_total_spin_trace(orb_v,1)
      print*, one_anhil_one_creat_inact_virt(iorb,vorb,state_target)
+     print*, ''
+    endif
+    if(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) .gt. 1.d-10)then
+     write(*,'(F11.8)'), one_anhil_one_creat_inact_virt(iorb,vorb,state_target)
+!    if(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) .lt. 1.d-2)then
+!     one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0
+!     print*, orb_i,orb_v
+!    endif
     endif
    enddo
   enddo
@@ -852,9 +906,6 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
     norm_bis = 0.d0
     do ispin = 1,2
      do i = 1, n_det_ref
-!     if(orb_a == 6 .and. orb_v == 12)then
-!       print*, 'i ref = ',i
-!     endif
       do j = 1, N_int
        psi_in_out(j,1,i) =  psi_ref(j,1,i) 
        psi_in_out(j,2,i) =  psi_ref(j,2,i) 
@@ -866,11 +917,6 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
        enddo
       else
        call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij)
-       if(orb_a == 6 .and. orb_v == 12)then
-         call debug_det(psi_ref(1,1,i),N_int)
-         call debug_det(psi_in_out(1,1,i),N_int)
-         print*, hij
-       endif
        do j = 1, n_states
          double precision ::  contrib
          psi_in_out_coef(i,j) = psi_ref_coef(i,j) * hij 
@@ -907,7 +953,6 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
      do state_target = 1, N_states
       energies_alpha_beta(state_target, ispin) = 0.d0
       if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
-!      call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
        call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
        energies_alpha_beta(state_target, ispin) +=  energies(state_target) 
       endif
@@ -915,7 +960,6 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
     enddo ! ispin 
    do state_target = 1, N_states
     if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
-!    one_creat_virt(aorb,vorb,state_target) =  energy_cas_dyall(state_target) - &
      one_creat_virt(aorb,vorb,state_target) =  energy_cas_dyall_no_exchange(state_target) - &
       ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) & 
      /( norm_bis(state_target,1) +  norm_bis(state_target,2) )
diff --git a/plugins/MRPT_Utils/excitations_cas.irp.f b/plugins/MRPT_Utils/excitations_cas.irp.f
index 768abe8c..4042d90b 100644
--- a/plugins/MRPT_Utils/excitations_cas.irp.f
+++ b/plugins/MRPT_Utils/excitations_cas.irp.f
@@ -25,6 +25,7 @@ subroutine apply_exc_to_psi(orb,hole_particle,spin_exc, &
   integer :: det_tmp(N_int), det_tmp_bis(N_int)
   double precision :: phase
   double precision :: norm_factor
+! print*, orb,hole_particle,spin_exc
 
   elec_num_tab_local = 0
   do i = 1, ndet
@@ -36,6 +37,7 @@ subroutine apply_exc_to_psi(orb,hole_particle,spin_exc, &
     exit
    endif
   enddo
+! print*, elec_num_tab_local(1),elec_num_tab_local(2)
   if(hole_particle == 1)then
    do i = 1, ndet
      call set_bit_to_integer(orb,psi_in_out(1,spin_exc,i),N_int) 
@@ -675,6 +677,7 @@ subroutine i_H_j_dyall_no_exchange(key_i,key_j,Nint,hij)
   integer                        :: n_occ_ab(2)
   logical                        :: has_mipi(Nint*bit_kind_size)
   double precision               :: mipi(Nint*bit_kind_size)
+  double precision               :: diag_H_mat_elem
   PROVIDE mo_bielec_integrals_in_map mo_integrals_map
   
   ASSERT (Nint > 0)
@@ -771,9 +774,11 @@ subroutine i_H_j_dyall_no_exchange(key_i,key_j,Nint,hij)
         
       endif
       hij = phase*(hij + mo_mono_elec_integral(m,p)  + fock_operator_active_from_core_inact(m,p) )
+!     hij = phase*(hij + mo_mono_elec_integral(m,p) ) 
       
     case (0)
       hij = diag_H_mat_elem_no_elec_check_no_exchange(key_i,Nint)
+!     hij = diag_H_mat_elem(key_i,Nint)
 !     hij = 0.d0
   end select
 end
@@ -799,7 +804,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
   ! alpha - alpha 
   do i = 1, elec_num_tab_local(1)
    iorb =  occ(i,1)
-   diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb)
+   diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb)   !+ fock_operator_active_from_core_inact(iorb,iorb)
    do j = i+1, elec_num_tab_local(1)
     jorb = occ(j,1)
     diag_H_mat_elem_no_elec_check_no_exchange +=  mo_bielec_integral_jj(jorb,iorb)
@@ -809,7 +814,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
   ! beta - beta   
   do i = 1, elec_num_tab_local(2)
    iorb =  occ(i,2)
-   diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb)
+   diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb)   !+ fock_operator_active_from_core_inact(iorb,iorb)
    do j = i+1, elec_num_tab_local(2)
     jorb = occ(j,2)
     diag_H_mat_elem_no_elec_check_no_exchange +=  mo_bielec_integral_jj(jorb,iorb)
@@ -835,7 +840,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
    do j = 1, n_core_inact_orb
     jorb = list_core_inact(j)
     diag_H_mat_elem_no_elec_check_no_exchange +=  2.d0 * mo_bielec_integral_jj(jorb,iorb) 
-    core_act_exchange(1) += - mo_bielec_integral_jj_exchange(jorb,iorb)
+!   core_act_exchange(1) += - mo_bielec_integral_jj_exchange(jorb,iorb)
 !   diag_H_mat_elem_no_elec_check_no_exchange += core_act_exchange(1) 
    enddo
   enddo 
@@ -846,7 +851,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
    do j = 1, n_core_inact_orb
     jorb = list_core_inact(j)
     diag_H_mat_elem_no_elec_check_no_exchange +=  2.d0 * mo_bielec_integral_jj(jorb,iorb) 
-    core_act_exchange(2) += - mo_bielec_integral_jj_exchange(jorb,iorb)
+!   core_act_exchange(2) += - mo_bielec_integral_jj_exchange(jorb,iorb)
 !   diag_H_mat_elem_no_elec_check_no_exchange += core_act_exchange(2) 
    enddo
   enddo 
@@ -884,3 +889,45 @@ subroutine u0_H_dyall_u0_no_exchange(energies,psi_in,psi_in_coef,ndet,dim_psi_in
  energies(state_target) = accu
  deallocate(psi_coef_tmp)
 end
+
+
+
+!subroutine u0_H_dyall_u0_no_exchange_bis(energies,psi_in,psi_in_active,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target)
+subroutine u0_H_dyall_u0_no_exchange_bis(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target)
+  use bitmasks
+ implicit none
+ integer, intent(in) :: N_states_in,ndet,dim_psi_in,dim_psi_coef,state_target
+!integer(bit_kind), intent(in) :: psi_in(N_int,2,dim_psi_in),psi_in_active(N_int,2,dim_psi_in)
+ integer(bit_kind), intent(in) :: psi_in(N_int,2,dim_psi_in)
+ double precision,  intent(in) :: psi_in_coef(dim_psi_coef,N_states_in)
+ double precision,  intent(out) :: energies(N_states_in)
+ 
+ integer :: i,j 
+ double precision :: hij,accu
+ energies = 0.d0
+ accu = 0.d0
+ double precision, allocatable :: psi_coef_tmp(:)
+ allocate(psi_coef_tmp(ndet))
+ 
+ do i = 1, ndet
+  psi_coef_tmp(i) = psi_in_coef(i,state_target)
+ enddo
+
+ double precision :: hij_bis,diag_H_mat_elem
+ do i = 1, ndet
+  if(psi_coef_tmp(i)==0.d0)cycle
+  do j = i+1, ndet
+   if(psi_coef_tmp(j)==0.d0)cycle
+!  call i_H_j_dyall_no_exchange(psi_in_active(1,1,i),psi_in_active(1,1,j),N_int,hij)
+   call i_H_j(psi_in(1,1,i),psi_in(1,1,j),N_int,hij)
+   accu += 2.d0 * psi_coef_tmp(i) * psi_coef_tmp(j) * hij
+  enddo
+ enddo
+ do i = 1, N_det
+  if(psi_coef_tmp(i)==0.d0)cycle
+  accu += psi_coef_tmp(i) * psi_coef_tmp(i) * diag_H_mat_elem(psi_in(1,1,i),N_int)
+ enddo
+ energies(state_target) = accu
+ deallocate(psi_coef_tmp)
+end
+
diff --git a/plugins/MRPT_Utils/fock_like_operators.irp.f b/plugins/MRPT_Utils/fock_like_operators.irp.f
index d4ce0661..f16aba26 100644
--- a/plugins/MRPT_Utils/fock_like_operators.irp.f
+++ b/plugins/MRPT_Utils/fock_like_operators.irp.f
@@ -197,7 +197,7 @@
     k_inact_core_orb = list_core_inact(k)
     coulomb  = get_mo_bielec_integral(k_inact_core_orb,iorb,k_inact_core_orb,jorb,mo_integrals_map)
     exchange = get_mo_bielec_integral(k_inact_core_orb,jorb,iorb,k_inact_core_orb,mo_integrals_map)
-    accu += 2.d0 * coulomb - exchange
+    accu += 2.d0 * coulomb  - exchange
    enddo
    fock_operator_active_from_core_inact(iorb,jorb) = accu
   enddo
diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f
index 91bb7a54..9699a1df 100644
--- a/plugins/MRPT_Utils/mrpt_dress.irp.f
+++ b/plugins/MRPT_Utils/mrpt_dress.irp.f
@@ -122,7 +122,7 @@ subroutine mrpt_dress(delta_ij_,  Ndet,i_generator,n_selected,det_buffer,Nint,ip
        enddo
       else 
        call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
-!      !!!!!!!!!!!!! SHIFTED BK 
+       !!!!!!!!!!!!! SHIFTED BK 
 !      double precision :: hjj
 !      call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj)
 !      delta_e(1) = CI_electronic_energy(1) - hjj
@@ -141,7 +141,11 @@ subroutine mrpt_dress(delta_ij_,  Ndet,i_generator,n_selected,det_buffer,Nint,ip
       do j = 1, idx_alpha(0)
        index_j = idx_alpha(j)
        !!!!!!!!!!!!!!!!!! WARNING TEST
-       if(index_j .ne. index_i)cycle
+       !!!!!!!!!!!!!!!!!! WARNING TEST
+!      if(index_j .ne. index_i)cycle
+       !!!!!!!!!!!!!!!!!! WARNING TEST
+       !!!!!!!!!!!!!!!!!! WARNING TEST
+       !!!!!!!!!!!!!!!!!! WARNING TEST
        do i_state=1,N_states
         ! standard dressing first order
          delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_inv_array(index_j,i_state)
diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f
index 9d0d1fc6..ec6bbe50 100644
--- a/plugins/MRPT_Utils/psi_active_prov.irp.f
+++ b/plugins/MRPT_Utils/psi_active_prov.irp.f
@@ -14,7 +14,7 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)]
    psi_active(j,1,i) = iand(psi_ref(j,1,i),cas_bitmask(j,1,1))
    psi_active(j,2,i) = iand(psi_ref(j,2,i),cas_bitmask(j,1,1))
   enddo
-  call debug_det(psi_active(1,1,i),N_int)
+! call debug_det(psi_active(1,1,i),N_int)
  enddo
 END_PROVIDER
 
@@ -330,6 +330,7 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
    i_part =  list_virt_reverse(p1)
    do i_state = 1, N_states
     delta_e_act(i_state) += one_creat_virt(i_hole,i_part,i_state)
+!   delta_e_act += 1.d12
    enddo
   else if (degree == 2)then
    do i_state = 1, N_states
diff --git a/plugins/loc_cele/loc_cele.irp.f b/plugins/loc_cele/loc_cele.irp.f
index 2dda522e..67e74f08 100644
--- a/plugins/loc_cele/loc_cele.irp.f
+++ b/plugins/loc_cele/loc_cele.irp.f
@@ -101,10 +101,12 @@
        cmoref = 0.d0
        irot = 0
 
-       irot(1,1) = 5 
-       irot(2,1) = 6 
-       cmoref(6,1,1)   = 1d0 
-       cmoref(26,2,1)   = 1d0 
+       irot(1,1) = 14 
+       irot(2,1) = 15
+!      cmoref(6,1,1)   = 1.d0 
+!      cmoref(26,2,1)   = 1.d0 
+       cmoref(36,1,1)   = 1.d0 
+       cmoref(56,2,1)   = 1.d0 
 
 !       !!! H2O
 !      irot(1,1) = 4 
diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f
index 411fe703..53e31647 100644
--- a/src/Determinants/H_apply.irp.f
+++ b/src/Determinants/H_apply.irp.f
@@ -195,6 +195,7 @@ subroutine copy_H_apply_buffer_to_wf
   !call remove_duplicates_in_psi_det(found_duplicates)
 end
 
+
 subroutine remove_duplicates_in_psi_det(found_duplicates)
   implicit none
   logical, intent(out) :: found_duplicates
@@ -270,6 +271,81 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
   deallocate (duplicate,bit_tmp)
 end
 
+subroutine remove_duplicates_in_psi_det_new(found_duplicates)
+  implicit none
+  logical, intent(out) :: found_duplicates
+  BEGIN_DOC
+! Removes duplicate determinants in the wave function.
+  END_DOC
+  integer                        :: i,j,k
+  integer(bit_kind), allocatable :: bit_tmp(:)
+  logical,allocatable            :: duplicate(:)
+
+  allocate (duplicate(N_det), bit_tmp(N_det))
+
+  do i=1,N_det
+    integer, external            :: det_search_key
+    !$DIR FORCEINLINE
+    bit_tmp(i) = det_search_key(psi_det_sorted_bit(1,1,i),N_int)
+    duplicate(i) = .False.
+  enddo
+
+  do i=1,N_det-1
+    if (duplicate(i)) then
+      cycle
+    endif
+    j = i+1
+    do while (bit_tmp(j)==bit_tmp(i))
+      if (duplicate(j)) then
+        j += 1
+        if (j > N_det) then
+          exit
+        else
+          cycle
+        endif
+      endif
+      duplicate(j) = .True.
+      do k=1,N_int
+        if ( (psi_det_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,j) ) &
+        .or. (psi_det_sorted_bit(k,2,i) /= psi_det_sorted_bit(k,2,j) ) ) then
+          duplicate(j) = .False.
+          exit
+        endif
+      enddo
+      j += 1
+      if (j > N_det) then
+        exit
+      endif
+    enddo
+  enddo
+
+  found_duplicates = .False.
+  do i=1,N_det
+    if (duplicate(i)) then
+      found_duplicates = .True.
+      exit
+    endif
+  enddo
+
+  if (found_duplicates) then
+    k=0
+    do i=1,N_det
+      if (.not.duplicate(i)) then
+        k += 1
+        psi_det(:,:,k) = psi_det_sorted_bit (:,:,i)
+        psi_coef(k,:)  = psi_coef_sorted_bit(i,:)
+      else
+       psi_det(:,:,k) = 0_bit_kind
+       psi_coef(k,:) = 0.d0
+      endif
+    enddo
+    N_det = k
+    call write_bool(output_determinants,found_duplicates,'Found duplicate determinants')
+    SOFT_TOUCH N_det psi_det psi_coef
+  endif
+  deallocate (duplicate,bit_tmp)
+end
+
 
 subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
   use bitmasks
diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f
index 6bafa287..56590d9c 100644
--- a/src/Determinants/density_matrix.irp.f
+++ b/src/Determinants/density_matrix.irp.f
@@ -134,7 +134,6 @@ END_PROVIDER
      !$OMP END CRITICAL
      deallocate(tmp_a,tmp_b)
      !$OMP END PARALLEL
-
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ]
diff --git a/src/Determinants/truncate_wf.irp.f b/src/Determinants/truncate_wf.irp.f
index aba16fa7..49b5e70a 100644
--- a/src/Determinants/truncate_wf.irp.f
+++ b/src/Determinants/truncate_wf.irp.f
@@ -1,8 +1,52 @@
 program s2_eig_restart
  implicit none
  read_wf = .True.
- call routine
+ call routine_2
 end
+
+subroutine routine_2
+ implicit none
+ integer :: i,j,k,l
+ use bitmasks
+ integer :: n_det_restart,degree
+ integer(bit_kind),allocatable :: psi_det_tmp(:,:,:)
+ double precision ,allocatable :: psi_coef_tmp(:,:),accu(:)
+ integer, allocatable :: index_restart(:)
+ allocate(index_restart(N_det))
+ print*, 'How many Slater determinants would ou like ?'
+ read(5,*)N_det_restart
+ do i = 1, N_det_restart
+  index_restart(i) = i
+ enddo
+ allocate (psi_det_tmp(N_int,2,N_det_restart),psi_coef_tmp(N_det_restart,N_states),accu(N_states))
+ accu = 0.d0
+ do i = 1, N_det_restart
+  do j = 1, N_int
+   psi_det_tmp(j,1,i) = psi_det(j,1,index_restart(i))
+   psi_det_tmp(j,2,i) = psi_det(j,2,index_restart(i))
+  enddo
+  do j = 1,N_states
+   psi_coef_tmp(i,j) = psi_coef(index_restart(i),j)
+   accu(j) += psi_coef_tmp(i,j) * psi_coef_tmp(i,j)
+  enddo
+ enddo
+ do j = 1, N_states
+  accu(j) = 1.d0/dsqrt(accu(j))
+ enddo
+ do j = 1,N_states
+  do i = 1, N_det_restart
+   psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j)
+  enddo
+ enddo
+ call save_wavefunction_general(N_det_restart,N_states,psi_det_tmp,N_det_restart,psi_coef_tmp)
+
+ deallocate (psi_det_tmp,psi_coef_tmp,accu,index_restart)
+
+
+
+end
+
+
 subroutine routine
  implicit none
   call make_s2_eigenfunction