Multistate dressing OK

ocaml/TaskServer.ml

@@ -830,10 +830,11 @@ let run ~port =
               let () = 
                 if debug_env then
                   begin
-                    Printf.sprintf "q:%d  r:%d  n:%d  : %s\n%!"
+                    Printf.sprintf "q:%d  r:%d  n:%d  c:%d : %s\n%!"
                     (Queuing_system.number_of_queued program_state.queue)
                     (Queuing_system.number_of_running program_state.queue)
                     (Queuing_system.number_of_tasks program_state.queue)
+                    (Queuing_system.number_of_clients program_state.queue)
                     (Message.to_string message)
                     |> debug
                   end

plugins/DavidsonDressed/diagonalize_CI.irp.f

@@ -55,66 +55,98 @@ END_PROVIDER
 
   if (diag_algorithm == "Davidson") then
     
-    allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)),&
-        eigenvectors_s2(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)
+        CI_eigenvectors_dressed(i,j) = psi_coef(i,j)
       enddo
     enddo
-    do mrcc_state=1,N_states
+    call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,&
-      do j=mrcc_state,min(N_states,N_det)
+        size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,&
-        do i=1,N_det
+        N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1)
-          eigenvectors(i,j) = psi_coef(i,j)
-        enddo
-      enddo
-      call davidson_diag_HS2(psi_det,eigenvectors, eigenvectors_s2,    &
-          size(eigenvectors,1),                                        &
-          eigenvalues,N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,&
-          mrcc_state)
-      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)
-    
     
   else if (diag_algorithm == "Lapack") then
     
     allocate (eigenvectors(size(H_matrix_dressed,1),N_det))
     allocate (eigenvalues(N_det))
 
-    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 j=mrcc_state,min(N_states,N_det)
-        do i=1,N_det
-          eigenvectors(i,j) = psi_coef(i,j)
-        enddo
-      enddo
-
     call lapack_diag(eigenvalues,eigenvectors,                         &
-        H_matrix_dressed(1,1,mrcc_state),size(H_matrix_dressed,1),N_det)
+        H_matrix_dressed,size(H_matrix_dressed,1),N_det)
-      CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
+     CI_electronic_energy_dressed(:) = 0.d0
-      CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
+     if (s2_eig) then
-    enddo
+       i_state = 0
-    do k=N_states+1,N_states_diag
+       allocate (s2_eigvalues(N_det))
-      CI_eigenvectors_dressed(1:N_det,k) = eigenvectors(1:N_det,k)
+       allocate(index_good_state_array(N_det),good_state_array(N_det))
-      CI_electronic_energy_dressed(k) = eigenvalues(k)
+       good_state_array = .False.
-    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))
 
