Cleaned Davidson

plugins/MRCC_Utils/davidson.irp.f

@@ -640,8 +640,10 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
   include 'constants.include.F'
   
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, S, y, h, lambda
-  if (N_st_diag > sze) then
+  if (N_st_diag*3 > sze) then
-     stop 'error in Davidson : N_st_diag > sze'
+     print *,  'error in Davidson :'
+     print *,  'Increase n_det_max_jacobi to ', N_st_diag*3
+     stop -1
   endif
 
   PROVIDE nuclear_repulsion
@@ -763,10 +765,11 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
           1.d0, U, size(U,1), S, size(S,1),                &
           0.d0, s_, size(s_,1))
 
-      ! Diagonalize S^2
+!      ! Diagonalize S^2
-      ! ---------------
+!      ! ---------------
-      call lapack_diag(s2,y,s_,size(s_,1),shift2)
+!
-
+!      call lapack_diag(s2,y,s_,size(s_,1),shift2)
+!
 !      ! Rotate H in the basis of eigenfunctions of s2
 !      ! ---------------------------------------------
 !
@@ -823,7 +826,7 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
       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.3d0)
+            state_ok(k) = (dabs(s2(k)-expected_s2) < 0.6d0)
           enddo
       else
         state_ok(k) = .True.
@@ -844,6 +847,8 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
         endif
       enddo
 
+      if (state_following) then
+
         ! Compute overlap with U_in
         ! -------------------------
         
@@ -878,6 +883,8 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
           endif
         enddo
         
+      endif
+
 
       ! Express eigenvectors of h in the determinant basis
       ! --------------------------------------------------
@@ -940,22 +947,18 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
 
     enddo
 
-    if (.not.converged) then
-      iter = itermax-1
-    endif
-
     ! Re-contract to u_in
     ! -----------
     
-    do k=1,N_st_diag
-      energies(k) = lambda(k)
-    enddo
-
     call dgemm('N','N', sze, N_st_diag, shift2,                    &
         1.d0, U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
 
   enddo
 
+  do k=1,N_st_diag
+    energies(k) = lambda(k)
+  enddo
+
   write_buffer = '===== '
   do i=1,N_st
     write_buffer = trim(write_buffer)//' ================ =========== ==========='
@@ -966,7 +969,7 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
 
   deallocate (                                                       &
       W, residual_norm,                                              &
-      U,                                                             &
+      U, overlap,                                                    &
       c, S,                                                       &
       h,                                                             &
       y, s_, s_tmp,                                                  &

src/Davidson/EZFIO.cfg

@@ -15,3 +15,16 @@ type: Strictly_positive_int
 doc: Number of micro-iterations before re-contracting
 default: 10
 interface: ezfio,provider,ocaml
+
+[state_following]
+type: logical
+doc: If true, the states are re-ordered to match the input states
+default: False
+interface: ezfio,provider,ocaml
+
+[disk_based_davidson]
+type: logical
+doc: If true, disk space is used to store the vectors
+default: False
+interface: ezfio,provider,ocaml
+

src/Davidson/diagonalization_hs2.irp.f

@@ -45,8 +45,11 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
   !$OMP END DO 
   !$OMP END PARALLEL
 
+  if (disk_based_davidson) then
+    call davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit)
+  else
     call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit)
-!  call davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit)
+  endif
   do i=1,N_st_diag
     s2_out(i) = S2_jj(i)
   enddo
@@ -99,7 +102,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
   integer                        :: k_pairs, kl
   
   integer                        :: iter2
-  double precision, allocatable  :: W(:,:),  U(:,:), S(:,:)
+  double precision, allocatable  :: W(:,:),  U(:,:), S(:,:), overlap(:,:)
   double precision, allocatable  :: y(:,:), h(:,:), lambda(:), s2(:)
   double precision, allocatable  :: c(:), s_(:,:), s_tmp(:,:)
   double precision               :: diag_h_mat_elem
