From 6e35f8f8f8735bd4a898fabbc6bf552f382e517a Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 12 Mar 2024 15:30:52 +0100 Subject: [PATCH] fixed n_states > 1 for TC --- .../cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 2 +- plugins/local/fci_tc_bi/diagonalize_ci.irp.f | 80 ++++++----- .../local/tc_bi_ortho/tc_h_eigvectors.irp.f | 27 ++-- src/iterations/summary_tc.irp.f | 125 ++++++++++++++++++ 4 files changed, 181 insertions(+), 53 deletions(-) create mode 100644 src/iterations/summary_tc.irp.f diff --git a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index 2a7273d3..59ea3f11 100644 --- a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -98,7 +98,7 @@ subroutine run_stochastic_cipsi call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection ! stop - call print_summary(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2) + call print_summary_tc(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2) call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) diff --git a/plugins/local/fci_tc_bi/diagonalize_ci.irp.f b/plugins/local/fci_tc_bi/diagonalize_ci.irp.f index a9ded70c..a5242b87 100644 --- a/plugins/local/fci_tc_bi/diagonalize_ci.irp.f +++ b/plugins/local/fci_tc_bi/diagonalize_ci.irp.f @@ -20,48 +20,44 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) PROVIDE mo_l_coef mo_r_coef - print*,'*****' - print*,'New wave function information' - print*,'N_det tc = ',N_det - do k = 1, N_states - print*,'************' - print*,'State ',k - pt2_plus = pt2_data % variance(k) - pt2_minus = pt2_data % pt2(k) - pt2_abs = pt2_plus - pt2_minus - pt2_tot = pt2_plus + pt2_minus -! error_pt2_minus = pt2_data_err % pt2(k) -! error_pt2_plus = pt2_data_err % variance(k) -! error_pt2_tot = dsqrt(error_pt2_minus**2+error_pt2_plus**2) -! error_pt2_abs = error_pt2_tot ! same variance because independent variables - - pt1_norm = pt2_data % overlap(k,k) - rpt2_tot = pt2_tot / (1.d0 + pt1_norm) - - - print*,'norm_ground_left_right_bi_orth = ',norm_ground_left_right_bi_orth(k) - print*,'eigval_right_tc = ',eigval_right_tc_bi_orth(k) - print*,'*****' - - if(print_pt2) then - print*,'*****' - print*,'previous wave function info' - print*,'norm(before) = ',norm - print*,'E(before) = ',E_tc - print*,'PT1 norm = ',dsqrt(pt1_norm) - print*,'PT2 = ',pt2_tot - print*,'rPT2 = ',rpt2_tot - print*,'|PT2| = ',pt2_abs - print*,'Positive PT2 = ',pt2_plus - print*,'Negative PT2 = ',pt2_minus - print*,'E(before) + PT2 = ',E_tc + pt2_tot/norm - print*,'E(before) +rPT2 = ',E_tc + rpt2_tot/norm - write(*,'(A28,X,I10,X,100(F16.8,X))')'Ndet,E,E+PT2,E+RPT2,|PT2|=',ndet,E_tc ,E_tc + pt2_tot/norm,E_tc + rpt2_tot/norm,pt2_minus, pt2_plus - print*,'*****' - endif - E_tc(k) = eigval_right_tc_bi_orth(k) - norm(k) = norm_ground_left_right_bi_orth(k) - enddo +! print*,'*****' +! print*,'New wave function information' +! print*,'N_det tc = ',N_det +! do k = 1, N_states +! print*,'************' +! print*,'State ',k +! pt2_plus = pt2_data % variance(k) +! pt2_minus = pt2_data % pt2(k) +! pt2_abs = pt2_plus - pt2_minus +! pt2_tot = pt2_plus + pt2_minus +! +! pt1_norm = pt2_data % overlap(k,k) +! rpt2_tot = pt2_tot / (1.d0 + pt1_norm) +! +! +! print*,'norm_ground_left_right_bi_orth = ',norm_ground_left_right_bi_orth(k) +! print*,'eigval_right_tc = ',eigval_right_tc_bi_orth(k) +! print*,'*****' +! +! if(print_pt2) then +! print*,'*****' +! print*,'previous wave function info' +! print*,'norm(before) = ',norm +! print*,'E(before) = ',E_tc +! print*,'PT1 norm = ',dsqrt(pt1_norm) +! print*,'PT2 = ',pt2_tot +! print*,'rPT2 = ',rpt2_tot +! print*,'|PT2| = ',pt2_abs +! print*,'Positive PT2 = ',pt2_plus +! print*,'Negative PT2 = ',pt2_minus +! print*,'E(before) + PT2 = ',E_tc + pt2_tot/norm +! print*,'E(before) +rPT2 = ',E_tc + rpt2_tot/norm +! write(*,'(A28,X,I10,X,100(F16.8,X))')'Ndet,E,E+PT2,E+RPT2,|PT2|=',ndet,E_tc ,E_tc + pt2_tot/norm,E_tc + rpt2_tot/norm,pt2_minus, pt2_plus +! print*,'*****' +! endif +! E_tc(k) = eigval_right_tc_bi_orth(k) +! norm(k) = norm_ground_left_right_bi_orth(k) +! enddo psi_energy(1:N_states) = eigval_right_tc_bi_orth(1:N_states) - nuclear_repulsion psi_s2(1:N_states) = s2_eigvec_tc_bi_orth(1:N_states) diff --git a/plugins/local/tc_bi_ortho/tc_h_eigvectors.irp.f b/plugins/local/tc_bi_ortho/tc_h_eigvectors.irp.f index c90c84c5..6bf3d99e 100644 --- a/plugins/local/tc_bi_ortho/tc_h_eigvectors.irp.f +++ b/plugins/local/tc_bi_ortho/tc_h_eigvectors.irp.f @@ -86,17 +86,20 @@ end endif call non_hrmt_real_diag(N_det, H_prime, leigvec_tc_bi_orth_tmp, reigvec_tc_bi_orth_tmp, n_real_tc_bi_orth_eigval_right, eigval_right_tmp) + if(N_states.gt.1)then + print*,'n_real_tc_bi_orth_eigval_right = ',n_real_tc_bi_orth_eigval_right + endif ! do i = 1, N_det ! call get_H_tc_s2_l0_r0(leigvec_tc_bi_orth_tmp(1,i),reigvec_tc_bi_orth_tmp(1,i),1,N_det,expect_e(i), s2_values_tmp(i)) ! enddo call get_H_tc_s2_l0_r0(leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,N_det,N_det,expect_e, s2_values_tmp) + allocate(index_good_state_array(N_det),good_state_array(N_det)) i_state = 0 good_state_array = .False. if(s2_eig) then - if(only_expected_s2) then do j = 1, N_det ! Select at least n_states states with S^2 values closed to "expected_s2" @@ -116,6 +119,9 @@ end good_state_array(j) = .True. enddo endif + if(N_states.gt.1)then + print*,'i_state = ',i_state + endif if(i_state .ne. 0) then ! Fill the first "i_state" states that have a correct S^2 value @@ -338,11 +344,6 @@ end TOUCH psi_r_coef_bi_ortho call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(buffer) deallocate(buffer) -! print*,'After diag' -! do i = 1, N_det! old version -! print*,'i',i,psi_l_coef_bi_ortho(i,1),psi_r_coef_bi_ortho(i,1) -! call debug_det(psi_det(1,1,i),N_int) -! enddo END_PROVIDER @@ -357,23 +358,29 @@ subroutine bi_normalize(u_l, u_r, n, ld, nstates) implicit none integer, intent(in) :: n, ld, nstates double precision, intent(inout) :: u_l(ld,nstates), u_r(ld,nstates) - integer :: i, j - double precision :: accu, tmp + integer :: i, j,j_loc + double precision :: accu, tmp, maxval_tmp do i = 1, nstates !!!! Normalization of right eigenvectors |Phi> accu = 0.d0 ! TODO: dot product lapack + maxval_tmp = 0.d0 do j = 1, n accu += u_r(j,i) * u_r(j,i) + if(dabs(u_r(j,i)).gt.maxval_tmp)then + maxval_tmp = dabs(u_r(j,i)) + j_loc = j + endif enddo accu = 1.d0/dsqrt(accu) print*,'accu_r = ',accu + print*,'j_loc = ',j_loc do j = 1, n u_r(j,i) *= accu enddo - tmp = u_r(1,i) / dabs(u_r(1,i)) + tmp = u_r(j_loc,i) / dabs(u_r(j_loc,i)) do j = 1, n u_r(j,i) *= tmp enddo @@ -390,7 +397,7 @@ subroutine bi_normalize(u_l, u_r, n, ld, nstates) else accu = 1.d0/dsqrt(-accu) endif - tmp = (u_l(1,i) * u_r(1,i) )/dabs(u_l(1,i) * u_r(1,i)) + tmp = (u_l(j_loc,i) * u_r(j_loc,i) )/dabs(u_l(j_loc,i) * u_r(j_loc,i)) do j = 1, n u_l(j,i) *= accu * tmp u_r(j,i) *= accu diff --git a/src/iterations/summary_tc.irp.f b/src/iterations/summary_tc.irp.f new file mode 100644 index 00000000..00c2ba38 --- /dev/null +++ b/src/iterations/summary_tc.irp.f @@ -0,0 +1,125 @@ +subroutine print_summary_tc(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s2_) + use selection_types + implicit none + BEGIN_DOC +! Print the extrapolated energy in the output + END_DOC + + integer, intent(in) :: n_det_, n_configuration_, n_st + double precision, intent(in) :: e_(n_st), s2_(n_st) + type(pt2_type) , intent(in) :: pt2_data, pt2_data_err + integer :: i, k + integer :: N_states_p + character*(9) :: pt2_string + character*(512) :: fmt + double precision, allocatable :: pt2_minus(:),pt2_plus(:),pt2_tot(:), pt2_abs(:),pt1_norm(:),rpt2_tot(:) + double precision, allocatable :: error_pt2_minus(:), error_pt2_plus(:), error_pt2_tot(:), error_pt2_abs(:) + + if (do_pt2) then + pt2_string = ' ' + else + pt2_string = '(approx)' + endif + + N_states_p = min(N_det_,n_st) + + allocate(pt2_minus(N_states_p),pt2_plus(N_states_p),pt2_tot(N_states_p), pt2_abs(N_states_p),pt1_norm(N_states_p),rpt2_tot(N_states_p)) + allocate(error_pt2_minus(N_states_p), error_pt2_plus(N_states_p), error_pt2_tot(N_states_p), error_pt2_abs(N_states_p)) + do k = 1, N_states_p + pt2_plus(k) = pt2_data % variance(k) + pt2_minus(k) = pt2_data % pt2(k) + pt2_abs(k) = pt2_plus(k) - pt2_minus(k) + pt2_tot(k) = pt2_plus(k) + pt2_minus(k) + pt1_norm(k) = pt2_data % overlap(k,k) + rpt2_tot(k) = pt2_tot(k) / (1.d0 + pt1_norm(k)) + error_pt2_minus(k) = pt2_data_err % pt2(k) + error_pt2_plus(k) = pt2_data_err % variance(k) + error_pt2_tot(k) = dsqrt(error_pt2_minus(k)**2+error_pt2_plus(k)**2) + error_pt2_abs(k) = error_pt2_tot(k) ! same variance because independent variables + enddo + k=1 + write(*,'(A40,X,I10,X,100(F16.8,X))')'Ndet,E,E+PT2,pt2_minus,pt2_plus,pt2_abs=',n_det_,e_(k),e_(k) + pt2_tot(k),e_(k) + rpt2_tot(k),pt2_minus(k), pt2_plus(k),pt2_abs(k) + + print *, '' + print '(A,I12)', 'Summary at N_det = ', N_det_ + print '(A)', '-----------------------------------' + print *, '' + + write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' + write(*,fmt) + write(fmt,*) '(13X,', N_states_p, '(6X,A7,1X,I6,10X))' + write(*,fmt) ('State',k, k=1,N_states_p) + write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' + write(*,fmt) + write(fmt,*) '(A13,', N_states_p, '(1X,F14.8,15X))' + write(*,fmt) '# E ', e_(1:N_states_p) + if (N_states_p > 1) then + write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1) + write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0 + endif + write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))' + write(*,fmt) '# PT2 '//pt2_string, (pt2_tot(k), error_pt2_tot(k), k=1,N_states_p) + write(*,fmt) '# rPT2'//pt2_string, (rpt2_tot(k), error_pt2_tot(k), k=1,N_states_p) + write(*,'(A)') '#' + write(*,fmt) '# E+PT2 ', (e_(k)+pt2_tot(k) ,error_pt2_tot(k), k=1,N_states_p) + write(*,fmt) '# E+rPT2 ', (e_(k)+rpt2_tot(k),error_pt2_tot(k), k=1,N_states_p) + if (N_states_p > 1) then + write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_tot(k)-e_(1)-pt2_tot(1)), & + dsqrt(error_pt2_tot(k)*error_pt2_tot(k)+error_pt2_tot(1)*error_pt2_tot(1)), k=1,N_states_p) + write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_tot(k)-e_(1)-pt2_tot(1))*27.211396641308d0, & + dsqrt(error_pt2_tot(k)*error_pt2_tot(k)+error_pt2_tot(1)*error_pt2_tot(1))*27.211396641308d0, k=1,N_states_p) + endif + write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' + write(*,fmt) + print *, '' + + print *, 'N_det = ', N_det_ + print *, 'N_states = ', n_st + if (s2_eig) then + print *, 'N_cfg = ', N_configuration_ + if (only_expected_s2) then + print *, 'N_csf = ', N_csf + endif + endif + print *, '' + + do k=1, N_states_p + print*,'* State ',k + print *, '< S^2 > = ', s2_(k) + print *, 'E = ', e_(k) + print *, 'PT norm = ', pt1_norm(k) + print *, 'PT2 = ', pt2_tot(k), ' +/- ', error_pt2_tot(k) + print *, 'rPT2 = ', rpt2_tot(k), ' +/- ', error_pt2_tot(k) + print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_tot(k) , ' +/- ', error_pt2_tot(k) + print *, 'E+rPT2'//pt2_string//' = ', e_(k)+rpt2_tot(k), ' +/- ', error_pt2_tot(k) + print *, 'Positive PT2 = ',pt2_plus(k),' +/- ',error_pt2_plus(k) + print *, 'Negative PT2 = ',pt2_minus(k),' +/- ',error_pt2_minus(k) + print *, 'Abs PT2 = ',pt2_abs(k), ' +/- ',error_pt2_abs(k) + print *, '' + enddo + + print *, '-----' + if(n_st.gt.1)then + print *, 'Variational Energy difference (au | eV)' + do i=2, N_states_p + print*,'Delta E = ', (e_(i) - e_(1)), & + (e_(i) - e_(1)) * 27.211396641308d0 + enddo + print *, '-----' + print*, 'Variational + perturbative Energy difference (au | eV)' + do i=2, N_states_p + print*,'Delta E = ', (e_(i)+ pt2_tot(i) - (e_(1) + pt2_tot(1))), & + (e_(i)+ pt2_tot(i) - (e_(1) + pt2_tot(1))) * 27.211396641308d0 + enddo + print *, '-----' + print*, 'Variational + renormalized perturbative Energy difference (au | eV)' + do i=2, N_states_p + print*,'Delta E = ', (e_(i)+ rpt2_tot(i) - (e_(1) + rpt2_tot(1))), & + (e_(i)+ rpt2_tot(i) - (e_(1) + rpt2_tot(1))) * 27.211396641308d0 + enddo + endif + +! call print_energy_components() + +end subroutine +