+       call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
+                  N_det,size(eigenvectors,1))
+       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
+           i_state +=1
+           index_good_state_array(i_state) = j
+           good_state_array(j) = .True.
+         endif
+         if(i_state.eq.N_states) then
+           exit
+         endif
+       enddo
+       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
+             CI_eigenvectors_dressed(i,j) = eigenvectors(i,index_good_state_array(j))
+           enddo
+           CI_electronic_energy_dressed(j) = eigenvalues(index_good_state_array(j))
+           CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j))
+         enddo
+         i_other_state = 0
+         do j = 1, N_det
+           if(good_state_array(j))cycle
+           i_other_state +=1
+           if(i_state+i_other_state.gt.n_states_diag)then
+             exit
+           endif
+           do i=1,N_det
+             CI_eigenvectors_dressed(i,i_state+i_other_state) = eigenvectors(i,j)
+           enddo
+           CI_electronic_energy_dressed(i_state+i_other_state) = eigenvalues(j)
+           CI_eigenvectors_s2_dressed(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
+         enddo
+       else
+         print*,''
+         print*,'!!!!!!!!   WARNING  !!!!!!!!!'
+         print*,'  Within the ',N_det,'determinants selected'
+         print*,'  and the ',N_states_diag,'states requested'
+         print*,'  We did not find any state 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_dressed'
+         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_dressed(i,j) = eigenvectors(i,j)
+           enddo
+           CI_electronic_energy_dressed(j) = eigenvalues(j)
+           CI_eigenvectors_s2_dressed(j) = s2_eigvalues(j)
+         enddo
+       endif
+       deallocate(index_good_state_array,good_state_array)
+       deallocate(s2_eigvalues)
+     else
+       call u_0_S2_u_0(CI_eigenvectors_s2_dressed,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
+           CI_eigenvectors_dressed(i,j) = eigenvectors(i,j)
+         enddo
+         CI_electronic_energy_dressed(j) = eigenvalues(j)
+       enddo
+     endif
     deallocate(eigenvectors,eigenvalues)
   endif
 
@@ -137,24 +169,23 @@ end
 
 
 
-BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
+BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ]
  implicit none
  BEGIN_DOC
  ! Dressed H with Delta_ij
  END_DOC
- integer                        :: i, j, ii,jj, dressing_state
+ integer                        :: i, j, k
- do dressing_state = 1,N_states
+
+  h_matrix_dressed(1:N_det,1:N_det) = h_matrix_all_dets(1:N_det,1:N_det)
+  do k=1,N_states
     do j=1,N_det
       do i=1,N_det
-       h_matrix_dressed(i,j,dressing_state) = h_matrix_all_dets(i,j) 
+        h_matrix_dressed(i,j) = h_matrix_dressed(i,j) +             &
+            0.5d0 * (dressing_column_h(i,k) * psi_coef(j,k) + &
+                     dressing_column_h(j,k) * psi_coef(i,k))
       enddo
     enddo
-   i = dressed_column_idx(dressing_state)
-   do j = 1, N_det
-     h_matrix_dressed(i,j,dressing_state) += dressing_column_h(j,dressing_state)
-     h_matrix_dressed(j,i,dressing_state) += dressing_column_h(j,dressing_state)
-   enddo
-   h_matrix_dressed(i,i,dressing_state) -= dressing_column_h(i,dressing_state)
   enddo
+
 END_PROVIDER
 

plugins/dress_zmq/dress_general.irp.f

@@ -1,5 +1,3 @@
-
-
 subroutine run_dressing(N_st,energy)
   implicit none
   

plugins/dress_zmq/dress_stoch_routines.irp.f

@@ -4,7 +4,7 @@ BEGIN_PROVIDER [ integer, fragment_first ]
 END_PROVIDER
 
 
-subroutine ZMQ_dress(E, dress, delta, delta_s2, relative_error)
+subroutine ZMQ_dress(E, dress, delta_out, delta_s2_out, relative_error)
   use f77_zmq
   
   implicit none
@@ -15,9 +15,11 @@ subroutine ZMQ_dress(E, dress, delta, delta_s2, relative_error)
   integer, external              :: omp_get_thread_num
   double precision, intent(in)   :: E(N_states), relative_error
   double precision, intent(out)  :: dress(N_states)
-  double precision, intent(out)  :: delta(N_states, N_det)
+  double precision, intent(out)  :: delta_out(N_states, N_det)
-  double precision, intent(out)  :: delta_s2(N_states, N_det)
+  double precision, intent(out)  :: delta_s2_out(N_states, N_det)
   
+  double precision, allocatable  :: delta(:,:)
+  double precision, allocatable  :: delta_s2(:,:)
   
   integer                        :: i, j, k, Ncp
   
@@ -27,6 +29,7 @@ subroutine ZMQ_dress(E, dress, delta, delta_s2, relative_error)
   double precision               :: state_average_weight_save(N_states)
   
   
+  allocate(delta(N_states,N_det), delta_s2(N_det,N_states))
   state_average_weight_save(:) = state_average_weight(:)
   do dress_stoch_istate=1,N_states
     SOFT_TOUCH dress_stoch_istate
@@ -108,6 +111,8 @@ subroutine ZMQ_dress(E, dress, delta, delta_s2, relative_error)
       call dress_slave_inproc(i)
     endif
     !$OMP END PARALLEL
+    delta_out(dress_stoch_istate,1:N_det) = delta(dress_stoch_istate,1:N_det)
+    delta_s2_out(dress_stoch_istate,1:N_det) = delta_s2_out(dress_stoch_istate,1:N_det)
     call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'dress')
     
     print *, '========== ================= ================= ================='
@@ -115,6 +120,7 @@ subroutine ZMQ_dress(E, dress, delta, delta_s2, relative_error)
   FREE dress_stoch_istate
   state_average_weight(:) = state_average_weight_save(:)
   TOUCH state_average_weight
+  deallocate(delta,delta_s2)
   
 end subroutine
 

plugins/dress_zmq/dressing_vector.irp.f

@@ -3,7 +3,7 @@
 &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
  implicit none
  BEGIN_DOC
- ! Null dressing vectors
+ ! \Delta_{state-specific}. \Psi
  END_DOC
 
  integer :: i,ii,k,j, l
@@ -14,18 +14,16 @@
  dressing_column_s(:,:) = 0.d0
 
  do k=1,N_states
-   l = dressed_column_idx(k)
-   f = 1.d0/psi_coef(l,k)
    do j = 1, n_det
-     dressing_column_h(j,k) = delta_ij(k,j,1) * f
+     dressing_column_h(j,k) = delta_ij(k,j,1) 
-     dressing_column_s(j,k) = delta_ij(k,j,2) * f
+     dressing_column_s(j,k) = delta_ij(k,j,2) 
    enddo