@@ -108,6 +111,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
   double precision               :: to_print(3,N_st)
   double precision               :: cpu, wall
   integer                        :: shift, shift2, itermax
+  double precision               :: r1, r2
+  logical                        :: state_ok(N_st_diag*davidson_sze_max)
   include 'constants.include.F'
   
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, S, y, h, lambda
@@ -137,7 +142,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
   write(iunit,'(A)') trim(write_buffer)
   write_buffer = ' Iter'
   do i=1,N_st
-    write_buffer = trim(write_buffer)//'      Energy          S^2      Residual'
+    write_buffer = trim(write_buffer)//'      Energy          S^2      Residual  '
   enddo
   write(iunit,'(A)') trim(write_buffer)
   write_buffer = '===== '
@@ -161,15 +166,16 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
       residual_norm(N_st_diag),                                      &
       c(N_st_diag*itermax),                                          &
       s2(N_st_diag*itermax),                                         &
+      overlap(N_st_diag*itermax, N_st_diag*itermax),                 &
       lambda(N_st_diag*itermax))
   
   h = 0.d0
-  s_ = 0.d0
-  s_tmp = 0.d0
   U = 0.d0
   W = 0.d0
   S = 0.d0
   y = 0.d0
+  s_ = 0.d0
+  s_tmp = 0.d0
 
 
   ASSERT (N_st > 0)
@@ -183,8 +189,8 @@ 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
+    u_in(k,k) = 10.d0
     do i=1,sze
       call random_number(r1)
       call random_number(r2)
@@ -233,8 +239,49 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
           0.d0, s_, size(s_,1))
 
 
+!      ! Diagonalize S^2
+!      ! ---------------
+!
+!      call lapack_diag(s2,y,s_,size(s_,1),shift2)
+!
+!
+!      ! Rotate H in the basis of eigenfunctions of s2
+!      ! ---------------------------------------------
+!
+!      call dgemm('N','N',shift2,shift2,shift2,                       &
+!          1.d0, h, size(h,1), y, size(y,1),                          &
+!          0.d0, s_tmp, size(s_tmp,1))
+!      
+!      call dgemm('T','N',shift2,shift2,shift2,                       &
+!          1.d0, y, size(y,1), s_tmp, size(s_tmp,1),                  &
+!          0.d0, h, size(h,1))
+!
+!      ! Damp interaction between different spin states
+!      ! ------------------------------------------------
+!
+!      do k=1,shift2
+!        do l=1,shift2
+!          if (dabs(s2(k) - s2(l)) > 1.d0) then
+!            h(k,l) = h(k,l)*(max(0.d0,1.d0 - dabs(s2(k) - s2(l))))
+!          endif
+!        enddo
+!      enddo
+!
+!      ! Rotate back H 
+!      ! -------------
+!
+!      call dgemm('N','T',shift2,shift2,shift2,                       &
+!          1.d0, h, size(h,1), y, size(y,1),                          &
+!          0.d0, s_tmp, size(s_tmp,1))
+!      
+!      call dgemm('N','N',shift2,shift2,shift2,                       &
+!          1.d0, y, size(y,1), s_tmp, size(s_tmp,1),                  &
+!          0.d0, h, size(h,1))
+
+      
       ! Diagonalize h
       ! -------------
+
       call lapack_diag(lambda,y,h,size(h,1),shift2)
       
       ! Compute S2 for each eigenvector
@@ -255,10 +302,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
       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
+      else
+        state_ok(k) = .True.
+      endif
+
       do k=1,shift2
         if (.not. state_ok(k)) then
           do l=k+1,shift2
@@ -273,6 +323,40 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
           enddo
         endif
       enddo
+
+      if (state_following) then
+
+        integer                        :: coord(2), order(N_st_diag)
+        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
+          coord = maxloc(overlap)
+          order( coord(2) )  = coord(1)
+          overlap(:,coord(2)) = -1.d0
+        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
 
 
@@ -290,11 +374,31 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
       ! -----------------------------------------
       
       do k=1,N_st_diag
