From 26bdbf7193b6f378b7cce2e88fd61fdaf4d365e9 Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 8 Feb 2023 10:55:03 +0100 Subject: [PATCH] modified cipsi_tc_bi_ortho/selection.irp.f --- .../pt2_stoch_routines.irp.f | 2 +- src/cipsi_tc_bi_ortho/selection.irp.f | 257 ++++++++++-------- src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 22 +- src/iterations_tc/EZFIO.cfg | 24 ++ src/iterations_tc/NEED | 0 src/iterations_tc/io.irp.f | 37 +++ src/iterations_tc/iterations.irp.f | 43 +++ src/iterations_tc/print_extrapolation.irp.f | 46 ++++ src/iterations_tc/print_summary.irp.f | 104 +++++++ .../{12.tc_scf.bats => 12.tc_bi_ortho.bats} | 0 10 files changed, 417 insertions(+), 118 deletions(-) create mode 100644 src/iterations_tc/EZFIO.cfg create mode 100644 src/iterations_tc/NEED create mode 100644 src/iterations_tc/io.irp.f create mode 100644 src/iterations_tc/iterations.irp.f create mode 100644 src/iterations_tc/print_extrapolation.irp.f create mode 100644 src/iterations_tc/print_summary.irp.f rename src/tc_bi_ortho/{12.tc_scf.bats => 12.tc_bi_ortho.bats} (100%) diff --git a/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f b/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f index e146efb1..027b74c5 100644 --- a/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f +++ b/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f @@ -162,7 +162,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) TOUCH state_average_weight pt2_stoch_istate selection_weight PROVIDE nproc pt2_F mo_two_e_integrals_in_map mo_one_e_integrals pt2_w - PROVIDE psi_selectors pt2_u pt2_J pt2_R + PROVIDE pt2_u pt2_J pt2_R call new_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') integer, external :: zmq_put_psi diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f index b293946a..13e6c510 100644 --- a/src/cipsi_tc_bi_ortho/selection.irp.f +++ b/src/cipsi_tc_bi_ortho/selection.irp.f @@ -1,49 +1,148 @@ use bitmasks -! --- +subroutine get_mask_phase(det1, pm, Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det1(Nint,2) + integer(bit_kind), intent(out) :: pm(Nint,2) + integer(bit_kind) :: tmp1, tmp2 + integer :: i + tmp1 = 0_8 + tmp2 = 0_8 + select case (Nint) -subroutine select_connected(i_generator, E0, pt2_data, b, subset, csubset) +BEGIN_TEMPLATE + case ($Nint) + do i=1,$Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do +SUBST [ Nint ] +1;; +2;; +3;; +4;; +END_TEMPLATE + case default + do i=1,Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do + end select +end subroutine + + +subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset) use bitmasks use selection_types - implicit none - integer, intent(in) :: i_generator, subset, csubset - double precision, intent(in) :: E0(N_states) + integer, intent(in) :: i_generator, subset, csubset type(selection_buffer), intent(inout) :: b - type(pt2_type), intent(inout) :: pt2_data + type(pt2_type), intent(inout) :: pt2_data + integer :: k,l + double precision, intent(in) :: E0(N_states) - integer :: k, l - integer(bit_kind) :: hole_mask(N_int,2), particle_mask(N_int,2) - double precision, allocatable :: fock_diag_tmp(:,:) + integer(bit_kind) :: hole_mask(N_int,2), particle_mask(N_int,2) + + double precision, allocatable :: fock_diag_tmp(:,:) allocate(fock_diag_tmp(2,mo_num+1)) call build_fock_tmp_tc(fock_diag_tmp, psi_det_generators(1,1,i_generator), N_int) - do k = 1, N_int - hole_mask(k,1) = iand(generators_bitmask(k,1,s_hole), psi_det_generators(k,1,i_generator)) - hole_mask(k,2) = iand(generators_bitmask(k,2,s_hole), psi_det_generators(k,2,i_generator)) + do k=1,N_int + hole_mask(k,1) = iand(generators_bitmask(k,1,s_hole), psi_det_generators(k,1,i_generator)) + hole_mask(k,2) = iand(generators_bitmask(k,2,s_hole), psi_det_generators(k,2,i_generator)) particle_mask(k,1) = iand(generators_bitmask(k,1,s_part), not(psi_det_generators(k,1,i_generator)) ) particle_mask(k,2) = iand(generators_bitmask(k,2,s_part), not(psi_det_generators(k,2,i_generator)) ) enddo - call select_singles_and_doubles(i_generator, hole_mask, particle_mask, fock_diag_tmp, E0, pt2_data, b, subset, csubset) - + call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,b,subset,csubset) deallocate(fock_diag_tmp) - end subroutine select_connected -! --- + + +double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2, Nint) + use bitmasks + implicit none + + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: phasemask(Nint,2) + integer, intent(in) :: s1, s2, h1, h2, p1, p2 + logical :: change + integer :: np + double precision, save :: res(0:1) = (/1d0, -1d0/) + + integer :: h1_int, h2_int + integer :: p1_int, p2_int + integer :: h1_bit, h2_bit + integer :: p1_bit, p2_bit + h1_int = shiftr(h1-1,bit_kind_shift)+1 + h1_bit = h1 - shiftl(h1_int-1,bit_kind_shift)-1 + + h2_int = shiftr(h2-1,bit_kind_shift)+1 + h2_bit = h2 - shiftl(h2_int-1,bit_kind_shift)-1 + + p1_int = shiftr(p1-1,bit_kind_shift)+1 + p1_bit = p1 - shiftl(p1_int-1,bit_kind_shift)-1 + + p2_int = shiftr(p2-1,bit_kind_shift)+1 + p2_bit = p2 - shiftl(p2_int-1,bit_kind_shift)-1 + + + ! Put the phasemask bits at position 0, and add them all + h1_bit = int(shiftr(phasemask(h1_int,s1),h1_bit)) + p1_bit = int(shiftr(phasemask(p1_int,s1),p1_bit)) + h2_bit = int(shiftr(phasemask(h2_int,s2),h2_bit)) + p2_bit = int(shiftr(phasemask(p2_int,s2),p2_bit)) + + np = h1_bit + p1_bit + h2_bit + p2_bit + + if(p1 < h1) np = np + 1 + if(p2 < h2) np = np + 1 + + if(s1 == s2 .and. max(h1, p1) > min(h2, p2)) np = np + 1 + get_phase_bi = res(iand(np,1)) +end + subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock_diag_tmp, E0, pt2_data, buf, subset, csubset) - - BEGIN_DOC - ! WARNING /!\ : It is assumed that the generators and selectors are psi_det_sorted_tc - END_DOC - use bitmasks use selection_types implicit none + BEGIN_DOC + ! WARNING /!\ : It is assumed that the generators and selectors are psi_det_sorted_tc + END_DOC integer, intent(in) :: i_generator, subset, csubset integer(bit_kind), intent(in) :: hole_mask(N_int,2), particle_mask(N_int,2) @@ -89,7 +188,8 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock monoAdo = .true. monoBdo = .true. - do k = 1, N_int + + do k=1,N_int hole (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1)) hole (k,2) = iand(psi_det_generators(k,2,i_generator), hole_mask(k,2)) particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), particle_mask(k,1)) @@ -102,35 +202,37 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock allocate( indices(N_det), exc_degree( max(N_det_alpha_unique, N_det_beta_unique) ) ) ! Pre-compute excitation degrees wrt alpha determinants - k = 1 - do i = 1, N_det_alpha_unique - call get_excitation_degree_spin(psi_det_alpha_unique(1,i), psi_det_generators(1,1,i_generator), exc_degree(i), N_int) + k=1 + do i=1,N_det_alpha_unique + call get_excitation_degree_spin(psi_det_alpha_unique(1,i), & + psi_det_generators(1,1,i_generator), exc_degree(i), N_int) enddo ! Iterate on 0SD beta, and find alphas 0SDTQ such that exc_degree <= 4 - do j = 1, N_det_beta_unique - call get_excitation_degree_spin(psi_det_beta_unique(1,j), psi_det_generators(1,2,i_generator), nt, N_int) + do j=1,N_det_beta_unique + call get_excitation_degree_spin(psi_det_beta_unique(1,j), & + psi_det_generators(1,2,i_generator), nt, N_int) if (nt > 2) cycle - do l_a = psi_bilinear_matrix_columns_loc(j), psi_bilinear_matrix_columns_loc(j+1)-1 + do l_a=psi_bilinear_matrix_columns_loc(j), psi_bilinear_matrix_columns_loc(j+1)-1 i = psi_bilinear_matrix_rows(l_a) if(nt + exc_degree(i) <= 4) then idx = psi_det_sorted_tc_order(psi_bilinear_matrix_order(l_a)) - if (psi_average_norm_contrib_sorted_tc(idx) > norm_thr) then +! if (psi_average_norm_contrib_sorted_tc(idx) > norm_thr) then indices(k) = idx k = k + 1 - endif +! endif endif enddo enddo ! Pre-compute excitation degrees wrt beta determinants - do i = 1, N_det_beta_unique + do i=1,N_det_beta_unique call get_excitation_degree_spin(psi_det_beta_unique(1,i), psi_det_generators(1,2,i_generator), exc_degree(i), N_int) enddo ! Iterate on 0S alpha, and find betas TQ such that exc_degree <= 4 ! Remove also contributions < 1.d-20) - do j = 1, N_det_alpha_unique + do j=1,N_det_alpha_unique call get_excitation_degree_spin(psi_det_alpha_unique(1,j), psi_det_generators(1,1,i_generator), nt, N_int) if (nt > 1) cycle do l_a = psi_bilinear_matrix_transp_rows_loc(j), psi_bilinear_matrix_transp_rows_loc(j+1)-1 @@ -140,10 +242,10 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock idx = psi_det_sorted_tc_order( & psi_bilinear_matrix_order( & psi_bilinear_matrix_transp_order(l_a))) - if(psi_average_norm_contrib_sorted_tc(idx) > norm_thr) then +! if(psi_average_norm_contrib_sorted_tc(idx) > norm_thr) then indices(k) = idx k = k + 1 - endif +! endif endif enddo enddo @@ -211,6 +313,9 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock allocate( mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num) ) maskInd = -1 + + + do s1 = 1, 2 do i1 = N_holes(s1), 1, -1 ! Generate low excitations first @@ -354,7 +459,6 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock fullinteresting(sze+1) = i endif enddo - allocate( fullminilist (N_int, 2, fullinteresting(0)), & minilist (N_int, 2, interesting(0)) ) @@ -579,16 +683,11 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int) if(nt == 4) then -! call get_d2 (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) -! call get_pm2(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i))) elseif(nt == 3) then -! call get_d1 (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) -! call get_pm1(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i))) else call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) -! call get_pm0(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i))) endif elseif(nt == 4) then call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) @@ -780,9 +879,19 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if(debug_tc_pt2 == 1)then !! Using the old version psi_h_alpha = 0.d0 alpha_h_psi = 0.d0 - do iii = 1, N_det + do iii = 1, N_det_selectors call htilde_mu_mat_bi_ortho_tot(psi_selectors(1,1,iii), det, N_int, i_h_alpha) call htilde_mu_mat_bi_ortho_tot(det, psi_selectors(1,1,iii), N_int, alpha_h_i) + call get_excitation_degree(psi_selectors(1,1,iii), det,degree,N_int) + if(degree == 0)then + print*,'problem !!!' + print*,'a determinant is already in the wave function !!' + print*,'it corresponds to the selector number ',iii + call debug_det(det,N_int) + stop + endif +! call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) +! call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function enddo @@ -791,7 +900,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version psi_h_alpha = 0.d0 alpha_h_psi = 0.d0 - do iii = 1, N_det ! old version + do iii = 1, N_det_selectors ! old version call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function @@ -881,70 +990,6 @@ end subroutine fill_buffer_double ! --- -subroutine get_mask_phase(det1, pm, Nint) - - use bitmasks - implicit none - - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: det1(Nint,2) - integer(bit_kind), intent(out) :: pm(Nint,2) - integer(bit_kind) :: tmp1, tmp2 - integer :: i - tmp1 = 0_8 - tmp2 = 0_8 - select case (Nint) - -BEGIN_TEMPLATE - case ($Nint) - do i=1,$Nint - pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) - pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) - pm(i,1) = ieor(pm(i,1), tmp1) - pm(i,2) = ieor(pm(i,2), tmp2) - if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) - if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) - end do -SUBST [ Nint ] -1;; -2;; -3;; -4;; -END_TEMPLATE - case default - do i=1,Nint - pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) - pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) - pm(i,1) = ieor(pm(i,1), tmp1) - pm(i,2) = ieor(pm(i,2), tmp2) - if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) - if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) - end do - end select - -end subroutine get_mask_phase - -! --- subroutine past_d1(bannedOrb, p) diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index d83c3689..c1e4af0c 100644 --- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -53,18 +53,18 @@ subroutine run_stochastic_cipsi ! call routine_save_right - if (N_det > N_det_max) then - psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det) - psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states) - N_det = N_det_max - soft_touch N_det psi_det psi_coef - if (s2_eig) then - call make_s2_eigenfunction - endif - print_pt2 = .False. - call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) +! if (N_det > N_det_max) then +! psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det) +! psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states) +! N_det = N_det_max +! soft_touch N_det psi_det psi_coef +! if (s2_eig) then +! call make_s2_eigenfunction +! endif +! print_pt2 = .False. +! call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) ! call routine_save_right - endif +! endif allocate(ept2(1000),pt1(1000),extrap_energy(100)) diff --git a/src/iterations_tc/EZFIO.cfg b/src/iterations_tc/EZFIO.cfg new file mode 100644 index 00000000..2a5e94a7 --- /dev/null +++ b/src/iterations_tc/EZFIO.cfg @@ -0,0 +1,24 @@ +[n_iter] +interface: ezfio +doc: Number of saved iterations +type:integer +default: 1 + +[n_det_iterations] +interface: ezfio, provider +doc: Number of determinants at each iteration +type: integer +size: (100) + +[energy_iterations] +interface: ezfio, provider +doc: The variational energy at each iteration +type: double precision +size: (determinants.n_states,100) + +[pt2_iterations] +interface: ezfio, provider +doc: The |PT2| correction at each iteration +type: double precision +size: (determinants.n_states,100) + diff --git a/src/iterations_tc/NEED b/src/iterations_tc/NEED new file mode 100644 index 00000000..e69de29b diff --git a/src/iterations_tc/io.irp.f b/src/iterations_tc/io.irp.f new file mode 100644 index 00000000..821f5e84 --- /dev/null +++ b/src/iterations_tc/io.irp.f @@ -0,0 +1,37 @@ +BEGIN_PROVIDER [ integer, n_iter ] + implicit none + BEGIN_DOC +! number of iterations + END_DOC + + logical :: has + PROVIDE ezfio_filename + if (mpi_master) then + + double precision :: zeros(N_states,100) + integer :: izeros(100) + zeros = 0.d0 + izeros = 0 + call ezfio_set_iterations_n_iter(0) + call ezfio_set_iterations_energy_iterations(zeros) + call ezfio_set_iterations_pt2_iterations(zeros) + call ezfio_set_iterations_n_det_iterations(izeros) + n_iter = 1 + endif + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( n_iter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read n_iter with MPI' + endif + IRP_ENDIF + + call write_time(6) + +END_PROVIDER + diff --git a/src/iterations_tc/iterations.irp.f b/src/iterations_tc/iterations.irp.f new file mode 100644 index 00000000..2f1cf0c1 --- /dev/null +++ b/src/iterations_tc/iterations.irp.f @@ -0,0 +1,43 @@ +BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter,N_states) ] + implicit none + BEGIN_DOC + ! Extrapolated energy, using E_var = f(PT2) where PT2=0 + END_DOC +! integer :: i + extrapolated_energy = 0.D0 +END_PROVIDER + + subroutine get_extrapolated_energy(Niter,ept2,pt1,extrap_energy) + implicit none + integer, intent(in) :: Niter + double precision, intent(in) :: ept2(Niter),pt1(Niter),extrap_energy(Niter) + call extrapolate_data(Niter,ept2,pt1,extrap_energy) + end + +subroutine save_iterations(e_, pt2_,n_) + implicit none + BEGIN_DOC +! Update the energy in the EZFIO file. + END_DOC + integer, intent(in) :: n_ + double precision, intent(in) :: e_(N_states), pt2_(N_states) + integer :: i + + if (N_iter == 101) then + do i=2,N_iter-1 + energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter) + pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter) + enddo + N_iter = N_iter-1 + TOUCH N_iter + endif + + energy_iterations(1:N_states,N_iter) = e_(1:N_states) + pt2_iterations(1:N_states,N_iter) = pt2_(1:N_states) + n_det_iterations(N_iter) = n_ + call ezfio_set_iterations_N_iter(N_iter) + call ezfio_set_iterations_energy_iterations(energy_iterations) + call ezfio_set_iterations_pt2_iterations(pt2_iterations) + call ezfio_set_iterations_n_det_iterations(n_det_iterations) +end + diff --git a/src/iterations_tc/print_extrapolation.irp.f b/src/iterations_tc/print_extrapolation.irp.f new file mode 100644 index 00000000..cb46fb67 --- /dev/null +++ b/src/iterations_tc/print_extrapolation.irp.f @@ -0,0 +1,46 @@ +subroutine print_extrapolated_energy + implicit none + BEGIN_DOC +! Print the extrapolated energy in the output + END_DOC + + integer :: i,k + + if (N_iter< 2) then + return + endif + write(*,'(A)') '' + write(*,'(A)') 'Extrapolated energies' + write(*,'(A)') '------------------------' + write(*,'(A)') '' + + print *, '' + print *, 'State ', 1 + print *, '' + write(*,*) '=========== ', '===================' + write(*,*) 'minimum PT2 ', 'Extrapolated energy' + write(*,*) '=========== ', '===================' + do k=2,min(N_iter,8) + write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter+1-k), extrapolated_energy(k,1) + enddo + write(*,*) '=========== ', '===================' + + do i=2, min(N_states,N_det) + print *, '' + print *, 'State ', i + print *, '' + write(*,*) '=========== ', '=================== ', '=================== ', '===================' + write(*,*) 'minimum PT2 ', 'Extrapolated energy ', ' Excitation (a.u) ', ' Excitation (eV) ' + write(*,*) '=========== ', '=================== ', '=================== ', '===================' + do k=2,min(N_iter,8) + write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), & + extrapolated_energy(k,i) - extrapolated_energy(k,1), & + (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0 + enddo + write(*,*) '=========== ', '=================== ', '=================== ', '===================' + enddo + + print *, '' + +end subroutine + diff --git a/src/iterations_tc/print_summary.irp.f b/src/iterations_tc/print_summary.irp.f new file mode 100644 index 00000000..8e6285e2 --- /dev/null +++ b/src/iterations_tc/print_summary.irp.f @@ -0,0 +1,104 @@ +subroutine print_summary(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 + + if (do_pt2) then + pt2_string = ' ' + else + pt2_string = '(approx)' + endif + + N_states_p = min(N_det_,n_st) + + 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_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p) + write(*,fmt) '# rPT2'//pt2_string, (pt2_data % rpt2(k), pt2_data_err % rpt2(k), k=1,N_states_p) + write(*,'(A)') '#' + write(*,fmt) '# E+PT2 ', (e_(k)+pt2_data % pt2(k),pt2_data_err % pt2(k), k=1,N_states_p) + write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % rpt2(k),pt2_data_err % rpt2(k), k=1,N_states_p) + if (N_states_p > 1) then + write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), & + dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p) + write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, & + dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(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 *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k) + print *, 'PT norm = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k)) + print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) + print *, 'rPT2 = ', pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k) + print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) + print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(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_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), & + (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0 + enddo + print *, '-----' + print*, 'Variational + renormalized perturbative Energy difference (au | eV)' + do i=2, N_states_p + print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), & + (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0 + enddo + endif + +! call print_energy_components() + +end subroutine + diff --git a/src/tc_bi_ortho/12.tc_scf.bats b/src/tc_bi_ortho/12.tc_bi_ortho.bats similarity index 100% rename from src/tc_bi_ortho/12.tc_scf.bats rename to src/tc_bi_ortho/12.tc_bi_ortho.bats