Merge pull request #325 from Ydrnan/dev-stable-add

state following

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -522,6 +522,84 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
         enddo
       endif
 
+      if (state_following) then
+         if (.not. only_expected_s2) then
+           print*,''
+           print*,'!!! State following only available with only_expected_s2 = .True. !!!'
+           STOP
+         endif
+       endif
+
+      if (state_following) then
+
+        integer ::  state(N_st), idx
+        double precision :: omax
+        logical :: used
+        logical, allocatable :: ok(:)
+        double precision, allocatable :: overlp(:,:)
+
+        allocate(overlp(shift2,N_st),ok(shift2))
+
+        overlp = 0d0
+        do j = 1, shift2-1, N_st_diag
+
+          ! Computes some states from the guess vectors
+          ! Psi(:,j:j+N_st_diag) = U y(:,j:j+N_st_diag) and put them
+          ! in U(1,shift2+1:shift2+1+N_st_diag) as temporary array
+          call dgemm('N','N', sze, N_st_diag, shift2,                    &
+          1.d0, U, size(U,1), y(1,j), size(y,1), 0.d0, U(1,shift2+1), size(U,1))
+
+          ! Overlap
+          do l = 1, N_st
+            do k = 1, N_st_diag
+              do i = 1, sze
+                overlp(k+j-1,l) += U(i,l) * U(i,shift2+k)
+              enddo
+            enddo
+          enddo
+
+        enddo
+
+        state = 0
+        do l = 1, N_st
+
+          omax = 0d0
+          idx = 0
+          do k = 1, shift2
+
+            ! Already used ?
+            used = .False.
+            do i = 1, N_st
+              if (state(i) == k) then
+                used = .True.
+              endif
+            enddo
+
+            ! Maximum overlap
+            if (dabs(overlp(k,l)) > omax .and. .not. used .and. state_ok(k)) then
+              omax = dabs(overlp(k,l))
+              idx = k
+            endif
+          enddo
+
+          state(l) = idx
+        enddo
+
+        ! tmp array before setting state_ok
+        ok = .False.
+        do l = 1, N_st
+          ok(state(l)) = .True.
+        enddo
+
+        do k = 1, shift2
+          if (.not. ok(k)) then
+            state_ok(k) = .False.
+          endif
+        enddo
+
+        deallocate(overlp,ok)
+      endif
+
       do k=1,shift2
         if (.not. state_ok(k)) then
           do l=k+1,shift2
@@ -537,46 +615,46 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
         endif
       enddo
 
-      if (state_following) then
+!      if (state_following) then
-
+!
-        overlap = -1.d0
+!        overlap = -1.d0
-        do k=1,shift2
+!        do k=1,shift2
-          do i=1,shift2
+!          do i=1,shift2
-            overlap(k,i) = dabs(y(k,i))
+!            overlap(k,i) = dabs(y(k,i))
-          enddo
+!          enddo
-        enddo
+!        enddo
-        do k=1,N_st
+!        do k=1,N_st
-          cmax = -1.d0
+!          cmax = -1.d0
-          do i=1,N_st
+!          do i=1,N_st
-            if (overlap(i,k) > cmax) then
+!            if (overlap(i,k) > cmax) then
-              cmax = overlap(i,k)
+!              cmax = overlap(i,k)
-              order(k) = i
+!              order(k) = i
-            endif
+!            endif
-          enddo
+!          enddo
-          do i=1,N_st_diag
+!          do i=1,N_st_diag
-            overlap(order(k),i) = -1.d0
+!            overlap(order(k),i) = -1.d0
-          enddo
+!          enddo
-        enddo
+!        enddo
-        overlap = y
+!        overlap = y
-        do k=1,N_st
+!        do k=1,N_st
-          l = order(k)
+!          l = order(k)
-          if (k /= l) then
+!          if (k /= l) then
-            y(1:shift2,k) = overlap(1:shift2,l)
+!            y(1:shift2,k) = overlap(1:shift2,l)
-          endif
+!          endif
-        enddo
+!        enddo
-        do k=1,N_st
+!        do k=1,N_st
-          overlap(k,1) = lambda(k)
+!          overlap(k,1) = lambda(k)
-          overlap(k,2) = s2(k)
+!          overlap(k,2) = s2(k)
-        enddo
+!        enddo
-        do k=1,N_st
+!        do k=1,N_st
-          l = order(k)
+!          l = order(k)
-          if (k /= l) then
+!          if (k /= l) then
-            lambda(k) = overlap(l,1)
+!            lambda(k) = overlap(l,1)
-            s2(k) = overlap(l,2)
+!            s2(k) = overlap(l,2)
-          endif
+!          endif
-        enddo
+!        enddo
-
+!
-      endif
+!      endif
 
 
       ! Express eigenvectors of h in the determinant basis

src/davidson/diagonalize_ci.irp.f

@@ -123,6 +123,7 @@ END_PROVIDER
         endif
 
      enddo
+
      if (N_states_diag > N_states_diag_save) then
        N_states_diag = N_states_diag_save
        TOUCH N_states_diag
@@ -133,24 +134,101 @@ END_PROVIDER
      print *,  'Diagonalization of H using Lapack'
      allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
      allocate (eigenvalues(N_det))
