From 4d5467218de414d76e81a5080f762519e49f54f4 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 31 Jul 2024 16:36:44 +0200
Subject: [PATCH] Fix state following when N_states_diag is too small

---
 .../diagonalization_hs2_dressed.irp.f         | 50 +++----------------
 src/davidson/diagonalize_ci.irp.f             |  5 +-
 2 files changed, 10 insertions(+), 45 deletions(-)

diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f
index fb04b29b..191e0021 100644
--- a/src/davidson/diagonalization_hs2_dressed.irp.f
+++ b/src/davidson/diagonalization_hs2_dressed.irp.f
@@ -594,6 +594,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
           state(l) = idx
         enddo
 
+        ! Check if all states are attributed. If not, exit and N_st_diag will be increased.
+        do l=1,N_st
+          if (state(l) == 0) then
+            return
+          endif
+        enddo
+
         ! tmp array before setting state_ok
         ok = .False.
         do l = 1, N_st
@@ -627,47 +634,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       ! Swapped eigenvectors
       prev_y = y
 
-!      if (state_following) then
-!
-!        overlap = -1.d0
-!        do k=1,shift2
-!          do i=1,shift2
-!            overlap(k,i) = dabs(y(k,i))
-!          enddo
-!        enddo
-!        do k=1,N_st
-!          cmax = -1.d0
-!          do i=1,N_st
-!            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
-!        overlap = y
-!        do k=1,N_st
-!          l = order(k)
-!          if (k /= l) then
-!            y(1:shift2,k) = overlap(1:shift2,l)
-!          endif
-!        enddo
-!        do k=1,N_st
-!          overlap(k,1) = lambda(k)
-!          overlap(k,2) = s2(k)
-!        enddo
-!        do k=1,N_st
-!          l = order(k)
-!          if (k /= l) then
-!            lambda(k) = overlap(l,1)
-!            s2(k) = overlap(l,2)
-!          endif
-!        enddo
-!
-!      endif
-
 
       ! Express eigenvectors of h in the determinant basis
       ! --------------------------------------------------
@@ -703,7 +669,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
 
 
       if ((itertot>1).and.(iter == 1)) then
-        !don't print
+        ! Don't print
         continue
       else
         write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,ES11.3))') iter-1, to_print(1:3,1:N_st)
diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f
index 59c8313a..6b852905 100644
--- a/src/davidson/diagonalize_ci.irp.f
+++ b/src/davidson/diagonalize_ci.irp.f
@@ -282,9 +282,8 @@ END_PROVIDER
          print*,'  Within the ',N_det,'determinants selected'
          print*,'  and the ',N_states_diag,'states requested'
          print*,'  We did not find only states with S^2 values close to ',expected_s2
-         print*,'  We will then set the first N_states eigenvectors of the H matrix'
-         print*,'  as the CI_eigenvectors'
-         print*,'  You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
+         print*,'  You should consider more states, or change s2_eig, or just enlarge the CI space'
+         print*,'!!!!!!!!   WARNING  !!!!!!!!!'
          print*,''
 
          do j=1,min(N_states_diag,N_det)