+        if (state_ok(k)) then
           do i=1,sze
             U(i,shift2+k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) )      &
                 * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz &
               )/max(H_jj(i) - lambda (k),1.d-2)
           enddo
+        else
+          ! Randomize components with bad <S2>
+            do i=1,sze-2,2
+              call random_number(r1)
+              call random_number(r2)
+              r1 = dsqrt(-2.d0*dlog(r1))
+              r2 = dtwo_pi*r2
+              U(i,shift2+k) = r1*dcos(r2)
+              U(i+1,shift2+k) = r1*dsin(r2)
+            enddo
+            do i=sze-2+1,sze
+              call random_number(r1)
+              call random_number(r2)
+              r1 = dsqrt(-2.d0*dlog(r1))
+              r2 = dtwo_pi*r2
+              U(i,shift2+k) = r1*dcos(r2)
+            enddo
+        endif
+
         if (k <= N_st) then
           residual_norm(k) = u_dot_u(U(1,shift2+k),sze)
           to_print(1,k) = lambda(k) + nuclear_repulsion
@@ -339,7 +443,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
 
   deallocate (                                                       &
       W, residual_norm,                                              &
-      U,                                                             &
+      U, overlap,                                                    &
       c, S,                                                          &
       h,                                                             &
       y, s_, s_tmp,                                                  &
@@ -393,7 +497,7 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
   integer                        :: k_pairs, kl
   
   integer                        :: iter2
-  double precision, pointer      :: W(:,:),  U(:,:), S(:,:)
+  double precision, pointer      :: W(:,:),  U(:,:), S(:,:), overlap(:,:)
   double precision, allocatable  :: y(:,:), h(:,:), lambda(:), s2(:)
   double precision, allocatable  :: c(:), s_(:,:), s_tmp(:,:)
   double precision               :: diag_h_mat_elem
@@ -401,6 +505,7 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
   character*(16384)              :: write_buffer
   double precision               :: to_print(3,N_st)
   double precision               :: cpu, wall
+  logical                        :: state_ok(N_st_diag*davidson_sze_max)
   integer                        :: shift, shift2, itermax
   include 'constants.include.F'
   
@@ -431,7 +536,7 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
   write(iunit,'(A)') trim(write_buffer)
   write_buffer = ' Iter'
   do i=1,N_st
-    write_buffer = trim(write_buffer)//'      Energy          S^2      Residual'
+    write_buffer = trim(write_buffer)//'      Energy          S^2      Residual  '
   enddo
   write(iunit,'(A)') trim(write_buffer)
   write_buffer = '===== '
@@ -470,18 +575,19 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
       y(N_st_diag*itermax,N_st_diag*itermax),                        &
       s_(N_st_diag*itermax,N_st_diag*itermax),                       &
       s_tmp(N_st_diag*itermax,N_st_diag*itermax),                    &
+      overlap(N_st_diag*itermax, N_st_diag*itermax),                 &
       residual_norm(N_st_diag),                                      &
       c(N_st_diag*itermax),                                          &
       s2(N_st_diag*itermax),                                         &
       lambda(N_st_diag*itermax))
   
   h = 0.d0
-  s_ = 0.d0
-  s_tmp = 0.d0
   U = 0.d0
   W = 0.d0
   S = 0.d0
   y = 0.d0
+  s_ = 0.d0
+  s_tmp = 0.d0
   
   
   ASSERT (N_st > 0)
@@ -497,6 +603,7 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
   
   double precision               :: r1, r2
   do k=N_st+1,N_st_diag
+      u_in(k,k) = 10.d0
       do i=1,sze
         call random_number(r1)
         r1 = dsqrt(-2.d0*dlog(r1))
@@ -546,6 +653,45 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
             0.d0, s_(shift+1,1), size(s_,1))
       enddo
 
