From 9175fb21c9dcbe931f89d96cf1297221693d5fde Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 12 Mar 2024 14:05:38 +0100 Subject: [PATCH 01/12] 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 From 6e35f8f8f8735bd4a898fabbc6bf552f382e517a Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 12 Mar 2024 15:30:52 +0100 Subject: [PATCH 02/12] 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 + From 0ef067337d9cffd2fba9b1bc29afe071c696f883 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 12 Mar 2024 16:37:16 +0100 Subject: [PATCH 03/12] Introducing cipsi_utils for CIPSI and TC-CIPSI --- plugins/local/cipsi_tc_bi_ortho/NEED | 1 + plugins/local/cipsi_tc_bi_ortho/cipsi.irp.f | 2 +- .../pt2_stoch_routines.irp.f | 869 +--------------- .../cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 2 +- plugins/local/fci_tc_bi/NEED | 1 + plugins/local/fci_tc_bi/selectors.irp.f | 2 +- src/cipsi/NEED | 1 + src/cipsi/pt2_stoch_routines.irp.f | 924 +----------------- src/generators_full_tc/README.rst | 9 + .../generators_full_tc}/generators.irp.f | 48 +- 10 files changed, 56 insertions(+), 1803 deletions(-) create mode 100644 src/generators_full_tc/README.rst rename {plugins/local/fci_tc_bi => src/generators_full_tc}/generators.irp.f (51%) diff --git a/plugins/local/cipsi_tc_bi_ortho/NEED b/plugins/local/cipsi_tc_bi_ortho/NEED index 8f05be69..d329326c 100644 --- a/plugins/local/cipsi_tc_bi_ortho/NEED +++ b/plugins/local/cipsi_tc_bi_ortho/NEED @@ -1,3 +1,4 @@ +cipsi_utils json mpi perturbation diff --git a/plugins/local/cipsi_tc_bi_ortho/cipsi.irp.f b/plugins/local/cipsi_tc_bi_ortho/cipsi.irp.f index fb907cb3..65e0790a 100644 --- a/plugins/local/cipsi_tc_bi_ortho/cipsi.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/cipsi.irp.f @@ -65,7 +65,7 @@ subroutine run_cipsi if (N_det > N_det_max) then psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(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) + psi_coef(1:N_det,1:N_states) = psi_coef_sorted_gen(1:N_det,1:N_states) N_det = N_det_max soft_touch N_det psi_det psi_coef if (s2_eig) then diff --git a/plugins/local/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f b/plugins/local/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f index 284b2bc8..6e1a6748 100644 --- a/plugins/local/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f @@ -1,868 +1,3 @@ -BEGIN_PROVIDER [ integer, pt2_stoch_istate ] - implicit none - BEGIN_DOC - ! State for stochatsic PT2 - END_DOC - pt2_stoch_istate = 1 -END_PROVIDER - - BEGIN_PROVIDER [ integer, pt2_F, (N_det_generators) ] -&BEGIN_PROVIDER [ integer, pt2_n_tasks_max ] - implicit none - logical, external :: testTeethBuilding - integer :: i,j - pt2_n_tasks_max = elec_alpha_num*elec_alpha_num + elec_alpha_num*elec_beta_num - n_core_orb*2 - pt2_n_tasks_max = min(pt2_n_tasks_max,1+N_det_generators/10000) - call write_int(6,pt2_n_tasks_max,'pt2_n_tasks_max') - - pt2_F(:) = max(int(sqrt(float(pt2_n_tasks_max))),1) - do i=1,pt2_n_0(1+pt2_N_teeth/4) - pt2_F(i) = pt2_n_tasks_max*pt2_min_parallel_tasks - enddo - do i=1+pt2_n_0(pt2_N_teeth-pt2_N_teeth/4), pt2_n_0(pt2_N_teeth-pt2_N_teeth/10) - pt2_F(i) = pt2_min_parallel_tasks - enddo - do i=1+pt2_n_0(pt2_N_teeth-pt2_N_teeth/10), N_det_generators - pt2_F(i) = 1 - enddo - -END_PROVIDER - - BEGIN_PROVIDER [ integer, pt2_N_teeth ] -&BEGIN_PROVIDER [ integer, pt2_minDetInFirstTeeth ] - implicit none - logical, external :: testTeethBuilding - - if(N_det_generators < 500) then - pt2_minDetInFirstTeeth = 1 - pt2_N_teeth = 1 - else - pt2_minDetInFirstTeeth = min(5, N_det_generators) - do pt2_N_teeth=100,2,-1 - if(testTeethBuilding(pt2_minDetInFirstTeeth, pt2_N_teeth)) exit - end do - end if - call write_int(6,pt2_N_teeth,'Number of comb teeth') -END_PROVIDER - - -logical function testTeethBuilding(minF, N) - implicit none - integer, intent(in) :: minF, N - integer :: n0, i - double precision :: u0, Wt, r - - double precision, allocatable :: tilde_w(:), tilde_cW(:) - integer, external :: dress_find_sample - - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - - rss = memory_of_double(2*N_det_generators+1) - call check_mem(rss,irp_here) - - allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators)) - - double precision :: norm2 - norm2 = 0.d0 - do i=N_det_generators,1,-1 - tilde_w(i) = psi_coef_sorted_tc_gen(i,pt2_stoch_istate) * & - psi_coef_sorted_tc_gen(i,pt2_stoch_istate) - norm2 = norm2 + tilde_w(i) - enddo - - f = 1.d0/norm2 - tilde_w(:) = tilde_w(:) * f - - tilde_cW(0) = -1.d0 - do i=1,N_det_generators - tilde_cW(i) = tilde_cW(i-1) + tilde_w(i) - enddo - tilde_cW(:) = tilde_cW(:) + 1.d0 - deallocate(tilde_w) - - n0 = 0 - testTeethBuilding = .false. - double precision :: f - integer :: minFN - minFN = N_det_generators - minF * N - f = 1.d0/dble(N) - do - u0 = tilde_cW(n0) - r = tilde_cW(n0 + minF) - Wt = (1d0 - u0) * f - if (dabs(Wt) <= 1.d-3) then - exit - endif - if(Wt >= r - u0) then - testTeethBuilding = .true. - exit - end if - n0 += 1 - if(n0 > minFN) then - exit - end if - end do - deallocate(tilde_cW) - -end function - - - -subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) - use f77_zmq - use selection_types - - implicit none - - integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull - integer, intent(in) :: N_in -! integer, intent(inout) :: N_in - double precision, intent(in) :: relative_error, E(N_states) - type(pt2_type), intent(inout) :: pt2_data, pt2_data_err -! - integer :: i, N - - double precision :: state_average_weight_save(N_states), w(N_states,4) - integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket - type(selection_buffer) :: b - - PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique - PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order - PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp_tc psi_det_sorted_tc - PROVIDE psi_det_hii selection_weight pseudo_sym - PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE excitation_beta_max excitation_alpha_max excitation_max - - if (h0_type == 'CFG') then - PROVIDE psi_configuration_hii det_to_configuration - endif - - if (N_det <= max(4,N_states) .or. pt2_N_teeth < 2) then - print*,'ZMQ_selection' - call ZMQ_selection(N_in, pt2_data) - else - print*,'else ZMQ_selection' - - N = max(N_in,1) * N_states - state_average_weight_save(:) = state_average_weight(:) - if (int(N,8)*2_8 > huge(1)) then - print *, irp_here, ': integer too large' - stop -1 - endif - call create_selection_buffer(N, N*2, b) - ASSERT (associated(b%det)) - ASSERT (associated(b%val)) - - do pt2_stoch_istate=1,N_states - state_average_weight(:) = 0.d0 - state_average_weight(pt2_stoch_istate) = 1.d0 - 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 pt2_u pt2_J pt2_R - call new_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') - - integer, external :: zmq_put_psi - integer, external :: zmq_put_N_det_generators - integer, external :: zmq_put_N_det_selectors - integer, external :: zmq_put_dvector - integer, external :: zmq_put_ivector - if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then - stop 'Unable to put psi on ZMQ server' - endif - if (zmq_put_N_det_generators(zmq_to_qp_run_socket, 1) == -1) then - stop 'Unable to put N_det_generators on ZMQ server' - endif - if (zmq_put_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) then - stop 'Unable to put N_det_selectors on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',pt2_e0_denominator,size(pt2_e0_denominator)) == -1) then - stop 'Unable to put energy on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then - stop 'Unable to put state_average_weight on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then - stop 'Unable to put selection_weight on ZMQ server' - endif - if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then - stop 'Unable to put pt2_stoch_istate on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',(/threshold_generators/),1) == -1) then - stop 'Unable to put threshold_generators on ZMQ server' - endif - - - integer, external :: add_task_to_taskserver - character(300000) :: task - - integer :: j,k,ipos,ifirst - ifirst=0 - - ipos=0 - do i=1,N_det_generators - if (pt2_F(i) > 1) then - ipos += 1 - endif - enddo - call write_int(6,sum(pt2_F),'Number of tasks') - call write_int(6,ipos,'Number of fragmented tasks') - - ipos=1 - do i= 1, N_det_generators - do j=1,pt2_F(pt2_J(i)) - write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') j, pt2_J(i), N_in - ipos += 30 - if (ipos > 300000-30) then - if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then - stop 'Unable to add task to task server' - endif - ipos=1 - if (ifirst == 0) then - ifirst=1 - if (zmq_set_running(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Failed in zmq_set_running' - endif - endif - endif - end do - enddo - if (ipos > 1) then - if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then - stop 'Unable to add task to task server' - endif - endif - - integer, external :: zmq_set_running - if (zmq_set_running(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Failed in zmq_set_running' - endif - - - double precision :: mem_collector, mem, rss - - call resident_memory(rss) - - mem_collector = 8.d0 * & ! bytes - ( 1.d0*pt2_n_tasks_max & ! task_id, index - + 0.635d0*N_det_generators & ! f,d - + pt2_n_tasks_max*pt2_type_size(N_states) & ! pt2_data_task - + N_det_generators*pt2_type_size(N_states) & ! pt2_data_I - + 4.d0*(pt2_N_teeth+1) & ! S, S2, T2, T3 - + 1.d0*(N_int*2.d0*N + N) & ! selection buffer - + 1.d0*(N_int*2.d0*N + N) & ! sort selection buffer - ) / 1024.d0**3 - - integer :: nproc_target, ii - nproc_target = nthreads_pt2 - ii = min(N_det, (elec_alpha_num*(mo_num-elec_alpha_num))**2) - - do - mem = mem_collector + & ! - nproc_target * 8.d0 * & ! bytes - ( 0.5d0*pt2_n_tasks_max & ! task_id - + 64.d0*pt2_n_tasks_max & ! task - + pt2_type_size(N_states)*pt2_n_tasks_max*N_states & ! pt2, variance, overlap - + 1.d0*pt2_n_tasks_max & ! i_generator, subset - + 1.d0*(N_int*2.d0*ii+ ii) & ! selection buffer - + 1.d0*(N_int*2.d0*ii+ ii) & ! sort selection buffer - + 2.0d0*(ii) & ! preinteresting, interesting, - ! prefullinteresting, fullinteresting - + 2.0d0*(N_int*2*ii) & ! minilist, fullminilist - + 1.0d0*(N_states*mo_num*mo_num) & ! mat - ) / 1024.d0**3 - - if (nproc_target == 0) then - call check_mem(mem,irp_here) - nproc_target = 1 - exit - endif - - if (mem+rss < qp_max_mem) then - exit - endif - - nproc_target = nproc_target - 1 - - enddo - call write_int(6,nproc_target,'Number of threads for PT2') - call write_double(6,mem,'Memory (Gb)') - - call omp_set_max_active_levels(1) - - - print '(A)', '========== ======================= ===================== ===================== ===========' - print '(A)', ' Samples Energy Variance Norm^2 Seconds' - print '(A)', '========== ======================= ===================== ===================== ===========' - - PROVIDE global_selection_buffer - - !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc_target+1) & - !$OMP PRIVATE(i) - i = omp_get_thread_num() - if (i==0) then - - call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, pt2_data, pt2_data_err, b, N) - pt2_data % rpt2(pt2_stoch_istate) = & - pt2_data % pt2(pt2_stoch_istate)/(1.d0+pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate)) - - !TODO : We should use here the correct formula for the error of X/Y - pt2_data_err % rpt2(pt2_stoch_istate) = & - pt2_data_err % pt2(pt2_stoch_istate)/(1.d0 + pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate)) - - else - call pt2_slave_inproc(i) - endif - !$OMP END PARALLEL - call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') - call omp_set_max_active_levels(8) - - print '(A)', '========== ======================= ===================== ===================== ===========' - - do k=1,N_states - pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) - enddo - SOFT_TOUCH pt2_overlap - - enddo - FREE pt2_stoch_istate - - ! Symmetrize overlap - do j=2,N_states - do i=1,j-1 - pt2_overlap(i,j) = 0.5d0 * (pt2_overlap(i,j) + pt2_overlap(j,i)) - pt2_overlap(j,i) = pt2_overlap(i,j) - enddo - enddo - - print *, 'Overlap of perturbed states:' - do k=1,N_states - print *, pt2_overlap(k,:) - enddo - print *, '-------' - - if (N_in > 0) then - b%cur = min(N_in,b%cur) - if (s2_eig) then - call make_selection_buffer_s2(b) - else - call remove_duplicates_in_selection_buffer(b) - endif - call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) - endif - call delete_selection_buffer(b) - - state_average_weight(:) = state_average_weight_save(:) - TOUCH state_average_weight - call update_pt2_and_variance_weights(pt2_data, N_states) - endif - - -end subroutine - - -subroutine pt2_slave_inproc(i) - implicit none - integer, intent(in) :: i - - PROVIDE global_selection_buffer - call run_pt2_slave(1,i,pt2_e0_denominator) +subroutine provide_for_zmq_pt2 + PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc psi_det_sorted_tc_order end - - -subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_err, b, N_) - use f77_zmq - use selection_types - use bitmasks - implicit none - - - integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - double precision, intent(in) :: relative_error, E - type(pt2_type), intent(inout) :: pt2_data, pt2_data_err - type(selection_buffer), intent(inout) :: b - integer, intent(in) :: N_ - - type(pt2_type), allocatable :: pt2_data_task(:) - type(pt2_type), allocatable :: pt2_data_I(:) - type(pt2_type), allocatable :: pt2_data_S(:) - type(pt2_type), allocatable :: pt2_data_S2(:) - type(pt2_type) :: pt2_data_teeth - integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - integer, external :: zmq_delete_tasks_async_send - integer, external :: zmq_delete_tasks_async_recv - integer, external :: zmq_abort - integer, external :: pt2_find_sample_lr - - PROVIDE pt2_stoch_istate - - integer :: more, n, i, p, c, t, n_tasks, U - integer, allocatable :: task_id(:) - integer, allocatable :: index(:) - - double precision :: v, x, x2, x3, avg, avg2, avg3(N_states), eqt, E0, v0, n0(N_states) - double precision :: eqta(N_states) - double precision :: time, time1, time0 - - integer, allocatable :: f(:) - logical, allocatable :: d(:) - logical :: do_exit, stop_now, sending - logical, external :: qp_stop - type(selection_buffer) :: b2 - - - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - - sending =.False. - - rss = memory_of_int(pt2_n_tasks_max*2+N_det_generators*2) - rss += memory_of_double(N_states*N_det_generators)*3.d0 - rss += memory_of_double(N_states*pt2_n_tasks_max)*3.d0 - rss += memory_of_double(pt2_N_teeth+1)*4.d0 - call check_mem(rss,irp_here) - - ! If an allocation is added here, the estimate of the memory should also be - ! updated in ZMQ_pt2 - allocate(task_id(pt2_n_tasks_max), index(pt2_n_tasks_max), f(N_det_generators)) - allocate(d(N_det_generators+1)) - allocate(pt2_data_task(pt2_n_tasks_max)) - allocate(pt2_data_I(N_det_generators)) - allocate(pt2_data_S(pt2_N_teeth+1)) - allocate(pt2_data_S2(pt2_N_teeth+1)) - - - - zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() - call create_selection_buffer(N_, N_*2, b2) - - - pt2_data % pt2(pt2_stoch_istate) = -huge(1.) - pt2_data_err % pt2(pt2_stoch_istate) = huge(1.) - pt2_data % variance(pt2_stoch_istate) = huge(1.) - pt2_data_err % variance(pt2_stoch_istate) = huge(1.) - pt2_data % overlap(:,pt2_stoch_istate) = 0.d0 - pt2_data_err % overlap(:,pt2_stoch_istate) = huge(1.) - n = 1 - t = 0 - U = 0 - do i=1,pt2_n_tasks_max - call pt2_alloc(pt2_data_task(i),N_states) - enddo - do i=1,pt2_N_teeth+1 - call pt2_alloc(pt2_data_S(i),N_states) - call pt2_alloc(pt2_data_S2(i),N_states) - enddo - do i=1,N_det_generators - call pt2_alloc(pt2_data_I(i),N_states) - enddo - f(:) = pt2_F(:) - d(:) = .false. - n_tasks = 0 - E0 = E - v0 = 0.d0 - n0(:) = 0.d0 - more = 1 - call wall_time(time0) - time1 = time0 - - do_exit = .false. - stop_now = .false. - do while (n <= N_det_generators) - if(f(pt2_J(n)) == 0) then - d(pt2_J(n)) = .true. - do while(d(U+1)) - U += 1 - end do - - ! Deterministic part - do while(t <= pt2_N_teeth) - if(U >= pt2_n_0(t+1)) then - t=t+1 - E0 = 0.d0 - v0 = 0.d0 - n0(:) = 0.d0 - do i=pt2_n_0(t),1,-1 - E0 += pt2_data_I(i) % pt2(pt2_stoch_istate) - v0 += pt2_data_I(i) % variance(pt2_stoch_istate) - n0(:) += pt2_data_I(i) % overlap(:,pt2_stoch_istate) - end do - else - exit - end if - end do - - ! Add Stochastic part - c = pt2_R(n) - if(c > 0) then - - call pt2_alloc(pt2_data_teeth,N_states) - do p=pt2_N_teeth, 1, -1 - v = pt2_u_0 + pt2_W_T * (pt2_u(c) + dble(p-1)) - i = pt2_find_sample_lr(v, pt2_cW,pt2_n_0(p),pt2_n_0(p+1)) - v = pt2_W_T / pt2_w(i) - call pt2_add ( pt2_data_teeth, v, pt2_data_I(i) ) - call pt2_add ( pt2_data_S(p), 1.d0, pt2_data_teeth ) - call pt2_add2( pt2_data_S2(p), 1.d0, pt2_data_teeth ) - enddo - call pt2_dealloc(pt2_data_teeth) - - avg = E0 + pt2_data_S(t) % pt2(pt2_stoch_istate) / dble(c) - avg2 = v0 + pt2_data_S(t) % variance(pt2_stoch_istate) / dble(c) - avg3(:) = n0(:) + pt2_data_S(t) % overlap(:,pt2_stoch_istate) / dble(c) - if ((avg /= 0.d0) .or. (n == N_det_generators) ) then - do_exit = .true. - endif - if (qp_stop()) then - stop_now = .True. - endif - pt2_data % pt2(pt2_stoch_istate) = avg - pt2_data % variance(pt2_stoch_istate) = avg2 - pt2_data % overlap(:,pt2_stoch_istate) = avg3(:) - call wall_time(time) - ! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969) - if(c > 2) then - eqt = dabs((pt2_data_S2(t) % pt2(pt2_stoch_istate) / c) - (pt2_data_S(t) % pt2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = sqrt(eqt / (dble(c) - 1.5d0)) - pt2_data_err % pt2(pt2_stoch_istate) = eqt - - eqt = dabs((pt2_data_S2(t) % variance(pt2_stoch_istate) / c) - (pt2_data_S(t) % variance(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = sqrt(eqt / (dble(c) - 1.5d0)) - pt2_data_err % variance(pt2_stoch_istate) = eqt - - eqta(:) = dabs((pt2_data_S2(t) % overlap(:,pt2_stoch_istate) / c) - (pt2_data_S(t) % overlap(:,pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqta(:) = sqrt(eqta(:) / (dble(c) - 1.5d0)) - pt2_data_err % overlap(:,pt2_stoch_istate) = eqta(:) - - - if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then - time1 = time - print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.4)', c, & - pt2_data % pt2(pt2_stoch_istate) +E, & - pt2_data_err % pt2(pt2_stoch_istate), & - pt2_data % variance(pt2_stoch_istate), & - pt2_data_err % variance(pt2_stoch_istate), & - pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate), & - pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), & - time-time0 - if (stop_now .or. ( & - (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & - (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then - if (zmq_abort(zmq_to_qp_run_socket) == -1) then - call sleep(10) - if (zmq_abort(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Error in sending abort signal (2)' - endif - endif - endif - endif - endif - end if - n += 1 - else if(more == 0) then - exit - else - call pull_pt2_results(zmq_socket_pull, index, pt2_data_task, task_id, n_tasks, b2) - if(n_tasks > pt2_n_tasks_max)then - print*,'PB !!!' - print*,'If you see this, send a bug report with the following content' - print*,irp_here - print*,'n_tasks,pt2_n_tasks_max = ',n_tasks,pt2_n_tasks_max - stop -1 - endif - if (zmq_delete_tasks_async_send(zmq_to_qp_run_socket,task_id,n_tasks,sending) == -1) then - stop 'PT2: Unable to delete tasks (send)' - endif - do i=1,n_tasks - if(index(i).gt.size(pt2_data_I,1).or.index(i).lt.1)then - print*,'PB !!!' - print*,'If you see this, send a bug report with the following content' - print*,irp_here - print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1) - stop -1 - endif - call pt2_add(pt2_data_I(index(i)),1.d0,pt2_data_task(i)) - f(index(i)) -= 1 - end do - do i=1, b2%cur - ! We assume the pulled buffer is sorted - if (b2%val(i) > b%mini) exit - call add_to_selection_buffer(b, b2%det(1,1,i), b2%val(i)) - end do - if (zmq_delete_tasks_async_recv(zmq_to_qp_run_socket,more,sending) == -1) then - stop 'PT2: Unable to delete tasks (recv)' - endif - end if - end do - do i=1,N_det_generators - call pt2_dealloc(pt2_data_I(i)) - enddo - do i=1,pt2_N_teeth+1 - call pt2_dealloc(pt2_data_S(i)) - call pt2_dealloc(pt2_data_S2(i)) - enddo - do i=1,pt2_n_tasks_max - call pt2_dealloc(pt2_data_task(i)) - enddo -!print *, 'deleting b2' - call delete_selection_buffer(b2) -!print *, 'sorting b' - call sort_selection_buffer(b) -!print *, 'done' - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - -end subroutine - - -integer function pt2_find_sample(v, w) - implicit none - double precision, intent(in) :: v, w(0:N_det_generators) - integer, external :: pt2_find_sample_lr - - pt2_find_sample = pt2_find_sample_lr(v, w, 0, N_det_generators) -end function - - -integer function pt2_find_sample_lr(v, w, l_in, r_in) - implicit none - double precision, intent(in) :: v, w(0:N_det_generators) - integer, intent(in) :: l_in,r_in - integer :: i,l,r - - l=l_in - r=r_in - - do while(r-l > 1) - i = shiftr(r+l,1) - if(w(i) < v) then - l = i - else - r = i - end if - end do - i = r - do r=i+1,N_det_generators - if (w(r) /= w(i)) then - exit - endif - enddo - pt2_find_sample_lr = r-1 -end function - - -BEGIN_PROVIDER [ integer, pt2_n_tasks ] - implicit none - BEGIN_DOC - ! Number of parallel tasks for the Monte Carlo - END_DOC - pt2_n_tasks = N_det_generators -END_PROVIDER - -BEGIN_PROVIDER[ double precision, pt2_u, (N_det_generators)] - implicit none - integer, allocatable :: seed(:) - integer :: m,i - call random_seed(size=m) - allocate(seed(m)) - do i=1,m - seed(i) = i - enddo - call random_seed(put=seed) - deallocate(seed) - - call RANDOM_NUMBER(pt2_u) - END_PROVIDER - - BEGIN_PROVIDER[ integer, pt2_J, (N_det_generators)] -&BEGIN_PROVIDER[ integer, pt2_R, (N_det_generators)] - implicit none - BEGIN_DOC -! pt2_J contains the list of generators after ordering them according to the -! Monte Carlo sampling. -! -! pt2_R(i) is the number of combs drawn when determinant i is computed. - END_DOC - integer :: N_c, N_j - integer :: U, t, i - double precision :: v - integer, external :: pt2_find_sample_lr - - logical, allocatable :: pt2_d(:) - integer :: m,l,r,k - integer :: ncache - integer, allocatable :: ii(:,:) - double precision :: dt - - ncache = min(N_det_generators,10000) - - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - rss = memory_of_int(ncache)*dble(pt2_N_teeth) + memory_of_int(N_det_generators) - call check_mem(rss,irp_here) - - allocate(ii(pt2_N_teeth,ncache),pt2_d(N_det_generators)) - - pt2_R(:) = 0 - pt2_d(:) = .false. - N_c = 0 - N_j = pt2_n_0(1) - do i=1,N_j - pt2_d(i) = .true. - pt2_J(i) = i - end do - - U = 0 - do while(N_j < pt2_n_tasks) - - if (N_c+ncache > N_det_generators) then - ncache = N_det_generators - N_c - endif - - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(dt,v,t,k) - do k=1, ncache - dt = pt2_u_0 - do t=1, pt2_N_teeth - v = dt + pt2_W_T *pt2_u(N_c+k) - dt = dt + pt2_W_T - ii(t,k) = pt2_find_sample_lr(v, pt2_cW,pt2_n_0(t),pt2_n_0(t+1)) - end do - enddo - !$OMP END PARALLEL DO - - do k=1,ncache - !ADD_COMB - N_c = N_c+1 - do t=1, pt2_N_teeth - i = ii(t,k) - if(.not. pt2_d(i)) then - N_j += 1 - pt2_J(N_j) = i - pt2_d(i) = .true. - end if - end do - - pt2_R(N_j) = N_c - - !FILL_TOOTH - do while(U < N_det_generators) - U += 1 - if(.not. pt2_d(U)) then - N_j += 1 - pt2_J(N_j) = U - pt2_d(U) = .true. - exit - end if - end do - if (N_j >= pt2_n_tasks) exit - end do - enddo - - if(N_det_generators > 1) then - pt2_R(N_det_generators-1) = 0 - pt2_R(N_det_generators) = N_c - end if - - deallocate(ii,pt2_d) - -END_PROVIDER - - - - BEGIN_PROVIDER [ double precision, pt2_w, (N_det_generators) ] -&BEGIN_PROVIDER [ double precision, pt2_cW, (0:N_det_generators) ] -&BEGIN_PROVIDER [ double precision, pt2_W_T ] -&BEGIN_PROVIDER [ double precision, pt2_u_0 ] -&BEGIN_PROVIDER [ integer, pt2_n_0, (pt2_N_teeth+1) ] - implicit none - integer :: i, t - double precision, allocatable :: tilde_w(:), tilde_cW(:) - double precision :: r, tooth_width - integer, external :: pt2_find_sample - - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - rss = memory_of_double(2*N_det_generators+1) - call check_mem(rss,irp_here) - - if (N_det_generators == 1) then - - pt2_w(1) = 1.d0 - pt2_cw(1) = 1.d0 - pt2_u_0 = 1.d0 - pt2_W_T = 0.d0 - pt2_n_0(1) = 0 - pt2_n_0(2) = 1 - - else - - allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators)) - - tilde_cW(0) = 0d0 - - do i=1,N_det_generators - tilde_w(i) = psi_coef_sorted_tc_gen(i,pt2_stoch_istate)**2 !+ 1.d-20 - enddo - - double precision :: norm2 - norm2 = 0.d0 - do i=N_det_generators,1,-1 - norm2 += tilde_w(i) - enddo - - tilde_w(:) = tilde_w(:) / norm2 - - tilde_cW(0) = -1.d0 - do i=1,N_det_generators - tilde_cW(i) = tilde_cW(i-1) + tilde_w(i) - enddo - tilde_cW(:) = tilde_cW(:) + 1.d0 - - pt2_n_0(1) = 0 - do - pt2_u_0 = tilde_cW(pt2_n_0(1)) - r = tilde_cW(pt2_n_0(1) + pt2_minDetInFirstTeeth) - pt2_W_T = (1d0 - pt2_u_0) / dble(pt2_N_teeth) - if(pt2_W_T >= r - pt2_u_0) then - exit - end if - pt2_n_0(1) += 1 - if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then - print *, "teeth building failed" - stop -1 - end if - end do - - do t=2, pt2_N_teeth - r = pt2_u_0 + pt2_W_T * dble(t-1) - pt2_n_0(t) = pt2_find_sample(r, tilde_cW) - end do - pt2_n_0(pt2_N_teeth+1) = N_det_generators - - pt2_w(:pt2_n_0(1)) = tilde_w(:pt2_n_0(1)) - do t=1, pt2_N_teeth - tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t)) - if (tooth_width == 0.d0) then - tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))) - endif - ASSERT(tooth_width > 0.d0) - do i=pt2_n_0(t)+1, pt2_n_0(t+1) - pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width - end do - end do - - pt2_cW(0) = 0d0 - do i=1,N_det_generators - pt2_cW(i) = pt2_cW(i-1) + pt2_w(i) - end do - pt2_n_0(pt2_N_teeth+1) = N_det_generators - - endif -END_PROVIDER - - - - - 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..2200373b 100644 --- a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -62,7 +62,7 @@ subroutine run_stochastic_cipsi ! if (N_det > N_det_max) then ! psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(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) +! psi_coef(1:N_det,1:N_states) = psi_coef_sorted_gen(1:N_det,1:N_states) ! N_det = N_det_max ! soft_touch N_det psi_det psi_coef ! if (s2_eig) then diff --git a/plugins/local/fci_tc_bi/NEED b/plugins/local/fci_tc_bi/NEED index 3bb9515a..8e9ae1c8 100644 --- a/plugins/local/fci_tc_bi/NEED +++ b/plugins/local/fci_tc_bi/NEED @@ -1,3 +1,4 @@ +generators_full_tc json tc_bi_ortho davidson_undressed diff --git a/plugins/local/fci_tc_bi/selectors.irp.f b/plugins/local/fci_tc_bi/selectors.irp.f index 7f93ae55..606660fd 100644 --- a/plugins/local/fci_tc_bi/selectors.irp.f +++ b/plugins/local/fci_tc_bi/selectors.irp.f @@ -40,7 +40,7 @@ END_PROVIDER enddo do k=1,N_states do i=1,N_det_selectors - psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k) + psi_selectors_coef(i,k) = psi_coef_sorted_gen(i,k) psi_selectors_coef_tc(i,1,k) = psi_l_coef_sorted_bi_ortho(i,k) psi_selectors_coef_tc(i,2,k) = psi_r_coef_sorted_bi_ortho(i,k) enddo diff --git a/src/cipsi/NEED b/src/cipsi/NEED index 89c128ec..ddd1e8cc 100644 --- a/src/cipsi/NEED +++ b/src/cipsi/NEED @@ -1,3 +1,4 @@ +cipsi_utils json perturbation zmq diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 3b048c14..228e0ef1 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -1,923 +1,3 @@ -BEGIN_PROVIDER [ integer, pt2_stoch_istate ] - implicit none - BEGIN_DOC - ! State for stochatsic PT2 - END_DOC - pt2_stoch_istate = 1 -END_PROVIDER - - BEGIN_PROVIDER [ integer, pt2_F, (N_det_generators) ] -&BEGIN_PROVIDER [ integer, pt2_n_tasks_max ] - implicit none - logical, external :: testTeethBuilding - integer :: i,j - pt2_n_tasks_max = elec_alpha_num*elec_alpha_num + elec_alpha_num*elec_beta_num - n_core_orb*2 - pt2_n_tasks_max = min(pt2_n_tasks_max,1+N_det_generators/10000) - call write_int(6,pt2_n_tasks_max,'pt2_n_tasks_max') - - pt2_F(:) = max(int(sqrt(float(pt2_n_tasks_max))),1) - do i=1,pt2_n_0(1+pt2_N_teeth/4) - pt2_F(i) = pt2_n_tasks_max*pt2_min_parallel_tasks - enddo - do i=1+pt2_n_0(pt2_N_teeth-pt2_N_teeth/4), pt2_n_0(pt2_N_teeth-pt2_N_teeth/10) - pt2_F(i) = pt2_min_parallel_tasks - enddo - do i=1+pt2_n_0(pt2_N_teeth-pt2_N_teeth/10), N_det_generators - pt2_F(i) = 1 - enddo - -END_PROVIDER - - BEGIN_PROVIDER [ integer, pt2_N_teeth ] -&BEGIN_PROVIDER [ integer, pt2_minDetInFirstTeeth ] - implicit none - logical, external :: testTeethBuilding - - if(N_det_generators < 1024) then - pt2_minDetInFirstTeeth = 1 - pt2_N_teeth = 1 - else - pt2_minDetInFirstTeeth = min(5, N_det_generators) - do pt2_N_teeth=100,2,-1 - if(testTeethBuilding(pt2_minDetInFirstTeeth, pt2_N_teeth)) exit - end do - end if - call write_int(6,pt2_N_teeth,'Number of comb teeth') -END_PROVIDER - - -logical function testTeethBuilding(minF, N) - implicit none - integer, intent(in) :: minF, N - integer :: n0, i - double precision :: u0, Wt, r - - double precision, allocatable :: tilde_w(:), tilde_cW(:) - integer, external :: dress_find_sample - - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - - rss = memory_of_double(2*N_det_generators+1) - call check_mem(rss,irp_here) - - allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators)) - - double precision :: norm2 - norm2 = 0.d0 - do i=N_det_generators,1,-1 - tilde_w(i) = psi_coef_sorted_gen(i,pt2_stoch_istate) * & - psi_coef_sorted_gen(i,pt2_stoch_istate) - norm2 = norm2 + tilde_w(i) - enddo - - f = 1.d0/norm2 - tilde_w(:) = tilde_w(:) * f - - tilde_cW(0) = -1.d0 - do i=1,N_det_generators - tilde_cW(i) = tilde_cW(i-1) + tilde_w(i) - enddo - tilde_cW(:) = tilde_cW(:) + 1.d0 - deallocate(tilde_w) - - n0 = 0 - testTeethBuilding = .false. - double precision :: f - integer :: minFN - minFN = N_det_generators - minF * N - f = 1.d0/dble(N) - do - u0 = tilde_cW(n0) - r = tilde_cW(n0 + minF) - Wt = (1d0 - u0) * f - if (dabs(Wt) <= 1.d-3) then - exit - endif - if(Wt >= r - u0) then - testTeethBuilding = .true. - exit - end if - n0 += 1 - if(n0 > minFN) then - exit - end if - end do - deallocate(tilde_cW) - -end function - - - -subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) - use f77_zmq - use selection_types - - implicit none - - integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull - integer, intent(in) :: N_in - double precision, intent(in) :: relative_error, E(N_states) - type(pt2_type), intent(inout) :: pt2_data, pt2_data_err -! - integer :: i, N - - double precision :: state_average_weight_save(N_states), w(N_states,4) - integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket - type(selection_buffer) :: b - - PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique - PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order - PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted - PROVIDE psi_det_hii selection_weight pseudo_sym - PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE excitation_beta_max excitation_alpha_max excitation_max - - if (h0_type == 'CFG') then - PROVIDE psi_configuration_hii det_to_configuration - endif - - if (N_det <= max(4,N_states) .or. pt2_N_teeth < 2) then - call ZMQ_selection(N_in, pt2_data) - else - - N = max(N_in,1) * N_states - state_average_weight_save(:) = state_average_weight(:) - if (int(N,8)*2_8 > huge(1)) then - print *, irp_here, ': integer too large' - stop -1 - endif - call create_selection_buffer(N, N*2, b) - ASSERT (associated(b%det)) - ASSERT (associated(b%val)) - - do pt2_stoch_istate=1,N_states - state_average_weight(:) = 0.d0 - state_average_weight(pt2_stoch_istate) = 1.d0 - 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 - call new_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') - - integer, external :: zmq_put_psi - integer, external :: zmq_put_N_det_generators - integer, external :: zmq_put_N_det_selectors - integer, external :: zmq_put_dvector - integer, external :: zmq_put_ivector - if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then - stop 'Unable to put psi on ZMQ server' - endif - if (zmq_put_N_det_generators(zmq_to_qp_run_socket, 1) == -1) then - stop 'Unable to put N_det_generators on ZMQ server' - endif - if (zmq_put_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) then - stop 'Unable to put N_det_selectors on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',pt2_e0_denominator,size(pt2_e0_denominator)) == -1) then - stop 'Unable to put energy on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then - stop 'Unable to put state_average_weight on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then - stop 'Unable to put selection_weight on ZMQ server' - endif - if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then - stop 'Unable to put pt2_stoch_istate on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',(/threshold_generators/),1) == -1) then - stop 'Unable to put threshold_generators on ZMQ server' - endif - - - integer, external :: add_task_to_taskserver - character(300000) :: task - - integer :: j,k,ipos,ifirst - ifirst=0 - - ipos=0 - do i=1,N_det_generators - if (pt2_F(i) > 1) then - ipos += 1 - endif - enddo - call write_int(6,sum(pt2_F),'Number of tasks') - call write_int(6,ipos,'Number of fragmented tasks') - - ipos=1 - do i= 1, N_det_generators - do j=1,pt2_F(pt2_J(i)) - write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') j, pt2_J(i), N_in - ipos += 30 - if (ipos > 300000-30) then - if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then - stop 'Unable to add task to task server' - endif - ipos=1 - if (ifirst == 0) then - ifirst=1 - if (zmq_set_running(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Failed in zmq_set_running' - endif - endif - endif - end do - enddo - if (ipos > 1) then - if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then - stop 'Unable to add task to task server' - endif - endif - - integer, external :: zmq_set_running - if (zmq_set_running(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Failed in zmq_set_running' - endif - - - double precision :: mem_collector, mem, rss - - call resident_memory(rss) - - mem_collector = 8.d0 * & ! bytes - ( 1.d0*pt2_n_tasks_max & ! task_id, index - + 0.635d0*N_det_generators & ! f,d - + pt2_n_tasks_max*pt2_type_size(N_states) & ! pt2_data_task - + N_det_generators*pt2_type_size(N_states) & ! pt2_data_I - + 4.d0*(pt2_N_teeth+1) & ! S, S2, T2, T3 - + 1.d0*(N_int*2.d0*N + N) & ! selection buffer - + 1.d0*(N_int*2.d0*N + N) & ! sort selection buffer - ) / 1024.d0**3 - - integer :: nproc_target, ii - nproc_target = nthreads_pt2 - ii = min(N_det, (elec_alpha_num*(mo_num-elec_alpha_num))**2) - - do - mem = mem_collector + & ! - nproc_target * 8.d0 * & ! bytes - ( 0.5d0*pt2_n_tasks_max & ! task_id - + 64.d0*pt2_n_tasks_max & ! task - + pt2_type_size(N_states)*pt2_n_tasks_max*N_states & ! pt2, variance, overlap - + 1.d0*pt2_n_tasks_max & ! i_generator, subset - + 1.d0*(N_int*2.d0*ii+ ii) & ! selection buffer - + 1.d0*(N_int*2.d0*ii+ ii) & ! sort selection buffer - + 2.0d0*(ii) & ! preinteresting, interesting, - ! prefullinteresting, fullinteresting - + 2.0d0*(N_int*2*ii) & ! minilist, fullminilist - + 1.0d0*(N_states*mo_num*mo_num) & ! mat - ) / 1024.d0**3 - - if (nproc_target == 0) then - call check_mem(mem,irp_here) - nproc_target = 1 - exit - endif - - if (mem+rss < qp_max_mem) then - exit - endif - - nproc_target = nproc_target - 1 - - enddo - call write_int(6,nproc_target,'Number of threads for PT2') - call write_double(6,mem,'Memory (Gb)') - - call set_multiple_levels_omp(.False.) - - - print '(A)', '========== ==================== ================ ================ ================ ============= ===========' - print '(A)', ' Samples Energy PT2 Variance Norm^2 Convergence Seconds' - print '(A)', '========== ==================== ================ ================ ================ ============= ===========' - - PROVIDE global_selection_buffer - - !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc_target+1) & - !$OMP PRIVATE(i) - i = omp_get_thread_num() - if (i==0) then - - call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, pt2_data, pt2_data_err, b, N) - pt2_data % rpt2(pt2_stoch_istate) = & - pt2_data % pt2(pt2_stoch_istate)/(1.d0+pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate)) - - !TODO : We should use here the correct formula for the error of X/Y - pt2_data_err % rpt2(pt2_stoch_istate) = & - pt2_data_err % pt2(pt2_stoch_istate)/(1.d0 + pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate)) - - else - call pt2_slave_inproc(i) - endif - !$OMP END PARALLEL - call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') - call set_multiple_levels_omp(.True.) - - print '(A)', '========== ==================== ================ ================ ================ ============= ===========' - - - do k=1,N_states - pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) - enddo - SOFT_TOUCH pt2_overlap - - enddo - FREE pt2_stoch_istate - - ! Symmetrize overlap - do j=2,N_states - do i=1,j-1 - pt2_overlap(i,j) = 0.5d0 * (pt2_overlap(i,j) + pt2_overlap(j,i)) - pt2_overlap(j,i) = pt2_overlap(i,j) - enddo - enddo - - print *, 'Overlap of perturbed states:' - do k=1,N_states - print *, pt2_overlap(k,:) - enddo - print *, '-------' - - if (N_in > 0) then - b%cur = min(N_in,b%cur) - if (s2_eig) then - call make_selection_buffer_s2(b) - else - call remove_duplicates_in_selection_buffer(b) - endif - call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) - endif - call delete_selection_buffer(b) - - state_average_weight(:) = state_average_weight_save(:) - TOUCH state_average_weight - call update_pt2_and_variance_weights(pt2_data, N_states) - endif - - -end subroutine - - -subroutine pt2_slave_inproc(i) - implicit none - integer, intent(in) :: i - - PROVIDE global_selection_buffer - call run_pt2_slave(1,i,pt2_e0_denominator) +subroutine provide_for_zmq_pt2 + PROVIDE psi_selectors_coef_transp psi_det_sorted psi_det_sorted_order end - - -subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_err, b, N_) - use f77_zmq - use selection_types - use bitmasks - implicit none - - - integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - double precision, intent(in) :: relative_error, E - type(pt2_type), intent(inout) :: pt2_data, pt2_data_err - type(selection_buffer), intent(inout) :: b - integer, intent(in) :: N_ - - type(pt2_type), allocatable :: pt2_data_task(:) - type(pt2_type), allocatable :: pt2_data_I(:) - type(pt2_type), allocatable :: pt2_data_S(:) - type(pt2_type), allocatable :: pt2_data_S2(:) - type(pt2_type) :: pt2_data_teeth - integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - integer, external :: zmq_delete_tasks_async_send - integer, external :: zmq_delete_tasks_async_recv - integer, external :: zmq_abort - integer, external :: pt2_find_sample_lr - - PROVIDE pt2_stoch_istate - - integer :: more, n, i, p, c, t, n_tasks, U - integer, allocatable :: task_id(:) - integer, allocatable :: index(:) - - double precision :: v, x, x2, x3, avg, avg2, avg3(N_states), eqt, E0, v0, n0(N_states) - double precision :: eqta(N_states) - double precision :: time, time1, time0 - - integer, allocatable :: f(:) - logical, allocatable :: d(:) - logical :: do_exit, stop_now, sending - logical, external :: qp_stop - type(selection_buffer) :: b2 - - - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - - character(len=20) :: format_str1, str_error1, format_str2, str_error2 - character(len=20) :: format_str3, str_error3, format_str4, str_error4 - character(len=20) :: format_value1, format_value2, format_value3, format_value4 - character(len=20) :: str_value1, str_value2, str_value3, str_value4 - character(len=20) :: str_conv - double precision :: value1, value2, value3, value4 - double precision :: error1, error2, error3, error4 - integer :: size1,size2,size3,size4 - - double precision :: conv_crit - - sending =.False. - - rss = memory_of_int(pt2_n_tasks_max*2+N_det_generators*2) - rss += memory_of_double(N_states*N_det_generators)*3.d0 - rss += memory_of_double(N_states*pt2_n_tasks_max)*3.d0 - rss += memory_of_double(pt2_N_teeth+1)*4.d0 - call check_mem(rss,irp_here) - - ! If an allocation is added here, the estimate of the memory should also be - ! updated in ZMQ_pt2 - allocate(task_id(pt2_n_tasks_max), index(pt2_n_tasks_max), f(N_det_generators)) - allocate(d(N_det_generators+1)) - allocate(pt2_data_task(pt2_n_tasks_max)) - allocate(pt2_data_I(N_det_generators)) - allocate(pt2_data_S(pt2_N_teeth+1)) - allocate(pt2_data_S2(pt2_N_teeth+1)) - - - - zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() - call create_selection_buffer(N_, N_*2, b2) - - - pt2_data % pt2(pt2_stoch_istate) = -huge(1.) - pt2_data_err % pt2(pt2_stoch_istate) = huge(1.) - pt2_data % variance(pt2_stoch_istate) = huge(1.) - pt2_data_err % variance(pt2_stoch_istate) = huge(1.) - pt2_data % overlap(:,pt2_stoch_istate) = 0.d0 - pt2_data_err % overlap(:,pt2_stoch_istate) = huge(1.) - n = 1 - t = 0 - U = 0 - do i=1,pt2_n_tasks_max - call pt2_alloc(pt2_data_task(i),N_states) - enddo - do i=1,pt2_N_teeth+1 - call pt2_alloc(pt2_data_S(i),N_states) - call pt2_alloc(pt2_data_S2(i),N_states) - enddo - do i=1,N_det_generators - call pt2_alloc(pt2_data_I(i),N_states) - enddo - f(:) = pt2_F(:) - d(:) = .false. - n_tasks = 0 - E0 = E - v0 = 0.d0 - n0(:) = 0.d0 - more = 1 - call wall_time(time0) - time1 = time0 - - do_exit = .false. - stop_now = .false. - do while (n <= N_det_generators) - if(f(pt2_J(n)) == 0) then - d(pt2_J(n)) = .true. - do while(d(U+1)) - U += 1 - end do - - ! Deterministic part - do while(t <= pt2_N_teeth) - if(U >= pt2_n_0(t+1)) then - t=t+1 - E0 = 0.d0 - v0 = 0.d0 - n0(:) = 0.d0 - do i=pt2_n_0(t),1,-1 - E0 += pt2_data_I(i) % pt2(pt2_stoch_istate) - v0 += pt2_data_I(i) % variance(pt2_stoch_istate) - n0(:) += pt2_data_I(i) % overlap(:,pt2_stoch_istate) - end do - else - exit - end if - end do - - ! Add Stochastic part - c = pt2_R(n) - if(c > 0) then - - call pt2_alloc(pt2_data_teeth,N_states) - do p=pt2_N_teeth, 1, -1 - v = pt2_u_0 + pt2_W_T * (pt2_u(c) + dble(p-1)) - i = pt2_find_sample_lr(v, pt2_cW,pt2_n_0(p),pt2_n_0(p+1)) - v = pt2_W_T / pt2_w(i) - call pt2_add ( pt2_data_teeth, v, pt2_data_I(i) ) - call pt2_add ( pt2_data_S(p), 1.d0, pt2_data_teeth ) - call pt2_add2( pt2_data_S2(p), 1.d0, pt2_data_teeth ) - enddo - call pt2_dealloc(pt2_data_teeth) - - avg = E0 + pt2_data_S(t) % pt2(pt2_stoch_istate) / dble(c) - avg2 = v0 + pt2_data_S(t) % variance(pt2_stoch_istate) / dble(c) - avg3(:) = n0(:) + pt2_data_S(t) % overlap(:,pt2_stoch_istate) / dble(c) - if ((avg /= 0.d0) .or. (n == N_det_generators) ) then - do_exit = .true. - endif - if (qp_stop()) then - stop_now = .True. - endif - pt2_data % pt2(pt2_stoch_istate) = avg - pt2_data % variance(pt2_stoch_istate) = avg2 - pt2_data % overlap(:,pt2_stoch_istate) = avg3(:) - call wall_time(time) - ! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969) - if(c > 2) then - eqt = dabs((pt2_data_S2(t) % pt2(pt2_stoch_istate) / c) - (pt2_data_S(t) % pt2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = dsqrt(eqt / (dble(c) - 1.5d0)) - pt2_data_err % pt2(pt2_stoch_istate) = eqt - - eqt = dabs((pt2_data_S2(t) % variance(pt2_stoch_istate) / c) - (pt2_data_S(t) % variance(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = dsqrt(eqt / (dble(c) - 1.5d0)) - pt2_data_err % variance(pt2_stoch_istate) = eqt - - eqta(:) = dabs((pt2_data_S2(t) % overlap(:,pt2_stoch_istate) / c) - (pt2_data_S(t) % overlap(:,pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqta(:) = dsqrt(eqta(:) / (dble(c) - 1.5d0)) - pt2_data_err % overlap(:,pt2_stoch_istate) = eqta(:) - - - if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then - time1 = time - - value1 = pt2_data % pt2(pt2_stoch_istate) + E - error1 = pt2_data_err % pt2(pt2_stoch_istate) - value2 = pt2_data % pt2(pt2_stoch_istate) - error2 = pt2_data_err % pt2(pt2_stoch_istate) - value3 = pt2_data % variance(pt2_stoch_istate) - error3 = pt2_data_err % variance(pt2_stoch_istate) - value4 = pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate) - error4 = pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate) - - ! Max size of the values (FX.Y) with X=size - size1 = 15 - size2 = 12 - size3 = 12 - size4 = 12 - - ! To generate the format: number(error) - call format_w_error(value1,error1,size1,8,format_value1,str_error1) - call format_w_error(value2,error2,size2,8,format_value2,str_error2) - call format_w_error(value3,error3,size3,8,format_value3,str_error3) - call format_w_error(value4,error4,size4,8,format_value4,str_error4) - - ! value > string with the right format - write(str_value1,'('//format_value1//')') value1 - write(str_value2,'('//format_value2//')') value2 - write(str_value3,'('//format_value3//')') value3 - write(str_value4,'('//format_value4//')') value4 - - ! Convergence criterion - conv_crit = dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & - (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) - write(str_conv,'(G10.3)') conv_crit - - write(*,'(I10,X,X,A20,X,A16,X,A16,X,A16,X,A12,X,F10.1)') c,& - adjustl(adjustr(str_value1)//'('//str_error1(1:1)//')'),& - adjustl(adjustr(str_value2)//'('//str_error2(1:1)//')'),& - adjustl(adjustr(str_value3)//'('//str_error3(1:1)//')'),& - adjustl(adjustr(str_value4)//'('//str_error4(1:1)//')'),& - adjustl(str_conv),& - time-time0 - - ! Old print - !print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1,ES16.6,ES16.6)', c, & - ! pt2_data % pt2(pt2_stoch_istate) +E, & - ! pt2_data_err % pt2(pt2_stoch_istate), & - ! pt2_data % variance(pt2_stoch_istate), & - ! pt2_data_err % variance(pt2_stoch_istate), & - ! pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate), & - ! pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), & - ! time-time0, & - ! pt2_data % pt2(pt2_stoch_istate), & - ! dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & - ! (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) - - if (stop_now .or. ( & - (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & - (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then - if (zmq_abort(zmq_to_qp_run_socket) == -1) then - call sleep(10) - if (zmq_abort(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Error in sending abort signal (2)' - endif - endif - endif - endif - endif - end if - n += 1 - else if(more == 0) then - exit - else - call pull_pt2_results(zmq_socket_pull, index, pt2_data_task, task_id, n_tasks, b2) - if(n_tasks > pt2_n_tasks_max)then - print*,'PB !!!' - print*,'If you see this, send a bug report with the following content' - print*,irp_here - print*,'n_tasks,pt2_n_tasks_max = ',n_tasks,pt2_n_tasks_max - stop -1 - endif - if (zmq_delete_tasks_async_send(zmq_to_qp_run_socket,task_id,n_tasks,sending) == -1) then - stop 'PT2: Unable to delete tasks (send)' - endif - do i=1,n_tasks - if(index(i).gt.size(pt2_data_I,1).or.index(i).lt.1)then - print*,'PB !!!' - print*,'If you see this, send a bug report with the following content' - print*,irp_here - print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1) - stop -1 - endif - call pt2_add(pt2_data_I(index(i)),1.d0,pt2_data_task(i)) - f(index(i)) -= 1 - end do - do i=1, b2%cur - ! We assume the pulled buffer is sorted - if (b2%val(i) > b%mini) exit - call add_to_selection_buffer(b, b2%det(1,1,i), b2%val(i)) - end do - if (zmq_delete_tasks_async_recv(zmq_to_qp_run_socket,more,sending) == -1) then - stop 'PT2: Unable to delete tasks (recv)' - endif - end if - end do - do i=1,N_det_generators - call pt2_dealloc(pt2_data_I(i)) - enddo - do i=1,pt2_N_teeth+1 - call pt2_dealloc(pt2_data_S(i)) - call pt2_dealloc(pt2_data_S2(i)) - enddo - do i=1,pt2_n_tasks_max - call pt2_dealloc(pt2_data_task(i)) - enddo -!print *, 'deleting b2' - call delete_selection_buffer(b2) -!print *, 'sorting b' - call sort_selection_buffer(b) -!print *, 'done' - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - -end subroutine - - -integer function pt2_find_sample(v, w) - implicit none - double precision, intent(in) :: v, w(0:N_det_generators) - integer, external :: pt2_find_sample_lr - - pt2_find_sample = pt2_find_sample_lr(v, w, 0, N_det_generators) -end function - - -integer function pt2_find_sample_lr(v, w, l_in, r_in) - implicit none - double precision, intent(in) :: v, w(0:N_det_generators) - integer, intent(in) :: l_in,r_in - integer :: i,l,r - - l=l_in - r=r_in - - do while(r-l > 1) - i = shiftr(r+l,1) - if(w(i) < v) then - l = i - else - r = i - end if - end do - i = r - do r=i+1,N_det_generators - if (w(r) /= w(i)) then - exit - endif - enddo - pt2_find_sample_lr = r-1 -end function - - -BEGIN_PROVIDER [ integer, pt2_n_tasks ] - implicit none - BEGIN_DOC - ! Number of parallel tasks for the Monte Carlo - END_DOC - pt2_n_tasks = N_det_generators -END_PROVIDER - -BEGIN_PROVIDER[ double precision, pt2_u, (N_det_generators)] - implicit none - integer, allocatable :: seed(:) - integer :: m,i - call random_seed(size=m) - allocate(seed(m)) - do i=1,m - seed(i) = i - enddo - call random_seed(put=seed) - deallocate(seed) - - call RANDOM_NUMBER(pt2_u) - END_PROVIDER - - BEGIN_PROVIDER[ integer, pt2_J, (N_det_generators)] -&BEGIN_PROVIDER[ integer, pt2_R, (N_det_generators)] - implicit none - BEGIN_DOC -! pt2_J contains the list of generators after ordering them according to the -! Monte Carlo sampling. -! -! pt2_R(i) is the number of combs drawn when determinant i is computed. - END_DOC - integer :: N_c, N_j - integer :: U, t, i - double precision :: v - integer, external :: pt2_find_sample_lr - - logical, allocatable :: pt2_d(:) - integer :: m,l,r,k - integer :: ncache - integer, allocatable :: ii(:,:) - double precision :: dt - - ncache = min(N_det_generators,10000) - - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - rss = memory_of_int(ncache)*dble(pt2_N_teeth) + memory_of_int(N_det_generators) - call check_mem(rss,irp_here) - - allocate(ii(pt2_N_teeth,ncache),pt2_d(N_det_generators)) - - pt2_R(:) = 0 - pt2_d(:) = .false. - N_c = 0 - N_j = pt2_n_0(1) - do i=1,N_j - pt2_d(i) = .true. - pt2_J(i) = i - end do - - U = 0 - do while(N_j < pt2_n_tasks) - - if (N_c+ncache > N_det_generators) then - ncache = N_det_generators - N_c - endif - - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(dt,v,t,k) - do k=1, ncache - dt = pt2_u_0 - do t=1, pt2_N_teeth - v = dt + pt2_W_T *pt2_u(N_c+k) - dt = dt + pt2_W_T - ii(t,k) = pt2_find_sample_lr(v, pt2_cW,pt2_n_0(t),pt2_n_0(t+1)) - end do - enddo - !$OMP END PARALLEL DO - - do k=1,ncache - !ADD_COMB - N_c = N_c+1 - do t=1, pt2_N_teeth - i = ii(t,k) - if(.not. pt2_d(i)) then - N_j += 1 - pt2_J(N_j) = i - pt2_d(i) = .true. - end if - end do - - pt2_R(N_j) = N_c - - !FILL_TOOTH - do while(U < N_det_generators) - U += 1 - if(.not. pt2_d(U)) then - N_j += 1 - pt2_J(N_j) = U - pt2_d(U) = .true. - exit - end if - end do - if (N_j >= pt2_n_tasks) exit - end do - enddo - - if(N_det_generators > 1) then - pt2_R(N_det_generators-1) = 0 - pt2_R(N_det_generators) = N_c - end if - - deallocate(ii,pt2_d) - -END_PROVIDER - - - - BEGIN_PROVIDER [ double precision, pt2_w, (N_det_generators) ] -&BEGIN_PROVIDER [ double precision, pt2_cW, (0:N_det_generators) ] -&BEGIN_PROVIDER [ double precision, pt2_W_T ] -&BEGIN_PROVIDER [ double precision, pt2_u_0 ] -&BEGIN_PROVIDER [ integer, pt2_n_0, (pt2_N_teeth+1) ] - implicit none - integer :: i, t - double precision, allocatable :: tilde_w(:), tilde_cW(:) - double precision :: r, tooth_width - integer, external :: pt2_find_sample - - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - rss = memory_of_double(2*N_det_generators+1) - call check_mem(rss,irp_here) - - if (N_det_generators == 1) then - - pt2_w(1) = 1.d0 - pt2_cw(1) = 1.d0 - pt2_u_0 = 1.d0 - pt2_W_T = 0.d0 - pt2_n_0(1) = 0 - pt2_n_0(2) = 1 - - else - - allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators)) - - tilde_cW(0) = 0d0 - - do i=1,N_det_generators - tilde_w(i) = psi_coef_sorted_gen(i,pt2_stoch_istate)**2 !+ 1.d-20 - enddo - - double precision :: norm2 - norm2 = 0.d0 - do i=N_det_generators,1,-1 - norm2 += tilde_w(i) - enddo - - tilde_w(:) = tilde_w(:) / norm2 - - tilde_cW(0) = -1.d0 - do i=1,N_det_generators - tilde_cW(i) = tilde_cW(i-1) + tilde_w(i) - enddo - tilde_cW(:) = tilde_cW(:) + 1.d0 - - pt2_n_0(1) = 0 - do - pt2_u_0 = tilde_cW(pt2_n_0(1)) - r = tilde_cW(pt2_n_0(1) + pt2_minDetInFirstTeeth) - pt2_W_T = (1d0 - pt2_u_0) / dble(pt2_N_teeth) - if(pt2_W_T >= r - pt2_u_0) then - exit - end if - pt2_n_0(1) += 1 - if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then - print *, "teeth building failed" - stop -1 - end if - end do - - do t=2, pt2_N_teeth - r = pt2_u_0 + pt2_W_T * dble(t-1) - pt2_n_0(t) = pt2_find_sample(r, tilde_cW) - end do - pt2_n_0(pt2_N_teeth+1) = N_det_generators - - pt2_w(:pt2_n_0(1)) = tilde_w(:pt2_n_0(1)) - do t=1, pt2_N_teeth - tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t)) - if (tooth_width == 0.d0) then - tooth_width = max(1.d-15,sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1)))) - endif - ASSERT(tooth_width > 0.d0) - do i=pt2_n_0(t)+1, pt2_n_0(t+1) - pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width - end do - end do - - pt2_cW(0) = 0d0 - do i=1,N_det_generators - pt2_cW(i) = pt2_cW(i-1) + pt2_w(i) - end do - pt2_n_0(pt2_N_teeth+1) = N_det_generators - - endif -END_PROVIDER - - - - - diff --git a/src/generators_full_tc/README.rst b/src/generators_full_tc/README.rst new file mode 100644 index 00000000..4e59ee3b --- /dev/null +++ b/src/generators_full_tc/README.rst @@ -0,0 +1,9 @@ +=============== +generators_full +=============== + +Module defining the generator determinants as all the determinants of the +variational space. + +This module is intended to be included in the :file:`NEED` file to define +a full set of generators. diff --git a/plugins/local/fci_tc_bi/generators.irp.f b/src/generators_full_tc/generators.irp.f similarity index 51% rename from plugins/local/fci_tc_bi/generators.irp.f rename to src/generators_full_tc/generators.irp.f index bf972423..a9da7dbc 100644 --- a/plugins/local/fci_tc_bi/generators.irp.f +++ b/src/generators_full_tc/generators.irp.f @@ -34,23 +34,49 @@ END_PROVIDER END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_gen, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc_gen, (psi_det_size,N_states) ] -&BEGIN_PROVIDER [ integer, psi_det_sorted_tc_gen_order, (psi_det_size) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen, (psi_det_size,N_states) ] +&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_order, (psi_det_size) ] implicit none BEGIN_DOC ! For Single reference wave functions, the generator is the ! Hartree-Fock determinant END_DOC - psi_det_sorted_tc_gen = psi_det_sorted_tc - psi_coef_sorted_tc_gen = psi_coef_sorted_tc - psi_det_sorted_tc_gen_order = psi_det_sorted_tc_order - integer :: i -! do i = 1,N_det -! print*,'i = ',i -! call debug_det(psi_det_sorted_tc(1,1,i),N_int) -! enddo + psi_det_sorted_gen = psi_det_sorted_tc + psi_coef_sorted_gen = psi_coef_sorted_tc + psi_det_sorted_gen_order = psi_det_sorted_tc_order END_PROVIDER +BEGIN_PROVIDER [integer, degree_max_generators] + implicit none + BEGIN_DOC +! Max degree of excitation (respect to HF) of the generators + END_DOC + integer :: i,degree + degree_max_generators = 0 + do i = 1, N_det_generators + call get_excitation_degree(HF_bitmask,psi_det_generators(1,1,i),degree,N_int) + if(degree .gt. degree_max_generators)then + degree_max_generators = degree + endif + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer, size_select_max] + implicit none + BEGIN_DOC + ! Size of the select_max array + END_DOC + size_select_max = 10000 +END_PROVIDER + +BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ] + implicit none + BEGIN_DOC + ! Memo to skip useless selectors + END_DOC + select_max = huge(1.d0) +END_PROVIDER + From 9a15fecd6a164375ead8684c6a836147b548f4fa Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 12 Mar 2024 16:42:08 +0100 Subject: [PATCH 04/12] Merging CIPSI and TC-CIPSI --- .../cipsi_tc_bi_ortho/zmq_selection.irp.f | 234 ------------------ .../zmq_selection.irp.f | 0 2 files changed, 234 deletions(-) delete mode 100644 plugins/local/cipsi_tc_bi_ortho/zmq_selection.irp.f rename src/{cipsi => cipsi_utils}/zmq_selection.irp.f (100%) diff --git a/plugins/local/cipsi_tc_bi_ortho/zmq_selection.irp.f b/plugins/local/cipsi_tc_bi_ortho/zmq_selection.irp.f deleted file mode 100644 index 22db643f..00000000 --- a/plugins/local/cipsi_tc_bi_ortho/zmq_selection.irp.f +++ /dev/null @@ -1,234 +0,0 @@ -subroutine ZMQ_selection(N_in, pt2_data) - use f77_zmq - use selection_types - - implicit none - - integer(ZMQ_PTR) :: zmq_to_qp_run_socket , zmq_socket_pull - integer, intent(in) :: N_in - type(selection_buffer) :: b - integer :: i, l, N - integer, external :: omp_get_thread_num - type(pt2_type), intent(inout) :: pt2_data - -! PROVIDE psi_det psi_coef N_det qp_max_mem N_states pt2_F s2_eig N_det_generators - - N = max(N_in,1) - N = min(N, (elec_alpha_num * (mo_num-elec_alpha_num))**2) - if (.True.) then - PROVIDE pt2_e0_denominator nproc - PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique - PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order - PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym - PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE excitation_beta_max excitation_alpha_max excitation_max - - call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection') - - integer, external :: zmq_put_psi - integer, external :: zmq_put_N_det_generators - integer, external :: zmq_put_N_det_selectors - integer, external :: zmq_put_dvector - - if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then - stop 'Unable to put psi on ZMQ server' - endif - if (zmq_put_N_det_generators(zmq_to_qp_run_socket, 1) == -1) then - stop 'Unable to put N_det_generators on ZMQ server' - endif - if (zmq_put_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) then - stop 'Unable to put N_det_selectors on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',pt2_e0_denominator,size(pt2_e0_denominator)) == -1) then - stop 'Unable to put energy on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then - stop 'Unable to put state_average_weight on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then - stop 'Unable to put selection_weight on ZMQ server' - endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',(/threshold_generators/),1) == -1) then - stop 'Unable to put threshold_generators on ZMQ server' - endif - call create_selection_buffer(N, N*2, b) - endif - - integer, external :: add_task_to_taskserver - character(len=100000) :: task - integer :: j,k,ipos - ipos=1 - task = ' ' - - do i= 1, N_det_generators - do j=1,pt2_F(i) - write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') j, i, N - ipos += 30 - if (ipos > 100000-30) then - if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then - stop 'Unable to add task to task server' - endif - ipos=1 - endif - end do - enddo - if (ipos > 1) then - if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then - stop 'Unable to add task to task server' - endif - endif - N = max(N_in,1) - - - ASSERT (associated(b%det)) - ASSERT (associated(b%val)) - - integer, external :: zmq_set_running - if (zmq_set_running(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Failed in zmq_set_running' - endif - - integer :: nproc_target - if (N_det < 3*nproc) then - nproc_target = N_det/4 - else - nproc_target = nproc - endif - double precision :: mem - mem = 8.d0 * N_det * (N_int * 2.d0 * 3.d0 + 3.d0 + 5.d0) / (1024.d0**3) - call write_double(6,mem,'Estimated memory/thread (Gb)') - if (qp_max_mem > 0) then - nproc_target = max(1,int(dble(qp_max_mem)/(0.1d0 + mem))) - nproc_target = min(nproc_target,nproc) - endif - - f(:) = 1.d0 - if (.not.do_pt2) then - double precision :: f(N_states), u_dot_u - do k=1,min(N_det,N_states) - f(k) = 1.d0 / u_dot_u(psi_selectors_coef(1,k), N_det_selectors) - enddo - endif - - !$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2_data) PRIVATE(i) NUM_THREADS(nproc_target+1) - i = omp_get_thread_num() - if (i==0) then - call selection_collector(zmq_socket_pull, b, N, pt2_data) - else - call selection_slave_inproc(i) - endif - !$OMP END PARALLEL - - call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'selection') - if (N_in > 0) then - if (s2_eig) then - call make_selection_buffer_s2(b) - endif - call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) - endif - call delete_selection_buffer(b) - - do k=1,N_states - pt2_data % pt2(k) = pt2_data % pt2(k) * f(k) - pt2_data % variance(k) = pt2_data % variance(k) * f(k) - do l=1,N_states - pt2_data % overlap(k,l) = pt2_data % overlap(k,l) * dsqrt(f(k)*f(l)) - pt2_data % overlap(l,k) = pt2_data % overlap(l,k) * dsqrt(f(k)*f(l)) - enddo - - pt2_data % rpt2(k) = & - pt2_data % pt2(k)/(1.d0 + pt2_data % overlap(k,k)) - enddo - - pt2_overlap(:,:) = pt2_data % overlap(:,:) - - print *, 'Overlap of perturbed states:' - do l=1,N_states - print *, pt2_overlap(l,:) - enddo - print *, '-------' - SOFT_TOUCH pt2_overlap - call update_pt2_and_variance_weights(pt2_data, N_states) - -end subroutine - - -subroutine selection_slave_inproc(i) - implicit none - integer, intent(in) :: i - - call run_selection_slave(1,i,pt2_e0_denominator) -end - -subroutine selection_collector(zmq_socket_pull, b, N, pt2_data) - use f77_zmq - use selection_types - use bitmasks - implicit none - - - integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - type(selection_buffer), intent(inout) :: b - integer, intent(in) :: N - type(pt2_type), intent(inout) :: pt2_data - type(pt2_type) :: pt2_data_tmp - - double precision :: pt2_mwen(N_states) - double precision :: variance_mwen(N_states) - double precision :: norm2_mwen(N_states) - integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - - integer(ZMQ_PTR), external :: new_zmq_pull_socket - - integer :: msg_size, rc, more - integer :: acc, i, j, robin, ntask - double precision, pointer :: val(:) - integer(bit_kind), pointer :: det(:,:,:) - integer, allocatable :: task_id(:) - type(selection_buffer) :: b2 - - - - - zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() - call create_selection_buffer(N, N*2, b2) - integer :: k - double precision :: rss - double precision, external :: memory_of_int - rss = memory_of_int(N_det_generators) - call check_mem(rss,irp_here) - allocate(task_id(N_det_generators)) - more = 1 - pt2_data % pt2(:) = 0d0 - pt2_data % variance(:) = 0.d0 - pt2_data % overlap(:,:) = 0.d0 - call pt2_alloc(pt2_data_tmp,N_states) - do while (more == 1) - call pull_selection_results(zmq_socket_pull, pt2_data_tmp, b2%val(1), b2%det(1,1,1), b2%cur, task_id, ntask) - - call pt2_add(pt2_data, 1.d0, pt2_data_tmp) - do i=1, b2%cur - call add_to_selection_buffer(b, b2%det(1,1,i), b2%val(i)) - if (b2%val(i) > b%mini) exit - end do - - do i=1, ntask - if(task_id(i) == 0) then - print *, "Error in collector" - endif - integer, external :: zmq_delete_task - if (zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more) == -1) then - stop 'Unable to delete task' - endif - end do - end do - call pt2_dealloc(pt2_data_tmp) - - - call delete_selection_buffer(b2) - call sort_selection_buffer(b) - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) -end subroutine - diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi_utils/zmq_selection.irp.f similarity index 100% rename from src/cipsi/zmq_selection.irp.f rename to src/cipsi_utils/zmq_selection.irp.f From 1769efddca34f996d0d3a169bea08030c6a9e0ed Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 12 Mar 2024 16:52:53 +0100 Subject: [PATCH 05/12] fixed the qp_test of tc_scf --- plugins/local/tc_scf/11.tc_scf.bats | 44 ++++++++++++++++------------- 1 file changed, 24 insertions(+), 20 deletions(-) diff --git a/plugins/local/tc_scf/11.tc_scf.bats b/plugins/local/tc_scf/11.tc_scf.bats index b81c2f4b..f5f2e3c1 100644 --- a/plugins/local/tc_scf/11.tc_scf.bats +++ b/plugins/local/tc_scf/11.tc_scf.bats @@ -10,16 +10,17 @@ function run_Ne() { qp create_ezfio -b cc-pcvdz Ne.xyz -o Ne_tc_scf qp run scf + qp set tc_keywords tc_integ_type numeric + qp set jastrow env_type Sum_Gauss qp set hamiltonian mu_erf 0.87 - qp set tc_keywords j1b_type 3 - qp set tc_keywords j1b_pen [1.5] - qp set tc_keywords bi_ortho True - qp set tc_keywords test_cycle_tc True + qp set jastrow j1e_type None + qp set jastrow env_coef "[1.]" + qp set jastrow env_expo "[1.5]" qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out eref=-128.552134 energy="$(qp get tc_scf bitc_energy)" - eq $energy $eref 1e-6 + eq $energy $eref 2e-4 } @@ -33,16 +34,17 @@ function run_C() { qp create_ezfio -b cc-pcvdz C.xyz -o C_tc_scf -m 3 qp run scf + qp set tc_keywords tc_integ_type numeric + qp set jastrow env_type Sum_Gauss qp set hamiltonian mu_erf 0.87 - qp set tc_keywords j1b_type 3 - qp set tc_keywords j1b_pen [1.5] - qp set tc_keywords bi_ortho True - qp set tc_keywords test_cycle_tc True + qp set jastrow j1e_type None + qp set jastrow env_coef "[1.]" + qp set jastrow env_expo "[1.5]" qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out eref=-37.691254356408791 energy="$(qp get tc_scf bitc_energy)" - eq $energy $eref 1e-6 + eq $energy $eref 2e-4 } @@ -57,16 +59,17 @@ function run_O() { qp create_ezfio -b cc-pcvdz O.xyz -o O_tc_scf -m 3 qp run scf + qp set tc_keywords tc_integ_type numeric + qp set jastrow env_type Sum_Gauss + qp set jastrow j1e_type None + qp set jastrow env_coef "[1.]" + qp set jastrow env_expo "[1.5]" qp set hamiltonian mu_erf 0.87 - qp set tc_keywords j1b_type 3 - qp set tc_keywords j1b_pen [1.5] - qp set tc_keywords bi_ortho True - qp set tc_keywords test_cycle_tc True qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out eref=-74.814687229354590 energy="$(qp get tc_scf bitc_energy)" - eq $energy $eref 1e-6 + eq $energy $eref 2e-4 } @@ -82,16 +85,17 @@ function run_ch2() { qp create_ezfio -b "C:cc-pcvdz|H:cc-pvdz" ch2.xyz -o ch2_tc_scf qp run scf + qp set tc_keywords tc_integ_type numeric + qp set jastrow env_type Sum_Gauss + qp set jastrow j1e_type None + qp set jastrow env_coef "[1., 1., 1.]" + qp set jastrow env_expo '[1.5,10000,10000]' qp set hamiltonian mu_erf 0.87 - qp set tc_keywords j1b_type 3 - qp set tc_keywords j1b_pen '[1.5,10000,10000]' - qp set tc_keywords bi_ortho True - qp set tc_keywords test_cycle_tc True qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out eref=-38.903247818077737 energy="$(qp get tc_scf bitc_energy)" - eq $energy $eref 1e-6 + eq $energy $eref 2e-4 } From a42c79ca34111ac449bdf2b18243ef38f9d4abe6 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 12 Mar 2024 17:09:58 +0100 Subject: [PATCH 06/12] The test works for fci_tc_bi but not for tc_bi_ortho --- plugins/local/tc_bi_ortho/31.tc_bi_ortho.bats | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/plugins/local/tc_bi_ortho/31.tc_bi_ortho.bats b/plugins/local/tc_bi_ortho/31.tc_bi_ortho.bats index 93bed2ab..33afcb92 100644 --- a/plugins/local/tc_bi_ortho/31.tc_bi_ortho.bats +++ b/plugins/local/tc_bi_ortho/31.tc_bi_ortho.bats @@ -14,7 +14,7 @@ function run_Ne() { qp run tc_bi_ortho | tee Ne_tc_scf.cisd_tc_bi_ortho.out eref=-128.77020441279302 energy=$(get_e Ne_tc_scf.cisd_tc_bi_ortho.out) - eq $energy $eref 1e-6 + eq $energy $eref 2e-4 } @@ -29,7 +29,7 @@ function run_C() { qp run tc_bi_ortho | tee C_tc_scf.cisd_tc_bi_ortho.out eref=-37.757536149952514 energy=$(get_e C_tc_scf.cisd_tc_bi_ortho.out) - eq $energy $eref 1e-6 + eq $energy $eref 2e-4 } @@ -43,7 +43,7 @@ function run_O() { qp run tc_bi_ortho | tee O_tc_scf.cisd_tc_bi_ortho.out eref=-74.908518517716161 energy=$(get_e O_tc_scf.cisd_tc_bi_ortho.out) - eq $energy $eref 1e-6 + eq $energy $eref 2e-4 } From f816773102c06547c1f8d3a5f5b492321b4fd84f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 12 Mar 2024 17:21:35 +0100 Subject: [PATCH 07/12] Refactor CIPSI / TC-CIPSI --- plugins/local/cipsi_tc_bi_ortho/energy.irp.f | 32 - .../local/cipsi_tc_bi_ortho/environment.irp.f | 14 - .../local/cipsi_tc_bi_ortho/lock_2rdm.irp.f | 0 .../cipsi_tc_bi_ortho/run_pt2_slave.irp.f | 546 ----------- .../run_selection_slave.irp.f | 261 +---- .../local/cipsi_tc_bi_ortho/selection.irp.f | 150 +-- .../cipsi_tc_bi_ortho/selection_buffer.irp.f | 424 --------- .../cipsi_tc_bi_ortho/selection_weight.irp.f | 134 --- .../local/cipsi_tc_bi_ortho/slave_cipsi.irp.f | 348 ------- src/cipsi/cipsi.irp.f | 11 +- src/cipsi/energy.irp.f | 9 - src/cipsi/lock_2rdm.irp.f | 0 src/cipsi/pt2_type.irp.f | 128 --- src/cipsi/run_selection_slave.irp.f | 259 +---- src/cipsi/selection.irp.f | 104 +- src/cipsi/selection_types.f90 | 25 - src/cipsi_utils/README.rst | 5 + src/{cipsi => cipsi_utils}/environment.irp.f | 0 src/cipsi_utils/pt2_stoch_routines.irp.f | 891 ++++++++++++++++++ .../cipsi_utils}/pt2_type.irp.f | 0 .../run_pt2_slave.irp.f | 0 src/cipsi_utils/run_selection_slave.irp.f | 257 +++++ .../selection_buffer.irp.f | 0 .../cipsi_utils}/selection_types.f90 | 0 .../selection_weight.irp.f | 0 src/{cipsi => cipsi_utils}/slave_cipsi.irp.f | 5 +- 26 files changed, 1303 insertions(+), 2300 deletions(-) delete mode 100644 plugins/local/cipsi_tc_bi_ortho/environment.irp.f delete mode 100644 plugins/local/cipsi_tc_bi_ortho/lock_2rdm.irp.f delete mode 100644 plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f delete mode 100644 plugins/local/cipsi_tc_bi_ortho/selection_buffer.irp.f delete mode 100644 plugins/local/cipsi_tc_bi_ortho/selection_weight.irp.f delete mode 100644 plugins/local/cipsi_tc_bi_ortho/slave_cipsi.irp.f delete mode 100644 src/cipsi/lock_2rdm.irp.f delete mode 100644 src/cipsi/pt2_type.irp.f delete mode 100644 src/cipsi/selection_types.f90 create mode 100644 src/cipsi_utils/README.rst rename src/{cipsi => cipsi_utils}/environment.irp.f (100%) create mode 100644 src/cipsi_utils/pt2_stoch_routines.irp.f rename {plugins/local/cipsi_tc_bi_ortho => src/cipsi_utils}/pt2_type.irp.f (100%) rename src/{cipsi => cipsi_utils}/run_pt2_slave.irp.f (100%) create mode 100644 src/cipsi_utils/run_selection_slave.irp.f rename src/{cipsi => cipsi_utils}/selection_buffer.irp.f (100%) rename {plugins/local/cipsi_tc_bi_ortho => src/cipsi_utils}/selection_types.f90 (100%) rename src/{cipsi => cipsi_utils}/selection_weight.irp.f (100%) rename src/{cipsi => cipsi_utils}/slave_cipsi.irp.f (98%) diff --git a/plugins/local/cipsi_tc_bi_ortho/energy.irp.f b/plugins/local/cipsi_tc_bi_ortho/energy.irp.f index 16f4528e..3698e5c2 100644 --- a/plugins/local/cipsi_tc_bi_ortho/energy.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/energy.irp.f @@ -15,37 +15,5 @@ BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ] pt2_E0_denominator = eigval_right_tc_bi_orth -! if (initialize_pt2_E0_denominator) then -! if (h0_type == "EN") then -! pt2_E0_denominator(1:N_states) = psi_energy(1:N_states) -! else if (h0_type == "HF") then -! do i=1,N_states -! j = maxloc(abs(psi_coef(:,i)),1) -! pt2_E0_denominator(i) = psi_det_hii(j) -! enddo -! else if (h0_type == "Barycentric") then -! pt2_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states) -! else if (h0_type == "CFG") then -! pt2_E0_denominator(1:N_states) = psi_energy(1:N_states) -! else -! print *, h0_type, ' not implemented' -! stop -! endif -! do i=1,N_states -! call write_double(6,pt2_E0_denominator(i)+nuclear_repulsion, 'PT2 Energy denominator') -! enddo -! else -! pt2_E0_denominator = -huge(1.d0) -! endif - -END_PROVIDER - - -BEGIN_PROVIDER [ double precision, pt2_overlap, (N_states, N_states) ] - implicit none - BEGIN_DOC - ! Overlap between the perturbed wave functions - END_DOC - pt2_overlap(1:N_states,1:N_states) = 0.d0 END_PROVIDER diff --git a/plugins/local/cipsi_tc_bi_ortho/environment.irp.f b/plugins/local/cipsi_tc_bi_ortho/environment.irp.f deleted file mode 100644 index 5c0e0820..00000000 --- a/plugins/local/cipsi_tc_bi_ortho/environment.irp.f +++ /dev/null @@ -1,14 +0,0 @@ -BEGIN_PROVIDER [ integer, nthreads_pt2 ] - implicit none - BEGIN_DOC - ! Number of threads for Davidson - END_DOC - nthreads_pt2 = nproc - character*(32) :: env - call getenv('QP_NTHREADS_PT2',env) - if (trim(env) /= '') then - read(env,*) nthreads_pt2 - call write_int(6,nthreads_pt2,'Target number of threads for PT2') - endif -END_PROVIDER - diff --git a/plugins/local/cipsi_tc_bi_ortho/lock_2rdm.irp.f b/plugins/local/cipsi_tc_bi_ortho/lock_2rdm.irp.f deleted file mode 100644 index e69de29b..00000000 diff --git a/plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f b/plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f deleted file mode 100644 index d4f45649..00000000 --- a/plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f +++ /dev/null @@ -1,546 +0,0 @@ - use omp_lib - use selection_types - use f77_zmq -BEGIN_PROVIDER [ integer(omp_lock_kind), global_selection_buffer_lock ] - use omp_lib - implicit none - BEGIN_DOC - ! Global buffer for the OpenMP selection - END_DOC - call omp_init_lock(global_selection_buffer_lock) -END_PROVIDER - -BEGIN_PROVIDER [ type(selection_buffer), global_selection_buffer ] - use omp_lib - implicit none - BEGIN_DOC - ! Global buffer for the OpenMP selection - END_DOC - call omp_set_lock(global_selection_buffer_lock) - call delete_selection_buffer(global_selection_buffer) - call create_selection_buffer(N_det_generators, 2*N_det_generators, & - global_selection_buffer) - call omp_unset_lock(global_selection_buffer_lock) -END_PROVIDER - - -subroutine run_pt2_slave(thread,iproc,energy) - use selection_types - use f77_zmq - implicit none - - double precision, intent(in) :: energy(N_states_diag) - integer, intent(in) :: thread, iproc - call run_pt2_slave_large(thread,iproc,energy) -! if (N_det > 100000 ) then -! call run_pt2_slave_large(thread,iproc,energy) -! else -! call run_pt2_slave_small(thread,iproc,energy) -! endif -end - -subroutine run_pt2_slave_small(thread,iproc,energy) - use selection_types - use f77_zmq - implicit none - - double precision, intent(in) :: energy(N_states_diag) - integer, intent(in) :: thread, iproc - integer :: rc, i - - integer :: worker_id, ctask, ltask - character*(512), allocatable :: task(:) - integer, allocatable :: task_id(:) - - integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - - integer(ZMQ_PTR), external :: new_zmq_push_socket - integer(ZMQ_PTR) :: zmq_socket_push - - type(selection_buffer) :: b - logical :: done, buffer_ready - - type(pt2_type), allocatable :: pt2_data(:) - integer :: n_tasks, k, N - integer, allocatable :: i_generator(:), subset(:) - - double precision, external :: memory_of_double, memory_of_int - integer :: bsize ! Size of selection buffers - - allocate(task_id(pt2_n_tasks_max), task(pt2_n_tasks_max)) - allocate(pt2_data(pt2_n_tasks_max), i_generator(pt2_n_tasks_max), subset(pt2_n_tasks_max)) - zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() - - integer, external :: connect_to_taskserver - if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - return - endif - - zmq_socket_push = new_zmq_push_socket(thread) - - b%N = 0 - buffer_ready = .False. - n_tasks = 1 - - done = .False. - do while (.not.done) - - n_tasks = max(1,n_tasks) - n_tasks = min(pt2_n_tasks_max,n_tasks) - - integer, external :: get_tasks_from_taskserver - if (get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task, n_tasks) == -1) then - exit - endif - done = task_id(n_tasks) == 0 - if (done) then - n_tasks = n_tasks-1 - endif - if (n_tasks == 0) exit - - do k=1,n_tasks - call sscanf_ddd(task(k), subset(k), i_generator(k), N) - enddo - if (b%N == 0) then - ! Only first time - bsize = min(N, (elec_alpha_num * (mo_num-elec_alpha_num))**2) - call create_selection_buffer(bsize, bsize*2, b) - buffer_ready = .True. - else - ASSERT (b%N == bsize) - endif - - double precision :: time0, time1 - call wall_time(time0) - do k=1,n_tasks - call pt2_alloc(pt2_data(k),N_states) - b%cur = 0 - call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k))) - enddo - call wall_time(time1) - - integer, external :: tasks_done_to_taskserver - if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then - done = .true. - endif - call sort_selection_buffer(b) - call push_pt2_results(zmq_socket_push, i_generator, pt2_data, b, task_id, n_tasks) - do k=1,n_tasks - call pt2_dealloc(pt2_data(k)) - enddo - b%cur=0 - -! ! Try to adjust n_tasks around nproc/2 seconds per job - n_tasks = min(2*n_tasks,int( dble(n_tasks * nproc/2) / (time1 - time0 + 1.d0))) - n_tasks = min(n_tasks, pt2_n_tasks_max) -! n_tasks = 1 - end do - - integer, external :: disconnect_from_taskserver - do i=1,300 - if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) /= -2) exit - call usleep(500) - print *, 'Retry disconnect...' - end do - - call end_zmq_push_socket(zmq_socket_push,thread) - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - if (buffer_ready) then - call delete_selection_buffer(b) - endif - deallocate(pt2_data) -end subroutine - - -subroutine run_pt2_slave_large(thread,iproc,energy) - use selection_types - use f77_zmq - implicit none - - double precision, intent(in) :: energy(N_states_diag) - integer, intent(in) :: thread, iproc - integer :: rc, i - - integer :: worker_id, ctask, ltask - character*(512) :: task - integer :: task_id(1) - - integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - - integer(ZMQ_PTR), external :: new_zmq_push_socket - integer(ZMQ_PTR) :: zmq_socket_push - - type(selection_buffer) :: b - logical :: done, buffer_ready - - type(pt2_type) :: pt2_data - integer :: n_tasks, k, N - integer :: i_generator, subset - integer :: ifirst - - integer :: bsize ! Size of selection buffers - logical :: sending - PROVIDE global_selection_buffer global_selection_buffer_lock - - - zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() - - integer, external :: connect_to_taskserver - if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - return - endif - - zmq_socket_push = new_zmq_push_socket(thread) - - ifirst = 0 - b%N = 0 - buffer_ready = .False. - n_tasks = 1 - - sending = .False. - done = .False. - do while (.not.done) - - integer, external :: get_tasks_from_taskserver - if (get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task, n_tasks) == -1) then - exit - endif - done = task_id(1) == 0 - if (done) then - n_tasks = n_tasks-1 - endif - if (n_tasks == 0) exit - - call sscanf_ddd(task, subset, i_generator, N) - if( pt2_F(i_generator) <= 0 .or. pt2_F(i_generator) > N_det ) then - print *, irp_here - stop 'bug in selection' - endif - if (b%N == 0) then - ! Only first time - bsize = min(N, (elec_alpha_num * (mo_num-elec_alpha_num))**2) - call create_selection_buffer(bsize, bsize*2, b) - buffer_ready = .True. - else - ASSERT (b%N == bsize) - endif - - double precision :: time0, time1 - call wall_time(time0) - call pt2_alloc(pt2_data,N_states) - b%cur = 0 - call select_connected(i_generator,energy,pt2_data,b,subset,pt2_F(i_generator)) - call wall_time(time1) - - integer, external :: tasks_done_to_taskserver - if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then - done = .true. - endif - call sort_selection_buffer(b) - call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending) - call omp_set_lock(global_selection_buffer_lock) - global_selection_buffer%mini = b%mini - call merge_selection_buffers(b,global_selection_buffer) - if (ifirst /= 0 ) then - b%cur=0 - else - ifirst = 1 - endif - call omp_unset_lock(global_selection_buffer_lock) - if ( iproc == 1 ) then - call omp_set_lock(global_selection_buffer_lock) - call push_pt2_results_async_send(zmq_socket_push, (/i_generator/), (/pt2_data/), global_selection_buffer, (/task_id/), 1,sending) - global_selection_buffer%cur = 0 - call omp_unset_lock(global_selection_buffer_lock) - else - call push_pt2_results_async_send(zmq_socket_push, (/i_generator/), (/pt2_data/), b, (/task_id/), 1,sending) - endif - - call pt2_dealloc(pt2_data) - end do - call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending) - - integer, external :: disconnect_from_taskserver - do i=1,300 - if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) /= -2) exit - call sleep(1) - print *, 'Retry disconnect...' - end do - - call end_zmq_push_socket(zmq_socket_push,thread) - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - if (buffer_ready) then - call delete_selection_buffer(b) - endif - FREE global_selection_buffer -end subroutine - - -subroutine push_pt2_results(zmq_socket_push, index, pt2_data, b, task_id, n_tasks) - use selection_types - use f77_zmq - implicit none - - integer(ZMQ_PTR), intent(in) :: zmq_socket_push - type(pt2_type), intent(in) :: pt2_data(n_tasks) - integer, intent(in) :: n_tasks, index(n_tasks), task_id(n_tasks) - type(selection_buffer), intent(inout) :: b - - logical :: sending - sending = .False. - call push_pt2_results_async_send(zmq_socket_push, index, pt2_data, b, task_id, n_tasks, sending) - call push_pt2_results_async_recv(zmq_socket_push, b%mini, sending) -end subroutine - - -subroutine push_pt2_results_async_send(zmq_socket_push, index, pt2_data, b, task_id, n_tasks, sending) - use selection_types - use f77_zmq - implicit none - - integer(ZMQ_PTR), intent(in) :: zmq_socket_push - type(pt2_type), intent(in) :: pt2_data(n_tasks) - integer, intent(in) :: n_tasks, index(n_tasks), task_id(n_tasks) - type(selection_buffer), intent(inout) :: b - logical, intent(inout) :: sending - integer :: rc, i - integer*8 :: rc8 - double precision, allocatable :: pt2_serialized(:,:) - - if (sending) then - print *, irp_here, ': sending is true' - stop -1 - endif - sending = .True. - - rc = f77_zmq_send( zmq_socket_push, n_tasks, 4, ZMQ_SNDMORE) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 1 - return - else if(rc /= 4) then - stop 'push' - endif - - - rc = f77_zmq_send( zmq_socket_push, index, 4*n_tasks, ZMQ_SNDMORE) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 2 - return - else if(rc /= 4*n_tasks) then - stop 'push' - endif - - - allocate(pt2_serialized (pt2_type_size(N_states),n_tasks) ) - do i=1,n_tasks - call pt2_serialize(pt2_data(i),N_states,pt2_serialized(1,i)) - enddo - - rc = f77_zmq_send( zmq_socket_push, pt2_serialized, size(pt2_serialized)*8, ZMQ_SNDMORE) - deallocate(pt2_serialized) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 3 - return - else if(rc /= size(pt2_serialized)*8) then - stop 'push' - endif - - - rc = f77_zmq_send( zmq_socket_push, task_id, n_tasks*4, ZMQ_SNDMORE) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 6 - return - else if(rc /= 4*n_tasks) then - stop 'push' - endif - - - if (b%cur == 0) then - - rc = f77_zmq_send( zmq_socket_push, b%cur, 4, 0) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 7 - return - else if(rc /= 4) then - stop 'push' - endif - - else - - rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 7 - return - else if(rc /= 4) then - stop 'push' - endif - - - rc8 = f77_zmq_send8( zmq_socket_push, b%val, 8_8*int(b%cur,8), ZMQ_SNDMORE) - if (rc8 == -1_8) then - print *, irp_here, ': error sending result' - stop 8 - return - else if(rc8 /= 8_8*int(b%cur,8)) then - stop 'push' - endif - - - rc8 = f77_zmq_send8( zmq_socket_push, b%det, int(bit_kind*N_int*2,8)*int(b%cur,8), 0) - if (rc8 == -1_8) then - print *, irp_here, ': error sending result' - stop 9 - return - else if(rc8 /= int(N_int*2*8,8)*int(b%cur,8)) then - stop 'push' - endif - - endif - -end subroutine - -subroutine push_pt2_results_async_recv(zmq_socket_push,mini,sending) - use selection_types - use f77_zmq - implicit none - - integer(ZMQ_PTR), intent(in) :: zmq_socket_push - double precision, intent(out) :: mini - logical, intent(inout) :: sending - integer :: rc - - if (.not.sending) return - -! Activate is zmq_socket_push is a REQ -IRP_IF ZMQ_PUSH -IRP_ELSE - character*(2) :: ok - rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 10 - return - else if ((rc /= 2).and.(ok(1:2) /= 'ok')) then - print *, irp_here//': error in receiving ok' - stop -1 - endif - rc = f77_zmq_recv( zmq_socket_push, mini, 8, 0) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 11 - return - else if (rc /= 8) then - print *, irp_here//': error in receiving mini' - stop 12 - endif -IRP_ENDIF - sending = .False. -end subroutine - - - -subroutine pull_pt2_results(zmq_socket_pull, index, pt2_data, task_id, n_tasks, b) - use selection_types - use f77_zmq - implicit none - integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - type(pt2_type), intent(inout) :: pt2_data(*) - type(selection_buffer), intent(inout) :: b - integer, intent(out) :: index(*) - integer, intent(out) :: n_tasks, task_id(*) - integer :: rc, rn, i - integer*8 :: rc8 - double precision, allocatable :: pt2_serialized(:,:) - - rc = f77_zmq_recv( zmq_socket_pull, n_tasks, 4, 0) - if (rc == -1) then - n_tasks = 1 - task_id(1) = 0 - else if(rc /= 4) then - stop 'pull' - endif - - rc = f77_zmq_recv( zmq_socket_pull, index, 4*n_tasks, 0) - if (rc == -1) then - n_tasks = 1 - task_id(1) = 0 - else if(rc /= 4*n_tasks) then - stop 'pull' - endif - - allocate(pt2_serialized (pt2_type_size(N_states),n_tasks) ) - rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized)*n_tasks, 0) - if (rc == -1) then - n_tasks = 1 - task_id(1) = 0 - else if(rc /= 8*size(pt2_serialized)) then - stop 'pull' - endif - - do i=1,n_tasks - call pt2_deserialize(pt2_data(i),N_states,pt2_serialized(1,i)) - enddo - deallocate(pt2_serialized) - - rc = f77_zmq_recv( zmq_socket_pull, task_id, n_tasks*4, 0) - if (rc == -1) then - n_tasks = 1 - task_id(1) = 0 - else if(rc /= 4*n_tasks) then - stop 'pull' - endif - - rc = f77_zmq_recv( zmq_socket_pull, b%cur, 4, 0) - if (rc == -1) then - n_tasks = 1 - task_id(1) = 0 - else if(rc /= 4) then - stop 'pull' - endif - - if (b%cur > 0) then - - rc8 = f77_zmq_recv8( zmq_socket_pull, b%val, 8_8*int(b%cur,8), 0) - if (rc8 == -1_8) then - n_tasks = 1 - task_id(1) = 0 - else if(rc8 /= 8_8*int(b%cur,8)) then - stop 'pull' - endif - - rc8 = f77_zmq_recv8( zmq_socket_pull, b%det, int(bit_kind*N_int*2,8)*int(b%cur,8), 0) - if (rc8 == -1_8) then - n_tasks = 1 - task_id(1) = 0 - else if(rc8 /= int(N_int*2*8,8)*int(b%cur,8)) then - stop 'pull' - endif - - endif - -! Activate is zmq_socket_pull is a REP -IRP_IF ZMQ_PUSH -IRP_ELSE - rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, ZMQ_SNDMORE) - if (rc == -1) then - n_tasks = 1 - task_id(1) = 0 - else if (rc /= 2) then - print *, irp_here//': error in sending ok' - stop -1 - endif - rc = f77_zmq_send( zmq_socket_pull, b%mini, 8, 0) -IRP_ENDIF - -end subroutine - diff --git a/plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f b/plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f index 39c83c4b..aaf2f31d 100644 --- a/plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f @@ -1,258 +1,5 @@ -subroutine run_selection_slave(thread, iproc, energy) - - use f77_zmq - use selection_types - - implicit none - - double precision, intent(in) :: energy(N_states) - integer, intent(in) :: thread, iproc - integer :: rc, i - - integer :: worker_id, task_id(1), ctask, ltask - character*(512) :: task - - integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - - integer(ZMQ_PTR), external :: new_zmq_push_socket - integer(ZMQ_PTR) :: zmq_socket_push - - type(selection_buffer) :: buf, buf2 - logical :: done, buffer_ready - type(pt2_type) :: pt2_data - - PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique - PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order - PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym - PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc weight_selection - - call pt2_alloc(pt2_data,N_states) - - zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() - - integer, external :: connect_to_taskserver - if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - return - endif - - zmq_socket_push = new_zmq_push_socket(thread) - - buf%N = 0 - buffer_ready = .False. - ctask = 1 - - do - integer, external :: get_task_from_taskserver - if (get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task) == -1) then - exit - endif - done = task_id(ctask) == 0 - if (done) then - ctask = ctask - 1 - else - integer :: i_generator, N, subset, bsize - call sscanf_ddd(task, subset, i_generator, N) - if(buf%N == 0) then - ! Only first time - call create_selection_buffer(N, N*2, buf) - buffer_ready = .True. - else - if (N /= buf%N) then - print *, 'N=', N - print *, 'buf%N=', buf%N - print *, 'bug in ', irp_here - stop '-1' - end if - end if - call select_connected(i_generator, energy, pt2_data, buf, subset, pt2_F(i_generator)) - endif - - integer, external :: task_done_to_taskserver - - if(done .or. ctask == size(task_id)) then - do i=1, ctask - if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i)) == -1) then - call usleep(100) - if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i)) == -1) then - ctask = 0 - done = .true. - exit - endif - endif - end do - if(ctask > 0) then - call sort_selection_buffer(buf) -! call merge_selection_buffers(buf,buf2) - call push_selection_results(zmq_socket_push, pt2_data, buf, task_id(1), ctask) - call pt2_dealloc(pt2_data) - call pt2_alloc(pt2_data,N_states) -! buf%mini = buf2%mini - buf%cur = 0 - end if - ctask = 0 - end if - - if(done) exit - ctask = ctask + 1 - end do - - if(ctask > 0) then - call sort_selection_buffer(buf) -! call merge_selection_buffers(buf,buf2) - call push_selection_results(zmq_socket_push, pt2_data, buf, task_id(1), ctask) -! buf%mini = buf2%mini - buf%cur = 0 - end if - ctask = 0 - call pt2_dealloc(pt2_data) - - integer, external :: disconnect_from_taskserver - if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then - continue - endif - - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - call end_zmq_push_socket(zmq_socket_push,thread) - if (buffer_ready) then - call delete_selection_buffer(buf) -! call delete_selection_buffer(buf2) - endif -end subroutine - - -subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntasks) - use f77_zmq - use selection_types - implicit none - - integer(ZMQ_PTR), intent(in) :: zmq_socket_push - type(pt2_type), intent(in) :: pt2_data - type(selection_buffer), intent(inout) :: b - integer, intent(in) :: ntasks, task_id(*) - integer :: rc - double precision, allocatable :: pt2_serialized(:) - - rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE) - if(rc /= 4) then - print *, 'f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE)' - endif - - - allocate(pt2_serialized (pt2_type_size(N_states)) ) - call pt2_serialize(pt2_data,N_states,pt2_serialized) - - rc = f77_zmq_send( zmq_socket_push, pt2_serialized, size(pt2_serialized)*8, ZMQ_SNDMORE) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 3 - return - else if(rc /= size(pt2_serialized)*8) then - stop 'push' - endif - deallocate(pt2_serialized) - - if (b%cur > 0) then - - rc = f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE) - if(rc /= 8*b%cur) then - print *, 'f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE)' - endif - - rc = f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE) - if(rc /= bit_kind*N_int*2*b%cur) then - print *, 'f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE)' - endif - - endif - - rc = f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE) - if(rc /= 4) then - print *, 'f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE)' - endif - - rc = f77_zmq_send( zmq_socket_push, task_id(1), ntasks*4, 0) - if(rc /= 4*ntasks) then - print *, 'f77_zmq_send( zmq_socket_push, task_id(1), ntasks*4, 0)' - endif - -! Activate is zmq_socket_push is a REQ -IRP_IF ZMQ_PUSH -IRP_ELSE - character*(2) :: ok - rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0) - if ((rc /= 2).and.(ok(1:2) /= 'ok')) then - print *, irp_here//': error in receiving ok' - stop -1 - endif -IRP_ENDIF - -end subroutine - - -subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_id, ntasks) - use f77_zmq - use selection_types - implicit none - integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - type(pt2_type), intent(inout) :: pt2_data - double precision, intent(out) :: val(*) - integer(bit_kind), intent(out) :: det(N_int, 2, *) - integer, intent(out) :: N, ntasks, task_id(*) - integer :: rc, rn, i - double precision, allocatable :: pt2_serialized(:) - - rc = f77_zmq_recv( zmq_socket_pull, N, 4, 0) - if(rc /= 4) then - print *, 'f77_zmq_recv( zmq_socket_pull, N, 4, 0)' - endif - - allocate(pt2_serialized (pt2_type_size(N_states)) ) - rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized), 0) - if (rc == -1) then - ntasks = 1 - task_id(1) = 0 - else if(rc /= 8*size(pt2_serialized)) then - stop 'pull' - endif - - call pt2_deserialize(pt2_data,N_states,pt2_serialized) - deallocate(pt2_serialized) - - if (N>0) then - rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, 0) - if(rc /= 8*N) then - print *, 'f77_zmq_recv( zmq_socket_pull, val(1), 8*N, 0)' - endif - - rc = f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, 0) - if(rc /= bit_kind*N_int*2*N) then - print *, 'f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, 0)' - endif - endif - - rc = f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0) - if(rc /= 4) then - print *, 'f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0)' - endif - - rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntasks*4, 0) - if(rc /= 4*ntasks) then - print *, 'f77_zmq_recv( zmq_socket_pull, task_id(1), ntasks*4, 0)' - endif - -! Activate is zmq_socket_pull is a REP -IRP_IF ZMQ_PUSH -IRP_ELSE - rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0) - if (rc /= 2) then - print *, irp_here//': error in sending ok' - stop -1 - endif -IRP_ENDIF -end subroutine - - +subroutine provide_for_selection_slave + PROVIDE psi_det_sorted_tc_order + PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc +end diff --git a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f index 06cf848b..9b8cc81e 100644 --- a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f @@ -76,6 +76,8 @@ subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset) double precision, allocatable :: fock_diag_tmp(:,:) + if (csubset == 0) return + 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) @@ -86,10 +88,13 @@ subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset) 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 + if ((subset == 1).and.(sum(hole_mask(:,2)) == 0_bit_kind)) then + ! No beta electron to excite + call select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,b) + endif 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 - +end subroutine double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2, Nint) @@ -136,7 +141,7 @@ double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2, Nint) end -subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock_diag_tmp, E0, pt2_data, buf, subset, csubset) +subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, fock_diag_tmp, E0, pt2_data, buf, subset, csubset) use bitmasks use selection_types implicit none @@ -151,8 +156,6 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock type(pt2_type), intent(inout) :: pt2_data type(selection_buffer), intent(inout) :: buf - double precision, parameter :: norm_thr = 1.d-16 - integer :: h1, h2, s1, s2, s3, i1, i2, ib, sp, k, i, j, nt, ii, sze integer :: maskInd integer :: N_holes(2), N_particles(2) @@ -170,6 +173,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock integer, allocatable :: preinteresting(:), prefullinteresting(:) integer, allocatable :: interesting(:), fullinteresting(:) integer, allocatable :: tmp_array(:) + integer, allocatable :: indices(:), exc_degree(:), iorder(:) integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) logical, allocatable :: banned(:,:,:), bannedOrb(:,:) @@ -178,15 +182,16 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique - PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order + PROVIDE psi_bilinear_matrix_rows psi_bilinear_matrix_order psi_bilinear_matrix_transp_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp_tc + PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc_order PROVIDE banned_excitation monoAdo = .true. monoBdo = .true. + if (csubset == 0) return do k=1,N_int hole (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1)) @@ -198,7 +203,11 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock call bitstring_to_list_ab(hole , hole_list , N_holes , N_int) call bitstring_to_list_ab(particle, particle_list, N_particles, N_int) - allocate( indices(N_det), exc_degree( max(N_det_alpha_unique, N_det_beta_unique) ) ) + ! Removed to avoid introducing determinants already presents in the wf + !double precision, parameter :: norm_thr = 1.d-16 + + allocate (indices(N_det), & + exc_degree(max(N_det_alpha_unique,N_det_beta_unique))) ! Pre-compute excitation degrees wrt alpha determinants k=1 @@ -214,73 +223,76 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock if (nt > 2) cycle 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 + 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 + ! Removed to avoid introducing determinants already presents in the wf + !if (psi_average_norm_contrib_sorted_tc(idx) > norm_thr) then indices(k) = idx - k = k + 1 -! endif + k=k+1 + !endif endif enddo enddo ! Pre-compute excitation degrees wrt beta determinants 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) + 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 - call get_excitation_degree_spin(psi_det_alpha_unique(1,j), psi_det_generators(1,1,i_generator), nt, N_int) + 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 + do l_a=psi_bilinear_matrix_transp_rows_loc(j), psi_bilinear_matrix_transp_rows_loc(j+1)-1 i = psi_bilinear_matrix_transp_columns(l_a) - if(exc_degree(i) < 3) cycle - if(nt + exc_degree(i) <= 4) then + if (exc_degree(i) < 3) cycle + if (nt + exc_degree(i) <= 4) then 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 + ! Removed to avoid introducing determinants already presents in the wf + !if(psi_average_norm_contrib_sorted_tc(idx) > norm_thr) then indices(k) = idx - k = k + 1 -! endif + k=k+1 + !endif endif enddo enddo deallocate(exc_degree) - nmax = k - 1 + nmax=k-1 call isort_noidx(indices,nmax) ! Start with 32 elements. Size will double along with the filtering. - allocate(preinteresting(0:32), prefullinteresting(0:32), interesting(0:32), fullinteresting(0:32)) + allocate(preinteresting(0:32), prefullinteresting(0:32), & + interesting(0:32), fullinteresting(0:32)) preinteresting(:) = 0 prefullinteresting(:) = 0 - do i = 1, N_int + do i=1,N_int negMask(i,1) = not(psi_det_generators(i,1,i_generator)) negMask(i,2) = not(psi_det_generators(i,2,i_generator)) - enddo - - do k = 1, nmax + end do + do k=1,nmax i = indices(k) mobMask(1,1) = iand(negMask(1,1), psi_det_sorted_tc(1,1,i)) mobMask(1,2) = iand(negMask(1,2), psi_det_sorted_tc(1,2,i)) nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) - do j = 2, N_int + do j=2,N_int mobMask(j,1) = iand(negMask(j,1), psi_det_sorted_tc(j,1,i)) mobMask(j,2) = iand(negMask(j,2), psi_det_sorted_tc(j,2,i)) nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) - enddo + end do if(nt <= 4) then if(i <= N_det_selectors) then sze = preinteresting(0) - if(sze+1 == size(preinteresting)) then - allocate(tmp_array(0:sze)) + if (sze+1 == size(preinteresting)) then + allocate (tmp_array(0:sze)) tmp_array(0:sze) = preinteresting(0:sze) deallocate(preinteresting) allocate(preinteresting(0:2*sze)) @@ -289,9 +301,9 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock endif preinteresting(0) = sze+1 preinteresting(sze+1) = i - elseif(nt <= 2) then + else if(nt <= 2) then sze = prefullinteresting(0) - if(sze+1 == size(prefullinteresting)) then + if (sze+1 == size(prefullinteresting)) then allocate (tmp_array(0:sze)) tmp_array(0:sze) = prefullinteresting(0:sze) deallocate(prefullinteresting) @@ -301,20 +313,16 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock endif prefullinteresting(0) = sze+1 prefullinteresting(sze+1) = i - endif - endif - - enddo + end if + end if + end do deallocate(indices) - allocate( banned(mo_num, mo_num,2), bannedOrb(mo_num, 2) ) - allocate( mat(N_states, mo_num, mo_num) ) - allocate( mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num) ) + allocate(banned(mo_num, mo_num,2), bannedOrb(mo_num, 2)) + allocate(mat(N_states, mo_num, mo_num)) + 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 @@ -347,17 +355,17 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock do ii = 1, preinteresting(0) i = preinteresting(ii) - select case(N_int) - case(1) + select case (N_int) + case (1) mobMask(1,1) = iand(negMask(1,1), psi_det_sorted_tc(1,1,i)) mobMask(1,2) = iand(negMask(1,2), psi_det_sorted_tc(1,2,i)) nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) - case(2) + case (2) mobMask(1:2,1) = iand(negMask(1:2,1), psi_det_sorted_tc(1:2,1,i)) mobMask(1:2,2) = iand(negMask(1:2,2), psi_det_sorted_tc(1:2,2,i)) nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) + & popcnt(mobMask(2, 1)) + popcnt(mobMask(2, 2)) - case(3) + case (3) mobMask(1:3,1) = iand(negMask(1:3,1), psi_det_sorted_tc(1:3,1,i)) mobMask(1:3,2) = iand(negMask(1:3,2), psi_det_sorted_tc(1:3,2,i)) nt = 0 @@ -370,8 +378,8 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock nt = nt+ popcnt(mobMask(j, 2)) if (nt > 4) exit endif - enddo - case(4) + end do + case (4) mobMask(1:4,1) = iand(negMask(1:4,1), psi_det_sorted_tc(1:4,1,i)) mobMask(1:4,2) = iand(negMask(1:4,2), psi_det_sorted_tc(1:4,2,i)) nt = 0 @@ -384,7 +392,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock nt = nt+ popcnt(mobMask(j, 2)) if (nt > 4) exit endif - enddo + end do case default mobMask(1:N_int,1) = iand(negMask(1:N_int,1), psi_det_sorted_tc(1:N_int,1,i)) mobMask(1:N_int,2) = iand(negMask(1:N_int,2), psi_det_sorted_tc(1:N_int,2,i)) @@ -398,12 +406,12 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock nt = nt+ popcnt(mobMask(j, 2)) if (nt > 4) exit endif - enddo + end do end select if(nt <= 4) then sze = interesting(0) - if(sze+1 == size(interesting)) then + if (sze+1 == size(interesting)) then allocate (tmp_array(0:sze)) tmp_array(0:sze) = interesting(0:sze) deallocate(interesting) @@ -425,8 +433,8 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock endif fullinteresting(0) = sze+1 fullinteresting(sze+1) = i - endif - endif + end if + end if enddo @@ -456,10 +464,10 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock endif fullinteresting(0) = sze+1 fullinteresting(sze+1) = i - endif - enddo - allocate( fullminilist (N_int, 2, fullinteresting(0)), & - minilist (N_int, 2, interesting(0)) ) + end if + end do + allocate (fullminilist (N_int, 2, fullinteresting(0)), & + minilist (N_int, 2, interesting(0)) ) do i = 1, fullinteresting(0) do k = 1, N_int @@ -517,7 +525,8 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting, mat_l, mat_r) call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, mat_l, mat_r) - endif + end if + enddo @@ -533,7 +542,8 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock deallocate(banned, bannedOrb,mat) deallocate(mat_l, mat_r) -end subroutine select_singles_and_doubles + +end subroutine ! --- @@ -924,13 +934,13 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d 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) print*,i_h_alpha,alpha_h_i - call debug_det(psi_selectors(1,1,iii),N_int) - enddo + call debug_det(psi_selectors(1,1,iii),N_int) + enddo ! print*,'psi_det ' ! do iii = 1, N_det! old version ! print*,'iii',iii,psi_l_coef_bi_ortho(iii,1),psi_r_coef_bi_ortho(iii,1) -! call debug_det(psi_det(1,1,iii),N_int) -! enddo +! call debug_det(psi_det(1,1,iii),N_int) +! enddo stop endif endif @@ -938,7 +948,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d psi_h_alpha = mat_l(istate, p1, p2) alpha_h_psi = mat_r(istate, p1, p2) endif - val = 4.d0 * psi_h_alpha * alpha_h_psi + val = 4.d0 * psi_h_alpha * alpha_h_psi tmp = dsqrt(delta_E * delta_E + val) ! if (delta_E < 0.d0) then ! tmp = -tmp @@ -946,21 +956,21 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d e_pert(istate) = 0.25 * val / delta_E ! e_pert(istate) = 0.5d0 * (tmp - delta_E) if(dsqrt(dabs(tmp)).gt.1.d-4.and.dabs(alpha_h_psi).gt.1.d-4)then - coef(istate) = e_pert(istate) / psi_h_alpha + coef(istate) = e_pert(istate) / psi_h_alpha else - coef(istate) = alpha_h_psi / delta_E + coef(istate) = alpha_h_psi / delta_E endif if(selection_tc == 1)then - if(e_pert(istate).lt.0.d0)then + if(e_pert(istate).lt.0.d0)then e_pert(istate)=0.d0 - else + else e_pert(istate)=-e_pert(istate) endif else if(selection_tc == -1)then if(e_pert(istate).gt.0.d0)e_pert(istate)=0.d0 endif - + ! if(selection_tc == 1 )then ! if(e_pert(istate).lt.0.d0)then ! e_pert(istate) = 0.d0 diff --git a/plugins/local/cipsi_tc_bi_ortho/selection_buffer.irp.f b/plugins/local/cipsi_tc_bi_ortho/selection_buffer.irp.f deleted file mode 100644 index 0bd51464..00000000 --- a/plugins/local/cipsi_tc_bi_ortho/selection_buffer.irp.f +++ /dev/null @@ -1,424 +0,0 @@ - -subroutine create_selection_buffer(N, size_in, res) - use selection_types - implicit none - BEGIN_DOC -! Allocates the memory for a selection buffer. -! The arrays have dimension size_in and the maximum number of elements is N - END_DOC - - integer, intent(in) :: N, size_in - type(selection_buffer), intent(out) :: res - - integer :: siz - siz = max(size_in,1) - - double precision :: rss - double precision, external :: memory_of_double - rss = memory_of_double(siz)*(N_int*2+1) - call check_mem(rss,irp_here) - - allocate(res%det(N_int, 2, siz), res%val(siz)) - - res%val(:) = 0d0 - res%det(:,:,:) = 0_8 - res%N = N - res%mini = 0d0 - res%cur = 0 -end subroutine - -subroutine delete_selection_buffer(b) - use selection_types - implicit none - type(selection_buffer), intent(inout) :: b - if (associated(b%det)) then - deallocate(b%det) - endif - if (associated(b%val)) then - deallocate(b%val) - endif - NULLIFY(b%det) - NULLIFY(b%val) - b%cur = 0 - b%mini = 0.d0 - b%N = 0 -end - - -subroutine add_to_selection_buffer(b, det, val) - use selection_types - implicit none - - type(selection_buffer), intent(inout) :: b - integer(bit_kind), intent(in) :: det(N_int, 2) - double precision, intent(in) :: val - integer :: i - - if(b%N > 0 .and. val <= b%mini) then - b%cur += 1 - b%det(1:N_int,1:2,b%cur) = det(1:N_int,1:2) - b%val(b%cur) = val - if(b%cur == size(b%val)) then - call sort_selection_buffer(b) - end if - end if -end subroutine - -subroutine merge_selection_buffers(b1, b2) - use selection_types - implicit none - BEGIN_DOC -! Merges the selection buffers b1 and b2 into b2 - END_DOC - type(selection_buffer), intent(inout) :: b1 - type(selection_buffer), intent(inout) :: b2 - integer(bit_kind), pointer :: detmp(:,:,:) - double precision, pointer :: val(:) - integer :: i, i1, i2, k, nmwen, sze - if (b1%cur == 0) return - do while (b1%val(b1%cur) > b2%mini) - b1%cur = b1%cur-1 - if (b1%cur == 0) then - return - endif - enddo - nmwen = min(b1%N, b1%cur+b2%cur) - double precision :: rss - double precision, external :: memory_of_double - sze = max(size(b1%val), size(b2%val)) - rss = memory_of_double(sze) + 2*N_int*memory_of_double(sze) - call check_mem(rss,irp_here) - allocate(val(sze), detmp(N_int, 2, sze)) - i1=1 - i2=1 - do i=1,nmwen - if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then - exit - else if (i1 > b1%cur) then - val(i) = b2%val(i2) - detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2) - detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2) - i2=i2+1 - else if (i2 > b2%cur) then - val(i) = b1%val(i1) - detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1) - detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1) - i1=i1+1 - else - if (b1%val(i1) <= b2%val(i2)) then - val(i) = b1%val(i1) - detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1) - detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1) - i1=i1+1 - else - val(i) = b2%val(i2) - detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2) - detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2) - i2=i2+1 - endif - endif - enddo - deallocate(b2%det, b2%val) - do i=nmwen+1,b2%N - val(i) = 0.d0 - detmp(1:N_int,1:2,i) = 0_bit_kind - enddo - b2%det => detmp - b2%val => val -! if(selection_tc == 1)then -! b2%mini = max(b2%mini,b2%val(b2%N)) -! else - b2%mini = min(b2%mini,b2%val(b2%N)) -! endif - b2%cur = nmwen -end - - -subroutine sort_selection_buffer(b) - use selection_types - implicit none - - type(selection_buffer), intent(inout) :: b - integer, allocatable :: iorder(:) - integer(bit_kind), pointer :: detmp(:,:,:) - integer :: i, nmwen - logical, external :: detEq - if (b%N == 0 .or. b%cur == 0) return - nmwen = min(b%N, b%cur) - - double precision :: rss - double precision, external :: memory_of_double, memory_of_int - rss = memory_of_int(b%cur) + 2*N_int*memory_of_double(size(b%det,3)) - call check_mem(rss,irp_here) - allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3))) - do i=1,b%cur - iorder(i) = i - end do - call dsort(b%val, iorder, b%cur) - do i=1, nmwen - detmp(1:N_int,1,i) = b%det(1:N_int,1,iorder(i)) - detmp(1:N_int,2,i) = b%det(1:N_int,2,iorder(i)) - end do - deallocate(b%det,iorder) - b%det => detmp -! if(selection_tc == 1)then -! b%mini = max(b%mini,b%val(b%N)) -! else - b%mini = min(b%mini,b%val(b%N)) -! endif - b%cur = nmwen -end subroutine - -subroutine make_selection_buffer_s2(b) - use selection_types - type(selection_buffer), intent(inout) :: b - - integer(bit_kind), allocatable :: o(:,:,:) - double precision, allocatable :: val(:) - - integer :: n_d - integer :: i,k,sze,n_alpha,j,n - logical :: dup - - ! Sort - integer, allocatable :: iorder(:) - integer*8, allocatable :: bit_tmp(:) - integer*8, external :: configuration_search_key - integer(bit_kind), allocatable :: tmp_array(:,:,:) - logical, allocatable :: duplicate(:) - - n_d = b%cur - double precision :: rss - double precision, external :: memory_of_double - rss = (4*N_int+4)*memory_of_double(n_d) - call check_mem(rss,irp_here) - allocate(o(N_int,2,n_d), iorder(n_d), duplicate(n_d), bit_tmp(n_d), & - tmp_array(N_int,2,n_d), val(n_d) ) - - do i=1,n_d - do k=1,N_int - o(k,1,i) = ieor(b%det(k,1,i), b%det(k,2,i)) - o(k,2,i) = iand(b%det(k,1,i), b%det(k,2,i)) - enddo - iorder(i) = i - bit_tmp(i) = configuration_search_key(o(1,1,i),N_int) - enddo - - deallocate(b%det) - - call i8sort(bit_tmp,iorder,n_d) - - do i=1,n_d - do k=1,N_int - tmp_array(k,1,i) = o(k,1,iorder(i)) - tmp_array(k,2,i) = o(k,2,iorder(i)) - enddo - val(i) = b%val(iorder(i)) - duplicate(i) = .False. - enddo - - ! Find duplicates - do i=1,n_d-1 - if (duplicate(i)) then - cycle - endif - j = i+1 - do while (bit_tmp(j)==bit_tmp(i)) - if (duplicate(j)) then - j+=1 - if (j>n_d) then - exit - endif - cycle - endif - dup = .True. - do k=1,N_int - if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) & - .or. (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then - dup = .False. - exit - endif - enddo - if (dup) then - val(i) = max(val(i), val(j)) - duplicate(j) = .True. - endif - j+=1 - if (j>n_d) then - exit - endif - enddo - enddo - - deallocate (b%val) - ! Copy filtered result - integer :: n_p - n_p=0 - do i=1,n_d - if (duplicate(i)) then - cycle - endif - n_p = n_p + 1 - do k=1,N_int - o(k,1,n_p) = tmp_array(k,1,i) - o(k,2,n_p) = tmp_array(k,2,i) - enddo - val(n_p) = val(i) - enddo - - ! Sort by importance - do i=1,n_p - iorder(i) = i - end do - call dsort(val,iorder,n_p) - do i=1,n_p - do k=1,N_int - tmp_array(k,1,i) = o(k,1,iorder(i)) - tmp_array(k,2,i) = o(k,2,iorder(i)) - enddo - enddo - do i=1,n_p - do k=1,N_int - o(k,1,i) = tmp_array(k,1,i) - o(k,2,i) = tmp_array(k,2,i) - enddo - enddo - - ! Create determinants - n_d = 0 - do i=1,n_p - call configuration_to_dets_size(o(1,1,i),sze,elec_alpha_num,N_int) - n_d = n_d + sze - if (n_d > b%cur) then -! if (n_d - b%cur > b%cur - n_d + sze) then -! n_d = n_d - sze -! endif - exit - endif - enddo - - rss = (4*N_int+2)*memory_of_double(n_d) - call check_mem(rss,irp_here) - allocate(b%det(N_int,2,2*n_d), b%val(2*n_d)) - k=1 - do i=1,n_p - n=n_d - call configuration_to_dets_size(o(1,1,i),n,elec_alpha_num,N_int) - call configuration_to_dets(o(1,1,i),b%det(1,1,k),n,elec_alpha_num,N_int) - do j=k,k+n-1 - b%val(j) = val(i) - enddo - k = k+n - if (k > n_d) exit - enddo - deallocate(o) - b%cur = n_d - b%N = n_d -end - - - - -subroutine remove_duplicates_in_selection_buffer(b) - use selection_types - type(selection_buffer), intent(inout) :: b - - integer(bit_kind), allocatable :: o(:,:,:) - double precision, allocatable :: val(:) - - integer :: n_d - integer :: i,k,sze,n_alpha,j,n - logical :: dup - - ! Sort - integer, allocatable :: iorder(:) - integer*8, allocatable :: bit_tmp(:) - integer*8, external :: det_search_key - integer(bit_kind), allocatable :: tmp_array(:,:,:) - logical, allocatable :: duplicate(:) - - n_d = b%cur - logical :: found_duplicates - double precision :: rss - double precision, external :: memory_of_double - rss = (4*N_int+4)*memory_of_double(n_d) - call check_mem(rss,irp_here) - - found_duplicates = .False. - allocate(iorder(n_d), duplicate(n_d), bit_tmp(n_d), & - tmp_array(N_int,2,n_d), val(n_d) ) - - do i=1,n_d - iorder(i) = i - bit_tmp(i) = det_search_key(b%det(1,1,i),N_int) - enddo - - call i8sort(bit_tmp,iorder,n_d) - - do i=1,n_d - do k=1,N_int - tmp_array(k,1,i) = b%det(k,1,iorder(i)) - tmp_array(k,2,i) = b%det(k,2,iorder(i)) - enddo - val(i) = b%val(iorder(i)) - duplicate(i) = .False. - enddo - - ! Find duplicates - do i=1,n_d-1 - if (duplicate(i)) then - cycle - endif - j = i+1 - do while (bit_tmp(j)==bit_tmp(i)) - if (duplicate(j)) then - j+=1 - if (j>n_d) then - exit - endif - cycle - endif - dup = .True. - do k=1,N_int - if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) & - .or. (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then - dup = .False. - exit - endif - enddo - if (dup) then - duplicate(j) = .True. - found_duplicates = .True. - endif - j+=1 - if (j>n_d) then - exit - endif - enddo - enddo - - if (found_duplicates) then - - ! Copy filtered result - integer :: n_p - n_p=0 - do i=1,n_d - if (duplicate(i)) then - cycle - endif - n_p = n_p + 1 - do k=1,N_int - b%det(k,1,n_p) = tmp_array(k,1,i) - b%det(k,2,n_p) = tmp_array(k,2,i) - enddo - val(n_p) = val(i) - enddo - b%cur=n_p - b%N=n_p - - endif - -end - - - diff --git a/plugins/local/cipsi_tc_bi_ortho/selection_weight.irp.f b/plugins/local/cipsi_tc_bi_ortho/selection_weight.irp.f deleted file mode 100644 index 3c09e59a..00000000 --- a/plugins/local/cipsi_tc_bi_ortho/selection_weight.irp.f +++ /dev/null @@ -1,134 +0,0 @@ -BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ] - implicit none - BEGIN_DOC - ! Weights adjusted along the selection to make the PT2 contributions - ! of each state coincide. - END_DOC - pt2_match_weight(:) = 1.d0 -END_PROVIDER - - - -BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ] - implicit none - BEGIN_DOC - ! Weights adjusted along the selection to make the variances - ! of each state coincide. - END_DOC - variance_match_weight(:) = 1.d0 -END_PROVIDER - - - -subroutine update_pt2_and_variance_weights(pt2_data, N_st) - implicit none - use selection_types - BEGIN_DOC -! Updates the PT2- and Variance- matching weights. - END_DOC - integer, intent(in) :: N_st - type(pt2_type), intent(in) :: pt2_data - double precision :: pt2(N_st) - double precision :: variance(N_st) - - double precision :: avg, element, dt, x - integer :: k - pt2(:) = pt2_data % pt2(:) - variance(:) = pt2_data % variance(:) - - avg = sum(pt2(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero - - dt = 8.d0 !* selection_factor - do k=1,N_st - element = exp(dt*(pt2(k)/avg - 1.d0)) - element = min(2.0d0 , element) - element = max(0.5d0 , element) - pt2_match_weight(k) *= element - enddo - - - avg = sum(variance(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero - - do k=1,N_st - element = exp(dt*(variance(k)/avg -1.d0)) - element = min(2.0d0 , element) - element = max(0.5d0 , element) - variance_match_weight(k) *= element - enddo - - if (N_det < 100) then - ! For tiny wave functions, weights are 1.d0 - pt2_match_weight(:) = 1.d0 - variance_match_weight(:) = 1.d0 - endif - - threshold_davidson_pt2 = min(1.d-6, & - max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) ) - - SOFT_TOUCH pt2_match_weight variance_match_weight threshold_davidson_pt2 -end - - - - -BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] - implicit none - BEGIN_DOC - ! Weights used in the selection criterion - END_DOC - select case (weight_selection) - - case (0) - print *, 'Using input weights in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * state_average_weight(1:N_states) - - case (1) - print *, 'Using 1/c_max^2 weight in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) - - case (2) - print *, 'Using pt2-matching weight in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) - print *, '# PT2 weight ', real(pt2_match_weight(:),4) - - case (3) - print *, 'Using variance-matching weight in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) - print *, '# var weight ', real(variance_match_weight(:),4) - - case (4) - print *, 'Using variance- and pt2-matching weights in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) - print *, '# PT2 weight ', real(pt2_match_weight(:),4) - print *, '# var weight ', real(variance_match_weight(:),4) - - case (5) - print *, 'Using variance-matching weight in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) - print *, '# var weight ', real(variance_match_weight(:),4) - - case (6) - print *, 'Using CI coefficient-based selection' - selection_weight(1:N_states) = c0_weight(1:N_states) - - case (7) - print *, 'Input weights multiplied by variance- and pt2-matching' - selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) * state_average_weight(1:N_states) - print *, '# PT2 weight ', real(pt2_match_weight(:),4) - print *, '# var weight ', real(variance_match_weight(:),4) - - case (8) - print *, 'Input weights multiplied by pt2-matching' - selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) * state_average_weight(1:N_states) - print *, '# PT2 weight ', real(pt2_match_weight(:),4) - - case (9) - print *, 'Input weights multiplied by variance-matching' - selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * state_average_weight(1:N_states) - print *, '# var weight ', real(variance_match_weight(:),4) - - end select - print *, '# Total weight ', real(selection_weight(:),4) - -END_PROVIDER - diff --git a/plugins/local/cipsi_tc_bi_ortho/slave_cipsi.irp.f b/plugins/local/cipsi_tc_bi_ortho/slave_cipsi.irp.f deleted file mode 100644 index 6343bf8b..00000000 --- a/plugins/local/cipsi_tc_bi_ortho/slave_cipsi.irp.f +++ /dev/null @@ -1,348 +0,0 @@ -subroutine run_slave_cipsi - - BEGIN_DOC - ! Helper program for distributed parallelism - END_DOC - - implicit none - - call omp_set_max_active_levels(1) - distributed_davidson = .False. - read_wf = .False. - SOFT_TOUCH read_wf distributed_davidson - call provide_everything - call switch_qp_run_to_master - call run_slave_main -end - -subroutine provide_everything - PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag - - PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context - PROVIDE psi_det psi_coef threshold_generators state_average_weight - PROVIDE N_det_selectors pt2_stoch_istate N_det selection_weight pseudo_sym -end - - -subroutine run_slave_main - - use f77_zmq - - implicit none - IRP_IF MPI - include 'mpif.h' - IRP_ENDIF - - integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - double precision :: energy(N_states) - character*(64) :: states(10) - character*(64) :: old_state - integer :: rc, i, ierr - double precision :: t0, t1 - - integer, external :: zmq_get_dvector, zmq_get_N_det_generators - integer, external :: zmq_get8_dvector - integer, external :: zmq_get_ivector - integer, external :: zmq_get_psi, zmq_get_N_det_selectors, zmq_get_psi_bilinear - integer, external :: zmq_get_psi_notouch - integer, external :: zmq_get_N_states_diag - - zmq_context = f77_zmq_ctx_new () - states(1) = 'selection' - states(2) = 'davidson' - states(3) = 'pt2' - old_state = 'Waiting' - - zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() - - PROVIDE psi_det psi_coef threshold_generators state_average_weight mpi_master - PROVIDE zmq_state N_det_selectors pt2_stoch_istate N_det pt2_e0_denominator - PROVIDE N_det_generators N_states N_states_diag pt2_e0_denominator mpi_rank - - IRP_IF MPI - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - do - - if (mpi_master) then - call wait_for_states(states,zmq_state,size(states)) - if (zmq_state(1:64) == old_state(1:64)) then - call usleep(200) - cycle - else - old_state(1:64) = zmq_state(1:64) - endif - print *, trim(zmq_state) - endif - - IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - IRP_IF MPI - call MPI_BCAST (zmq_state, 128, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - print *, irp_here, 'error in broadcast of zmq_state' - endif - IRP_ENDIF - - if(zmq_state(1:7) == 'Stopped') then - exit - endif - - - if (zmq_state(1:9) == 'selection') then - - ! Selection - ! --------- - - call wall_time(t0) - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_psi') - IRP_ENDIF - if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_dvector threshold_generators') - IRP_ENDIF - if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_generators',(/threshold_generators/),1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_dvector energy') - IRP_ENDIF - if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_N_det_generators') - IRP_ENDIF - if (zmq_get_N_det_generators (zmq_to_qp_run_socket, 1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_N_det_selectors') - IRP_ENDIF - if (zmq_get_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_dvector state_average_weight') - IRP_ENDIF - if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_dvector selection_weight') - IRP_ENDIF - if (zmq_get_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) cycle - pt2_e0_denominator(1:N_states) = energy(1:N_states) - TOUCH pt2_e0_denominator state_average_weight threshold_generators selection_weight psi_det psi_coef - - if (mpi_master) then - print *, 'N_det', N_det - print *, 'N_det_generators', N_det_generators - print *, 'N_det_selectors', N_det_selectors - print *, 'pt2_e0_denominator', pt2_e0_denominator - print *, 'pt2_stoch_istate', pt2_stoch_istate - print *, 'state_average_weight', state_average_weight - print *, 'selection_weight', selection_weight - endif - call wall_time(t1) - call write_double(6,(t1-t0),'Broadcast time') - - IRP_IF MPI_DEBUG - call mpi_print('Entering OpenMP section') - IRP_ENDIF - !$OMP PARALLEL PRIVATE(i) - i = omp_get_thread_num() - call run_selection_slave(0,i,energy) - !$OMP END PARALLEL - print *, mpi_rank, ': Selection done' - IRP_IF MPI - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - print *, irp_here, 'error in barrier' - endif - IRP_ENDIF - call mpi_print('----------') - - else if (zmq_state(1:8) == 'davidson') then - - ! Davidson - ! -------- - - call wall_time(t0) - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_N_states_diag') - IRP_ENDIF - if (zmq_get_N_states_diag(zmq_to_qp_run_socket,1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_psi') - IRP_ENDIF - if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle - - call wall_time(t1) - call write_double(6,(t1-t0),'Broadcast time') - - !--- - call omp_set_max_active_levels(8) - call davidson_slave_tcp(0) - call omp_set_max_active_levels(1) - print *, mpi_rank, ': Davidson done' - !--- - - IRP_IF MPI - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - print *, irp_here, 'error in barrier' - endif - IRP_ENDIF - call mpi_print('----------') - - else if (zmq_state(1:3) == 'pt2') then - - ! PT2 - ! --- - - IRP_IF MPI - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - print *, irp_here, 'error in barrier' - endif - IRP_ENDIF - call wall_time(t0) - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_psi') - IRP_ENDIF - if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_N_det_generators') - IRP_ENDIF - if (zmq_get_N_det_generators (zmq_to_qp_run_socket, 1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_N_det_selectors') - IRP_ENDIF - if (zmq_get_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_dvector threshold_generators') - IRP_ENDIF - if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_generators',(/threshold_generators/),1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_dvector energy') - IRP_ENDIF - if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_ivector pt2_stoch_istate') - IRP_ENDIF - if (zmq_get_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_dvector state_average_weight') - IRP_ENDIF - if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle - IRP_IF MPI_DEBUG - call mpi_print('zmq_get_dvector selection_weight') - IRP_ENDIF - if (zmq_get_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) cycle - pt2_e0_denominator(1:N_states) = energy(1:N_states) - SOFT_TOUCH pt2_e0_denominator state_average_weight pt2_stoch_istate threshold_generators selection_weight psi_det psi_coef N_det_generators N_det_selectors - - - call wall_time(t1) - call write_double(6,(t1-t0),'Broadcast time') - IRP_IF MPI - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - print *, irp_here, 'error in barrier' - endif - IRP_ENDIF - - - IRP_IF MPI_DEBUG - call mpi_print('Entering OpenMP section') - IRP_ENDIF - if (.true.) then - integer :: nproc_target, ii - double precision :: mem_collector, mem, rss - - call resident_memory(rss) - - nproc_target = nthreads_pt2 - ii = min(N_det, (elec_alpha_num*(mo_num-elec_alpha_num))**2) - - do - mem = rss + & ! - nproc_target * 8.d0 * & ! bytes - ( 0.5d0*pt2_n_tasks_max & ! task_id - + 64.d0*pt2_n_tasks_max & ! task - + 3.d0*pt2_n_tasks_max*N_states & ! pt2, variance, norm - + 1.d0*pt2_n_tasks_max & ! i_generator, subset - + 3.d0*(N_int*2.d0*ii+ ii) & ! selection buffer - + 1.d0*(N_int*2.d0*ii+ ii) & ! sort selection buffer - + 2.0d0*(ii) & ! preinteresting, interesting, - ! prefullinteresting, fullinteresting - + 2.0d0*(N_int*2*ii) & ! minilist, fullminilist - + 1.0d0*(N_states*mo_num*mo_num) & ! mat - ) / 1024.d0**3 - - if (nproc_target == 0) then - call check_mem(mem,irp_here) - nproc_target = 1 - exit - endif - - if (mem+rss < qp_max_mem) then - exit - endif - - nproc_target = nproc_target - 1 - - enddo - - if (N_det > 100000) then - - if (mpi_master) then - print *, 'N_det', N_det - print *, 'N_det_generators', N_det_generators - print *, 'N_det_selectors', N_det_selectors - print *, 'pt2_e0_denominator', pt2_e0_denominator - print *, 'pt2_stoch_istate', pt2_stoch_istate - print *, 'state_average_weight', state_average_weight - print *, 'selection_weight', selection_weight - print *, 'Number of threads', nproc_target - endif - - if (h0_type == 'CFG') then - PROVIDE det_to_configuration - endif - - PROVIDE global_selection_buffer pt2_N_teeth pt2_F N_det_generators - PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique - PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order - PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted_tc - - PROVIDE psi_det_hii selection_weight pseudo_sym pt2_min_parallel_tasks - - if (mpi_master) then - print *, 'Running PT2' - endif - !$OMP PARALLEL PRIVATE(i) NUM_THREADS(nproc_target+1) - i = omp_get_thread_num() - call run_pt2_slave(0,i,pt2_e0_denominator) - !$OMP END PARALLEL - FREE state_average_weight - print *, mpi_rank, ': PT2 done' - print *, '-------' - - endif - endif - - IRP_IF MPI - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - print *, irp_here, 'error in barrier' - endif - IRP_ENDIF - call mpi_print('----------') - - endif - - end do - IRP_IF MPI - call MPI_finalize(ierr) - IRP_ENDIF -end - - - diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index cf770049..446e8d87 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -1,10 +1,13 @@ subroutine run_cipsi - implicit none - use selection_types + BEGIN_DOC -! Selected Full Configuration Interaction with deterministic selection and -! stochastic PT2. + ! Selected Full Configuration Interaction with deterministic selection and + ! stochastic PT2. END_DOC + + use selection_types + + implicit none integer :: i,j,k type(pt2_type) :: pt2_data, pt2_data_err double precision, allocatable :: zeros(:) diff --git a/src/cipsi/energy.irp.f b/src/cipsi/energy.irp.f index 1f7cf122..4b496c11 100644 --- a/src/cipsi/energy.irp.f +++ b/src/cipsi/energy.irp.f @@ -36,12 +36,3 @@ BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ] endif END_PROVIDER - -BEGIN_PROVIDER [ double precision, pt2_overlap, (N_states, N_states) ] - implicit none - BEGIN_DOC - ! Overlap between the perturbed wave functions - END_DOC - pt2_overlap(1:N_states,1:N_states) = 0.d0 -END_PROVIDER - diff --git a/src/cipsi/lock_2rdm.irp.f b/src/cipsi/lock_2rdm.irp.f deleted file mode 100644 index e69de29b..00000000 diff --git a/src/cipsi/pt2_type.irp.f b/src/cipsi/pt2_type.irp.f deleted file mode 100644 index ee90d421..00000000 --- a/src/cipsi/pt2_type.irp.f +++ /dev/null @@ -1,128 +0,0 @@ -subroutine pt2_alloc(pt2_data,N) - implicit none - use selection_types - type(pt2_type), intent(inout) :: pt2_data - integer, intent(in) :: N - integer :: k - - allocate(pt2_data % pt2(N) & - ,pt2_data % variance(N) & - ,pt2_data % rpt2(N) & - ,pt2_data % overlap(N,N) & - ) - - pt2_data % pt2(:) = 0.d0 - pt2_data % variance(:) = 0.d0 - pt2_data % rpt2(:) = 0.d0 - pt2_data % overlap(:,:) = 0.d0 - -end subroutine - -subroutine pt2_dealloc(pt2_data) - implicit none - use selection_types - type(pt2_type), intent(inout) :: pt2_data - deallocate(pt2_data % pt2 & - ,pt2_data % variance & - ,pt2_data % rpt2 & - ,pt2_data % overlap & - ) -end subroutine - -subroutine pt2_add(p1, w, p2) - implicit none - use selection_types - BEGIN_DOC -! p1 += w * p2 - END_DOC - type(pt2_type), intent(inout) :: p1 - double precision, intent(in) :: w - type(pt2_type), intent(in) :: p2 - - if (w == 1.d0) then - - p1 % pt2(:) = p1 % pt2(:) + p2 % pt2(:) - p1 % rpt2(:) = p1 % rpt2(:) + p2 % rpt2(:) - p1 % variance(:) = p1 % variance(:) + p2 % variance(:) - p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:) - - else - - p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) - p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) - p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) - p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) - - endif - -end subroutine - - -subroutine pt2_add2(p1, w, p2) - implicit none - use selection_types - BEGIN_DOC -! p1 += w * p2**2 - END_DOC - type(pt2_type), intent(inout) :: p1 - double precision, intent(in) :: w - type(pt2_type), intent(in) :: p2 - - if (w == 1.d0) then - - p1 % pt2(:) = p1 % pt2(:) + p2 % pt2(:) * p2 % pt2(:) - p1 % rpt2(:) = p1 % rpt2(:) + p2 % rpt2(:) * p2 % rpt2(:) - p1 % variance(:) = p1 % variance(:) + p2 % variance(:) * p2 % variance(:) - p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:) * p2 % overlap(:,:) - - else - - p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) * p2 % pt2(:) - p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) * p2 % rpt2(:) - p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) * p2 % variance(:) - p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) * p2 % overlap(:,:) - - endif - -end subroutine - - -subroutine pt2_serialize(pt2_data, n, x) - implicit none - use selection_types - type(pt2_type), intent(in) :: pt2_data - integer, intent(in) :: n - double precision, intent(out) :: x(*) - - integer :: i,k,n2 - - n2 = n*n - x(1:n) = pt2_data % pt2(1:n) - k=n - x(k+1:k+n) = pt2_data % rpt2(1:n) - k=k+n - x(k+1:k+n) = pt2_data % variance(1:n) - k=k+n - x(k+1:k+n2) = reshape(pt2_data % overlap(1:n,1:n), (/ n2 /)) - -end - -subroutine pt2_deserialize(pt2_data, n, x) - implicit none - use selection_types - type(pt2_type), intent(inout) :: pt2_data - integer, intent(in) :: n - double precision, intent(in) :: x(*) - - integer :: i,k,n2 - - n2 = n*n - pt2_data % pt2(1:n) = x(1:n) - k=n - pt2_data % rpt2(1:n) = x(k+1:k+n) - k=k+n - pt2_data % variance(1:n) = x(k+1:k+n) - k=k+n - pt2_data % overlap(1:n,1:n) = reshape(x(k+1:k+n2), (/ n, n /)) - -end diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index 87ebca40..38a8f362 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -1,256 +1,5 @@ -subroutine run_selection_slave(thread,iproc,energy) - use f77_zmq - use selection_types - implicit none - - double precision, intent(in) :: energy(N_states) - integer, intent(in) :: thread, iproc - integer :: rc, i - - integer :: worker_id, task_id(1), ctask, ltask - character*(512) :: task - - integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - - integer(ZMQ_PTR), external :: new_zmq_push_socket - integer(ZMQ_PTR) :: zmq_socket_push - - type(selection_buffer) :: buf, buf2 - logical :: done, buffer_ready - type(pt2_type) :: pt2_data - - PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique - PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order - PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym - PROVIDE psi_selectors_coef_transp psi_det_sorted weight_selection - - call pt2_alloc(pt2_data,N_states) - - zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() - - integer, external :: connect_to_taskserver - if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - return - endif - - zmq_socket_push = new_zmq_push_socket(thread) - - buf%N = 0 - buffer_ready = .False. - ctask = 1 - - do - integer, external :: get_task_from_taskserver - if (get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task) == -1) then - exit - endif - done = task_id(ctask) == 0 - if (done) then - ctask = ctask - 1 - else - integer :: i_generator, N, subset, bsize - call sscanf_ddd(task, subset, i_generator, N) - if(buf%N == 0) then - ! Only first time - call create_selection_buffer(N, N*2, buf) - buffer_ready = .True. - else - if (N /= buf%N) then - print *, 'N=', N - print *, 'buf%N=', buf%N - print *, 'bug in ', irp_here - stop '-1' - end if - end if - call select_connected(i_generator, energy, pt2_data, buf, subset, pt2_F(i_generator)) - endif - - integer, external :: task_done_to_taskserver - - if(done .or. ctask == size(task_id)) then - do i=1, ctask - if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i)) == -1) then - call usleep(100) - if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i)) == -1) then - ctask = 0 - done = .true. - exit - endif - endif - end do - if(ctask > 0) then - call sort_selection_buffer(buf) -! call merge_selection_buffers(buf,buf2) - call push_selection_results(zmq_socket_push, pt2_data, buf, task_id(1), ctask) - call pt2_dealloc(pt2_data) - call pt2_alloc(pt2_data,N_states) -! buf%mini = buf2%mini - buf%cur = 0 - end if - ctask = 0 - end if - - if(done) exit - ctask = ctask + 1 - end do - - if(ctask > 0) then - call sort_selection_buffer(buf) -! call merge_selection_buffers(buf,buf2) - call push_selection_results(zmq_socket_push, pt2_data, buf, task_id(1), ctask) -! buf%mini = buf2%mini - buf%cur = 0 - end if - ctask = 0 - call pt2_dealloc(pt2_data) - - integer, external :: disconnect_from_taskserver - if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then - continue - endif - - call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) - call end_zmq_push_socket(zmq_socket_push,thread) - if (buffer_ready) then - call delete_selection_buffer(buf) -! call delete_selection_buffer(buf2) - endif -end subroutine - - -subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntasks) - use f77_zmq - use selection_types - implicit none - - integer(ZMQ_PTR), intent(in) :: zmq_socket_push - type(pt2_type), intent(in) :: pt2_data - type(selection_buffer), intent(inout) :: b - integer, intent(in) :: ntasks, task_id(*) - integer :: rc - double precision, allocatable :: pt2_serialized(:) - - rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE) - if(rc /= 4) then - print *, 'f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE)' - endif - - - allocate(pt2_serialized (pt2_type_size(N_states)) ) - call pt2_serialize(pt2_data,N_states,pt2_serialized) - - rc = f77_zmq_send( zmq_socket_push, pt2_serialized, size(pt2_serialized)*8, ZMQ_SNDMORE) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 3 - return - else if(rc /= size(pt2_serialized)*8) then - stop 'push' - endif - deallocate(pt2_serialized) - - if (b%cur > 0) then - - rc = f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE) - if(rc /= 8*b%cur) then - print *, 'f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE)' - endif - - rc = f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE) - if(rc /= bit_kind*N_int*2*b%cur) then - print *, 'f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE)' - endif - - endif - - rc = f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE) - if(rc /= 4) then - print *, 'f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE)' - endif - - rc = f77_zmq_send( zmq_socket_push, task_id(1), ntasks*4, 0) - if(rc /= 4*ntasks) then - print *, 'f77_zmq_send( zmq_socket_push, task_id(1), ntasks*4, 0)' - endif - -! Activate is zmq_socket_push is a REQ -IRP_IF ZMQ_PUSH -IRP_ELSE - character*(2) :: ok - rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0) - if ((rc /= 2).and.(ok(1:2) /= 'ok')) then - print *, irp_here//': error in receiving ok' - stop -1 - endif -IRP_ENDIF - -end subroutine - - -subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_id, ntasks) - use f77_zmq - use selection_types - implicit none - integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - type(pt2_type), intent(inout) :: pt2_data - double precision, intent(out) :: val(*) - integer(bit_kind), intent(out) :: det(N_int, 2, *) - integer, intent(out) :: N, ntasks, task_id(*) - integer :: rc, rn, i - double precision, allocatable :: pt2_serialized(:) - - rc = f77_zmq_recv( zmq_socket_pull, N, 4, 0) - if(rc /= 4) then - print *, 'f77_zmq_recv( zmq_socket_pull, N, 4, 0)' - endif - - allocate(pt2_serialized (pt2_type_size(N_states)) ) - rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized), 0) - if (rc == -1) then - ntasks = 1 - task_id(1) = 0 - else if(rc /= 8*size(pt2_serialized)) then - stop 'pull' - endif - - call pt2_deserialize(pt2_data,N_states,pt2_serialized) - deallocate(pt2_serialized) - - if (N>0) then - rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, 0) - if(rc /= 8*N) then - print *, 'f77_zmq_recv( zmq_socket_pull, val(1), 8*N, 0)' - endif - - rc = f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, 0) - if(rc /= bit_kind*N_int*2*N) then - print *, 'f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, 0)' - endif - endif - - rc = f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0) - if(rc /= 4) then - print *, 'f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0)' - endif - - rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntasks*4, 0) - if(rc /= 4*ntasks) then - print *, 'f77_zmq_recv( zmq_socket_pull, task_id(1), ntasks*4, 0)' - endif - -! Activate is zmq_socket_pull is a REP -IRP_IF ZMQ_PUSH -IRP_ELSE - rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0) - if (rc /= 2) then - print *, irp_here//': error in sending ok' - stop -1 - endif -IRP_ENDIF -end subroutine - - +subroutine provide_for_selection_slave + PROVIDE psi_det_sorted_order + PROVIDE psi_selectors_coef_transp psi_det_sorted +end diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index ae84f84e..50749272 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -141,12 +141,12 @@ double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2, Nint) end -subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,buf,subset,csubset) +subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, fock_diag_tmp, E0, pt2_data, buf, subset, csubset) use bitmasks use selection_types implicit none BEGIN_DOC -! WARNING /!\ : It is assumed that the generators and selectors are psi_det_sorted + ! WARNING /!\ : It is assumed that the generators and selectors are psi_det_sorted END_DOC integer, intent(in) :: i_generator, subset, csubset @@ -156,28 +156,35 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d type(pt2_type), intent(inout) :: pt2_data type(selection_buffer), intent(inout) :: buf - integer :: h1,h2,s1,s2,s3,i1,i2,ib,sp,k,i,j,nt,ii,sze - integer(bit_kind) :: hole(N_int,2), particle(N_int,2), mask(N_int, 2), pmask(N_int, 2) - logical :: fullMatch, ok + integer :: h1, h2, s1, s2, s3, i1, i2, ib, sp, k, i, j, nt, ii, sze + integer :: maskInd + integer :: N_holes(2), N_particles(2) + integer :: hole_list(N_int*bit_kind_size,2) + integer :: particle_list(N_int*bit_kind_size,2) + integer :: l_a, nmax, idx + integer :: nb_count, maskInd_save + integer(bit_kind) :: hole(N_int,2), particle(N_int,2), mask(N_int, 2), pmask(N_int, 2) + integer(bit_kind) :: mobMask(N_int, 2), negMask(N_int, 2) + logical :: fullMatch, ok + logical :: monoAdo, monoBdo + logical :: monoBdo_save + logical :: found - integer(bit_kind) :: mobMask(N_int, 2), negMask(N_int, 2) - integer,allocatable :: preinteresting(:), prefullinteresting(:) - integer,allocatable :: interesting(:), fullinteresting(:) - integer,allocatable :: tmp_array(:) - integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) - logical, allocatable :: banned(:,:,:), bannedOrb(:,:) - double precision, allocatable :: coef_fullminilist_rev(:,:) + integer, allocatable :: preinteresting(:), prefullinteresting(:) + integer, allocatable :: interesting(:), fullinteresting(:) + integer, allocatable :: tmp_array(:) + integer, allocatable :: indices(:), exc_degree(:), iorder(:) + integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) + logical, allocatable :: banned(:,:,:), bannedOrb(:,:) + double precision, allocatable :: coef_fullminilist_rev(:,:) + double precision, allocatable :: mat(:,:,:) - double precision, allocatable :: mat(:,:,:) - - logical :: monoAdo, monoBdo - integer :: maskInd PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique - PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order + PROVIDE psi_bilinear_matrix_rows psi_bilinear_matrix_order psi_bilinear_matrix_transp_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp + PROVIDE psi_selectors_coef_transp psi_det_sorted_order PROVIDE banned_excitation monoAdo = .true. @@ -192,17 +199,9 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), particle_mask(k,2)) enddo - - integer :: N_holes(2), N_particles(2) - integer :: hole_list(N_int*bit_kind_size,2) - integer :: particle_list(N_int*bit_kind_size,2) - call bitstring_to_list_ab(hole , hole_list , N_holes , N_int) call bitstring_to_list_ab(particle, particle_list, N_particles, N_int) - integer :: l_a, nmax, idx - integer, allocatable :: indices(:), exc_degree(:), iorder(:) - ! Removed to avoid introducing determinants already presents in the wf !double precision, parameter :: norm_thr = 1.d-16 @@ -320,22 +319,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d allocate(banned(mo_num, mo_num,2), bannedOrb(mo_num, 2)) - allocate (mat(N_states, mo_num, mo_num)) + allocate(mat(N_states, mo_num, mo_num)) maskInd = -1 - integer :: nb_count, maskInd_save - logical :: monoBdo_save - logical :: found - do s1=1,2 - do i1=N_holes(s1),1,-1 ! Generate low excitations first + do s1 = 1, 2 + do i1 = N_holes(s1), 1, -1 ! Generate low excitations first found = .False. monoBdo_save = monoBdo maskInd_save = maskInd - do s2=s1,2 + do s2 = s1, 2 ib = 1 if(s1 == s2) ib = i1+1 - do i2=N_holes(s2),ib,-1 + do i2 = N_holes(s2), ib, -1 maskInd = maskInd + 1 if(mod(maskInd, csubset) == (subset-1)) then found = .True. @@ -349,14 +345,14 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d maskInd = maskInd_save h1 = hole_list(i1,s1) - call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int) + call apply_hole(psi_det_generators(1,1,i_generator), s1, h1, pmask, ok, N_int) negMask = not(pmask) interesting(0) = 0 fullinteresting(0) = 0 - do ii=1,preinteresting(0) + do ii = 1, preinteresting(0) i = preinteresting(ii) select case (N_int) case (1) @@ -372,7 +368,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d mobMask(1:3,1) = iand(negMask(1:3,1), psi_det_sorted(1:3,1,i)) mobMask(1:3,2) = iand(negMask(1:3,2), psi_det_sorted(1:3,2,i)) nt = 0 - do j=3,1,-1 + do j = 3, 1, -1 if (mobMask(j,1) /= 0_bit_kind) then nt = nt+ popcnt(mobMask(j, 1)) if (nt > 4) exit @@ -386,7 +382,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d mobMask(1:4,1) = iand(negMask(1:4,1), psi_det_sorted(1:4,1,i)) mobMask(1:4,2) = iand(negMask(1:4,2), psi_det_sorted(1:4,2,i)) nt = 0 - do j=4,1,-1 + do j = 4, 1, -1 if (mobMask(j,1) /= 0_bit_kind) then nt = nt+ popcnt(mobMask(j, 1)) if (nt > 4) exit @@ -400,7 +396,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d mobMask(1:N_int,1) = iand(negMask(1:N_int,1), psi_det_sorted(1:N_int,1,i)) mobMask(1:N_int,2) = iand(negMask(1:N_int,2), psi_det_sorted(1:N_int,2,i)) nt = 0 - do j=N_int,1,-1 + do j = N_int, 1, -1 if (mobMask(j,1) /= 0_bit_kind) then nt = nt+ popcnt(mobMask(j, 1)) if (nt > 4) exit @@ -441,7 +437,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d end do - do ii=1,prefullinteresting(0) + do ii = 1, prefullinteresting(0) i = prefullinteresting(ii) nt = 0 mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i)) @@ -480,40 +476,38 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d minilist(:,:,i) = psi_det_sorted(:,:,interesting(i)) enddo - do s2=s1,2 + do s2 = s1, 2 sp = s1 - if(s1 /= s2) then - sp = 3 - endif + if(s1 /= s2) sp = 3 ib = 1 if(s1 == s2) ib = i1+1 monoAdo = .true. - do i2=N_holes(s2),ib,-1 ! Generate low excitations first + do i2 = N_holes(s2), ib, -1 ! Generate low excitations first h2 = hole_list(i2,s2) call apply_hole(pmask, s2,h2, mask, ok, N_int) banned(:,:,1) = banned_excitation(:,:) banned(:,:,2) = banned_excitation(:,:) - do j=1,mo_num + do j = 1, mo_num bannedOrb(j, 1) = .true. bannedOrb(j, 2) = .true. enddo - do s3=1,2 - do i=1,N_particles(s3) + do s3 = 1, 2 + do i = 1, N_particles(s3) bannedOrb(particle_list(i,s3), s3) = .false. enddo enddo if(s1 /= s2) then if(monoBdo) then bannedOrb(h1,s1) = .false. - end if + endif if(monoAdo) then bannedOrb(h2,s2) = .false. monoAdo = .false. - end if - end if + endif + endif maskInd = maskInd + 1 if(mod(maskInd, csubset) == (subset-1)) then @@ -522,12 +516,18 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d if(fullMatch) cycle call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) + call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) end if + + enddo + if(s1 /= s2) monoBdo = .false. enddo - deallocate(fullminilist,minilist) + + deallocate(fullminilist, minilist) + enddo enddo deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) diff --git a/src/cipsi/selection_types.f90 b/src/cipsi/selection_types.f90 deleted file mode 100644 index 58ce0e03..00000000 --- a/src/cipsi/selection_types.f90 +++ /dev/null @@ -1,25 +0,0 @@ -module selection_types - type selection_buffer - integer :: N, cur - integer(8) , pointer :: det(:,:,:) - double precision, pointer :: val(:) - double precision :: mini - endtype - - type pt2_type - double precision, allocatable :: pt2(:) - double precision, allocatable :: rpt2(:) - double precision, allocatable :: variance(:) - double precision, allocatable :: overlap(:,:) - endtype - - contains - - integer function pt2_type_size(N) - implicit none - integer, intent(in) :: N - pt2_type_size = (3*n + n*n) - end function - -end module - diff --git a/src/cipsi_utils/README.rst b/src/cipsi_utils/README.rst new file mode 100644 index 00000000..8e98e3ac --- /dev/null +++ b/src/cipsi_utils/README.rst @@ -0,0 +1,5 @@ +=========== +cipsi_utils +=========== + +Common functions for CIPSI and TC-CIPSI diff --git a/src/cipsi/environment.irp.f b/src/cipsi_utils/environment.irp.f similarity index 100% rename from src/cipsi/environment.irp.f rename to src/cipsi_utils/environment.irp.f diff --git a/src/cipsi_utils/pt2_stoch_routines.irp.f b/src/cipsi_utils/pt2_stoch_routines.irp.f new file mode 100644 index 00000000..f067d0be --- /dev/null +++ b/src/cipsi_utils/pt2_stoch_routines.irp.f @@ -0,0 +1,891 @@ +BEGIN_PROVIDER [ integer, pt2_stoch_istate ] + implicit none + BEGIN_DOC + ! State for stochatsic PT2 + END_DOC + pt2_stoch_istate = 1 +END_PROVIDER + + BEGIN_PROVIDER [ integer, pt2_F, (N_det_generators) ] +&BEGIN_PROVIDER [ integer, pt2_n_tasks_max ] + implicit none + logical, external :: testTeethBuilding + integer :: i,j + pt2_n_tasks_max = elec_alpha_num*elec_alpha_num + elec_alpha_num*elec_beta_num - n_core_orb*2 + pt2_n_tasks_max = min(pt2_n_tasks_max,1+N_det_generators/10000) + call write_int(6,pt2_n_tasks_max,'pt2_n_tasks_max') + + pt2_F(:) = max(int(sqrt(float(pt2_n_tasks_max))),1) + do i=1,pt2_n_0(1+pt2_N_teeth/4) + pt2_F(i) = pt2_n_tasks_max*pt2_min_parallel_tasks + enddo + do i=1+pt2_n_0(pt2_N_teeth-pt2_N_teeth/4), pt2_n_0(pt2_N_teeth-pt2_N_teeth/10) + pt2_F(i) = pt2_min_parallel_tasks + enddo + do i=1+pt2_n_0(pt2_N_teeth-pt2_N_teeth/10), N_det_generators + pt2_F(i) = 1 + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ integer, pt2_N_teeth ] +&BEGIN_PROVIDER [ integer, pt2_minDetInFirstTeeth ] + implicit none + logical, external :: testTeethBuilding + + if(N_det_generators < 1024) then + pt2_minDetInFirstTeeth = 1 + pt2_N_teeth = 1 + else + pt2_minDetInFirstTeeth = min(5, N_det_generators) + do pt2_N_teeth=100,2,-1 + if(testTeethBuilding(pt2_minDetInFirstTeeth, pt2_N_teeth)) exit + end do + end if + call write_int(6,pt2_N_teeth,'Number of comb teeth') +END_PROVIDER + + +logical function testTeethBuilding(minF, N) + implicit none + integer, intent(in) :: minF, N + integer :: n0, i + double precision :: u0, Wt, r + + double precision, allocatable :: tilde_w(:), tilde_cW(:) + integer, external :: dress_find_sample + + double precision :: rss + double precision, external :: memory_of_double, memory_of_int + + rss = memory_of_double(2*N_det_generators+1) + call check_mem(rss,irp_here) + + allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators)) + + double precision :: norm2 + norm2 = 0.d0 + do i=N_det_generators,1,-1 + tilde_w(i) = psi_coef_sorted_gen(i,pt2_stoch_istate) * & + psi_coef_sorted_gen(i,pt2_stoch_istate) + norm2 = norm2 + tilde_w(i) + enddo + + f = 1.d0/norm2 + tilde_w(:) = tilde_w(:) * f + + tilde_cW(0) = -1.d0 + do i=1,N_det_generators + tilde_cW(i) = tilde_cW(i-1) + tilde_w(i) + enddo + tilde_cW(:) = tilde_cW(:) + 1.d0 + deallocate(tilde_w) + + n0 = 0 + testTeethBuilding = .false. + double precision :: f + integer :: minFN + minFN = N_det_generators - minF * N + f = 1.d0/dble(N) + do + u0 = tilde_cW(n0) + r = tilde_cW(n0 + minF) + Wt = (1d0 - u0) * f + if (dabs(Wt) <= 1.d-3) then + exit + endif + if(Wt >= r - u0) then + testTeethBuilding = .true. + exit + end if + n0 += 1 + if(n0 > minFN) then + exit + end if + end do + deallocate(tilde_cW) + +end function + + +!subroutine provide_for_zmq_pt2 +! PROVIDE psi_det_sorted_order psi_selectors_coef_transp psi_det_sorted +!end + +subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) + use f77_zmq + use selection_types + + implicit none + + integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull + integer, intent(in) :: N_in + double precision, intent(in) :: relative_error, E(N_states) + type(pt2_type), intent(inout) :: pt2_data, pt2_data_err +! + integer :: i, N + + double precision :: state_average_weight_save(N_states), w(N_states,4) + integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket + type(selection_buffer) :: b + + PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique + PROVIDE psi_bilinear_matrix_rows psi_bilinear_matrix_order + PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns + PROVIDE psi_bilinear_matrix_transp_order + PROVIDE psi_det_hii selection_weight pseudo_sym + PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max + + call provide_for_zmq_pt2 + + if (h0_type == 'CFG') then + PROVIDE psi_configuration_hii det_to_configuration + endif + + if (N_det <= max(4,N_states) .or. pt2_N_teeth < 2) then + call ZMQ_selection(N_in, pt2_data) + else + + N = max(N_in,1) * N_states + state_average_weight_save(:) = state_average_weight(:) + if (int(N,8)*2_8 > huge(1)) then + print *, irp_here, ': integer too large' + stop -1 + endif + call create_selection_buffer(N, N*2, b) + ASSERT (associated(b%det)) + ASSERT (associated(b%val)) + + do pt2_stoch_istate=1,N_states + state_average_weight(:) = 0.d0 + state_average_weight(pt2_stoch_istate) = 1.d0 + 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 + call new_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') + + integer, external :: zmq_put_psi + integer, external :: zmq_put_N_det_generators + integer, external :: zmq_put_N_det_selectors + integer, external :: zmq_put_dvector + integer, external :: zmq_put_ivector + if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then + stop 'Unable to put psi on ZMQ server' + endif + if (zmq_put_N_det_generators(zmq_to_qp_run_socket, 1) == -1) then + stop 'Unable to put N_det_generators on ZMQ server' + endif + if (zmq_put_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) then + stop 'Unable to put N_det_selectors on ZMQ server' + endif + if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',pt2_e0_denominator,size(pt2_e0_denominator)) == -1) then + stop 'Unable to put energy on ZMQ server' + endif + if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then + stop 'Unable to put state_average_weight on ZMQ server' + endif + if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then + stop 'Unable to put selection_weight on ZMQ server' + endif + if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then + stop 'Unable to put pt2_stoch_istate on ZMQ server' + endif + if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',(/threshold_generators/),1) == -1) then + stop 'Unable to put threshold_generators on ZMQ server' + endif + + + integer, external :: add_task_to_taskserver + character(300000) :: task + + integer :: j,k,ipos,ifirst + ifirst=0 + + ipos=0 + do i=1,N_det_generators + if (pt2_F(i) > 1) then + ipos += 1 + endif + enddo + call write_int(6,sum(pt2_F),'Number of tasks') + call write_int(6,ipos,'Number of fragmented tasks') + + ipos=1 + do i= 1, N_det_generators + do j=1,pt2_F(pt2_J(i)) + write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') j, pt2_J(i), N_in + ipos += 30 + if (ipos > 300000-30) then + if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then + stop 'Unable to add task to task server' + endif + ipos=1 + if (ifirst == 0) then + ifirst=1 + if (zmq_set_running(zmq_to_qp_run_socket) == -1) then + print *, irp_here, ': Failed in zmq_set_running' + endif + endif + endif + end do + enddo + if (ipos > 1) then + if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then + stop 'Unable to add task to task server' + endif + endif + + integer, external :: zmq_set_running + if (zmq_set_running(zmq_to_qp_run_socket) == -1) then + print *, irp_here, ': Failed in zmq_set_running' + endif + + + double precision :: mem_collector, mem, rss + + call resident_memory(rss) + + mem_collector = 8.d0 * & ! bytes + ( 1.d0*pt2_n_tasks_max & ! task_id, index + + 0.635d0*N_det_generators & ! f,d + + pt2_n_tasks_max*pt2_type_size(N_states) & ! pt2_data_task + + N_det_generators*pt2_type_size(N_states) & ! pt2_data_I + + 4.d0*(pt2_N_teeth+1) & ! S, S2, T2, T3 + + 1.d0*(N_int*2.d0*N + N) & ! selection buffer + + 1.d0*(N_int*2.d0*N + N) & ! sort selection buffer + ) / 1024.d0**3 + + integer :: nproc_target, ii + nproc_target = nthreads_pt2 + ii = min(N_det, (elec_alpha_num*(mo_num-elec_alpha_num))**2) + + do + mem = mem_collector + & ! + nproc_target * 8.d0 * & ! bytes + ( 0.5d0*pt2_n_tasks_max & ! task_id + + 64.d0*pt2_n_tasks_max & ! task + + pt2_type_size(N_states)*pt2_n_tasks_max*N_states & ! pt2, variance, overlap + + 1.d0*pt2_n_tasks_max & ! i_generator, subset + + 1.d0*(N_int*2.d0*ii+ ii) & ! selection buffer + + 1.d0*(N_int*2.d0*ii+ ii) & ! sort selection buffer + + 2.0d0*(ii) & ! preinteresting, interesting, + ! prefullinteresting, fullinteresting + + 2.0d0*(N_int*2*ii) & ! minilist, fullminilist + + 1.0d0*(N_states*mo_num*mo_num) & ! mat + ) / 1024.d0**3 + + if (nproc_target == 0) then + call check_mem(mem,irp_here) + nproc_target = 1 + exit + endif + + if (mem+rss < qp_max_mem) then + exit + endif + + nproc_target = nproc_target - 1 + + enddo + call write_int(6,nproc_target,'Number of threads for PT2') + call write_double(6,mem,'Memory (Gb)') + + call set_multiple_levels_omp(.False.) + + + print '(A)', '========== ==================== ================ ================ ================ ============= ===========' + print '(A)', ' Samples Energy PT2 Variance Norm^2 Convergence Seconds' + print '(A)', '========== ==================== ================ ================ ================ ============= ===========' + + PROVIDE global_selection_buffer + + !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc_target+1) & + !$OMP PRIVATE(i) + i = omp_get_thread_num() + if (i==0) then + + call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, pt2_data, pt2_data_err, b, N) + pt2_data % rpt2(pt2_stoch_istate) = & + pt2_data % pt2(pt2_stoch_istate)/(1.d0+pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate)) + + !TODO : We should use here the correct formula for the error of X/Y + pt2_data_err % rpt2(pt2_stoch_istate) = & + pt2_data_err % pt2(pt2_stoch_istate)/(1.d0 + pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate)) + + else + call pt2_slave_inproc(i) + endif + !$OMP END PARALLEL + call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') + call set_multiple_levels_omp(.True.) + + print '(A)', '========== ==================== ================ ================ ================ ============= ===========' + + + do k=1,N_states + pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) + enddo + SOFT_TOUCH pt2_overlap + + enddo + FREE pt2_stoch_istate + + ! Symmetrize overlap + do j=2,N_states + do i=1,j-1 + pt2_overlap(i,j) = 0.5d0 * (pt2_overlap(i,j) + pt2_overlap(j,i)) + pt2_overlap(j,i) = pt2_overlap(i,j) + enddo + enddo + + print *, 'Overlap of perturbed states:' + do k=1,N_states + print *, pt2_overlap(k,:) + enddo + print *, '-------' + + if (N_in > 0) then + b%cur = min(N_in,b%cur) + if (s2_eig) then + call make_selection_buffer_s2(b) + else + call remove_duplicates_in_selection_buffer(b) + endif + call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) + endif + call delete_selection_buffer(b) + + state_average_weight(:) = state_average_weight_save(:) + TOUCH state_average_weight + call update_pt2_and_variance_weights(pt2_data, N_states) + endif + + +end subroutine + + +subroutine pt2_slave_inproc(i) + implicit none + integer, intent(in) :: i + + PROVIDE global_selection_buffer + call run_pt2_slave(1,i,pt2_e0_denominator) +end + + +subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_err, b, N_) + use f77_zmq + use selection_types + use bitmasks + implicit none + + + integer(ZMQ_PTR), intent(in) :: zmq_socket_pull + double precision, intent(in) :: relative_error, E + type(pt2_type), intent(inout) :: pt2_data, pt2_data_err + type(selection_buffer), intent(inout) :: b + integer, intent(in) :: N_ + + type(pt2_type), allocatable :: pt2_data_task(:) + type(pt2_type), allocatable :: pt2_data_I(:) + type(pt2_type), allocatable :: pt2_data_S(:) + type(pt2_type), allocatable :: pt2_data_S2(:) + type(pt2_type) :: pt2_data_teeth + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + integer, external :: zmq_delete_tasks_async_send + integer, external :: zmq_delete_tasks_async_recv + integer, external :: zmq_abort + integer, external :: pt2_find_sample_lr + + PROVIDE pt2_stoch_istate + + integer :: more, n, i, p, c, t, n_tasks, U + integer, allocatable :: task_id(:) + integer, allocatable :: index(:) + + double precision :: v, x, x2, x3, avg, avg2, avg3(N_states), eqt, E0, v0, n0(N_states) + double precision :: eqta(N_states) + double precision :: time, time1, time0 + + integer, allocatable :: f(:) + logical, allocatable :: d(:) + logical :: do_exit, stop_now, sending + logical, external :: qp_stop + type(selection_buffer) :: b2 + + + double precision :: rss + double precision, external :: memory_of_double, memory_of_int + + character(len=20) :: format_str1, str_error1, format_str2, str_error2 + character(len=20) :: format_str3, str_error3, format_str4, str_error4 + character(len=20) :: format_value1, format_value2, format_value3, format_value4 + character(len=20) :: str_value1, str_value2, str_value3, str_value4 + character(len=20) :: str_conv + double precision :: value1, value2, value3, value4 + double precision :: error1, error2, error3, error4 + integer :: size1,size2,size3,size4 + + double precision :: conv_crit + + sending =.False. + + rss = memory_of_int(pt2_n_tasks_max*2+N_det_generators*2) + rss += memory_of_double(N_states*N_det_generators)*3.d0 + rss += memory_of_double(N_states*pt2_n_tasks_max)*3.d0 + rss += memory_of_double(pt2_N_teeth+1)*4.d0 + call check_mem(rss,irp_here) + + ! If an allocation is added here, the estimate of the memory should also be + ! updated in ZMQ_pt2 + allocate(task_id(pt2_n_tasks_max), index(pt2_n_tasks_max), f(N_det_generators)) + allocate(d(N_det_generators+1)) + allocate(pt2_data_task(pt2_n_tasks_max)) + allocate(pt2_data_I(N_det_generators)) + allocate(pt2_data_S(pt2_N_teeth+1)) + allocate(pt2_data_S2(pt2_N_teeth+1)) + + + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + call create_selection_buffer(N_, N_*2, b2) + + + pt2_data % pt2(pt2_stoch_istate) = -huge(1.) + pt2_data_err % pt2(pt2_stoch_istate) = huge(1.) + pt2_data % variance(pt2_stoch_istate) = huge(1.) + pt2_data_err % variance(pt2_stoch_istate) = huge(1.) + pt2_data % overlap(:,pt2_stoch_istate) = 0.d0 + pt2_data_err % overlap(:,pt2_stoch_istate) = huge(1.) + n = 1 + t = 0 + U = 0 + do i=1,pt2_n_tasks_max + call pt2_alloc(pt2_data_task(i),N_states) + enddo + do i=1,pt2_N_teeth+1 + call pt2_alloc(pt2_data_S(i),N_states) + call pt2_alloc(pt2_data_S2(i),N_states) + enddo + do i=1,N_det_generators + call pt2_alloc(pt2_data_I(i),N_states) + enddo + f(:) = pt2_F(:) + d(:) = .false. + n_tasks = 0 + E0 = E + v0 = 0.d0 + n0(:) = 0.d0 + more = 1 + call wall_time(time0) + time1 = time0 + + do_exit = .false. + stop_now = .false. + do while (n <= N_det_generators) + if(f(pt2_J(n)) == 0) then + d(pt2_J(n)) = .true. + do while(d(U+1)) + U += 1 + end do + + ! Deterministic part + do while(t <= pt2_N_teeth) + if(U >= pt2_n_0(t+1)) then + t=t+1 + E0 = 0.d0 + v0 = 0.d0 + n0(:) = 0.d0 + do i=pt2_n_0(t),1,-1 + E0 += pt2_data_I(i) % pt2(pt2_stoch_istate) + v0 += pt2_data_I(i) % variance(pt2_stoch_istate) + n0(:) += pt2_data_I(i) % overlap(:,pt2_stoch_istate) + end do + else + exit + end if + end do + + ! Add Stochastic part + c = pt2_R(n) + if(c > 0) then + + call pt2_alloc(pt2_data_teeth,N_states) + do p=pt2_N_teeth, 1, -1 + v = pt2_u_0 + pt2_W_T * (pt2_u(c) + dble(p-1)) + i = pt2_find_sample_lr(v, pt2_cW,pt2_n_0(p),pt2_n_0(p+1)) + v = pt2_W_T / pt2_w(i) + call pt2_add ( pt2_data_teeth, v, pt2_data_I(i) ) + call pt2_add ( pt2_data_S(p), 1.d0, pt2_data_teeth ) + call pt2_add2( pt2_data_S2(p), 1.d0, pt2_data_teeth ) + enddo + call pt2_dealloc(pt2_data_teeth) + + avg = E0 + pt2_data_S(t) % pt2(pt2_stoch_istate) / dble(c) + avg2 = v0 + pt2_data_S(t) % variance(pt2_stoch_istate) / dble(c) + avg3(:) = n0(:) + pt2_data_S(t) % overlap(:,pt2_stoch_istate) / dble(c) + if ((avg /= 0.d0) .or. (n == N_det_generators) ) then + do_exit = .true. + endif + if (qp_stop()) then + stop_now = .True. + endif + pt2_data % pt2(pt2_stoch_istate) = avg + pt2_data % variance(pt2_stoch_istate) = avg2 + pt2_data % overlap(:,pt2_stoch_istate) = avg3(:) + call wall_time(time) + ! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969) + if(c > 2) then + eqt = dabs((pt2_data_S2(t) % pt2(pt2_stoch_istate) / c) - (pt2_data_S(t) % pt2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability + eqt = sqrt(eqt / (dble(c) - 1.5d0)) + pt2_data_err % pt2(pt2_stoch_istate) = eqt + + eqt = dabs((pt2_data_S2(t) % variance(pt2_stoch_istate) / c) - (pt2_data_S(t) % variance(pt2_stoch_istate)/c)**2) ! dabs for numerical stability + eqt = sqrt(eqt / (dble(c) - 1.5d0)) + pt2_data_err % variance(pt2_stoch_istate) = eqt + + eqta(:) = dabs((pt2_data_S2(t) % overlap(:,pt2_stoch_istate) / c) - (pt2_data_S(t) % overlap(:,pt2_stoch_istate)/c)**2) ! dabs for numerical stability + eqta(:) = sqrt(eqta(:) / (dble(c) - 1.5d0)) + pt2_data_err % overlap(:,pt2_stoch_istate) = eqta(:) + + + if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then + time1 = time + print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.4)', c, & + pt2_data % pt2(pt2_stoch_istate) +E, & + pt2_data_err % pt2(pt2_stoch_istate), & + pt2_data % variance(pt2_stoch_istate), & + pt2_data_err % variance(pt2_stoch_istate), & + pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate), & + pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), & + time-time0 + if (stop_now .or. ( & + (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & + (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then + if (zmq_abort(zmq_to_qp_run_socket) == -1) then + call sleep(10) + if (zmq_abort(zmq_to_qp_run_socket) == -1) then + print *, irp_here, ': Error in sending abort signal (2)' + endif + endif + endif + endif + endif + end if + n += 1 + else if(more == 0) then + exit + else + call pull_pt2_results(zmq_socket_pull, index, pt2_data_task, task_id, n_tasks, b2) + if(n_tasks > pt2_n_tasks_max)then + print*,'PB !!!' + print*,'If you see this, send a bug report with the following content' + print*,irp_here + print*,'n_tasks,pt2_n_tasks_max = ',n_tasks,pt2_n_tasks_max + stop -1 + endif + if (zmq_delete_tasks_async_send(zmq_to_qp_run_socket,task_id,n_tasks,sending) == -1) then + stop 'PT2: Unable to delete tasks (send)' + endif + do i=1,n_tasks + if(index(i).gt.size(pt2_data_I,1).or.index(i).lt.1)then + print*,'PB !!!' + print*,'If you see this, send a bug report with the following content' + print*,irp_here + print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1) + stop -1 + endif + call pt2_add(pt2_data_I(index(i)),1.d0,pt2_data_task(i)) + f(index(i)) -= 1 + end do + do i=1, b2%cur + ! We assume the pulled buffer is sorted + if (b2%val(i) > b%mini) exit + call add_to_selection_buffer(b, b2%det(1,1,i), b2%val(i)) + end do + if (zmq_delete_tasks_async_recv(zmq_to_qp_run_socket,more,sending) == -1) then + stop 'PT2: Unable to delete tasks (recv)' + endif + end if + end do + do i=1,N_det_generators + call pt2_dealloc(pt2_data_I(i)) + enddo + do i=1,pt2_N_teeth+1 + call pt2_dealloc(pt2_data_S(i)) + call pt2_dealloc(pt2_data_S2(i)) + enddo + do i=1,pt2_n_tasks_max + call pt2_dealloc(pt2_data_task(i)) + enddo +!print *, 'deleting b2' + call delete_selection_buffer(b2) +!print *, 'sorting b' + call sort_selection_buffer(b) +!print *, 'done' + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + +end subroutine + + +integer function pt2_find_sample(v, w) + implicit none + double precision, intent(in) :: v, w(0:N_det_generators) + integer, external :: pt2_find_sample_lr + + pt2_find_sample = pt2_find_sample_lr(v, w, 0, N_det_generators) +end function + + +integer function pt2_find_sample_lr(v, w, l_in, r_in) + implicit none + double precision, intent(in) :: v, w(0:N_det_generators) + integer, intent(in) :: l_in,r_in + integer :: i,l,r + + l=l_in + r=r_in + + do while(r-l > 1) + i = shiftr(r+l,1) + if(w(i) < v) then + l = i + else + r = i + end if + end do + i = r + do r=i+1,N_det_generators + if (w(r) /= w(i)) then + exit + endif + enddo + pt2_find_sample_lr = r-1 +end function + + +BEGIN_PROVIDER [ integer, pt2_n_tasks ] + implicit none + BEGIN_DOC + ! Number of parallel tasks for the Monte Carlo + END_DOC + pt2_n_tasks = N_det_generators +END_PROVIDER + +BEGIN_PROVIDER[ double precision, pt2_u, (N_det_generators)] + implicit none + integer, allocatable :: seed(:) + integer :: m,i + call random_seed(size=m) + allocate(seed(m)) + do i=1,m + seed(i) = i + enddo + call random_seed(put=seed) + deallocate(seed) + + call RANDOM_NUMBER(pt2_u) + END_PROVIDER + + BEGIN_PROVIDER[ integer, pt2_J, (N_det_generators)] +&BEGIN_PROVIDER[ integer, pt2_R, (N_det_generators)] + implicit none + BEGIN_DOC +! pt2_J contains the list of generators after ordering them according to the +! Monte Carlo sampling. +! +! pt2_R(i) is the number of combs drawn when determinant i is computed. + END_DOC + integer :: N_c, N_j + integer :: U, t, i + double precision :: v + integer, external :: pt2_find_sample_lr + + logical, allocatable :: pt2_d(:) + integer :: m,l,r,k + integer :: ncache + integer, allocatable :: ii(:,:) + double precision :: dt + + ncache = min(N_det_generators,10000) + + double precision :: rss + double precision, external :: memory_of_double, memory_of_int + rss = memory_of_int(ncache)*dble(pt2_N_teeth) + memory_of_int(N_det_generators) + call check_mem(rss,irp_here) + + allocate(ii(pt2_N_teeth,ncache),pt2_d(N_det_generators)) + + pt2_R(:) = 0 + pt2_d(:) = .false. + N_c = 0 + N_j = pt2_n_0(1) + do i=1,N_j + pt2_d(i) = .true. + pt2_J(i) = i + end do + + U = 0 + do while(N_j < pt2_n_tasks) + + if (N_c+ncache > N_det_generators) then + ncache = N_det_generators - N_c + endif + + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(dt,v,t,k) + do k=1, ncache + dt = pt2_u_0 + do t=1, pt2_N_teeth + v = dt + pt2_W_T *pt2_u(N_c+k) + dt = dt + pt2_W_T + ii(t,k) = pt2_find_sample_lr(v, pt2_cW,pt2_n_0(t),pt2_n_0(t+1)) + end do + enddo + !$OMP END PARALLEL DO + + do k=1,ncache + !ADD_COMB + N_c = N_c+1 + do t=1, pt2_N_teeth + i = ii(t,k) + if(.not. pt2_d(i)) then + N_j += 1 + pt2_J(N_j) = i + pt2_d(i) = .true. + end if + end do + + pt2_R(N_j) = N_c + + !FILL_TOOTH + do while(U < N_det_generators) + U += 1 + if(.not. pt2_d(U)) then + N_j += 1 + pt2_J(N_j) = U + pt2_d(U) = .true. + exit + end if + end do + if (N_j >= pt2_n_tasks) exit + end do + enddo + + if(N_det_generators > 1) then + pt2_R(N_det_generators-1) = 0 + pt2_R(N_det_generators) = N_c + end if + + deallocate(ii,pt2_d) + +END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, pt2_w, (N_det_generators) ] +&BEGIN_PROVIDER [ double precision, pt2_cW, (0:N_det_generators) ] +&BEGIN_PROVIDER [ double precision, pt2_W_T ] +&BEGIN_PROVIDER [ double precision, pt2_u_0 ] +&BEGIN_PROVIDER [ integer, pt2_n_0, (pt2_N_teeth+1) ] + implicit none + integer :: i, t + double precision, allocatable :: tilde_w(:), tilde_cW(:) + double precision :: r, tooth_width + integer, external :: pt2_find_sample + + double precision :: rss + double precision, external :: memory_of_double, memory_of_int + rss = memory_of_double(2*N_det_generators+1) + call check_mem(rss,irp_here) + + if (N_det_generators == 1) then + + pt2_w(1) = 1.d0 + pt2_cw(1) = 1.d0 + pt2_u_0 = 1.d0 + pt2_W_T = 0.d0 + pt2_n_0(1) = 0 + pt2_n_0(2) = 1 + + else + + allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators)) + + tilde_cW(0) = 0d0 + + do i=1,N_det_generators + tilde_w(i) = psi_coef_sorted_gen(i,pt2_stoch_istate)**2 !+ 1.d-20 + enddo + + double precision :: norm2 + norm2 = 0.d0 + do i=N_det_generators,1,-1 + norm2 += tilde_w(i) + enddo + + tilde_w(:) = tilde_w(:) / norm2 + + tilde_cW(0) = -1.d0 + do i=1,N_det_generators + tilde_cW(i) = tilde_cW(i-1) + tilde_w(i) + enddo + tilde_cW(:) = tilde_cW(:) + 1.d0 + + pt2_n_0(1) = 0 + do + pt2_u_0 = tilde_cW(pt2_n_0(1)) + r = tilde_cW(pt2_n_0(1) + pt2_minDetInFirstTeeth) + pt2_W_T = (1d0 - pt2_u_0) / dble(pt2_N_teeth) + if(pt2_W_T >= r - pt2_u_0) then + exit + end if + pt2_n_0(1) += 1 + if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then + print *, "teeth building failed" + stop -1 + end if + end do + + do t=2, pt2_N_teeth + r = pt2_u_0 + pt2_W_T * dble(t-1) + pt2_n_0(t) = pt2_find_sample(r, tilde_cW) + end do + pt2_n_0(pt2_N_teeth+1) = N_det_generators + + pt2_w(:pt2_n_0(1)) = tilde_w(:pt2_n_0(1)) + do t=1, pt2_N_teeth + tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t)) + if (tooth_width == 0.d0) then + tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))) + endif + ASSERT(tooth_width > 0.d0) + do i=pt2_n_0(t)+1, pt2_n_0(t+1) + pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width + end do + end do + + pt2_cW(0) = 0d0 + do i=1,N_det_generators + pt2_cW(i) = pt2_cW(i-1) + pt2_w(i) + end do + pt2_n_0(pt2_N_teeth+1) = N_det_generators + + endif +END_PROVIDER + + + + + +BEGIN_PROVIDER [ double precision, pt2_overlap, (N_states, N_states) ] + implicit none + BEGIN_DOC + ! Overlap between the perturbed wave functions + END_DOC + pt2_overlap(1:N_states,1:N_states) = 0.d0 +END_PROVIDER + + diff --git a/plugins/local/cipsi_tc_bi_ortho/pt2_type.irp.f b/src/cipsi_utils/pt2_type.irp.f similarity index 100% rename from plugins/local/cipsi_tc_bi_ortho/pt2_type.irp.f rename to src/cipsi_utils/pt2_type.irp.f diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi_utils/run_pt2_slave.irp.f similarity index 100% rename from src/cipsi/run_pt2_slave.irp.f rename to src/cipsi_utils/run_pt2_slave.irp.f diff --git a/src/cipsi_utils/run_selection_slave.irp.f b/src/cipsi_utils/run_selection_slave.irp.f new file mode 100644 index 00000000..783bed0f --- /dev/null +++ b/src/cipsi_utils/run_selection_slave.irp.f @@ -0,0 +1,257 @@ +subroutine run_selection_slave(thread,iproc,energy) + use f77_zmq + use selection_types + implicit none + + double precision, intent(in) :: energy(N_states) + integer, intent(in) :: thread, iproc + integer :: rc, i + + integer :: worker_id, task_id(1), ctask, ltask + character*(512) :: task + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_push_socket + integer(ZMQ_PTR) :: zmq_socket_push + + type(selection_buffer) :: buf, buf2 + logical :: done, buffer_ready + type(pt2_type) :: pt2_data + + PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique + PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns + PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym + PROVIDE psi_bilinear_matrix_rows psi_bilinear_matrix_order weight_selection + + call provide_for_selection_slave + + call pt2_alloc(pt2_data,N_states) + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + + integer, external :: connect_to_taskserver + if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + return + endif + + zmq_socket_push = new_zmq_push_socket(thread) + + buf%N = 0 + buffer_ready = .False. + ctask = 1 + + do + integer, external :: get_task_from_taskserver + if (get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task) == -1) then + exit + endif + done = task_id(ctask) == 0 + if (done) then + ctask = ctask - 1 + else + integer :: i_generator, N, subset, bsize + call sscanf_ddd(task, subset, i_generator, N) + if(buf%N == 0) then + ! Only first time + call create_selection_buffer(N, N*2, buf) + buffer_ready = .True. + else + if (N /= buf%N) then + print *, 'N=', N + print *, 'buf%N=', buf%N + print *, 'bug in ', irp_here + stop '-1' + end if + end if + call select_connected(i_generator, energy, pt2_data, buf, subset, pt2_F(i_generator)) + endif + + integer, external :: task_done_to_taskserver + + if(done .or. ctask == size(task_id)) then + do i=1, ctask + if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i)) == -1) then + call usleep(100) + if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i)) == -1) then + ctask = 0 + done = .true. + exit + endif + endif + end do + if(ctask > 0) then + call sort_selection_buffer(buf) +! call merge_selection_buffers(buf,buf2) + call push_selection_results(zmq_socket_push, pt2_data, buf, task_id(1), ctask) + call pt2_dealloc(pt2_data) + call pt2_alloc(pt2_data,N_states) +! buf%mini = buf2%mini + buf%cur = 0 + end if + ctask = 0 + end if + + if(done) exit + ctask = ctask + 1 + end do + + if(ctask > 0) then + call sort_selection_buffer(buf) +! call merge_selection_buffers(buf,buf2) + call push_selection_results(zmq_socket_push, pt2_data, buf, task_id(1), ctask) +! buf%mini = buf2%mini + buf%cur = 0 + end if + ctask = 0 + call pt2_dealloc(pt2_data) + + integer, external :: disconnect_from_taskserver + if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then + continue + endif + + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_push_socket(zmq_socket_push,thread) + if (buffer_ready) then + call delete_selection_buffer(buf) +! call delete_selection_buffer(buf2) + endif +end subroutine + + +subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntasks) + use f77_zmq + use selection_types + implicit none + + integer(ZMQ_PTR), intent(in) :: zmq_socket_push + type(pt2_type), intent(in) :: pt2_data + type(selection_buffer), intent(inout) :: b + integer, intent(in) :: ntasks, task_id(*) + integer :: rc + double precision, allocatable :: pt2_serialized(:) + + rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE) + if(rc /= 4) then + print *, 'f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE)' + endif + + + allocate(pt2_serialized (pt2_type_size(N_states)) ) + call pt2_serialize(pt2_data,N_states,pt2_serialized) + + rc = f77_zmq_send( zmq_socket_push, pt2_serialized, size(pt2_serialized)*8, ZMQ_SNDMORE) + if (rc == -1) then + print *, irp_here, ': error sending result' + stop 3 + return + else if(rc /= size(pt2_serialized)*8) then + stop 'push' + endif + deallocate(pt2_serialized) + + if (b%cur > 0) then + + rc = f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE) + if(rc /= 8*b%cur) then + print *, 'f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE)' + endif + + rc = f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE) + if(rc /= bit_kind*N_int*2*b%cur) then + print *, 'f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE)' + endif + + endif + + rc = f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE) + if(rc /= 4) then + print *, 'f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE)' + endif + + rc = f77_zmq_send( zmq_socket_push, task_id(1), ntasks*4, 0) + if(rc /= 4*ntasks) then + print *, 'f77_zmq_send( zmq_socket_push, task_id(1), ntasks*4, 0)' + endif + +! Activate is zmq_socket_push is a REQ +IRP_IF ZMQ_PUSH +IRP_ELSE + character*(2) :: ok + rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0) + if ((rc /= 2).and.(ok(1:2) /= 'ok')) then + print *, irp_here//': error in receiving ok' + stop -1 + endif +IRP_ENDIF + +end subroutine + + +subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_id, ntasks) + use f77_zmq + use selection_types + implicit none + integer(ZMQ_PTR), intent(in) :: zmq_socket_pull + type(pt2_type), intent(inout) :: pt2_data + double precision, intent(out) :: val(*) + integer(bit_kind), intent(out) :: det(N_int, 2, *) + integer, intent(out) :: N, ntasks, task_id(*) + integer :: rc, rn, i + double precision, allocatable :: pt2_serialized(:) + + rc = f77_zmq_recv( zmq_socket_pull, N, 4, 0) + if(rc /= 4) then + print *, 'f77_zmq_recv( zmq_socket_pull, N, 4, 0)' + endif + + allocate(pt2_serialized (pt2_type_size(N_states)) ) + rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized), 0) + if (rc == -1) then + ntasks = 1 + task_id(1) = 0 + else if(rc /= 8*size(pt2_serialized)) then + stop 'pull' + endif + + call pt2_deserialize(pt2_data,N_states,pt2_serialized) + deallocate(pt2_serialized) + + if (N>0) then + rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, 0) + if(rc /= 8*N) then + print *, 'f77_zmq_recv( zmq_socket_pull, val(1), 8*N, 0)' + endif + + rc = f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, 0) + if(rc /= bit_kind*N_int*2*N) then + print *, 'f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, 0)' + endif + endif + + rc = f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0) + if(rc /= 4) then + print *, 'f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0)' + endif + + rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntasks*4, 0) + if(rc /= 4*ntasks) then + print *, 'f77_zmq_recv( zmq_socket_pull, task_id(1), ntasks*4, 0)' + endif + +! Activate is zmq_socket_pull is a REP +IRP_IF ZMQ_PUSH +IRP_ELSE + rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0) + if (rc /= 2) then + print *, irp_here//': error in sending ok' + stop -1 + endif +IRP_ENDIF +end subroutine + + + diff --git a/src/cipsi/selection_buffer.irp.f b/src/cipsi_utils/selection_buffer.irp.f similarity index 100% rename from src/cipsi/selection_buffer.irp.f rename to src/cipsi_utils/selection_buffer.irp.f diff --git a/plugins/local/cipsi_tc_bi_ortho/selection_types.f90 b/src/cipsi_utils/selection_types.f90 similarity index 100% rename from plugins/local/cipsi_tc_bi_ortho/selection_types.f90 rename to src/cipsi_utils/selection_types.f90 diff --git a/src/cipsi/selection_weight.irp.f b/src/cipsi_utils/selection_weight.irp.f similarity index 100% rename from src/cipsi/selection_weight.irp.f rename to src/cipsi_utils/selection_weight.irp.f diff --git a/src/cipsi/slave_cipsi.irp.f b/src/cipsi_utils/slave_cipsi.irp.f similarity index 98% rename from src/cipsi/slave_cipsi.irp.f rename to src/cipsi_utils/slave_cipsi.irp.f index ddfc050e..8be48f40 100644 --- a/src/cipsi/slave_cipsi.irp.f +++ b/src/cipsi_utils/slave_cipsi.irp.f @@ -303,10 +303,11 @@ subroutine run_slave_main PROVIDE global_selection_buffer pt2_N_teeth pt2_F N_det_generators PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique - PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order + PROVIDE psi_bilinear_matrix_rows psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted + PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp PROVIDE psi_det_hii selection_weight pseudo_sym pt2_min_parallel_tasks + call provide_for_zmq_pt2 if (mpi_master) then print *, 'Running PT2' From 6b7f2411b17c87368cbe56a03aad157819fcd1aa Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 12 Mar 2024 17:31:49 +0100 Subject: [PATCH 08/12] Add NEED in cipsi_utils --- src/cipsi_utils/NEED | 1 + 1 file changed, 1 insertion(+) create mode 100644 src/cipsi_utils/NEED diff --git a/src/cipsi_utils/NEED b/src/cipsi_utils/NEED new file mode 100644 index 00000000..d3d4d2c7 --- /dev/null +++ b/src/cipsi_utils/NEED @@ -0,0 +1 @@ +determinants From 37588e520766f303acaecd26b1dc16484b69f80f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 12 Mar 2024 17:32:38 +0100 Subject: [PATCH 09/12] Add NEED in generators_full_tc --- src/generators_full_tc/NEED | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 src/generators_full_tc/NEED diff --git a/src/generators_full_tc/NEED b/src/generators_full_tc/NEED new file mode 100644 index 00000000..0cf7d3aa --- /dev/null +++ b/src/generators_full_tc/NEED @@ -0,0 +1,2 @@ +determinants +hartree_fock From 0618372b29284e16aeb3dd0cfc9b62377571a03d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 12 Mar 2024 17:38:30 +0100 Subject: [PATCH 10/12] Commented out select_singles in TC --- plugins/local/cipsi_tc_bi_ortho/selection.irp.f | 8 ++++---- src/.gitignore | 1 + 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f index 9b8cc81e..b1c02102 100644 --- a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f @@ -88,10 +88,10 @@ subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset) 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 - if ((subset == 1).and.(sum(hole_mask(:,2)) == 0_bit_kind)) then - ! No beta electron to excite - call select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,b) - endif +! if ((subset == 1).and.(sum(hole_mask(:,2)) == 0_bit_kind)) then +! ! No beta electron to excite +! call select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,b) +! endif 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 diff --git a/src/.gitignore b/src/.gitignore index 6353c21a..abc6a4c0 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1,5 +1,6 @@ * !README.rst +!NEED !*/ */* !*/*.* From fdc418d72a12eb307a0cf875225794fbd37dde11 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 12 Mar 2024 17:45:50 +0100 Subject: [PATCH 11/12] fixed print in TC --- plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 1 + plugins/local/fci_tc_bi/diagonalize_ci.irp.f | 6 ++++-- 2 files changed, 5 insertions(+), 2 deletions(-) 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 8863b7bc..721564e6 100644 --- a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -91,6 +91,7 @@ subroutine run_stochastic_cipsi to_select = max(N_states_diag, to_select) E_denom = E_tc ! TC Energy of the current wave function + print*,'E_tc = ',E_tc call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) call pt2_alloc(pt2_data, N_states) diff --git a/plugins/local/fci_tc_bi/diagonalize_ci.irp.f b/plugins/local/fci_tc_bi/diagonalize_ci.irp.f index a5242b87..5fcce5eb 100644 --- a/plugins/local/fci_tc_bi/diagonalize_ci.irp.f +++ b/plugins/local/fci_tc_bi/diagonalize_ci.irp.f @@ -55,9 +55,11 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) ! 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 + do k = 1, N_states + 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) From a56488e3a865dccc98d7984dd2cc4a7be1885539 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 12 Mar 2024 18:23:09 +0100 Subject: [PATCH 12/12] fci_tc_bi_ortho works for multi state ninja --- .../cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 57 ++++--------------- plugins/local/fci_tc_bi/diagonalize_ci.irp.f | 42 +------------- .../local/tc_bi_ortho/psi_det_tc_sorted.irp.f | 8 ++- src/cipsi/pt2_stoch_routines.irp.f | 2 +- src/cipsi_utils/slave_cipsi.irp.f | 2 +- 5 files changed, 20 insertions(+), 91 deletions(-) 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 721564e6..99a8de7e 100644 --- a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -11,15 +11,13 @@ subroutine run_stochastic_cipsi implicit none integer :: i, j, k, ndet integer :: to_select - logical :: print_pt2 logical :: has type(pt2_type) :: pt2_data, pt2_data_err double precision :: rss - double precision :: correlation_energy_ratio, E_denom, E_tc, norm + double precision :: correlation_energy_ratio double precision :: hf_energy_ref double precision :: relative_error - double precision, allocatable :: ept2(:), pt1(:), extrap_energy(:) - double precision, allocatable :: zeros(:) + double precision, allocatable :: zeros(:),E_tc(:), norm(:) logical, external :: qp_stop double precision, external :: memory_of_double @@ -32,14 +30,13 @@ subroutine run_stochastic_cipsi write(*,*) i, Fock_matrix_tc_mo_tot(i,i) enddo - N_iter = 1 threshold_generators = 1.d0 SOFT_TOUCH threshold_generators rss = memory_of_double(N_states)*4.d0 call check_mem(rss, irp_here) - allocate(zeros(N_states)) + allocate(zeros(N_states),E_tc(N_states), norm(N_states)) call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data_err, N_states) @@ -55,8 +52,7 @@ subroutine run_stochastic_cipsi ! 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 diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm) ! if (N_det > N_det_max) then @@ -67,19 +63,16 @@ subroutine run_stochastic_cipsi ! 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 diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm) ! call routine_save_right ! endif - allocate(ept2(1000),pt1(1000),extrap_energy(100)) correlation_energy_ratio = 0.d0 ! thresh_it_dav = 5.d-5 ! soft_touch thresh_it_dav - print_pt2 = .True. do while( (N_det < N_det_max) .and. & (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max)) @@ -90,13 +83,12 @@ subroutine run_stochastic_cipsi to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) - E_denom = E_tc ! TC Energy of the current wave function print*,'E_tc = ',E_tc call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data_err, N_states) - call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection + call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection ! stop call print_summary_tc(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2) @@ -117,48 +109,19 @@ subroutine run_stochastic_cipsi PROVIDE psi_det PROVIDE psi_det_sorted_tc - ept2(N_iter-1) = E_tc + nuclear_repulsion + (pt2_data % pt2(1))/norm - pt1(N_iter-1) = dsqrt(pt2_data % overlap(1,1)) - call diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) + call diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm) ! stop if (qp_stop()) exit enddo -! print*,'data to extrapolate ' -! do i = 2, N_iter -! print*,'iteration ',i -! print*,'pt1,Ept2',pt1(i),ept2(i) -! call get_extrapolated_energy(i-1,ept2(i),pt1(i),extrap_energy(i)) -! do j = 2, i -! print*,'j,e,energy',j,extrap_energy(j) -! enddo -! enddo - -! thresh_it_dav = 5.d-6 -! soft_touch thresh_it_dav call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data_err, N_states) call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection - call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) -! if (.not.qp_stop()) then -! if (N_det < N_det_max) then -! thresh_it_dav = 5.d-7 -! soft_touch thresh_it_dav -! call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) -! endif -! -! call pt2_dealloc(pt2_data) -! call pt2_dealloc(pt2_data_err) -! call pt2_alloc(pt2_data, N_states) -! call pt2_alloc(pt2_data_err, N_states) -! call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error, 0) ! Stochastic PT2 -! call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) -! endif -! call pt2_dealloc(pt2_data) -! call pt2_dealloc(pt2_data_err) -! call routine_save_right + call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm) + call pt2_dealloc(pt2_data) + call pt2_dealloc(pt2_data_err) end diff --git a/plugins/local/fci_tc_bi/diagonalize_ci.irp.f b/plugins/local/fci_tc_bi/diagonalize_ci.irp.f index 5fcce5eb..85518116 100644 --- a/plugins/local/fci_tc_bi/diagonalize_ci.irp.f +++ b/plugins/local/fci_tc_bi/diagonalize_ci.irp.f @@ -1,7 +1,7 @@ ! --- -subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) +subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm ) BEGIN_DOC ! Replace the coefficients of the CI states by the coefficients of the @@ -12,50 +12,10 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) implicit none integer, intent(inout) :: ndet ! number of determinants from before 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,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 -! 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 -! enddo do k = 1, N_states E_tc(k) = eigval_right_tc_bi_orth(k) norm(k) = norm_ground_left_right_bi_orth(k) diff --git a/plugins/local/tc_bi_ortho/psi_det_tc_sorted.irp.f b/plugins/local/tc_bi_ortho/psi_det_tc_sorted.irp.f index 5dad91ca..eef99de8 100644 --- a/plugins/local/tc_bi_ortho/psi_det_tc_sorted.irp.f +++ b/plugins/local/tc_bi_ortho/psi_det_tc_sorted.irp.f @@ -11,10 +11,16 @@ BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_tc, (psi_det_size) ] psi_average_norm_contrib_tc(:) = 0.d0 do k=1,N_states do i=1,N_det - psi_average_norm_contrib_tc(i) = & +! print*,dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k)),psi_l_coef_bi_ortho(i,k),psi_r_coef_bi_ortho(i,k) + psi_average_norm_contrib_tc(i) += & dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k))*state_average_weight(k) enddo enddo +! print*,'***' +! do i = 1, N_det +! print*,psi_average_norm_contrib_tc(i) +! enddo + print*,'sum(psi_average_norm_contrib_tc(1:N_det))',sum(psi_average_norm_contrib_tc(1:N_det)) f = 1.d0/sum(psi_average_norm_contrib_tc(1:N_det)) do i=1,N_det psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i)*f diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 228e0ef1..bd5943da 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -1,3 +1,3 @@ subroutine provide_for_zmq_pt2 - PROVIDE psi_selectors_coef_transp psi_det_sorted psi_det_sorted_order + PROVIDE psi_selectors_coef_transp psi_det_sorted psi_det_sorted_order psi_det_hii end diff --git a/src/cipsi_utils/slave_cipsi.irp.f b/src/cipsi_utils/slave_cipsi.irp.f index 8be48f40..3e778270 100644 --- a/src/cipsi_utils/slave_cipsi.irp.f +++ b/src/cipsi_utils/slave_cipsi.irp.f @@ -306,7 +306,7 @@ subroutine run_slave_main PROVIDE psi_bilinear_matrix_rows psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp - PROVIDE psi_det_hii selection_weight pseudo_sym pt2_min_parallel_tasks + PROVIDE selection_weight pseudo_sym pt2_min_parallel_tasks call provide_for_zmq_pt2 if (mpi_master) then