From 9175fb21c9dcbe931f89d96cf1297221693d5fde Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 12 Mar 2024 14:05:38 +0100 Subject: [PATCH] modifs in json and diagonalize_ci for fci tc bi --- .../local/cipsi_tc_bi_ortho/selection.irp.f | 7 +- .../cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 1 - .../cipsi_tc_bi_ortho/write_cipsi_json.irp.f | 29 +++- plugins/local/fci_tc_bi/diagonalize_ci.irp.f | 124 ++++++------------ .../local/tc_bi_ortho/tc_h_eigvectors.irp.f | 18 +-- 5 files changed, 82 insertions(+), 97 deletions(-) diff --git a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f index 06cf848b..a01d4131 100644 --- a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f @@ -980,8 +980,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d psi_h_alpha = mat_l(istate, p1, p2) pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate) - pt2_data % variance(istate) = pt2_data % variance(istate) + dabs(e_pert(istate)) - pt2_data % pt2(istate) = pt2_data % pt2(istate) + e_pert(istate) + if(e_pert(istate).gt.0.d0)then! accumulate the positive part of the pt2 + pt2_data % variance(istate) = pt2_data % variance(istate) + e_pert(istate) + else ! accumulate the negative part of the pt2 + pt2_data % pt2(istate) = pt2_data % pt2(istate) + e_pert(istate) + endif select case (weight_selection) case(5) 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 66d82964..2a7273d3 100644 --- a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -57,7 +57,6 @@ subroutine run_stochastic_cipsi ! endif print_pt2 = .False. call diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) -! call routine_save_right ! if (N_det > N_det_max) then diff --git a/plugins/local/cipsi_tc_bi_ortho/write_cipsi_json.irp.f b/plugins/local/cipsi_tc_bi_ortho/write_cipsi_json.irp.f index 98a402a2..f8c95d38 100644 --- a/plugins/local/cipsi_tc_bi_ortho/write_cipsi_json.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/write_cipsi_json.irp.f @@ -9,6 +9,8 @@ subroutine write_cipsi_json(pt2_data, pt2_data_err) call lock_io character*(64), allocatable :: fmtk(:) + double precision:: pt2_minus,pt2_plus,pt2_tot, pt2_abs + double precision :: error_pt2_minus, error_pt2_plus, error_pt2_tot, error_pt2_abs integer :: N_states_p, N_iter_p N_states_p = min(N_states,N_det) N_iter_p = min(N_iter,8) @@ -26,15 +28,34 @@ subroutine write_cipsi_json(pt2_data, pt2_data_err) endif write(json_unit, json_array_open_fmt) 'states' do k=1,N_states_p + 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 write(json_unit, json_dict_uopen_fmt) write(json_unit, json_real_fmt) 'energy', psi_energy_with_nucl_rep(k) write(json_unit, json_real_fmt) 's2', psi_s2(k) - write(json_unit, json_real_fmt) 'pt2', pt2_data % pt2(k) - write(json_unit, json_real_fmt) 'pt2_err', pt2_data_err % pt2(k) + + write(json_unit, json_real_fmt) 'pt2', pt2_tot + write(json_unit, json_real_fmt) 'pt2_err', error_pt2_tot + + write(json_unit, json_real_fmt) 'pt2_minus', pt2_minus + write(json_unit, json_real_fmt) 'pt2_minus_err', error_pt2_minus + + write(json_unit, json_real_fmt) 'pt2_abs', pt2_abs + write(json_unit, json_real_fmt) 'pt2_abs_err', error_pt2_abs + + write(json_unit, json_real_fmt) 'pt2_plus', pt2_plus + write(json_unit, json_real_fmt) 'pt2_plus_err', error_pt2_plus + write(json_unit, json_real_fmt) 'rpt2', pt2_data % rpt2(k) write(json_unit, json_real_fmt) 'rpt2_err', pt2_data_err % rpt2(k) - write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k) - write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k) +! write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k) +! write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k) write(json_unit, json_array_open_fmt) 'ex_energy' do i=2,N_iter_p write(json_unit, fmtk(i)) extrapolated_energy(i,k) diff --git a/plugins/local/fci_tc_bi/diagonalize_ci.irp.f b/plugins/local/fci_tc_bi/diagonalize_ci.irp.f index 6c8f3431..a9ded70c 100644 --- a/plugins/local/fci_tc_bi/diagonalize_ci.irp.f +++ b/plugins/local/fci_tc_bi/diagonalize_ci.irp.f @@ -11,49 +11,61 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) use selection_types implicit none integer, intent(inout) :: ndet ! number of determinants from before - double precision, intent(inout) :: E_tc, norm ! E and norm from previous wave function + double precision, intent(inout) :: E_tc(N_states), norm(N_states) ! E and norm from previous wave function type(pt2_type) , intent(in) :: pt2_data ! PT2 from previous wave function logical, intent(in) :: print_pt2 - integer :: i, j - double precision :: pt2_tmp, pt1_norm, rpt2_tmp, abs_pt2 + integer :: i, j,k + double precision:: pt2_minus,pt2_plus,pt2_tot, pt2_abs,pt1_norm,rpt2_tot + double precision :: error_pt2_minus, error_pt2_plus, error_pt2_tot, error_pt2_abs PROVIDE mo_l_coef mo_r_coef - pt2_tmp = pt2_data % pt2(1) - abs_pt2 = pt2_data % variance(1) - pt1_norm = pt2_data % overlap(1,1) - rpt2_tmp = pt2_tmp/(1.d0 + pt1_norm) - print*,'*****' print*,'New wave function information' print*,'N_det tc = ',N_det - print*,'norm_ground_left_right_bi_orth = ',norm_ground_left_right_bi_orth - print*,'eigval_right_tc = ',eigval_right_tc_bi_orth(1) - print*,'Ndet, E_tc = ',N_det,eigval_right_tc_bi_orth(1) - 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_tmp - print*,'rPT2 = ',rpt2_tmp - print*,'|PT2| = ',abs_pt2 - print*,'Positive PT2 = ',(pt2_tmp + abs_pt2)*0.5d0 - print*,'Negative PT2 = ',(pt2_tmp - abs_pt2)*0.5d0 - print*,'E(before) + PT2 = ',E_tc + pt2_tmp/norm - print*,'E(before) +rPT2 = ',E_tc + rpt2_tmp/norm - write(*,'(A28,X,I10,X,100(F16.8,X))')'Ndet,E,E+PT2,E+RPT2,|PT2|=',ndet,E_tc ,E_tc + pt2_tmp/norm,E_tc + rpt2_tmp/norm,abs_pt2 - print*,'*****' - endif + 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 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) - E_tc = eigval_right_tc_bi_orth(1) - norm = norm_ground_left_right_bi_orth ndet = N_det do j = 1, N_states do i = 1, N_det @@ -71,53 +83,3 @@ end ! --- -subroutine print_CI_dressed(ndet, E_tc, norm, pt2_data, print_pt2) - - BEGIN_DOC - ! Replace the coefficients of the CI states by the coefficients of the - ! eigenstates of the CI matrix - END_DOC - - use selection_types - implicit none - integer, intent(inout) :: ndet ! number of determinants from before - double precision, intent(inout) :: E_tc,norm ! E and norm from previous wave function - type(pt2_type) , intent(in) :: pt2_data ! PT2 from previous wave function - logical, intent(in) :: print_pt2 - integer :: i, j - - print*,'*****' - print*,'New wave function information' - print*,'N_det tc = ',N_det - print*,'norm_ground_left_right_bi_orth = ',norm_ground_left_right_bi_orth - print*,'eigval_right_tc = ',eigval_right_tc_bi_orth(1) - print*,'Ndet, E_tc = ',N_det,eigval_right_tc_bi_orth(1) - print*,'*****' - - if(print_pt2) then - print*,'*****' - print*,'previous wave function info' - print*,'norm(before) = ',norm - print*,'E(before) = ',E_tc - print*,'PT1 norm = ',dsqrt(pt2_data % overlap(1,1)) - print*,'E(before) + PT2 = ',E_tc + (pt2_data % pt2(1))/norm - print*,'PT2 = ',pt2_data % pt2(1) - print*,'Ndet, E_tc, E+PT2 = ',ndet,E_tc,E_tc + (pt2_data % pt2(1))/norm,dsqrt(pt2_data % overlap(1,1)) - print*,'*****' - endif - - E_tc = eigval_right_tc_bi_orth(1) - norm = norm_ground_left_right_bi_orth - ndet = N_det - - do j = 1, N_states - do i = 1, N_det - psi_coef(i,j) = reigvec_tc_bi_orth(i,j) - enddo - enddo - SOFT_TOUCH eigval_left_tc_bi_orth eigval_right_tc_bi_orth leigvec_tc_bi_orth norm_ground_left_right_bi_orth psi_coef reigvec_tc_bi_orth - -end - -! --- - 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 75f3dfbe..c90c84c5 100644 --- a/plugins/local/tc_bi_ortho/tc_h_eigvectors.irp.f +++ b/plugins/local/tc_bi_ortho/tc_h_eigvectors.irp.f @@ -45,12 +45,12 @@ end ! --- - BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth, (N_states) ] -&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth , (N_states) ] -&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth , (N_det,N_states)] -&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth , (N_det,N_states)] -&BEGIN_PROVIDER [double precision, s2_eigvec_tc_bi_orth , (N_states) ] -&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ] + BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth , (N_states) ] +&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth , (N_states) ] +&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth , (N_det,N_states)] +&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth , (N_det,N_states)] +&BEGIN_PROVIDER [double precision, s2_eigvec_tc_bi_orth , (N_states) ] +&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth , (N_states) ] BEGIN_DOC ! eigenvalues, right and left eigenvectors of the transcorrelated Hamiltonian on the BI-ORTHO basis @@ -309,13 +309,13 @@ end deallocate(Stmp) print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ', leigvec_tc_bi_orth(1,1), reigvec_tc_bi_orth(1,1) + norm_ground_left_right_bi_orth = 0.d0 do i = 1, N_states - norm_ground_left_right_bi_orth = 0.d0 do j = 1, N_det - norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,i) * reigvec_tc_bi_orth(j,i) + norm_ground_left_right_bi_orth(i) += leigvec_tc_bi_orth(j,i) * reigvec_tc_bi_orth(j,i) enddo print*,' state ', i - print*,' norm l/r = ', norm_ground_left_right_bi_orth + print*,' norm l/r = ', norm_ground_left_right_bi_orth(i) print*,' = ', s2_eigvec_tc_bi_orth(i) enddo