+
      if (s2_eig) then
+
        double precision, parameter :: alpha = 0.1d0
        allocate (H_prime(N_det,N_det) )
+
        H_prime(1:N_det,1:N_det) = H_matrix_all_dets(1:N_det,1:N_det) +  &
          alpha * S2_matrix_all_dets(1:N_det,1:N_det)
+
        do j=1,N_det
          H_prime(j,j) = H_prime(j,j) - alpha*expected_s2
        enddo
+
        call lapack_diag(eigenvalues,eigenvectors,H_prime,size(H_prime,1),N_det)
        call nullify_small_elements(N_det,N_det,eigenvectors,size(eigenvectors,1),1.d-12)
+
        CI_electronic_energy(:) = 0.d0
        i_state = 0
+
        allocate (s2_eigvalues(N_det))
        allocate(index_good_state_array(N_det),good_state_array(N_det))
+
        good_state_array = .False.
        call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
          N_det,size(eigenvectors,1))
-       if (only_expected_s2) then
+
+       if (state_following) then
+         if (.not. only_expected_s2) then
+           print*,''
+           print*,'!!! State following only available with only_expected_s2 = .True. !!!'
+           STOP
+         endif
+         if (N_det < N_states) then
+           print*,''
+           print*,'!!! State following requires at least N_states determinants to be activated !!!'
+           STOP
+         endif
+       endif
+
+       if (state_following .and. only_expected_s2) then
+
+         integer :: state(N_states), idx,l
+         double precision :: omax
+         double precision, allocatable :: overlp(:)
+         logical :: used
+         logical, allocatable :: ok(:)
+
+         allocate(overlp(N_det), ok(N_det))
+
+         i_state = 0
+         state = 0
+         do l = 1, N_states
+
+           ! Overlap wrt each state
+           overlp = 0d0
+           do k = 1, N_det
+             do i = 1, N_det
+               overlp(k) = overlp(k) + psi_coef(i,l) * eigenvectors(i,k)
+             enddo
+           enddo
+
+           ! Idx of the state with the maximum overlap not already "used"
+           omax = 0d0
+           idx = 0
+           do k = 1, N_det
+
+             ! Already used ?
+             used = .False.
+             do i = 1, N_states
+               if (state(i) == k) then
+                 used = .True.
+               endif
+             enddo
+
+             ! Maximum overlap
+             if (dabs(overlp(k)) > omax .and. .not. used) then
+               if (dabs(s2_eigvalues(k)-expected_s2) > 0.5d0) cycle
+               omax = dabs(overlp(k))
+               idx = k
+             endif
+           enddo
+
+           state(l) = idx
+           i_state +=1
+         enddo
+
+         deallocate(overlp, ok)
+
+         do i = 1, i_state
+           index_good_state_array(i) = state(i)
+           good_state_array(i) = .True.
+         enddo
+
+       else if (only_expected_s2) then
+
          do j=1,N_det
            ! Select at least n_states states with S^2 values closed to "expected_s2"
            if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
@@ -158,17 +236,23 @@ END_PROVIDER
              index_good_state_array(i_state) = j
              good_state_array(j) = .True.
            endif
+
            if(i_state.eq.N_states) then
              exit
            endif
          enddo
+
        else
+
          do j=1,N_det
            index_good_state_array(j) = j
            good_state_array(j) = .True.
          enddo
+
        endif
+
        if(i_state .ne.0)then
+
          ! Fill the first "i_state" states that have a correct S^2 value
          do j = 1, i_state
            do i=1,N_det
@@ -177,6 +261,7 @@ END_PROVIDER
            CI_electronic_energy(j) = eigenvalues(index_good_state_array(j))
            CI_s2(j) = s2_eigvalues(index_good_state_array(j))
          enddo
+
          i_other_state = 0
          do j = 1, N_det
            if(good_state_array(j))cycle
@@ -201,6 +286,7 @@ END_PROVIDER
          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*,''
+
          do j=1,min(N_states_diag,N_det)
            do i=1,N_det
              CI_eigenvectors(i,j) = eigenvectors(i,j)
@@ -209,14 +295,18 @@ END_PROVIDER
            CI_s2(j) = s2_eigvalues(j)
          enddo
        endif
+
        deallocate(index_good_state_array,good_state_array)
        deallocate(s2_eigvalues)
+
      else
+
        call lapack_diag(eigenvalues,eigenvectors,                    &
            H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
        CI_electronic_energy(:) = 0.d0
        call u_0_S2_u_0(CI_s2,eigenvectors,N_det,psi_det,N_int,       &
            min(N_det,N_states_diag),size(eigenvectors,1))
+
        ! Select the "N_states_diag" states of lowest energy
        do j=1,min(N_det,N_states_diag)
          do i=1,N_det
@@ -224,7 +314,9 @@ END_PROVIDER
          enddo
          CI_electronic_energy(j) = eigenvalues(j)
        enddo
+
      endif
+
      do k=1,N_states_diag
        CI_electronic_energy(k) = 0.d0
        do j=1,N_det
@@ -235,6 +327,7 @@ END_PROVIDER
          enddo
        enddo
      enddo
+
      deallocate(eigenvectors,eigenvalues)
    endif