From 6be57e3c01f1ba71d22495fcfaed52448792163e Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 12 Apr 2023 17:10:06 +0200 Subject: [PATCH] fixed bug in S2 for TC davidson --- .../diagonalization_hs2_dressed.irp.f | 4 ++-- src/tc_bi_ortho/dav_h_tc_s2.irp.f | 16 +++++++++++---- src/tc_bi_ortho/tc_h_eigvectors.irp.f | 20 +++++++++++++++++-- src/tc_bi_ortho/test_s2_tc.irp.f | 7 +++++-- 4 files changed, 37 insertions(+), 10 deletions(-) diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index 8117f320..ac71d1d4 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -465,8 +465,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ integer :: lwork, info double precision, allocatable :: work(:) -! y = h - y = h_p + y = h +! y = h_p lwork = -1 allocate(work(1)) call dsygv(1,'V','U',shift2,y,size(y,1), & diff --git a/src/tc_bi_ortho/dav_h_tc_s2.irp.f b/src/tc_bi_ortho/dav_h_tc_s2.irp.f index c0ea054a..ea9cacff 100644 --- a/src/tc_bi_ortho/dav_h_tc_s2.irp.f +++ b/src/tc_bi_ortho/dav_h_tc_s2.irp.f @@ -1,7 +1,7 @@ ! --- -subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N_st_diag_in, converged, hcalc) +subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N_st_diag_in, n_it_max_dav, converged, hcalc) use mmap_module @@ -21,11 +21,17 @@ subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N ! Initial guess vectors are not necessarily orthonormal ! ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output) + ! + ! !!! WARNING !!! IT SEEMS THAT IF THE NUMBER OF MACRO ITERATIONS EXCEEDS n_it_max_dav, + ! + ! THE RECONTRACTION IS WRONG. YOU SHOULD CONSIDER CALLING MULTIPLE TIME THE ROUTINE + ! + ! SEE FOR INSTANCE IN tc_bi_ortho/tc_h_eigvectors.irp.f END_DOC implicit none - integer, intent(in) :: sze, N_st, N_st_diag_in + integer, intent(in) :: sze, N_st, N_st_diag_in, n_it_max_dav double precision, intent(in) :: H_jj(sze) logical, intent(inout) :: converged double precision, intent(inout) :: u_in(sze,N_st_diag_in) @@ -246,7 +252,9 @@ subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N itertot = 0 - do while (.not.converged) +! do while (.not.converged.or.itertot.le.n_it_max_dav) + integer :: iiii + do iiii = 1, n_it_max_dav itertot = itertot + 1 if(itertot == 8) then @@ -522,7 +530,7 @@ subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N enddo endif enddo - + if(converged)exit enddo ! loop over while ! --- diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f index 71dad8d6..91775cf1 100644 --- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f +++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f @@ -208,7 +208,11 @@ end call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i)) enddo !!!! Preparing the left-eigenvector + print*,'---------------------------------' + print*,'---------------------------------' print*,'Computing the left-eigenvector ' + print*,'---------------------------------' + print*,'---------------------------------' vec_tmp = 0.d0 do istate = 1, N_states vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate) @@ -218,12 +222,21 @@ end enddo ! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav) ! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_dagger_u_0_opt) - call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt) + integer :: n_it_max + n_it_max = 1 + converged = .False. + do while (.not.converged) + call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt) + enddo do istate = 1, N_states leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) enddo + print*,'---------------------------------' + print*,'---------------------------------' print*,'Computing the right-eigenvector ' + print*,'---------------------------------' + print*,'---------------------------------' !!!! Preparing the right-eigenvector vec_tmp = 0.d0 do istate = 1, N_states @@ -234,7 +247,10 @@ end enddo ! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav) ! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_u_0_opt) - call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt) + converged = .False. + do while (.not.converged) + call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt) + enddo do istate = 1, N_states reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) enddo diff --git a/src/tc_bi_ortho/test_s2_tc.irp.f b/src/tc_bi_ortho/test_s2_tc.irp.f index 4229fef1..4debe2e2 100644 --- a/src/tc_bi_ortho/test_s2_tc.irp.f +++ b/src/tc_bi_ortho/test_s2_tc.irp.f @@ -105,7 +105,9 @@ subroutine routine_test_s2_davidson do istate = 1, N_states leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) enddo - call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt) + integer :: n_it_max + n_it_max = 1 + call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt) double precision, allocatable :: v_0_new(:,:),s_0_new(:,:) integer :: sze,N_st logical :: do_right @@ -138,7 +140,8 @@ subroutine routine_test_s2_davidson do istate = 1, N_states leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) enddo - call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_u_0_opt) + n_it_max = 1 + call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt) sze = N_det N_st = 1 do_right = .True.