-   tmp = u_dot_v(dressing_column_h(1,k), psi_coef(1,k), N_det) &
+!   tmp = u_dot_v(dressing_column_h(1,k), psi_coef(1,k), N_det) &
-     - dressing_column_h(l,k) * psi_coef(l,k)
+!     - dressing_column_h(l,k) * psi_coef(l,k)
-   dressing_column_h(l,k) -= tmp * f
+!   dressing_column_h(l,k) -= tmp * f
-   tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) &
+!   tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) &
-     - dressing_column_s(l,k) * psi_coef(l,k)
+!     - dressing_column_s(l,k) * psi_coef(l,k)
-   dressing_column_s(l,k) -= tmp * f
+!   dressing_column_s(l,k) -= tmp * f
  enddo
 
 END_PROVIDER

plugins/dress_zmq/run_dress_slave.irp.f

@@ -62,9 +62,9 @@ subroutine run_dress_slave(thread,iproc,energy)
       exit
     end if
   end do
-  call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
+  call disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id)
-  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
   call end_zmq_push_socket(zmq_socket_push,thread)
+  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
 end subroutine
 
 

plugins/shiftedbk/shifted_bk.irp.f

@@ -1,4 +1,3 @@
-
 program shifted_bk
   implicit none
   BEGIN_DOC

src/Davidson/diagonalization_hs2_dressed.irp.f

@@ -38,7 +38,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
   double precision, allocatable  :: H_jj(:), S2_jj(:)
   
   double precision, external     :: diag_H_mat_elem, diag_S_mat_elem
-  integer                        :: i
+  integer                        :: i,k
   ASSERT (N_st > 0)
   ASSERT (sze > 0)
   ASSERT (Nint > 0)
@@ -58,7 +58,11 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
   !$OMP END PARALLEL
 
   if (dressing_state > 0) then
-    H_jj(dressed_column_idx(dressing_state)) += dressing_column_h(dressed_column_idx(dressing_state),dressing_state)
+    do k=1,N_st
+      do i=1,sze
+        H_jj(i) += u_in(i,k) * dressing_column_h(i,k)
+      enddo
+    enddo
   endif
 
   call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state)
@@ -150,17 +154,17 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
   do i=1,N_st
     write_buffer = trim(write_buffer)//' ================ =========== ==========='
   enddo
-  write(6,'(A)') write_buffer(1:6+41*N_states)
+  write(6,'(A)') write_buffer(1:6+41*N_st)
   write_buffer = 'Iter'
   do i=1,N_st
     write_buffer = trim(write_buffer)//'       Energy          S^2       Residual '
   enddo
-  write(6,'(A)') write_buffer(1:6+41*N_states)
+  write(6,'(A)') write_buffer(1:6+41*N_st)
   write_buffer = '====='
   do i=1,N_st
     write_buffer = trim(write_buffer)//' ================ =========== ==========='
   enddo
-  write(6,'(A)') write_buffer(1:6+41*N_states)
+  write(6,'(A)') write_buffer(1:6+41*N_st)
   
 
   allocate(                                                          &
@@ -242,17 +246,35 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       
       if (dressing_state > 0) then
 
-        l = dressed_column_idx(dressing_state)
+        call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
-        do istate=1,N_st_diag
+          psi_coef, size(psi_coef,1), &
-          do i=1,sze
+          U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
-              W(i,shift+istate) += dressing_column_h(i,dressing_state) * U(l,shift+istate)
+
-              S(i,shift+istate) += dressing_column_s(i,dressing_state) * U(l,shift+istate)
+        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
-              W(l,shift+istate) += dressing_column_h(i,dressing_state) * U(i,shift+istate)
+          dressing_column_h, size(dressing_column_h,1), s_tmp, size(s_tmp,1), &
-              S(l,shift+istate) += dressing_column_s(i,dressing_state) * U(i,shift+istate)
+          1.d0, W(1,shift+1), size(W,1))
-          enddo
+
-          W(l,shift+istate) -= dressing_column_h(l,dressing_state) * U(l,shift+istate)
+        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
-          S(l,shift+istate) -= dressing_column_s(l,dressing_state) * U(l,shift+istate)
+          dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), &
-        enddo
+          1.d0, S(1,shift+1), size(S,1))
+
+
+        call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
+          dressing_column_h, size(dressing_column_h,1), &
+          U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+
+        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
+          psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
+          1.d0, W(1,shift+1), size(W,1))
+
+        call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
+          dressing_column_s, size(dressing_column_s,1), &
+          U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+
+        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
+          psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
+          1.d0, S(1,shift+1), size(S,1))
+
       endif
 
       ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>