diff --git a/ocaml/TaskServer.ml b/ocaml/TaskServer.ml index 170e011a..ca295971 100644 --- a/ocaml/TaskServer.ml +++ b/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 diff --git a/plugins/DavidsonDressed/diagonalize_CI.irp.f b/plugins/DavidsonDressed/diagonalize_CI.irp.f index 3d1c1118..2ee540e2 100644 --- a/plugins/DavidsonDressed/diagonalize_CI.irp.f +++ b/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 - 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_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 davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,& + size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,& + N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1) 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) - 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)) + H_matrix_dressed,size(H_matrix_dressed,1),N_det) + CI_electronic_energy_dressed(:) = 0.d0 + if (s2_eig) then + 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)) + 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 - do dressing_state = 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) - 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 + integer :: i, j, k + + 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) = 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 + enddo + END_PROVIDER diff --git a/plugins/dress_zmq/dress_general.irp.f b/plugins/dress_zmq/dress_general.irp.f index 0bf7e715..068d811e 100644 --- a/plugins/dress_zmq/dress_general.irp.f +++ b/plugins/dress_zmq/dress_general.irp.f @@ -1,5 +1,3 @@ - - subroutine run_dressing(N_st,energy) implicit none diff --git a/plugins/dress_zmq/dress_stoch_routines.irp.f b/plugins/dress_zmq/dress_stoch_routines.irp.f index 06b5d538..f1406b7b 100644 --- a/plugins/dress_zmq/dress_stoch_routines.irp.f +++ b/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_s2(N_states, N_det) + double precision, intent(out) :: delta_out(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 diff --git a/plugins/dress_zmq/dressing_vector.irp.f b/plugins/dress_zmq/dressing_vector.irp.f index 494e7c4b..5a8fee3b 100644 --- a/plugins/dress_zmq/dressing_vector.irp.f +++ b/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_s(j,k) = delta_ij(k,j,2) * f + dressing_column_h(j,k) = delta_ij(k,j,1) + 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) & - - dressing_column_h(l,k) * psi_coef(l,k) - dressing_column_h(l,k) -= tmp * f - 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) -= tmp * f +! 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) -= tmp * f +! 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) -= tmp * f enddo END_PROVIDER diff --git a/plugins/dress_zmq/run_dress_slave.irp.f b/plugins/dress_zmq/run_dress_slave.irp.f index b0896c00..0311d2ed 100644 --- a/plugins/dress_zmq/run_dress_slave.irp.f +++ b/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 end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) call end_zmq_push_socket(zmq_socket_push,thread) + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) end subroutine diff --git a/plugins/shiftedbk/shifted_bk.irp.f b/plugins/shiftedbk/shifted_bk.irp.f index 4c0408d8..d2ab57ad 100644 --- a/plugins/shiftedbk/shifted_bk.irp.f +++ b/plugins/shiftedbk/shifted_bk.irp.f @@ -1,4 +1,3 @@ - program shifted_bk implicit none BEGIN_DOC diff --git a/src/Davidson/diagonalization_hs2_dressed.irp.f b/src/Davidson/diagonalization_hs2_dressed.irp.f index 4e853f32..59e9c9fe 100644 --- a/src/Davidson/diagonalization_hs2_dressed.irp.f +++ b/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) - do istate=1,N_st_diag - do i=1,sze - 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) - W(l,shift+istate) += dressing_column_h(i,dressing_state) * U(i,shift+istate) - S(l,shift+istate) += dressing_column_s(i,dressing_state) * U(i,shift+istate) - enddo - W(l,shift+istate) -= dressing_column_h(l,dressing_state) * U(l,shift+istate) - S(l,shift+istate) -= dressing_column_s(l,dressing_state) * U(l,shift+istate) - enddo + call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & + psi_coef, size(psi_coef,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, & + dressing_column_h, size(dressing_column_h,1), s_tmp, size(s_tmp,1), & + 1.d0, W(1,shift+1), size(W,1)) + + call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, & + dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), & + 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 = =