+!      ! Diagonalize S^2
+!      ! ---------------
+!
+!      call lapack_diag(s2,y,s_,size(s_,1),shift2)
+!
+!
+!      ! Rotate H in the basis of eigenfunctions of s2
+!      ! ---------------------------------------------
+!
+!      call dgemm('N','N',shift2,shift2,shift2,                       &
+!          1.d0, h, size(h,1), y, size(y,1),                          &
+!          0.d0, s_tmp, size(s_tmp,1))
+!      
+!      call dgemm('T','N',shift2,shift2,shift2,                       &
+!          1.d0, y, size(y,1), s_tmp, size(s_tmp,1),                  &
+!          0.d0, h, size(h,1))
+!
+!      ! Damp interaction between different spin states
+!      ! ------------------------------------------------
+!
+!      do k=1,shift2
+!        do l=1,shift2
+!          if (dabs(s2(k) - s2(l)) > 1.d0) then
+!            h(k,l) = h(k,l)*(max(0.d0,1.d0 - dabs(s2(k) - s2(l))))
+!          endif
+!        enddo
+!      enddo
+!
+!      ! Rotate back H 
+!      ! -------------
+!
+!      call dgemm('N','T',shift2,shift2,shift2,                       &
+!          1.d0, h, size(h,1), y, size(y,1),                          &
+!          0.d0, s_tmp, size(s_tmp,1))
+!      
+!      call dgemm('N','N',shift2,shift2,shift2,                       &
+!          1.d0, y, size(y,1), s_tmp, size(s_tmp,1),                  &
+!          0.d0, h, size(h,1))
+
 
       ! Diagonalize h
       ! -------------
@@ -568,11 +714,15 @@ subroutine davidson_diag_hjj_sjj_mmap(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
+      else
+        state_ok(k) = .True.
+      endif
+
       do k=1,shift2
         if (.not. state_ok(k)) then
           do l=k+1,shift2
@@ -587,17 +737,40 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
           enddo
         endif
       enddo
-        ! Randomize components with bad <S2>
+
-        if (.not. state_ok(k)) then
+      if (state_following) then
+
+        integer                        :: coord(2), order(N_st_diag)
+        overlap = -1.d0
+        do k=1,shift2
           do i=1,shift2
-              call random_number(r1)
+            overlap(k,i) = dabs(y(k,i))
-              call random_number(r2)
-              r1 = dsqrt(-2.d0*dlog(r1))
-              r2 = dtwo_pi*r2
-              y(i,k) = r1*dcos(r2)
-              lambda(k) = 1.d0
           enddo
+        enddo
+        do k=1,N_st
+          coord = maxloc(overlap)
+          order( coord(2) )  = coord(1)
+          overlap(:,coord(2)) = -1.d0
+        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
 
 
@@ -615,11 +788,31 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
       ! -----------------------------------------
       
       do k=1,N_st_diag
+        if (state_ok(k)) then
           do i=1,sze
             U(i,shift2+k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) )      &
                 * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz &
               )/max(H_jj(i) - lambda (k),1.d-2)
           enddo
+        else
+          ! Randomize components with bad <S2>
+            do i=1,sze-2,2
+              call random_number(r1)
+              call random_number(r2)
+              r1 = dsqrt(-2.d0*dlog(r1))
+              r2 = dtwo_pi*r2
+              U(i,shift2+k) = r1*dcos(r2)
+              U(i+1,shift2+k) = r1*dsin(r2)
+            enddo
+            do i=sze-2+1,sze
+              call random_number(r1)
+              call random_number(r2)
+              r1 = dsqrt(-2.d0*dlog(r1))
+              r2 = dtwo_pi*r2
+              U(i,shift2+k) = r1*dcos(r2)
+            enddo
+        endif
+
         if (k <= N_st) then
           residual_norm(k) = u_dot_u(U(1,shift2+k),sze)
           to_print(1,k) = lambda(k) + nuclear_repulsion
@@ -676,7 +869,7 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
 
   deallocate (                                                       &
       residual_norm,                                                 &
-      c,                                                             &
+      c, overlap,                                                    &
       h,                                                             &
       y, s_, s_tmp,                                                  &
       lambda                                                         &