diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f
index b3b2f427..3b05aaeb 100644
--- a/plugins/MRCC_Utils/mrcc_utils.irp.f
+++ b/plugins/MRCC_Utils/mrcc_utils.irp.f
@@ -149,26 +149,31 @@ END_PROVIDER
      
      allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)), &
      eigenvalues(size(CI_electronic_energy_dressed,1)))
+     do j=1,min(N_states,N_det)
+       do i=1,N_det
+         eigenvectors(i,j) = psi_coef(i,j)
+       enddo
+     enddo
      do mrcc_state=1,N_states
-      do i=1,N_det
-        eigenvectors(i,1) = psi_coef(i,mrcc_state)
+      do j=mrcc_state,min(N_states,N_det)
+        do i=1,N_det
+          eigenvectors(i,j) = psi_coef(i,j)
+        enddo
       enddo
       call davidson_diag_mrcc_HS2(psi_det,eigenvectors,&
             size(eigenvectors,1), &
-            eigenvalues,N_det,1,N_states_diag,N_int, &
+            eigenvalues,N_det,N_states,N_states_diag,N_int, &
             output_determinants,mrcc_state)
-      CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,1)
-      CI_electronic_energy_dressed(mrcc_state) = eigenvalues(1)
-      if (mrcc_state == 1) then
-        do k=N_states+1,N_states_diag
-          CI_eigenvectors_dressed(1:N_det,k) = eigenvectors(1:N_det,k)
-          CI_electronic_energy_dressed(k) = eigenvalues(k)
-        enddo
-      endif
-     enddo
-     call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
+      CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
+      CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
+   enddo
+   do k=N_states+1,N_states_diag
+     CI_eigenvectors_dressed(1:N_det,k) = eigenvectors(1:N_det,k)
+     CI_electronic_energy_dressed(k) = eigenvalues(k)
+   enddo
+   call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
           N_states_diag,size(CI_eigenvectors_dressed,1))
-     deallocate (eigenvectors,eigenvalues)
+   deallocate (eigenvectors,eigenvalues)
 
      
    else if (diag_algorithm == "Lapack") then
@@ -716,7 +721,7 @@ END_PROVIDER
     factor = 1.d0
     resold = huge(1.d0)
 
-    do k=0,hh_nex
+    do k=0,hh_nex*hh_nex
       !$OMP PARALLEL default(shared) private(cx, i, a_col, a_coll)
       
       !$OMP DO
@@ -751,15 +756,15 @@ END_PROVIDER
         X(a_col) = X_new(a_col)
       end do
       if (res > resold) then
-        factor = -factor * 0.5d0
+        factor = factor * 0.5d0
       endif
       resold = res
       
-      if(mod(k, 100) == 0) then
+      if(iand(k, 4095) == 0) then
         print *, "res ", k, res
       end if
       
-      if(res < 1d-9) exit
+      if(res < 1d-12) exit
     end do
     
     norm = 0.d0
diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f
index 1e89cc2c..d9607e6a 100644
--- a/plugins/mrcepa0/mrcepa0_general.irp.f
+++ b/plugins/mrcepa0/mrcepa0_general.irp.f
@@ -37,9 +37,9 @@ subroutine run(N_st,energy)
     lambda = 1.d0
     do while (delta_E > thresh_mrcc)
       iteration += 1
-      print *,  '===========================' 
-      print *,  'MRCEPA0 Iteration', iteration
-      print *,  '===========================' 
+      print *,  '===============================================' 
+      print *,  'MRCEPA0 Iteration', iteration, '/', n_it_mrcc_max
+      print *,  '===============================================' 
       print *,  ''
       E_old = sum(ci_energy_dressed(1:N_states))
       do i=1,N_st
@@ -48,6 +48,8 @@ subroutine run(N_st,energy)
       call diagonalize_ci_dressed(lambda)
       E_new = sum(ci_energy_dressed(1:N_states))
       delta_E = (E_new - E_old)/dble(N_states)
+      print *,  ''
+      call write_double(6,thresh_mrcc,"thresh_mrcc")
       call write_double(6,delta_E,"delta_E")
       delta_E = dabs(delta_E)
       call save_wavefunction
diff --git a/src/Davidson/diagonalization_hs2.irp.f b/src/Davidson/diagonalization_hs2.irp.f
index 2e05df48..8dc6e00d 100644
--- a/src/Davidson/diagonalization_hs2.irp.f
+++ b/src/Davidson/diagonalization_hs2.irp.f
@@ -326,32 +326,32 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
 
       if (state_following) then
 
-        integer                        :: coord(2), order(N_st_diag)
+        integer                        :: order(N_st_diag)
+        double precision               :: cmax
+
         overlap = -1.d0
-        do i=1,shift2
-          do k=1,shift2
+        do k=1,shift2
+          do i=1,shift2
             overlap(k,i) = dabs(y(k,i))
           enddo
         enddo
         do k=1,N_st
-          coord = maxloc(overlap)
-          order( coord(2) )  = coord(1)
-          do i=1,shift2
-            overlap(coord(1),i) = -1.d0
+          cmax = -1.d0
+          do i=1,N_st_diag
+            if (overlap(i,k) > cmax) then
+              cmax = overlap(i,k) 
+              order(k) = i
+            endif
+          enddo
+          do i=1,N_st_diag
+            overlap(order(k),i) = -1.d0
           enddo
         enddo
-        print *,  order(1:N_st)
-        do i=1,shift2
-         do k=1,shift2
-           overlap(k,i) = y(k,i)
-         enddo
-        enddo
+        overlap = y
         do k=1,N_st
           l = order(k)
           if (k /= l) then
-            do i=1,shift2
-              y(i,k) = overlap(i,l)
-            enddo
+            y(1:shift2,k) = overlap(1:shift2,l)
           endif
         enddo
         do k=1,N_st