diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f
index d8b0a2c3..8e0af39e 100644
--- a/plugins/MRCC_Utils/davidson.irp.f
+++ b/plugins/MRCC_Utils/davidson.irp.f
@@ -678,14 +678,6 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
   integer, external :: align_double
   sze_8 = align_double(sze)
 
-  double precision :: delta
-
-  if (s2_eig) then
-    delta = 1.d0
-  else
-    delta = 0.d0
-  endif
-  
   itermax = min(davidson_sze_max, sze/N_st_diag)
   allocate(                                                          &
       W(sze_8,N_st_diag*itermax),                                    &
@@ -722,24 +714,17 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
   converged = .False.
   
   double precision               :: r1, r2
-  do k=N_st+1,N_st_diag-2,2
+  do k=N_st+1,N_st_diag
       do i=1,sze
         call random_number(r1)
         call random_number(r2)
         r1 = dsqrt(-2.d0*dlog(r1))
         r2 = dtwo_pi*r2
         u_in(i,k) = r1*dcos(r2)
-        u_in(i,k+1) = r1*dsin(r2)
       enddo
   enddo
-  do k=N_st_diag-1,N_st_diag
-      do i=1,sze
-        call random_number(r1)
-        call random_number(r2)
-        r1 = dsqrt(-2.d0*dlog(r1))
-        r2 = dtwo_pi*r2
-        u_in(i,k) = r1*dcos(r2)
-      enddo
+  do k=1,N_st_diag
+      call normalize(u_in(1,k),sze)
   enddo
 
   
@@ -796,14 +781,36 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
         s2(k) = s_(k,k) + S_z2_Sz
       enddo
 
+    if (s2_eig) then
+        logical :: state_ok(N_st_diag*davidson_sze_max)
+        do k=1,shift2
+          state_ok(k) = (dabs(s2(k)-expected_s2) < 0.6d0)
+        enddo
+        do k=1,shift2
+          if (.not. state_ok(k)) then
+            do l=k+1,shift2
+              if (state_ok(l)) then
+                call dswap(shift2, y(1,k), 1, y(1,l), 1)
+                call dswap(1, s2(k), 1, s2(l), 1)
+                call dswap(1, lambda(k), 1, lambda(l), 1)
+                state_ok(k) = .True.
+                state_ok(l) = .False.
+                exit
+              endif
+            enddo
+          endif
+        enddo
+      endif
+
+
       ! Compute overlap with U_in
       ! -------------------------
       
-      integer :: coord(2), order(N_st)
+      integer :: coord(2), order(N_st_diag)
       overlap = -1.d0
-      do k=1,N_st
+      do k=1,shift2 
         do i=1,shift2
-          overlap(i,k) = dabs(y(i,k))
+          overlap(k,i) = dabs(y(k,i))
         enddo
       enddo
       do k=1,N_st
diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f
index 8ad922cf..f28ccf25 100644
--- a/plugins/MRCC_Utils/mrcc_utils.irp.f
+++ b/plugins/MRCC_Utils/mrcc_utils.irp.f
@@ -924,6 +924,9 @@ END_PROVIDER
 
      norm = norm*f
      print *,  'norm of |T Psi_0> = ', dsqrt(norm)
+     if (dsqrt(norm) > 1.d0) then
+       stop 'Error : Norm of the SD larger than the norm of the reference.'
+     endif
 
      do i=1,N_det_ref
        norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
diff --git a/plugins/mrcepa0/mrcc.irp.f b/plugins/mrcepa0/mrcc.irp.f
index 91592e62..a28d4be3 100644
--- a/plugins/mrcepa0/mrcc.irp.f
+++ b/plugins/mrcepa0/mrcc.irp.f
@@ -8,8 +8,16 @@ program mrsc2sub
   
   read_wf = .True.
   SOFT_TOUCH read_wf
-  call print_cas_coefs
   call set_generators_bitmasks_as_holes_and_particles
+  if (.True.) then
+    integer :: i,j
+    do j=1,N_states
+      do i=1,N_det
+        psi_coef(i,j) = CI_eigenvectors(i,j)
+      enddo
+    enddo
+    TOUCH psi_coef 
+  endif
   call run(N_states,energy)
   if(do_pt2_end)then
     call run_pt2(N_states,energy)
diff --git a/plugins/mrcepa0/mrcepa0.irp.f b/plugins/mrcepa0/mrcepa0.irp.f
index 34d3dec5..aeacbb39 100644
--- a/plugins/mrcepa0/mrcepa0.irp.f
+++ b/plugins/mrcepa0/mrcepa0.irp.f
@@ -8,8 +8,18 @@ program mrcepa0
   
   read_wf = .True.
   SOFT_TOUCH read_wf
-  call print_cas_coefs
   call set_generators_bitmasks_as_holes_and_particles
+  if (.True.) then
+    integer :: i,j
+    do j=1,N_states
+      do i=1,N_det
+        psi_coef(i,j) = CI_eigenvectors(i,j)
+      enddo
+    enddo
+    TOUCH psi_coef 
+  endif
+  call print_cas_coefs
+  
   call run(N_states,energy)
   if(do_pt2_end)then
     call run_pt2(N_states,energy)
diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f
index 25fe5f53..09c35e52 100644
--- a/plugins/mrcepa0/mrcepa0_general.irp.f
+++ b/plugins/mrcepa0/mrcepa0_general.irp.f
@@ -17,7 +17,6 @@ subroutine run(N_st,energy)
   double precision, allocatable :: lambda(:)
   allocate (lambda(N_states))
 
-  
   thresh_mrcc = thresh_dressed_ci
   n_it_mrcc_max = n_it_max_dressed_ci
 
diff --git a/plugins/mrcepa0/mrsc2.irp.f b/plugins/mrcepa0/mrsc2.irp.f
index d0f44a33..948b1b5c 100644
--- a/plugins/mrcepa0/mrsc2.irp.f
+++ b/plugins/mrcepa0/mrsc2.irp.f
@@ -7,8 +7,16 @@ program mrsc2
   mrmode = 2
   read_wf = .True.
   SOFT_TOUCH read_wf
-  call print_cas_coefs
   call set_generators_bitmasks_as_holes_and_particles
+  if (.True.) then
+    integer :: i,j
+    do j=1,N_states
+      do i=1,N_det
+        psi_coef(i,j) = CI_eigenvectors(i,j)
+      enddo
+    enddo
+    TOUCH psi_coef 
+  endif
   call run(N_states,energy)
   if(do_pt2_end)then
     call run_pt2(N_states,energy)
diff --git a/src/Davidson/diagonalization_hs2.irp.f b/src/Davidson/diagonalization_hs2.irp.f
index fddac471..102dcfb8 100644
--- a/src/Davidson/diagonalization_hs2.irp.f
+++ b/src/Davidson/diagonalization_hs2.irp.f
@@ -183,24 +183,16 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
   converged = .False.
   
   double precision               :: r1, r2
-  do k=N_st+1,N_st_diag-2,2
+  do k=N_st+1,N_st_diag
       do i=1,sze
         call random_number(r1)
-        call random_number(r2)
         r1 = dsqrt(-2.d0*dlog(r1))
         r2 = dtwo_pi*r2
         u_in(i,k) = r1*dcos(r2)
-        u_in(i,k+1) = r1*dsin(r2)
       enddo
   enddo
-  do k=N_st_diag-1,N_st_diag
-      do i=1,sze
-        call random_number(r1)
-        call random_number(r2)
-        r1 = dsqrt(-2.d0*dlog(r1))
-        r2 = dtwo_pi*r2
-        u_in(i,k) = r1*dcos(r2)
-      enddo
+  do k=1,N_st_diag
+    call normalize(u_in(1,k),sze)
   enddo