diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f
index 66f4975a..49aff533 100644
--- a/plugins/MRCC_Utils/davidson.irp.f
+++ b/plugins/MRCC_Utils/davidson.irp.f
@@ -364,6 +364,38 @@ end
 
 
 
+subroutine u0_H_u_0_mrcc(e_0,u_0,n,keys_tmp,Nint,istate)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Computes e_0 = <u_0|H|u_0>/<u_0|u_0>
+  !
+  ! n : number of determinants
+  !
+  END_DOC
+  integer, intent(in)            :: n,Nint
+  double precision, intent(out)  :: e_0
+  double precision, intent(in)   :: u_0(n)
+  integer(bit_kind),intent(in)   :: keys_tmp(Nint,2,n)
+  integer,intent(in)             :: istate
+  
+  double precision               :: H_jj(n)
+  double precision               :: v_0(n)
+  double precision               :: u_dot_u,u_dot_v,diag_H_mat_elem
+  integer :: i,j
+  do i = 1, n
+   H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint)
+  enddo
+
+  do i=1,N_det_ref
+    H_jj(idx_ref(i)) +=  delta_ii(istate,i)
+  enddo
+  
+  call H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
+  e_0 = u_dot_v(v_0,u_0,n)/u_dot_u(u_0,n)
+end
+
+
 subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
   use bitmasks
   implicit none
diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f
index 1d12d569..b6dd2bd1 100644
--- a/plugins/MRCC_Utils/mrcc_utils.irp.f
+++ b/plugins/MRCC_Utils/mrcc_utils.irp.f
@@ -119,9 +119,23 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ]
    implicit none
    BEGIN_DOC
-   ! Eigenvectors/values of the CI matrix
+   ! Eigenvectors/values of the dressed CI matrix
    END_DOC
-   integer                        :: i,j
+   double precision               :: ovrlp,u_dot_v
+   integer                        :: i_good_state
+   integer, allocatable           :: index_good_state_array(:)
+   logical, allocatable           :: good_state_array(:)
+   double precision, allocatable  :: s2_values_tmp(:)
+   integer                        :: i_other_state
+   double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:)
+   integer                        :: i_state
+   double precision               :: s2,e_0
+   integer                        :: i,j,k
+   double precision, allocatable  :: s2_eigvalues(:)
+   double precision, allocatable  :: e_array(:)
+   integer, allocatable           :: iorder(:)
+
+   integer, parameter :: mrcc_state = 1
    
    do j=1,N_states_diag
      do i=1,N_det
@@ -131,14 +145,17 @@ END_PROVIDER
    
    if (diag_algorithm == "Davidson") then
      
-     integer                        :: istate
-     istate = 1
-     call davidson_diag_mrcc(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed,&
-         size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_determinants,istate)
+     do i_state=1,N_states
+        call davidson_diag_mrcc(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed,&
+          size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_determinants,mrcc_state)
+     enddo
+     do j=1,N_states_diag
+       call get_s2_u0(psi_det,CI_eigenvectors_dressed(1,j),N_det,size(CI_eigenvectors_dressed,1),CI_eigenvectors_s2_dressed(j))
+     enddo
+
      
    else if (diag_algorithm == "Lapack") then
      
-     double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:)
      allocate (eigenvectors(size(H_matrix_dressed,1),N_det))
      allocate (eigenvalues(N_det))
      call lapack_diag(eigenvalues,eigenvectors,                      &
@@ -147,8 +164,6 @@ END_PROVIDER
      do i=1,N_det
        CI_eigenvectors_dressed(i,1) = eigenvectors(i,1)
      enddo
-     integer                        :: i_state
-     double precision               :: s2
      i_state = 0
      if (s2_eig) then
        do j=1,N_det
@@ -179,6 +194,98 @@ END_PROVIDER
      deallocate(eigenvectors,eigenvalues)
    endif
    
+   if(diagonalize_s2.and.n_states_diag > 1.and. n_det >= n_states_diag)then
+     ! Diagonalizing S^2 within the "n_states_diag" states found
+     allocate(s2_eigvalues(N_states_diag))
+     call diagonalize_s2_betweenstates(psi_det,CI_eigenvectors_dressed,n_det,size(psi_det,3),size(CI_eigenvectors_dressed,1),min(n_states_diag,n_det),s2_eigvalues)
+
+     do j = 1, N_states_diag
+       do i = 1, N_det
+         psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
+       enddo
+     enddo
+
+     if(s2_eig)then
+
+       ! Browsing the "n_states_diag" states and getting the lowest in energy "n_states" ones that have the S^2 value
+       ! closer to the "expected_s2" set as input
+
+       allocate(index_good_state_array(N_det),good_state_array(N_det))
+       good_state_array = .False.
+       i_state = 0
+       do j = 1, N_states_diag
+         if(dabs(s2_eigvalues(j)-expected_s2).le.0.3d0)then
+           good_state_array(j) = .True.
+           i_state +=1
+           index_good_state_array(i_state) = j
+         endif
+       enddo
+       ! Sorting the i_state good states by energy
+       allocate(e_array(i_state),iorder(i_state))
+       do j = 1, i_state
+         do i = 1, N_det
+           CI_eigenvectors_dressed(i,j) = psi_coef(i,index_good_state_array(j))
+         enddo
+         CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j))
+         call u0_H_u_0_mrcc(e_0,CI_eigenvectors_dressed(1,j),n_det,psi_det,N_int,mrcc_state)
+         CI_electronic_energy_dressed(j) = e_0
+         e_array(j) = e_0
+         iorder(j) = j
+       enddo
+       call dsort(e_array,iorder,i_state)
+       do j = 1, i_state
+         CI_electronic_energy_dressed(j) = e_array(j)
+         CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(iorder(j)))
+         do i = 1, N_det
+           CI_eigenvectors_dressed(i,j) = psi_coef(i,index_good_state_array(iorder(j)))
+         enddo
+         !    call u0_H_u_0_mrcc(e_0,CI_eigenvectors_dressed(1,j),n_det,psi_det,N_int,mrcc_state)
+         !    print*,'e    = ',CI_electronic_energy_dressed(j)
+         !    print*,'<e>  = ',e_0
+         !    call get_s2_u0(psi_det,CI_eigenvectors_dressed(1,j),N_det,size(CI_eigenvectors_dressed,1),s2)
+         !    print*,'s^2  = ',CI_eigenvectors_s2_dressed(j)
+         !    print*,'<s^2>= ',s2
+       enddo
+       deallocate(e_array,iorder)
+
+       ! Then setting the other states without any specific energy order
+       i_other_state = 0
+       do j = 1, N_states_diag
+         if(good_state_array(j))cycle
+         i_other_state +=1
+         do i = 1, N_det
+           CI_eigenvectors_dressed(i,i_state + i_other_state) = psi_coef(i,j)
+         enddo
+         CI_eigenvectors_s2_dressed(i_state + i_other_state) = s2_eigvalues(j)
+         call u0_H_u_0_mrcc(e_0,CI_eigenvectors_dressed(1,i_state + i_other_state),n_det,psi_det,N_int,mrcc_state)
+         CI_electronic_energy_dressed(i_state + i_other_state) = e_0
+       enddo
+       deallocate(index_good_state_array,good_state_array)
+
+     else
+
+       ! Sorting the N_states_diag by energy, whatever the S^2 value is
+
+       allocate(e_array(n_states_diag),iorder(n_states_diag))
+       do j = 1, N_states_diag
+         call u0_H_u_0_mrcc(e_0,CI_eigenvectors_dressed(1,j),n_det,psi_det,N_int,mrcc_state)
+         e_array(j) = e_0
+         iorder(j) = j
+       enddo
+       call dsort(e_array,iorder,n_states_diag)
+       do j = 1, N_states_diag
+         CI_electronic_energy_dressed(j) = e_array(j)
+         do i = 1, N_det
+           CI_eigenvectors_dressed(i,j) = psi_coef(i,iorder(j))
+         enddo
+         CI_eigenvectors_s2_dressed(j) = s2_eigvalues(iorder(j))
+       enddo
+       deallocate(e_array,iorder)
+     endif
+     deallocate(s2_eigvalues)
+
+   endif
+
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]