From dc0d668f3829fee928b645f469a6c43f90608a8d Mon Sep 17 00:00:00 2001 From: amandadumi Date: Tue, 16 Jun 2020 10:57:24 -0400 Subject: [PATCH 01/27] changing molden 'Atoms' label to match coordinate units --- src/tools/molden.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f index f70ed0de..417b25ad 100644 --- a/src/tools/molden.irp.f +++ b/src/tools/molden.irp.f @@ -17,7 +17,7 @@ program molden write(i_unit_output,'(A)') '[Molden Format]' - write(i_unit_output,'(A)') '[Atoms] AU' + write(i_unit_output,'(A)') '[Atoms] Angs' do i = 1, nucl_num write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') & trim(element_name(int(nucl_charge(i)))), & From e0f17d516b91ba2ada7032791e228725ce9ae7ce Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Aug 2020 00:10:46 +0200 Subject: [PATCH 02/27] Introduced pt2_type --- src/cipsi/cipsi.irp.f | 57 +++++++++++++++++------------ src/cipsi/pt2_stoch_routines.irp.f | 28 +++++++------- src/cipsi/selection.irp.f | 14 +++++-- src/cipsi/selection_types.f90 | 6 +++ src/cipsi/stochastic_cipsi.irp.f | 59 +++++++++++++++++++----------- src/cipsi/zmq_selection.irp.f | 25 ++++++------- src/fci/pt2.irp.f | 22 +++++++---- 7 files changed, 126 insertions(+), 85 deletions(-) diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 413a9d9a..060722a1 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -1,18 +1,20 @@ subroutine run_cipsi implicit none + use selection_types BEGIN_DOC ! Selected Full Configuration Interaction with deterministic selection and ! stochastic PT2. END_DOC integer :: i,j,k - double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:) + type(pt2_type) :: pt2_data + double precision, allocatable :: rpt2(:), zeros(:) integer :: to_select logical, external :: qp_stop double precision :: threshold_generators_save double precision :: rss double precision, external :: memory_of_double - PROVIDE H_apply_buffer_allocated + PROVIDE H_apply_buffer_allocated N_iter = 1 threshold_generators = 1.d0 @@ -21,7 +23,8 @@ subroutine run_cipsi rss = memory_of_double(N_states)*4.d0 call check_mem(rss,irp_here) - allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states)) + allocate (zeros(N_states), rpt2(N_states)) + allocate (pt2_data % pt2(N_states), pt2_data % norm2(N_states), pt2_data % variance(N_states)) double precision :: hf_energy_ref logical :: has @@ -30,10 +33,10 @@ subroutine run_cipsi relative_error=PT2_relative_error zeros = 0.d0 - pt2 = -huge(1.e0) rpt2 = -huge(1.e0) - norm = 0.d0 - variance = huge(1.e0) + pt2_data % pt2 = -huge(1.e0) + pt2_data % norm2 = 0.d0 + pt2_data % variance = huge(1.e0) if (s2_eig) then call make_s2_eigenfunction @@ -67,8 +70,8 @@ subroutine run_cipsi do while ( & (N_det < N_det_max) .and. & - (maxval(abs(rpt2(1:N_states))) > pt2_max) .and. & - (maxval(abs(variance(1:N_states))) > variance_max) .and. & + (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. & + (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and. & (correlation_energy_ratio <= correlation_energy_ratio_max) & ) write(*,'(A)') '--------------------------------------------------------------------------------' @@ -77,22 +80,21 @@ subroutine run_cipsi to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) if (do_pt2) then - pt2 = 0.d0 - variance = 0.d0 - norm = 0.d0 + pt2_data % pt2 = 0.d0 + pt2_data % variance = 0.d0 + pt2_data % norm2 = 0.d0 threshold_generators_save = threshold_generators threshold_generators = 1.d0 SOFT_TOUCH threshold_generators - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & - norm, 0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,error, 0) ! Stochastic PT2 threshold_generators = threshold_generators_save SOFT_TOUCH threshold_generators else - call ZMQ_selection(to_select, pt2, variance, norm) + call ZMQ_selection(to_select, pt2_data) endif do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & @@ -100,7 +102,11 @@ subroutine run_cipsi correlation_energy_ratio = min(1.d0,correlation_energy_ratio) call write_double(6,correlation_energy_ratio, 'Correlation ratio') - call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep, & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, rpt2) @@ -121,7 +127,7 @@ subroutine run_cipsi call diagonalize_CI call save_wavefunction call save_energy(psi_energy_with_nucl_rep, zeros) - if (qp_stop()) exit + if (qp_stop()) exit enddo if (.not.qp_stop()) then @@ -132,13 +138,12 @@ subroutine run_cipsi endif if (do_pt2) then - pt2(:) = 0.d0 - variance(:) = 0.d0 - norm(:) = 0.d0 + pt2_data % pt2(:) = 0.d0 + pt2_data % variance(:) = 0.d0 + pt2_data % norm2(:) = 0.d0 threshold_generators = 1d0 SOFT_TOUCH threshold_generators - call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, & - norm,0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, error, 0) ! Stochastic PT2 SOFT_TOUCH threshold_generators endif print *, 'N_det = ', N_det @@ -148,11 +153,15 @@ subroutine run_cipsi do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo call save_energy(psi_energy_with_nucl_rep, rpt2) - call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep(1:N_states), & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call print_extrapolated_energy() endif diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 52831c16..ca4b3ba7 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -107,7 +107,7 @@ end function -subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) +subroutine ZMQ_pt2(E, pt2_data, relative_error, error, N_in) use f77_zmq use selection_types @@ -117,10 +117,9 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) integer, intent(in) :: N_in ! integer, intent(inout) :: N_in double precision, intent(in) :: relative_error, E(N_states) - double precision, intent(out) :: pt2(N_states),error(N_states) - double precision, intent(out) :: variance(N_states),norm2(N_states) - - + double precision, intent(out) :: error(N_states) + type(pt2_type), intent(inout) :: pt2_data +! integer :: i, N double precision :: state_average_weight_save(N_states), w(N_states,4) @@ -138,10 +137,10 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) endif if (N_det <= max(4,N_states) .or. pt2_N_teeth < 2) then - pt2=0.d0 - variance=0.d0 - norm2=0.d0 - call ZMQ_selection(N_in, pt2, variance, norm2) + pt2_data % pt2=0.d0 + pt2_data % variance=0.d0 + pt2_data % norm2=0.d0 + call ZMQ_selection(N_in, pt2_data % pt2, pt2_data % variance, pt2_data % norm2) error(:) = 0.d0 else @@ -304,10 +303,11 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) if (i==0) then call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, w(1,1), w(1,2), w(1,3), w(1,4), b, N) - pt2(pt2_stoch_istate) = w(pt2_stoch_istate,1) + pt2_data % pt2(pt2_stoch_istate) = w(pt2_stoch_istate,1) error(pt2_stoch_istate) = w(pt2_stoch_istate,2) - variance(pt2_stoch_istate) = w(pt2_stoch_istate,3) - norm2(pt2_stoch_istate) = w(pt2_stoch_istate,4) + pt2_data % variance(pt2_stoch_istate) = w(pt2_stoch_istate,3) + pt2_data % norm2(pt2_stoch_istate) = w(pt2_stoch_istate,4) + !TODO SEGV else call pt2_slave_inproc(i) @@ -335,10 +335,10 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) TOUCH state_average_weight endif do k=N_det+1,N_states - pt2(k) = 0.d0 + pt2_data % pt2(k) = 0.d0 enddo - call update_pt2_and_variance_weights(pt2, variance, norm2, N_states) + call update_pt2_and_variance_weights(pt2_data, N_states) end subroutine diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 0ca08a0f..288a155c 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -19,15 +19,17 @@ BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ] variance_match_weight(:) = 1.d0 END_PROVIDER -subroutine update_pt2_and_variance_weights(pt2, variance, norm2, N_st) +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 - double precision, intent(in) :: pt2(N_st) - double precision, intent(in) :: variance(N_st) - double precision, intent(in) :: norm2(N_st) + type(pt2_type), intent(in) :: pt2_data + double precision :: pt2(N_st) + double precision :: variance(N_st) + double precision :: norm2(N_st) double precision :: avg, rpt2(N_st), element, dt, x integer :: k @@ -35,6 +37,10 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm2, N_st) integer, parameter :: i_itermax = 1 double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:) + pt2(:) = pt2_data % pt2(:) + variance(:) = pt2_data % variance(:) + norm2(:) = pt2_data % norm2(:) + if (i_iter == 0) then allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax)) memo_pt2(:,:) = 1.d0 diff --git a/src/cipsi/selection_types.f90 b/src/cipsi/selection_types.f90 index 29e48524..fd76215a 100644 --- a/src/cipsi/selection_types.f90 +++ b/src/cipsi/selection_types.f90 @@ -5,5 +5,11 @@ module selection_types double precision, pointer :: val(:) double precision :: mini endtype + + type pt2_type + double precision, allocatable :: pt2(:) + double precision, allocatable :: variance(:) + double precision, allocatable :: norm2(:) + endtype end module diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index b1534eae..449fd61a 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -1,12 +1,15 @@ subroutine run_stochastic_cipsi + use selection_types implicit none BEGIN_DOC ! Selected Full Configuration Interaction with Stochastic selection and PT2. END_DOC integer :: i,j,k - double precision, allocatable :: pt2(:), variance(:), norm2(:), rpt2(:), zeros(:) + double precision, allocatable :: rpt2(:), zeros(:) integer :: to_select - logical, external :: qp_stop + type(pt2_type) :: pt2_data + logical, external :: qp_stop + double precision :: rss double precision, external :: memory_of_double @@ -19,7 +22,11 @@ subroutine run_stochastic_cipsi rss = memory_of_double(N_states)*4.d0 call check_mem(rss,irp_here) - allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm2(N_states), variance(N_states)) + allocate (zeros(N_states), rpt2(N_states)) + allocate( pt2_data % pt2(N_states) ) + allocate( pt2_data % variance(N_states) ) + allocate( pt2_data % norm2(N_states) ) + double precision :: hf_energy_ref logical :: has @@ -28,10 +35,10 @@ subroutine run_stochastic_cipsi relative_error=PT2_relative_error zeros = 0.d0 - pt2 = -huge(1.e0) + pt2_data % pt2 = -huge(1.e0) rpt2 = -huge(1.e0) - norm2 = 0.d0 - variance = huge(1.e0) + pt2_data % norm2 = 0.d0 + pt2_data % variance = huge(1.e0) if (s2_eig) then call make_s2_eigenfunction @@ -65,8 +72,8 @@ subroutine run_stochastic_cipsi do while ( & (N_det < N_det_max) .and. & - (maxval(abs(rpt2(1:N_states))) > pt2_max) .and. & - (maxval(abs(variance(1:N_states))) > variance_max) .and. & + (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. & + (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and. & (correlation_energy_ratio <= correlation_energy_ratio_max) & ) write(*,'(A)') '--------------------------------------------------------------------------------' @@ -75,14 +82,15 @@ subroutine run_stochastic_cipsi to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) - pt2 = 0.d0 - variance = 0.d0 - norm2 = 0.d0 - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & - norm2, to_select) ! Stochastic PT2 and selection + + pt2_data % pt2 = 0.d0 + pt2_data % variance = 0.d0 + pt2_data % norm2 = 0.d0 + call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,error, & + to_select) ! Stochastic PT2 and selection do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm2(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & @@ -90,7 +98,11 @@ subroutine run_stochastic_cipsi correlation_energy_ratio = min(1.d0,correlation_energy_ratio) call write_double(6,correlation_energy_ratio, 'Correlation ratio') - call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm2,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep, & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, rpt2) @@ -121,18 +133,21 @@ subroutine run_stochastic_cipsi call save_energy(psi_energy_with_nucl_rep, zeros) endif - pt2(:) = 0.d0 - variance(:) = 0.d0 - norm2(:) = 0.d0 - call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, & - norm2,0) ! Stochastic PT2 + pt2_data % pt2(:) = 0.d0 + pt2_data % variance(:) = 0.d0 + pt2_data % norm2(:) = 0.d0 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, error, 0) ! Stochastic PT2 do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm2(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo call save_energy(psi_energy_with_nucl_rep, rpt2) - call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm2,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep, & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call print_extrapolated_energy() endif diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 22919117..19e8b33b 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -1,4 +1,4 @@ -subroutine ZMQ_selection(N_in, pt2, variance, norm2) +subroutine ZMQ_selection(N_in, pt2_data) use f77_zmq use selection_types @@ -9,9 +9,7 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm2) type(selection_buffer) :: b integer :: i, N integer, external :: omp_get_thread_num - double precision, intent(out) :: pt2(N_states) - double precision, intent(out) :: variance(N_states) - double precision, intent(out) :: norm2(N_states) + 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 @@ -112,19 +110,20 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm2) enddo endif - !$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2, variance, norm2) PRIVATE(i) NUM_THREADS(nproc_target+1) + !$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, variance, norm2) + call selection_collector(zmq_socket_pull, b, N, & + pt2_data % pt2, pt2_data % variance, pt2_data % norm2) else call selection_slave_inproc(i) endif !$OMP END PARALLEL call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'selection') do i=N_det+1,N_states - pt2(i) = 0.d0 - variance(i) = 0.d0 - norm2(i) = 0.d0 + pt2_data % pt2(i) = 0.d0 + pt2_data % variance(i) = 0.d0 + pt2_data % norm2(i) = 0.d0 enddo if (N_in > 0) then if (s2_eig) then @@ -134,12 +133,12 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm2) endif call delete_selection_buffer(b) do k=1,N_states - pt2(k) = pt2(k) * f(k) - variance(k) = variance(k) * f(k) - norm2(k) = norm2(k) * f(k) + pt2_data % pt2(k) = pt2_data % pt2(k) * f(k) + pt2_data % variance(k) = pt2_data % variance(k) * f(k) + pt2_data % norm2(k) = pt2_data % norm2(k) * f(k) enddo - call update_pt2_and_variance_weights(pt2, variance, norm2, N_states) + call update_pt2_and_variance_weights(pt2_data, N_states) end subroutine diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f index 83994781..8c076ade 100644 --- a/src/fci/pt2.irp.f +++ b/src/fci/pt2.irp.f @@ -28,35 +28,41 @@ end subroutine run implicit none + use selection_types integer :: i,j,k logical, external :: detEq - double precision :: pt2(N_states) + type(pt2_type) :: pt2_data integer :: degree integer :: n_det_before, to_select double precision :: threshold_davidson_in double precision :: E_CI_before(N_states), relative_error, error(N_states), variance(N_states), norm2(N_states), rpt2(N_states) - pt2(:) = 0.d0 + allocate( pt2_data % pt2(N_states) ) + allocate( pt2_data % variance(N_states) ) + allocate( pt2_data % norm2(N_states) ) E_CI_before(:) = psi_energy(:) + nuclear_repulsion relative_error=PT2_relative_error if (do_pt2) then - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & - norm2,0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,error,0) ! Stochastic PT2 else - call ZMQ_selection(0, pt2, variance, norm2) + call ZMQ_selection(0, pt2_data) endif do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm2(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo - call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm2,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep(1:N_states), & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) - call save_energy(E_CI_before,pt2) + call save_energy(E_CI_before,pt2_data % pt2) end From 7bde6f74514ae0034dfacdf68c6708eb434d8566 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Aug 2020 15:39:01 +0200 Subject: [PATCH 03/27] Factorized pt2 data --- src/cipsi/cipsi.irp.f | 15 ++---- src/cipsi/pt2_stoch_routines.irp.f | 15 ++---- src/cipsi/pt2_type.irp.f | 74 ++++++++++++++++++++++++++++++ src/cipsi/selection_types.f90 | 9 ++++ src/cipsi/stochastic_cipsi.irp.f | 24 +++------- src/cipsi/zmq_selection.irp.f | 49 ++++++++++---------- src/fci/pt2.irp.f | 18 ++------ src/iterations/print_summary.irp.f | 44 +++++++++--------- 8 files changed, 148 insertions(+), 100 deletions(-) create mode 100644 src/cipsi/pt2_type.irp.f diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 060722a1..df699ad4 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -64,7 +64,6 @@ subroutine run_cipsi endif double precision :: correlation_energy_ratio - double precision :: error(N_states) correlation_energy_ratio = 0.d0 @@ -86,7 +85,7 @@ subroutine run_cipsi threshold_generators_save = threshold_generators threshold_generators = 1.d0 SOFT_TOUCH threshold_generators - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,error, 0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error, 0) ! Stochastic PT2 threshold_generators = threshold_generators_save SOFT_TOUCH threshold_generators else @@ -103,10 +102,7 @@ subroutine run_cipsi call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep, & - pt2_data % pt2, error, & - pt2_data % variance, & - pt2_data % norm2, & - N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, N_det,N_occ_pattern,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, rpt2) @@ -143,7 +139,7 @@ subroutine run_cipsi pt2_data % norm2(:) = 0.d0 threshold_generators = 1d0 SOFT_TOUCH threshold_generators - call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, error, 0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2 SOFT_TOUCH threshold_generators endif print *, 'N_det = ', N_det @@ -158,10 +154,7 @@ subroutine run_cipsi call save_energy(psi_energy_with_nucl_rep, rpt2) call print_summary(psi_energy_with_nucl_rep(1:N_states), & - pt2_data % pt2, error, & - pt2_data % variance, & - pt2_data % norm2, & - N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, N_det,N_occ_pattern,N_states,psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call print_extrapolated_energy() endif diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index ca4b3ba7..78555dc9 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -107,7 +107,7 @@ end function -subroutine ZMQ_pt2(E, pt2_data, relative_error, error, N_in) +subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in) use f77_zmq use selection_types @@ -117,7 +117,6 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, error, N_in) integer, intent(in) :: N_in ! integer, intent(inout) :: N_in double precision, intent(in) :: relative_error, E(N_states) - double precision, intent(out) :: error(N_states) type(pt2_type), intent(inout) :: pt2_data ! integer :: i, N @@ -137,11 +136,7 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, error, N_in) endif if (N_det <= max(4,N_states) .or. pt2_N_teeth < 2) then - pt2_data % pt2=0.d0 - pt2_data % variance=0.d0 - pt2_data % norm2=0.d0 - call ZMQ_selection(N_in, pt2_data % pt2, pt2_data % variance, pt2_data % norm2) - error(:) = 0.d0 + call ZMQ_selection(N_in, pt2_data) else N = max(N_in,1) * N_states @@ -304,10 +299,9 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, error, N_in) call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, w(1,1), w(1,2), w(1,3), w(1,4), b, N) pt2_data % pt2(pt2_stoch_istate) = w(pt2_stoch_istate,1) - error(pt2_stoch_istate) = w(pt2_stoch_istate,2) + pt2_data % pt2_err(pt2_stoch_istate) = w(pt2_stoch_istate,2) pt2_data % variance(pt2_stoch_istate) = w(pt2_stoch_istate,3) pt2_data % norm2(pt2_stoch_istate) = w(pt2_stoch_istate,4) - !TODO SEGV else call pt2_slave_inproc(i) @@ -334,9 +328,6 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, error, N_in) state_average_weight(:) = state_average_weight_save(:) TOUCH state_average_weight endif - do k=N_det+1,N_states - pt2_data % pt2(k) = 0.d0 - enddo call update_pt2_and_variance_weights(pt2_data, N_states) diff --git a/src/cipsi/pt2_type.irp.f b/src/cipsi/pt2_type.irp.f new file mode 100644 index 00000000..8f8bb225 --- /dev/null +++ b/src/cipsi/pt2_type.irp.f @@ -0,0 +1,74 @@ +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 % pt2_err(N) & + ,pt2_data % variance(N) & + ,pt2_data % variance_err(N) & + ,pt2_data % norm2(N) & + ,pt2_data % norm2_err(N) & + ,pt2_data % rpt2(N) & + ,pt2_data % rpt2_err(N) & + ,pt2_data % overlap(N,N) & + ,pt2_data % overlap_err(N,N) & + ) + + pt2_data % pt2(:) = 0.d0 + pt2_data % pt2_err(:) = 0.d0 + pt2_data % variance(:) = 0.d0 + pt2_data % variance_err(:) = 0.d0 + pt2_data % norm2(:) = 0.d0 + pt2_data % norm2_err(:) = 0.d0 + pt2_data % rpt2(:) = 0.d0 + pt2_data % rpt2_err(:) = 0.d0 + pt2_data % overlap(:,:) = 0.d0 + pt2_data % overlap_err(:,:) = 0.d0 + + do k=1,N + pt2_data % overlap(k,k) = 1.d0 + enddo +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 % pt2_err & + ,pt2_data % variance & + ,pt2_data % variance_err& + ,pt2_data % norm2 & + ,pt2_data % norm2_err & + ,pt2_data % rpt2 & + ,pt2_data % rpt2_err & + ,pt2_data % overlap & + ,pt2_data % overlap_err & + ) +end subroutine + +subroutine pt2_add(p1, p2) + implicit none + use selection_types + BEGIN_DOC +! p1 += p2 + END_DOC + type(pt2_type), intent(inout) :: p1 + type(pt2_type), intent(in) :: p2 + + p1 % pt2(:) += p2 % pt2(:) + p1 % pt2_err(:) += p2 % pt2_err(:) + p1 % rpt2(:) += p2 % rpt2(:) + p1 % rpt2_err(:) += p2 % rpt2_err(:) + p1 % variance(:) += p2 % variance(:) + p1 % variance_err(:) += p2 % variance_err(:) + p1 % norm2(:) += p2 % norm2(:) + p1 % norm2_err(:) += p2 % norm2_err(:) + p1 % overlap(:,:) += p2 % overlap(:,:) + p1 % overlap_err(:,:) += p2 % overlap_err(:,:) +end subroutine + + diff --git a/src/cipsi/selection_types.f90 b/src/cipsi/selection_types.f90 index fd76215a..1905f2d3 100644 --- a/src/cipsi/selection_types.f90 +++ b/src/cipsi/selection_types.f90 @@ -8,8 +8,17 @@ module selection_types type pt2_type double precision, allocatable :: pt2(:) + double precision, allocatable :: pt2_err(:) + double precision, allocatable :: rpt2(:) + double precision, allocatable :: rpt2_err(:) double precision, allocatable :: variance(:) + double precision, allocatable :: variance_err(:) double precision, allocatable :: norm2(:) + double precision, allocatable :: norm2_err(:) + double precision, allocatable :: overlap(:,:) + double precision, allocatable :: overlap_err(:,:) endtype + + end module diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 449fd61a..00495e2b 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -23,10 +23,7 @@ subroutine run_stochastic_cipsi call check_mem(rss,irp_here) allocate (zeros(N_states), rpt2(N_states)) - allocate( pt2_data % pt2(N_states) ) - allocate( pt2_data % variance(N_states) ) - allocate( pt2_data % norm2(N_states) ) - + call pt2_alloc(pt2_data, N_states) double precision :: hf_energy_ref logical :: has @@ -35,8 +32,8 @@ subroutine run_stochastic_cipsi relative_error=PT2_relative_error zeros = 0.d0 - pt2_data % pt2 = -huge(1.e0) rpt2 = -huge(1.e0) + pt2_data % pt2 = -huge(1.e0) pt2_data % norm2 = 0.d0 pt2_data % variance = huge(1.e0) @@ -66,7 +63,6 @@ subroutine run_stochastic_cipsi endif double precision :: correlation_energy_ratio - double precision :: error(N_states) correlation_energy_ratio = 0.d0 @@ -86,8 +82,7 @@ subroutine run_stochastic_cipsi pt2_data % pt2 = 0.d0 pt2_data % variance = 0.d0 pt2_data % norm2 = 0.d0 - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,error, & - to_select) ! Stochastic PT2 and selection + call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,to_select) ! Stochastic PT2 and selection do k=1,N_states rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) @@ -99,10 +94,7 @@ subroutine run_stochastic_cipsi call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep, & - pt2_data % pt2, error, & - pt2_data % variance, & - pt2_data % norm2, & - N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, N_det,N_occ_pattern,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, rpt2) @@ -136,7 +128,7 @@ subroutine run_stochastic_cipsi pt2_data % pt2(:) = 0.d0 pt2_data % variance(:) = 0.d0 pt2_data % norm2(:) = 0.d0 - call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, error, 0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2 do k=1,N_states rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) @@ -144,12 +136,10 @@ subroutine run_stochastic_cipsi call save_energy(psi_energy_with_nucl_rep, rpt2) call print_summary(psi_energy_with_nucl_rep, & - pt2_data % pt2, error, & - pt2_data % variance, & - pt2_data % norm2, & - N_det,N_occ_pattern,N_states,psi_s2) + pt2_data , N_det, N_occ_pattern, N_states, psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call print_extrapolated_energy() endif + call pt2_dealloc(pt2_data) end diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 19e8b33b..754dd80a 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -104,7 +104,7 @@ subroutine ZMQ_selection(N_in, pt2_data) f(:) = 1.d0 if (.not.do_pt2) then - double precision :: f(N_states), u_dot_u + 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 @@ -113,18 +113,13 @@ subroutine ZMQ_selection(N_in, pt2_data) !$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 % pt2, pt2_data % variance, pt2_data % norm2) + 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') - do i=N_det+1,N_states - pt2_data % pt2(i) = 0.d0 - pt2_data % variance(i) = 0.d0 - pt2_data % norm2(i) = 0.d0 - enddo if (N_in > 0) then if (s2_eig) then call make_selection_buffer_s2(b) @@ -132,11 +127,14 @@ subroutine ZMQ_selection(N_in, pt2_data) 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) - pt2_data % norm2(k) = pt2_data % norm2(k) * f(k) - enddo + + if (.not.do_pt2) then + 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) + pt2_data % norm2(k) = pt2_data % norm2(k) * f(k) + enddo + endif call update_pt2_and_variance_weights(pt2_data, N_states) @@ -150,7 +148,7 @@ subroutine selection_slave_inproc(i) call run_selection_slave(1,i,pt2_e0_denominator) end -subroutine selection_collector(zmq_socket_pull, b, N, pt2, variance, norm2) +subroutine selection_collector(zmq_socket_pull, b, N, pt2_data) use f77_zmq use selection_types use bitmasks @@ -160,12 +158,11 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2, variance, norm2) integer(ZMQ_PTR), intent(in) :: zmq_socket_pull type(selection_buffer), intent(inout) :: b integer, intent(in) :: N - double precision, intent(out) :: pt2(N_states) - double precision, intent(out) :: variance(N_states) - double precision, intent(out) :: norm2(N_states) - double precision :: pt2_mwen(N_states) - double precision :: variance_mwen(N_states) - double precision :: norm2_mwen(N_states) + type(pt2_type), intent(inout) :: pt2_data + + 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 @@ -189,18 +186,18 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2, variance, norm2) call check_mem(rss,irp_here) allocate(task_id(N_det_generators)) more = 1 - pt2(:) = 0d0 - variance(:) = 0.d0 - norm2(:) = 0.d0 + pt2_data % pt2(:) = 0d0 + pt2_data % variance(:) = 0.d0 + pt2_data % norm2(:) = 0.d0 pt2_mwen(:) = 0.d0 variance_mwen(:) = 0.d0 norm2_mwen(:) = 0.d0 do while (more == 1) call pull_selection_results(zmq_socket_pull, pt2_mwen, variance_mwen, norm2_mwen, b2%val(1), b2%det(1,1,1), b2%cur, task_id, ntask) - pt2(:) += pt2_mwen(:) - variance(:) += variance_mwen(:) - norm2(:) += norm2_mwen(:) + pt2_data % pt2(:) += pt2_mwen(:) + pt2_data % variance(:) += variance_mwen(:) + pt2_data % norm2(:) += norm2_mwen(:) 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 diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f index 8c076ade..c34705f5 100644 --- a/src/fci/pt2.irp.f +++ b/src/fci/pt2.irp.f @@ -37,32 +37,24 @@ subroutine run integer :: n_det_before, to_select double precision :: threshold_davidson_in - double precision :: E_CI_before(N_states), relative_error, error(N_states), variance(N_states), norm2(N_states), rpt2(N_states) + double precision :: E_CI_before(N_states), relative_error - allocate( pt2_data % pt2(N_states) ) - allocate( pt2_data % variance(N_states) ) - allocate( pt2_data % norm2(N_states) ) + call pt2_alloc(pt2_data, N_states) E_CI_before(:) = psi_energy(:) + nuclear_repulsion relative_error=PT2_relative_error if (do_pt2) then - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,error,0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2 else call ZMQ_selection(0, pt2_data) endif - do k=1,N_states - rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) - enddo - call print_summary(psi_energy_with_nucl_rep(1:N_states), & - pt2_data % pt2, error, & - pt2_data % variance, & - pt2_data % norm2, & - N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, N_det,N_occ_pattern,N_states,psi_s2) call save_energy(E_CI_before,pt2_data % pt2) + call pt2_dealloc(pt2_data) end diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index ad87bc8e..3b00b783 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -1,11 +1,13 @@ -subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_st,s2_) +subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,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_occ_pattern_, n_st - double precision, intent(in) :: e_(n_st), pt2_(n_st), variance_(n_st), norm_(n_st), error_(n_st), s2_(n_st) + double precision, intent(in) :: e_(n_st), s2_(n_st) + type(pt2_type) , intent(in) :: pt2_data integer :: i, k integer :: N_states_p character*(9) :: pt2_string @@ -21,7 +23,7 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_ N_states_p = min(N_det_,n_st) do i=1,N_states_p - f(i) = 1.d0/(1.d0+norm_(i)) + f(i) = 1.d0/(1.d0+pt2_data % norm2(i)) enddo print *, '' @@ -42,16 +44,16 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_ 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_(k), error_(k), k=1,N_states_p) - write(*,fmt) '# rPT2'//pt2_string, (pt2_(k)*f(k), error_(k)*f(k), k=1,N_states_p) + write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data % pt2_err(k), k=1,N_states_p) + write(*,fmt) '# rPT2'//pt2_string, (pt2_data % pt2(k)*f(k), pt2_data % pt2_err(k)*f(k), k=1,N_states_p) write(*,'(A)') '#' - write(*,fmt) '# E+PT2 ', (e_(k)+pt2_(k),error_(k), k=1,N_states_p) - write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_(k)*f(k),error_(k)*f(k), k=1,N_states_p) + write(*,fmt) '# E+PT2 ', (e_(k)+pt2_data % pt2(k),pt2_data % pt2_err(k), k=1,N_states_p) + write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % pt2(k)*f(k),pt2_data % pt2_err(k)*f(k), k=1,N_states_p) if (N_states_p > 1) then - write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_(k)-e_(1)-pt2_(1)), & - dsqrt(error_(k)*error_(k)+error_(1)*error_(1)), k=1,N_states_p) - write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_(k)-e_(1)-pt2_(1))*27.211396641308d0, & - dsqrt(error_(k)*error_(k)+error_(1)*error_(1))*27.211396641308d0, k=1,N_states_p) + write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), & + dsqrt(pt2_data % pt2_err(k)*pt2_data % pt2_err(k)+pt2_data % pt2_err(1)*pt2_data % pt2_err(1)), k=1,N_states_p) + write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, & + dsqrt(pt2_data % pt2_err(k)*pt2_data % pt2_err(k)+pt2_data % pt2_err(1)*pt2_data % pt2_err(1))*27.211396641308d0, k=1,N_states_p) endif write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' write(*,fmt) @@ -68,12 +70,12 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_ print*,'* State ',k print *, '< S^2 > = ', s2_(k) print *, 'E = ', e_(k) - print *, 'Variance = ', variance_(k) - print *, 'PT norm = ', dsqrt(norm_(k)) - print *, 'PT2 = ', pt2_(k) - print *, 'rPT2 = ', pt2_(k)*f(k) - print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_(k), ' +/- ', error_(k) - print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_(k)*f(k), ' +/- ', error_(k)*f(k) + print *, 'Variance = ', pt2_data % variance(k) + print *, 'PT norm = ', dsqrt(pt2_data % norm2(k)) + print *, 'PT2 = ', pt2_data % pt2(k) + print *, 'rPT2 = ', pt2_data % pt2(k)*f(k) + print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data % pt2_err(k) + print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % pt2(k)*f(k), ' +/- ', pt2_data % pt2_err(k)*f(k) print *, '' enddo @@ -87,14 +89,14 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_ print *, '-----' print*, 'Variational + perturbative Energy difference (au | eV)' do i=2, N_states_p - print*,'Delta E = ', (e_(i)+ pt2_(i) - (e_(1) + pt2_(1))), & - (e_(i)+ pt2_(i) - (e_(1) + pt2_(1))) * 27.211396641308d0 + print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), & + (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0 enddo print *, '-----' print*, 'Variational + renormalized perturbative Energy difference (au | eV)' do i=2, N_states_p - print*,'Delta E = ', (e_(i)+ pt2_(i)*f(i) - (e_(1) + pt2_(1)*f(1))), & - (e_(i)+ pt2_(i)*f(i) - (e_(1) + pt2_(1)*f(1))) * 27.211396641308d0 + print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i)*f(i) - (e_(1) + pt2_data % pt2(1)*f(1))), & + (e_(i)+ pt2_data % pt2(i)*f(i) - (e_(1) + pt2_data % pt2(1)*f(1))) * 27.211396641308d0 enddo endif From cf2d78fced92a660ff145b33d0a79b5f70318976 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Aug 2020 16:05:53 +0200 Subject: [PATCH 04/27] Cleaned rPT2 --- src/cipsi/cipsi.irp.f | 29 ++++++++--------------- src/cipsi/pt2_stoch_routines.irp.f | 38 ++++++++++++++++-------------- src/cipsi/selection.irp.f | 7 +----- src/cipsi/stochastic_cipsi.irp.f | 26 +++++++------------- src/cipsi/zmq_selection.irp.f | 9 +++---- 5 files changed, 45 insertions(+), 64 deletions(-) diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index df699ad4..44b87394 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -7,7 +7,7 @@ subroutine run_cipsi END_DOC integer :: i,j,k type(pt2_type) :: pt2_data - double precision, allocatable :: rpt2(:), zeros(:) + double precision, allocatable :: zeros(:) integer :: to_select logical, external :: qp_stop @@ -23,8 +23,8 @@ subroutine run_cipsi rss = memory_of_double(N_states)*4.d0 call check_mem(rss,irp_here) - allocate (zeros(N_states), rpt2(N_states)) - allocate (pt2_data % pt2(N_states), pt2_data % norm2(N_states), pt2_data % variance(N_states)) + allocate (zeros(N_states)) + call pt2_alloc(pt2_data, N_states) double precision :: hf_energy_ref logical :: has @@ -33,8 +33,8 @@ subroutine run_cipsi relative_error=PT2_relative_error zeros = 0.d0 - rpt2 = -huge(1.e0) - pt2_data % pt2 = -huge(1.e0) + pt2_data % pt2 = -huge(1.e0) + pt2_data % rpt2 = -huge(1.e0) pt2_data % norm2 = 0.d0 pt2_data % variance = huge(1.e0) @@ -92,21 +92,17 @@ subroutine run_cipsi call ZMQ_selection(to_select, pt2_data) endif - do k=1,N_states - rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) - enddo - correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & - (psi_energy_with_nucl_rep(1) + rpt2(1) - hf_energy_ref) + (psi_energy_with_nucl_rep(1) + pt2_data % rpt2(1) - hf_energy_ref) correlation_energy_ratio = min(1.d0,correlation_energy_ratio) call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep, & pt2_data, N_det,N_occ_pattern,N_states,psi_s2) - call save_energy(psi_energy_with_nucl_rep, rpt2) + call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) - call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) + call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call print_extrapolated_energy() N_iter += 1 @@ -147,15 +143,10 @@ subroutine run_cipsi print *, 'N_states = ', N_states print*, 'correlation_ratio = ', correlation_energy_ratio - - do k=1,N_states - rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) - enddo - - call save_energy(psi_energy_with_nucl_rep, rpt2) + call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) call print_summary(psi_energy_with_nucl_rep(1:N_states), & pt2_data, N_det,N_occ_pattern,N_states,psi_s2) - call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) + call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call print_extrapolated_energy() endif diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 78555dc9..92b598a9 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -292,16 +292,19 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in) 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, w(1,1), w(1,2), w(1,3), w(1,4), b, N) - pt2_data % pt2(pt2_stoch_istate) = w(pt2_stoch_istate,1) - pt2_data % pt2_err(pt2_stoch_istate) = w(pt2_stoch_istate,2) - pt2_data % variance(pt2_stoch_istate) = w(pt2_stoch_istate,3) - pt2_data % norm2(pt2_stoch_istate) = w(pt2_stoch_istate,4) + call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, pt2_data, b, N) + pt2_data % rpt2(pt2_stoch_istate) = & + pt2_data % pt2(pt2_stoch_istate)/(1.d0 + pt2_data % norm2(pt2_stoch_istate)) + + !TODO : We should use here the correct formula for the error of X/Y + pt2_data % rpt2_err(pt2_stoch_istate) = & + pt2_data % pt2_err(pt2_stoch_istate)/(1.d0 + pt2_data % norm2(pt2_stoch_istate)) else call pt2_slave_inproc(i) @@ -343,7 +346,7 @@ subroutine pt2_slave_inproc(i) end -subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error, variance, norm2, b, N_) +subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) use f77_zmq use selection_types use bitmasks @@ -352,8 +355,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error, varianc integer(ZMQ_PTR), intent(in) :: zmq_socket_pull double precision, intent(in) :: relative_error, E - double precision, intent(out) :: pt2(N_states), error(N_states) - double precision, intent(out) :: variance(N_states), norm2(N_states) + type(pt2_type), intent(inout) :: pt2_data type(selection_buffer), intent(inout) :: b integer, intent(in) :: N_ @@ -409,10 +411,10 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error, varianc call create_selection_buffer(N_, N_*2, b2) - pt2(:) = -huge(1.) - error(:) = huge(1.) - variance(:) = huge(1.) - norm2(:) = 0.d0 + pt2_data % pt2(pt2_stoch_istate) = -huge(1.) + pt2_data % pt2_err(pt2_stoch_istate) = huge(1.) + pt2_data % variance(pt2_stoch_istate) = huge(1.) + pt2_data % norm2(pt2_stoch_istate) = 0.d0 S(:) = 0d0 S2(:) = 0d0 T2(:) = 0d0 @@ -486,21 +488,21 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error, varianc if (qp_stop()) then stop_now = .True. endif - pt2(pt2_stoch_istate) = avg - variance(pt2_stoch_istate) = avg2 - norm2(pt2_stoch_istate) = avg3 + pt2_data % pt2(pt2_stoch_istate) = avg + pt2_data % variance(pt2_stoch_istate) = avg2 + pt2_data % norm2(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((S2(t) / c) - (S(t)/c)**2) ! dabs for numerical stability eqt = sqrt(eqt / (dble(c) - 1.5d0)) - error(pt2_stoch_istate) = eqt + pt2_data % pt2_err(pt2_stoch_istate) = eqt if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, F14.10, 2X, F14.10, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3, time-time0, '' if (stop_now .or. ( & - (do_exit .and. (dabs(error(pt2_stoch_istate)) / & - (1.d-20 + dabs(pt2(pt2_stoch_istate)) ) <= relative_error))) ) then + (do_exit .and. (dabs(pt2_data % pt2_err(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 diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 288a155c..cf52aaff 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -31,7 +31,7 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) double precision :: variance(N_st) double precision :: norm2(N_st) - double precision :: avg, rpt2(N_st), element, dt, x + double precision :: avg, element, dt, x integer :: k integer, save :: i_iter=0 integer, parameter :: i_itermax = 1 @@ -54,11 +54,6 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) dt = 2.0d0 - do k=1,N_st - ! rPT2 - rpt2(k) = pt2(k)/(1.d0 + norm2(k)) - enddo - avg = sum(pt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero do k=1,N_st element = exp(dt*(pt2(k)/avg -1.d0)) diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 00495e2b..e80bd856 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -5,7 +5,7 @@ subroutine run_stochastic_cipsi ! Selected Full Configuration Interaction with Stochastic selection and PT2. END_DOC integer :: i,j,k - double precision, allocatable :: rpt2(:), zeros(:) + double precision, allocatable :: zeros(:) integer :: to_select type(pt2_type) :: pt2_data logical, external :: qp_stop @@ -22,7 +22,7 @@ subroutine run_stochastic_cipsi rss = memory_of_double(N_states)*4.d0 call check_mem(rss,irp_here) - allocate (zeros(N_states), rpt2(N_states)) + allocate (zeros(N_states)) call pt2_alloc(pt2_data, N_states) double precision :: hf_energy_ref @@ -32,8 +32,8 @@ subroutine run_stochastic_cipsi relative_error=PT2_relative_error zeros = 0.d0 - rpt2 = -huge(1.e0) - pt2_data % pt2 = -huge(1.e0) + pt2_data % pt2 = -huge(1.e0) + pt2_data % rpt2 = -huge(1.e0) pt2_data % norm2 = 0.d0 pt2_data % variance = huge(1.e0) @@ -84,21 +84,17 @@ subroutine run_stochastic_cipsi pt2_data % norm2 = 0.d0 call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,to_select) ! Stochastic PT2 and selection - do k=1,N_states - rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) - enddo - correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & - (psi_energy_with_nucl_rep(1) + rpt2(1) - hf_energy_ref) + (psi_energy_with_nucl_rep(1) + pt2_data % rpt2(1) - hf_energy_ref) correlation_energy_ratio = min(1.d0,correlation_energy_ratio) call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep, & pt2_data, N_det,N_occ_pattern,N_states,psi_s2) - call save_energy(psi_energy_with_nucl_rep, rpt2) + call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) - call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) + call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call print_extrapolated_energy() N_iter += 1 @@ -130,14 +126,10 @@ subroutine run_stochastic_cipsi pt2_data % norm2(:) = 0.d0 call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2 - do k=1,N_states - rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) - enddo - - call save_energy(psi_energy_with_nucl_rep, rpt2) + call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) call print_summary(psi_energy_with_nucl_rep, & pt2_data , N_det, N_occ_pattern, N_states, psi_s2) - call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) + call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call print_extrapolated_energy() endif call pt2_dealloc(pt2_data) diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 754dd80a..006e6578 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -128,13 +128,14 @@ subroutine ZMQ_selection(N_in, pt2_data) endif call delete_selection_buffer(b) - if (.not.do_pt2) then - do k=1,N_states + 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) pt2_data % norm2(k) = pt2_data % norm2(k) * f(k) - enddo - endif + + pt2_data % rpt2(k) = & + pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) + enddo call update_pt2_and_variance_weights(pt2_data, N_states) From 3ec31857f9f89426f31c9fce9a91e0d4565cab93 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Aug 2020 23:57:11 +0200 Subject: [PATCH 05/27] Try to fix travis NaN --- src/hartree_fock/10.hf.bats | 1 + 1 file changed, 1 insertion(+) diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats index 6f13e595..54845685 100644 --- a/src/hartree_fock/10.hf.bats +++ b/src/hartree_fock/10.hf.bats @@ -11,6 +11,7 @@ function run() { qp edit --check qp reset --mos qp set scf_utils n_it_scf_max 50 + qp set ao_one_e_ints lin_dep_cutoff 1.e-50 qp run scf # qp set_frozen_core energy="$(ezfio get hartree_fock energy)" From 32dd686f9696d70ec4f1233e216e75bc4efb333e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 29 Aug 2020 00:22:48 +0200 Subject: [PATCH 06/27] Introduced error bars in variance and norm --- src/cipsi/pt2_stoch_routines.irp.f | 52 +++++++++++++++++++----------- src/iterations/print_summary.irp.f | 12 +++---- 2 files changed, 39 insertions(+), 25 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 92b598a9..04df5e16 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -360,9 +360,9 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) integer, intent(in) :: N_ - double precision, allocatable :: eI(:,:), eI_task(:,:), S(:), S2(:) - double precision, allocatable :: vI(:,:), vI_task(:,:), T2(:) - double precision, allocatable :: nI(:,:), nI_task(:,:), T3(:) + double precision, allocatable :: eI(:,:), eI_task(:,:), Se(:), Se2(:) + double precision, allocatable :: vI(:,:), vI_task(:,:), Sv(:), Sv2(:) + double precision, allocatable :: nI(:,:), nI_task(:,:), Sn(:), Sn2(:) 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 @@ -402,8 +402,9 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) allocate(eI(N_states, N_det_generators), eI_task(N_states, pt2_n_tasks_max)) allocate(vI(N_states, N_det_generators), vI_task(N_states, pt2_n_tasks_max)) allocate(nI(N_states, N_det_generators), nI_task(N_states, pt2_n_tasks_max)) - allocate(S(pt2_N_teeth+1), S2(pt2_N_teeth+1)) - allocate(T2(pt2_N_teeth+1), T3(pt2_N_teeth+1)) + allocate(Se(pt2_N_teeth+1), Se2(pt2_N_teeth+1)) + allocate(Sv(pt2_N_teeth+1), Sv2(pt2_N_teeth+1)) + allocate(Sn(pt2_N_teeth+1), Sn2(pt2_N_teeth+1)) @@ -415,10 +416,12 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) pt2_data % pt2_err(pt2_stoch_istate) = huge(1.) pt2_data % variance(pt2_stoch_istate) = huge(1.) pt2_data % norm2(pt2_stoch_istate) = 0.d0 - S(:) = 0d0 - S2(:) = 0d0 - T2(:) = 0d0 - T3(:) = 0d0 + Se(:) = 0d0 + Sv(:) = 0d0 + Sn(:) = 0d0 + Se2(:) = 0d0 + Sv2(:) = 0d0 + Sn2(:) = 0d0 n = 1 t = 0 U = 0 @@ -474,14 +477,16 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) x += eI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i) x2 += vI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i) x3 += nI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i) - S(p) += x - S2(p) += x*x - T2(p) += x2 - T3(p) += x3 + Se(p) += x + Sv(p) += x2 + Sn(p) += x3 + Se2(p) += x*x + Sv2(p) += x2*x2 + Sn2(p) += x3*3 end do - avg = E0 + S(t) / dble(c) - avg2 = v0 + T2(t) / dble(c) - avg3 = n0 + T3(t) / dble(c) + avg = E0 + Se(t) / dble(c) + avg2 = v0 + Sv(t) / dble(c) + avg3 = n0 + Sn(t) / dble(c) if ((avg /= 0.d0) .or. (n == N_det_generators) ) then do_exit = .true. endif @@ -494,9 +499,18 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) call wall_time(time) ! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969) if(c > 2) then - eqt = dabs((S2(t) / c) - (S(t)/c)**2) ! dabs for numerical stability + eqt = dabs((Se2(t) / c) - (Se(t)/c)**2) ! dabs for numerical stability eqt = sqrt(eqt / (dble(c) - 1.5d0)) pt2_data % pt2_err(pt2_stoch_istate) = eqt + + eqt = dabs((Sv2(t) / c) - (Sv(t)/c)**2) ! dabs for numerical stability + eqt = sqrt(eqt / (dble(c) - 1.5d0)) + pt2_data % variance_err(pt2_stoch_istate) = eqt + + eqt = dabs((Sn2(t) / c) - (Sn(t)/c)**2) ! dabs for numerical stability + eqt = sqrt(eqt / (dble(c) - 1.5d0)) + pt2_data % norm2_err(pt2_stoch_istate) = eqt + if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, F14.10, 2X, F14.10, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3, time-time0, '' @@ -520,7 +534,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) call pull_pt2_results(zmq_socket_pull, index, eI_task, vI_task, nI_task, task_id, n_tasks, b2) if(n_tasks > pt2_n_tasks_max)then print*,'PB !!!' - print*,'If you see this, send an email to Anthony scemama with the following content' + 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 @@ -531,7 +545,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) do i=1,n_tasks if(index(i).gt.size(eI,2).or.index(i).lt.1)then print*,'PB !!!' - print*,'If you see this, send an email to Anthony scemama with the following content' + print*,'If you see this, send a bug report with the following content' print*,irp_here print*,'i,index(i),size(ei,2) = ',i,index(i),size(ei,2) stop -1 diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 3b00b783..32c07ce5 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -70,12 +70,12 @@ subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,n_st,s2_) print*,'* State ',k print *, '< S^2 > = ', s2_(k) print *, 'E = ', e_(k) - print *, 'Variance = ', pt2_data % variance(k) - print *, 'PT norm = ', dsqrt(pt2_data % norm2(k)) - print *, 'PT2 = ', pt2_data % pt2(k) - print *, 'rPT2 = ', pt2_data % pt2(k)*f(k) - print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data % pt2_err(k) - print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % pt2(k)*f(k), ' +/- ', pt2_data % pt2_err(k)*f(k) + print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data % variance_err(k) + print *, 'PT norm = ', dsqrt(pt2_data % norm2(k)), ' +/- ', 0.5d0*dsqrt(pt2_data % norm2(k)) * pt2_data % norm2_err(k) / pt2_data % norm2(k) + print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data % pt2_err(k) + print *, 'rPT2 = ', pt2_data % pt2(k)*f(k), ' +/- ', pt2_data % rpt2_err(k) + print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data % pt2_err(k) + print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % pt2(k)*f(k), ' +/- ', pt2_data % pt2_err(k)*f(k) print *, '' enddo From 622aee8bf5d7e1be0424b14b27a98b0290e6b4f5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 29 Aug 2020 01:15:48 +0200 Subject: [PATCH 07/27] Introduced overlap --- src/cipsi/pert_rdm_providers.irp.f | 12 ++-- src/cipsi/pt2_stoch_routines.irp.f | 25 +++---- src/cipsi/run_pt2_slave.irp.f | 109 +++++++++-------------------- src/cipsi/selection.irp.f | 64 ++++++++++------- src/cipsi/selection_types.f90 | 7 ++ 5 files changed, 96 insertions(+), 121 deletions(-) diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f index e2917261..82ef1e63 100644 --- a/src/cipsi/pert_rdm_providers.irp.f +++ b/src/cipsi/pert_rdm_providers.irp.f @@ -31,7 +31,7 @@ BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert END_PROVIDER -subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf, psi_det_connection, psi_coef_connection_reverse, n_det_connection) +subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, psi_det_connection, psi_coef_connection_reverse, n_det_connection) use bitmasks use selection_types implicit none @@ -44,9 +44,7 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num) double precision, intent(in) :: fock_diag_tmp(mo_num) double precision, intent(in) :: E0(N_states) - double precision, intent(inout) :: pt2(N_states) - double precision, intent(inout) :: variance(N_states) - double precision, intent(inout) :: norm(N_states) + type(pt2_type), intent(inout) :: pt2_data type(selection_buffer), intent(inout) :: buf logical :: ok integer :: s1, s2, p1, p2, ib, j, istate @@ -152,9 +150,9 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo e_pert = 0.5d0 * (tmp - delta_E) coef(istate) = e_pert / alpha_h_psi print*,e_pert,coef,alpha_h_psi - pt2(istate) = pt2(istate) + e_pert - variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi - norm(istate) = norm(istate) + coef(istate) * coef(istate) + pt2_data % pt2(istate) += e_pert + pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi + pt2_data % norm2(istate) = coef(istate) * coef(istate) if (weight_selection /= 5) then ! Energy selection diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 04df5e16..5acd4752 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -243,7 +243,7 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in) ( 1.d0*pt2_n_tasks_max & ! task_id, index + 0.635d0*N_det_generators & ! f,d + 3.d0*N_det_generators*N_states & ! eI, vI, nI - + 3.d0*pt2_n_tasks_max*N_states & ! eI_task, vI_task, nI_task + + pt2_n_tasks_max*pt2_type_size(N_states)/8 & ! pt2_data_task + 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 @@ -359,10 +359,10 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) type(selection_buffer), intent(inout) :: b integer, intent(in) :: N_ - - double precision, allocatable :: eI(:,:), eI_task(:,:), Se(:), Se2(:) - double precision, allocatable :: vI(:,:), vI_task(:,:), Sv(:), Sv2(:) - double precision, allocatable :: nI(:,:), nI_task(:,:), Sn(:), Sn2(:) + type(pt2_type), allocatable :: pt2_data_task(:) + double precision, allocatable :: eI(:,:), Se(:), Se2(:) + double precision, allocatable :: vI(:,:), Sv(:), Sv2(:) + double precision, allocatable :: nI(:,:), Sn(:), Sn2(:) 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 @@ -399,9 +399,10 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) ! 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(eI(N_states, N_det_generators), eI_task(N_states, pt2_n_tasks_max)) - allocate(vI(N_states, N_det_generators), vI_task(N_states, pt2_n_tasks_max)) - allocate(nI(N_states, N_det_generators), nI_task(N_states, pt2_n_tasks_max)) + allocate(pt2_data_task(pt2_n_tasks_max)) + allocate(eI(N_states, N_det_generators)) + allocate(vI(N_states, N_det_generators)) + allocate(nI(N_states, N_det_generators)) allocate(Se(pt2_N_teeth+1), Se2(pt2_N_teeth+1)) allocate(Sv(pt2_N_teeth+1), Sv2(pt2_N_teeth+1)) allocate(Sn(pt2_N_teeth+1), Sn2(pt2_N_teeth+1)) @@ -531,7 +532,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) else if(more == 0) then exit else - call pull_pt2_results(zmq_socket_pull, index, eI_task, vI_task, nI_task, task_id, n_tasks, b2) + 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' @@ -550,9 +551,9 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) print*,'i,index(i),size(ei,2) = ',i,index(i),size(ei,2) stop -1 endif - eI(1:N_states, index(i)) += eI_task(1:N_states,i) - vI(1:N_states, index(i)) += vI_task(1:N_states,i) - nI(1:N_states, index(i)) += nI_task(1:N_states,i) + eI(1:N_states, index(i)) += pt2_data_task(i) % pt2(1:N_states) + vI(1:N_states, index(i)) += pt2_data_task(i) % variance(1:N_states) + nI(1:N_states, index(i)) += pt2_data_task(i) % norm2(1:N_states) f(index(i)) -= 1 end do do i=1, b2%cur diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi/run_pt2_slave.irp.f index 7df98a87..ae0ac2a0 100644 --- a/src/cipsi/run_pt2_slave.irp.f +++ b/src/cipsi/run_pt2_slave.irp.f @@ -1,8 +1,8 @@ - use omp_lib + use omp_lib use selection_types use f77_zmq BEGIN_PROVIDER [ integer(omp_lock_kind), global_selection_buffer_lock ] - use omp_lib + use omp_lib implicit none BEGIN_DOC ! Global buffer for the OpenMP selection @@ -11,7 +11,7 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), global_selection_buffer_lock ] END_PROVIDER BEGIN_PROVIDER [ type(selection_buffer), global_selection_buffer ] - use omp_lib + use omp_lib implicit none BEGIN_DOC ! Global buffer for the OpenMP selection @@ -61,7 +61,7 @@ subroutine run_pt2_slave_small(thread,iproc,energy) type(selection_buffer) :: b logical :: done, buffer_ready - double precision,allocatable :: pt2(:,:), variance(:,:), norm(:,:) + type(pt2_type), allocatable :: pt2_data(:) integer :: n_tasks, k, N integer, allocatable :: i_generator(:), subset(:) @@ -70,10 +70,7 @@ subroutine run_pt2_slave_small(thread,iproc,energy) ! logical :: sending allocate(task_id(pt2_n_tasks_max), task(pt2_n_tasks_max)) - allocate(pt2(N_states,pt2_n_tasks_max), i_generator(pt2_n_tasks_max), subset(pt2_n_tasks_max)) - allocate(variance(N_states,pt2_n_tasks_max)) - allocate(norm(N_states,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 @@ -120,13 +117,11 @@ subroutine run_pt2_slave_small(thread,iproc,energy) double precision :: time0, time1 call wall_time(time0) do k=1,n_tasks - pt2(:,k) = 0.d0 - variance(:,k) = 0.d0 - norm(:,k) = 0.d0 + call pt2_alloc(pt2_data(k),N_states) b%cur = 0 !double precision :: time2 !call wall_time(time2) - call select_connected(i_generator(k),energy,pt2(1,k),variance(1,k),norm(1,k),b,subset(k),pt2_F(i_generator(k))) + call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k))) !call wall_time(time1) !print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1)) enddo @@ -138,7 +133,10 @@ subroutine run_pt2_slave_small(thread,iproc,energy) done = .true. endif call sort_selection_buffer(b) - call push_pt2_results(zmq_socket_push, i_generator, pt2, variance, norm, b, task_id, n_tasks) + 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 @@ -183,7 +181,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy) type(selection_buffer) :: b logical :: done, buffer_ready - double precision,allocatable :: pt2(:,:), variance(:,:), norm(:,:) + type(pt2_type), allocatable :: pt2_data(:) integer :: n_tasks, k, N integer, allocatable :: i_generator(:), subset(:) @@ -193,9 +191,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy) allocate(task_id(pt2_n_tasks_max), task(pt2_n_tasks_max)) - allocate(pt2(N_states,pt2_n_tasks_max), i_generator(pt2_n_tasks_max), subset(pt2_n_tasks_max)) - allocate(variance(N_states,pt2_n_tasks_max)) - allocate(norm(N_states,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() @@ -243,13 +239,11 @@ subroutine run_pt2_slave_large(thread,iproc,energy) double precision :: time0, time1 call wall_time(time0) do k=1,n_tasks - pt2(:,k) = 0.d0 - variance(:,k) = 0.d0 - norm(:,k) = 0.d0 + call pt2_alloc(pt2_data(k),N_states) b%cur = 0 !double precision :: time2 !call wall_time(time2) - call select_connected(i_generator(k),energy,pt2(1,k),variance(1,k),norm(1,k),b,subset(k),pt2_F(i_generator(k))) + call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k))) !call wall_time(time1) !print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1)) enddo @@ -269,13 +263,16 @@ subroutine run_pt2_slave_large(thread,iproc,energy) 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, variance, norm, global_selection_buffer, task_id, n_tasks,sending) + call push_pt2_results_async_send(zmq_socket_push, i_generator, pt2_data, global_selection_buffer, task_id, n_tasks,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, variance, norm, b, task_id, n_tasks,sending) + call push_pt2_results_async_send(zmq_socket_push, i_generator, pt2_data, b, task_id, n_tasks,sending) endif + do k=1,n_tasks + call pt2_dealloc(pt2_data(k)) + enddo ! ! 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 = 1 @@ -298,34 +295,30 @@ subroutine run_pt2_slave_large(thread,iproc,energy) end subroutine -subroutine push_pt2_results(zmq_socket_push, index, pt2, variance, norm, b, task_id, n_tasks) +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 - double precision, intent(in) :: pt2(N_states,n_tasks) - double precision, intent(in) :: variance(N_states,n_tasks) - double precision, intent(in) :: norm(N_states,n_tasks) + 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, variance, norm, b, task_id, n_tasks, sending) + 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, variance, norm, b, task_id, n_tasks, sending) +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 - double precision, intent(in) :: pt2(N_states,n_tasks) - double precision, intent(in) :: variance(N_states,n_tasks) - double precision, intent(in) :: norm(N_states,n_tasks) + 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 @@ -358,32 +351,12 @@ subroutine push_pt2_results_async_send(zmq_socket_push, index, pt2, variance, no endif - rc = f77_zmq_send( zmq_socket_push, pt2, 8*N_states*n_tasks, ZMQ_SNDMORE) + rc = f77_zmq_send( zmq_socket_push, pt2_data, pt2_type_size(N_states)*n_tasks, ZMQ_SNDMORE) if (rc == -1) then print *, irp_here, ': error sending result' stop 3 return - else if(rc /= 8*N_states*n_tasks) then - stop 'push' - endif - - - rc = f77_zmq_send( zmq_socket_push, variance, 8*N_states*n_tasks, ZMQ_SNDMORE) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 4 - return - else if(rc /= 8*N_states*n_tasks) then - stop 'push' - endif - - - rc = f77_zmq_send( zmq_socket_push, norm, 8*N_states*n_tasks, ZMQ_SNDMORE) - if (rc == -1) then - print *, irp_here, ': error sending result' - stop 5 - return - else if(rc /= 8*N_states*n_tasks) then + else if(rc /= pt2_type_size(N_states)*n_tasks) then stop 'push' endif @@ -475,7 +448,7 @@ IRP_ELSE stop 11 return else if (rc /= 8) then - print *, irp_here//': error in receiving mini' + print *, irp_here//': error in receiving mini' stop 12 endif IRP_ENDIF @@ -484,14 +457,12 @@ end subroutine -subroutine pull_pt2_results(zmq_socket_pull, index, pt2, variance, norm, task_id, n_tasks, b) +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 - double precision, intent(inout) :: pt2(N_states,*) - double precision, intent(inout) :: variance(N_states,*) - double precision, intent(inout) :: norm(N_states,*) + type(pt2_type), intent(inout) :: pt2_data(n_tasks) type(selection_buffer), intent(inout) :: b integer, intent(out) :: index(*) integer, intent(out) :: n_tasks, task_id(*) @@ -514,27 +485,11 @@ subroutine pull_pt2_results(zmq_socket_pull, index, pt2, variance, norm, task_id stop 'pull' endif - rc = f77_zmq_recv( zmq_socket_pull, pt2, N_states*8*n_tasks, 0) + rc = f77_zmq_recv( zmq_socket_pull, pt2_data, pt2_type_size(N_states)*n_tasks, 0) if (rc == -1) then n_tasks = 1 task_id(1) = 0 - else if(rc /= 8*N_states*n_tasks) then - stop 'pull' - endif - - rc = f77_zmq_recv( zmq_socket_pull, variance, N_states*8*n_tasks, 0) - if (rc == -1) then - n_tasks = 1 - task_id(1) = 0 - else if(rc /= 8*N_states*n_tasks) then - stop 'pull' - endif - - rc = f77_zmq_recv( zmq_socket_pull, norm, N_states*8*n_tasks, 0) - if (rc == -1) then - n_tasks = 1 - task_id(1) = 0 - else if(rc /= 8*N_states*n_tasks) then + else if(rc /= pt2_type_size(N_states)*n_tasks) then stop 'pull' endif diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index cf52aaff..e68e4b79 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -180,15 +180,13 @@ subroutine get_mask_phase(det1, pm, Nint) end subroutine -subroutine select_connected(i_generator,E0,pt2,variance,norm2,b,subset,csubset) +subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset) use bitmasks use selection_types implicit none integer, intent(in) :: i_generator, subset, csubset type(selection_buffer), intent(inout) :: b - double precision, intent(inout) :: pt2(N_states) - double precision, intent(inout) :: variance(N_states) - double precision, intent(inout) :: norm2(N_states) + type(pt2_type), intent(inout) :: pt2_data integer :: k,l double precision, intent(in) :: E0(N_states) @@ -206,7 +204,7 @@ subroutine select_connected(i_generator,E0,pt2,variance,norm2,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 - call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,variance,norm2,b,subset,csubset) + call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,b,subset,csubset) deallocate(fock_diag_tmp) end subroutine @@ -255,7 +253,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,variance,norm2,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 @@ -267,9 +265,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d integer(bit_kind), intent(in) :: hole_mask(N_int,2), particle_mask(N_int,2) double precision, intent(in) :: fock_diag_tmp(mo_num) double precision, intent(in) :: E0(N_states) - double precision, intent(inout) :: pt2(N_states) - double precision, intent(inout) :: variance(N_states) - double precision, intent(inout) :: norm2(N_states) + 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 @@ -645,9 +641,9 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) if(.not.pert_2rdm)then - call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm2, mat, buf) + call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) else - call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm2, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) + call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) endif end if enddo @@ -665,7 +661,7 @@ end subroutine -subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm2, mat, buf) +subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) use bitmasks use selection_types implicit none @@ -673,16 +669,14 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d integer, intent(in) :: i_generator, sp, h1, h2 double precision, intent(in) :: mat(N_states, mo_num, mo_num) logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num) - double precision, intent(in) :: fock_diag_tmp(mo_num) + double precision, intent(in) :: fock_diag_tmp(mo_num) double precision, intent(in) :: E0(N_states) - double precision, intent(inout) :: pt2(N_states) - double precision, intent(inout) :: variance(N_states) - double precision, intent(inout) :: norm2(N_states) + type(pt2_type), intent(inout) :: pt2_data type(selection_buffer), intent(inout) :: buf logical :: ok - integer :: s1, s2, p1, p2, ib, j, istate + integer :: s1, s2, p1, p2, ib, j, istate, jstate integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert, delta_E, val, Hii, w, tmp, alpha_h_psi, coef + double precision :: e_pert, delta_E, val, Hii, w, tmp, alpha_h_psi, coef(N_states) double precision, external :: diag_H_mat_elem_fock double precision :: E_shift @@ -782,16 +776,36 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d endif e_pert = 0.5d0 * (tmp - delta_E) if (dabs(alpha_h_psi) > 1.d-4) then - coef = e_pert / alpha_h_psi + coef(istate) = e_pert / alpha_h_psi else - coef = alpha_h_psi / delta_E + coef(istate) = alpha_h_psi / delta_E endif - pt2(istate) = pt2(istate) + e_pert - variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi - norm2(istate) = norm2(istate) + coef * coef + enddo + + do istate=1,N_states + do jstate=1,N_states + pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate) + enddo + enddo + + ! Gram-Schmidt using input overlap matrix +! do istate=1,N_states +! do jstate=1,istate-1 +! coef(istate) = coef(istate) - pt2_overlap(jstate,istate)/(pt2_overlap(jstate,jstate)) * coef(jstate) +! enddo +! enddo + + do istate=1,N_states + alpha_h_psi = mat(istate, p1, p2) + e_pert = coef(istate) * alpha_h_psi + + pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi + pt2_data % norm2(istate) += coef(istate) * coef(istate) + pt2_data % pt2(istate) += e_pert !!!DEBUG -! pt2(istate) = pt2(istate) - e_pert + alpha_h_psi**2/delta_E +! delta_E = E0(istate) - Hii + E_shift +! pt2_data % pt2(istate) = pt2_data % pt2(istate) + alpha_h_psi**2/delta_E ! ! integer :: k ! double precision :: alpha_h_psi_2,hij @@ -815,7 +829,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d w = w - alpha_h_psi * alpha_h_psi * selection_weight(istate) case(6) - w = w - coef * coef * selection_weight(istate) + w = w - coef(istate) * coef(istate) * selection_weight(istate) case default ! Energy selection diff --git a/src/cipsi/selection_types.f90 b/src/cipsi/selection_types.f90 index 1905f2d3..52b84cf1 100644 --- a/src/cipsi/selection_types.f90 +++ b/src/cipsi/selection_types.f90 @@ -19,6 +19,13 @@ module selection_types double precision, allocatable :: overlap_err(:,:) endtype + contains + + integer function pt2_type_size(N) + implicit none + integer, intent(in) :: N + pt2_type_size = 8*(8*n + 2*n*n) + end function end module From 061e7100ca6c17ac737fadf557cfba2f1733e148 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 29 Aug 2020 11:28:59 +0200 Subject: [PATCH 08/27] Removed eI --- src/cipsi/pt2_stoch_routines.irp.f | 105 ++++++++++++++++------------- src/cipsi/pt2_type.irp.f | 50 ++++++++++---- src/cipsi/run_pt2_slave.irp.f | 2 +- 3 files changed, 97 insertions(+), 60 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 5acd4752..fa1438dc 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -242,8 +242,8 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in) mem_collector = 8.d0 * & ! bytes ( 1.d0*pt2_n_tasks_max & ! task_id, index + 0.635d0*N_det_generators & ! f,d - + 3.d0*N_det_generators*N_states & ! eI, vI, nI + pt2_n_tasks_max*pt2_type_size(N_states)/8 & ! pt2_data_task + + N_det_generators*pt2_type_size(N_states)/8 & ! 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 @@ -360,9 +360,10 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) integer, intent(in) :: N_ type(pt2_type), allocatable :: pt2_data_task(:) - double precision, allocatable :: eI(:,:), Se(:), Se2(:) - double precision, allocatable :: vI(:,:), Sv(:), Sv2(:) - double precision, allocatable :: nI(:,:), Sn(:), Sn2(:) + 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 @@ -370,6 +371,8 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) 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(:) @@ -400,12 +403,9 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) 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(eI(N_states, N_det_generators)) - allocate(vI(N_states, N_det_generators)) - allocate(nI(N_states, N_det_generators)) - allocate(Se(pt2_N_teeth+1), Se2(pt2_N_teeth+1)) - allocate(Sv(pt2_N_teeth+1), Sv2(pt2_N_teeth+1)) - allocate(Sn(pt2_N_teeth+1), Sn2(pt2_N_teeth+1)) + allocate(pt2_data_I(N_det_generators)) + allocate(pt2_data_S(pt2_N_teeth+1)) + allocate(pt2_data_S2(pt2_N_teeth+1)) @@ -417,18 +417,19 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) pt2_data % pt2_err(pt2_stoch_istate) = huge(1.) pt2_data % variance(pt2_stoch_istate) = huge(1.) pt2_data % norm2(pt2_stoch_istate) = 0.d0 - Se(:) = 0d0 - Sv(:) = 0d0 - Sn(:) = 0d0 - Se2(:) = 0d0 - Sv2(:) = 0d0 - Sn2(:) = 0d0 n = 1 t = 0 U = 0 - eI(:,:) = 0d0 - vI(:,:) = 0d0 - nI(:,:) = 0d0 + 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 @@ -456,9 +457,9 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) v0 = 0.d0 n0 = 0.d0 do i=pt2_n_0(t),1,-1 - E0 += eI(pt2_stoch_istate, i) - v0 += vI(pt2_stoch_istate, i) - n0 += nI(pt2_stoch_istate, i) + E0 += pt2_data_I(i) % pt2(pt2_stoch_istate) + v0 += pt2_data_I(i) % variance(pt2_stoch_istate) + n0 += pt2_data_I(i) % norm2(pt2_stoch_istate) end do else exit @@ -468,26 +469,19 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) ! Add Stochastic part c = pt2_R(n) if(c > 0) then -!print *, 'c>0' - x = 0d0 - x2 = 0d0 - x3 = 0d0 + + 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)) - x += eI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i) - x2 += vI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i) - x3 += nI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i) - Se(p) += x - Sv(p) += x2 - Sn(p) += x3 - Se2(p) += x*x - Sv2(p) += x2*x2 - Sn2(p) += x3*3 + call pt2_add ( pt2_data_teeth, pt2_W_T / pt2_w(i), 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 ) end do - avg = E0 + Se(t) / dble(c) - avg2 = v0 + Sv(t) / dble(c) - avg3 = n0 + Sn(t) / dble(c) + 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) % norm2(pt2_stoch_istate) / dble(c) if ((avg /= 0.d0) .or. (n == N_det_generators) ) then do_exit = .true. endif @@ -500,18 +494,19 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) call wall_time(time) ! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969) if(c > 2) then - eqt = dabs((Se2(t) / c) - (Se(t)/c)**2) ! dabs for numerical stability + 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 % pt2_err(pt2_stoch_istate) = eqt - eqt = dabs((Sv2(t) / c) - (Sv(t)/c)**2) ! dabs for numerical stability + 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 % variance_err(pt2_stoch_istate) = eqt - eqt = dabs((Sn2(t) / c) - (Sn(t)/c)**2) ! dabs for numerical stability + eqt = dabs((pt2_data_S2(t) % norm2(pt2_stoch_istate) / c) - (pt2_data_S(t) % norm2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability eqt = sqrt(eqt / (dble(c) - 1.5d0)) pt2_data % norm2_err(pt2_stoch_istate) = eqt + if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, F14.10, 2X, F14.10, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3, time-time0, '' @@ -544,16 +539,22 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) stop 'PT2: Unable to delete tasks (send)' endif do i=1,n_tasks - if(index(i).gt.size(eI,2).or.index(i).lt.1)then + 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*,'If you see this, send a bug report with the following content' print*,irp_here - print*,'i,index(i),size(ei,2) = ',i,index(i),size(ei,2) + print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1) stop -1 endif - eI(1:N_states, index(i)) += pt2_data_task(i) % pt2(1:N_states) - vI(1:N_states, index(i)) += pt2_data_task(i) % variance(1:N_states) - nI(1:N_states, index(i)) += pt2_data_task(i) % norm2(1:N_states) +! print *, pt2_data_I(index(i))%pt2 +! print *, pt2_data_I(index(i))%variance +! print *, pt2_data_I(index(i))%norm2 +! print *, '' +! print *, pt2_data_task(i)%pt2 +! print *, pt2_data_task(i)%variance +! print *, pt2_data_task(i)%norm2 +! print *, '' + 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 @@ -566,6 +567,16 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) 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' diff --git a/src/cipsi/pt2_type.irp.f b/src/cipsi/pt2_type.irp.f index 8f8bb225..af8cf6a7 100644 --- a/src/cipsi/pt2_type.irp.f +++ b/src/cipsi/pt2_type.irp.f @@ -50,25 +50,51 @@ subroutine pt2_dealloc(pt2_data) ) end subroutine -subroutine pt2_add(p1, p2) +subroutine pt2_add(p1, w, p2) implicit none use selection_types BEGIN_DOC -! p1 += p2 +! p1 += w * p2 END_DOC type(pt2_type), intent(inout) :: p1 + double precision, intent(in) :: w type(pt2_type), intent(in) :: p2 - p1 % pt2(:) += p2 % pt2(:) - p1 % pt2_err(:) += p2 % pt2_err(:) - p1 % rpt2(:) += p2 % rpt2(:) - p1 % rpt2_err(:) += p2 % rpt2_err(:) - p1 % variance(:) += p2 % variance(:) - p1 % variance_err(:) += p2 % variance_err(:) - p1 % norm2(:) += p2 % norm2(:) - p1 % norm2_err(:) += p2 % norm2_err(:) - p1 % overlap(:,:) += p2 % overlap(:,:) - p1 % overlap_err(:,:) += p2 % overlap_err(:,:) + p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) + p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:) + p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) + p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:) + p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) + p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:) + p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) + p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:) + p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) + p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) + +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 + + p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) * p2 % pt2(:) + p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:) * p2 % pt2_err(:) + p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) * p2 % rpt2(:) + p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:) * p2 % rpt2_err(:) + p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) * p2 % variance(:) + p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:) * p2 % variance_err(:) + p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) * p2 % norm2(:) + p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:) * p2 % norm2_err(:) + p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) * p2 % overlap(:,:) + p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) * p2 % overlap_err(:,:) + end subroutine diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi/run_pt2_slave.irp.f index ae0ac2a0..1859ca88 100644 --- a/src/cipsi/run_pt2_slave.irp.f +++ b/src/cipsi/run_pt2_slave.irp.f @@ -462,7 +462,7 @@ subroutine pull_pt2_results(zmq_socket_pull, index, pt2_data, task_id, n_tasks, use f77_zmq implicit none integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - type(pt2_type), intent(inout) :: pt2_data(n_tasks) + type(pt2_type), intent(inout) :: pt2_data(*) type(selection_buffer), intent(inout) :: b integer, intent(out) :: index(*) integer, intent(out) :: n_tasks, task_id(*) From 93fc49000cdf1c059204790562abf1e86207ddc8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 30 Aug 2020 22:16:39 +0200 Subject: [PATCH 09/27] Fixed PT2 --- src/cipsi/pt2_stoch_routines.irp.f | 14 +--- src/cipsi/pt2_type.irp.f | 123 +++++++++++++++++++++++----- src/cipsi/run_pt2_slave.irp.f | 60 +++++++------- src/cipsi/run_selection_slave.irp.f | 89 +++++++++----------- src/cipsi/selection_types.f90 | 2 +- src/cipsi/zmq_selection.irp.f | 22 +++-- 6 files changed, 185 insertions(+), 125 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index fa1438dc..4b781dd8 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -474,10 +474,12 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) 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)) - call pt2_add ( pt2_data_teeth, pt2_W_T / pt2_w(i), pt2_data_I(i) ) + 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 ) - end do + 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) @@ -546,14 +548,6 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1) stop -1 endif -! print *, pt2_data_I(index(i))%pt2 -! print *, pt2_data_I(index(i))%variance -! print *, pt2_data_I(index(i))%norm2 -! print *, '' -! print *, pt2_data_task(i)%pt2 -! print *, pt2_data_task(i)%variance -! print *, pt2_data_task(i)%norm2 -! print *, '' call pt2_add(pt2_data_I(index(i)),1.d0,pt2_data_task(i)) f(index(i)) -= 1 end do diff --git a/src/cipsi/pt2_type.irp.f b/src/cipsi/pt2_type.irp.f index af8cf6a7..e6f31799 100644 --- a/src/cipsi/pt2_type.irp.f +++ b/src/cipsi/pt2_type.irp.f @@ -60,16 +60,33 @@ subroutine pt2_add(p1, w, p2) double precision, intent(in) :: w type(pt2_type), intent(in) :: p2 - p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) - p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:) - p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) - p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:) - p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) - p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:) - p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) - p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:) - p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) - p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) + if (w == 1.d0) then + + p1 % pt2(:) = p1 % pt2(:) + p2 % pt2(:) + p1 % pt2_err(:) = p1 % pt2_err(:) + p2 % pt2_err(:) + p1 % rpt2(:) = p1 % rpt2(:) + p2 % rpt2(:) + p1 % rpt2_err(:) = p1 % rpt2_err(:) + p2 % rpt2_err(:) + p1 % variance(:) = p1 % variance(:) + p2 % variance(:) + p1 % variance_err(:) = p1 % variance_err(:) + p2 % variance_err(:) + p1 % norm2(:) = p1 % norm2(:) + p2 % norm2(:) + p1 % norm2_err(:) = p1 % norm2_err(:) + p2 % norm2_err(:) + p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:) + p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + p2 % overlap_err(:,:) + + else + + p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) + p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:) + p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) + p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:) + p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) + p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:) + p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) + p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:) + p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) + p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) + + endif end subroutine @@ -84,17 +101,83 @@ subroutine pt2_add2(p1, w, p2) double precision, intent(in) :: w type(pt2_type), intent(in) :: p2 - p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) * p2 % pt2(:) - p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:) * p2 % pt2_err(:) - p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) * p2 % rpt2(:) - p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:) * p2 % rpt2_err(:) - p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) * p2 % variance(:) - p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:) * p2 % variance_err(:) - p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) * p2 % norm2(:) - p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:) * p2 % norm2_err(:) - p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) * p2 % overlap(:,:) - p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) * p2 % overlap_err(:,:) + if (w == 1.d0) then + + p1 % pt2(:) = p1 % pt2(:) + p2 % pt2(:) * p2 % pt2(:) + p1 % pt2_err(:) = p1 % pt2_err(:) + p2 % pt2_err(:) * p2 % pt2_err(:) + p1 % rpt2(:) = p1 % rpt2(:) + p2 % rpt2(:) * p2 % rpt2(:) + p1 % rpt2_err(:) = p1 % rpt2_err(:) + p2 % rpt2_err(:) * p2 % rpt2_err(:) + p1 % variance(:) = p1 % variance(:) + p2 % variance(:) * p2 % variance(:) + p1 % variance_err(:) = p1 % variance_err(:) + p2 % variance_err(:) * p2 % variance_err(:) + p1 % norm2(:) = p1 % norm2(:) + p2 % norm2(:) * p2 % norm2(:) + p1 % norm2_err(:) = p1 % norm2_err(:) + p2 % norm2_err(:) * p2 % norm2_err(:) + p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:) * p2 % overlap(:,:) + p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + p2 % overlap_err(:,:) * p2 % overlap_err(:,:) + + else + + p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) * p2 % pt2(:) + p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:) * p2 % pt2_err(:) + p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) * p2 % rpt2(:) + p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:) * p2 % rpt2_err(:) + p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) * p2 % variance(:) + p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:) * p2 % variance_err(:) + p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) * p2 % norm2(:) + p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:) * p2 % norm2_err(:) + p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) * p2 % overlap(:,:) + p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) * p2 % overlap_err(:,:) + + 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) + x(n+1:2*n) = pt2_data % pt2_err(1:n) + x(2*n+1:3*n) = pt2_data % rpt2(1:n) + x(3*n+1:4*n) = pt2_data % rpt2_err(1:n) + x(4*n+1:5*n) = pt2_data % variance(1:n) + x(5*n+1:6*n) = pt2_data % variance_err(1:n) + x(6*n+1:7*n) = pt2_data % norm2(1:n) + x(7*n+1:8*n) = pt2_data % norm2_err(1:n) + k=8*n + x(k+1:k+n2) = reshape(pt2_data % overlap(1:n,1:n), (/ n2 /)) + k=8*n+n2 + x(k+1:k+n2) = reshape(pt2_data % overlap_err(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) + pt2_data % pt2_err(1:n) = x(n+1:2*n) + pt2_data % rpt2(1:n) = x(2*n+1:3*n) + pt2_data % rpt2_err(1:n) = x(3*n+1:4*n) + pt2_data % variance(1:n) = x(4*n+1:5*n) + pt2_data % variance_err(1:n) = x(5*n+1:6*n) + pt2_data % norm2(1:n) = x(6*n+1:7*n) + pt2_data % norm2_err(1:n) = x(7*n+1:8*n) + k=8*n + pt2_data % overlap(1:n,1:n) = reshape(x(k+1:k+n2), (/ n, n /)) + k=8*n+n2 + pt2_data % overlap_err(1:n,1:n) = reshape(x(k+1:k+n2), (/ n, n /)) + +end diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi/run_pt2_slave.irp.f index 1859ca88..3c72dac0 100644 --- a/src/cipsi/run_pt2_slave.irp.f +++ b/src/cipsi/run_pt2_slave.irp.f @@ -141,6 +141,7 @@ subroutine run_pt2_slave_small(thread,iproc,energy) ! ! 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 @@ -169,8 +170,8 @@ subroutine run_pt2_slave_large(thread,iproc,energy) integer :: rc, i integer :: worker_id, ctask, ltask - character*(512), allocatable :: task(:) - integer, allocatable :: task_id(:) + character*(512) :: task + integer :: task_id integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket integer(ZMQ_PTR) :: zmq_to_qp_run_socket @@ -181,18 +182,15 @@ subroutine run_pt2_slave_large(thread,iproc,energy) type(selection_buffer) :: b logical :: done, buffer_ready - type(pt2_type), allocatable :: pt2_data(:) + type(pt2_type) :: pt2_data integer :: n_tasks, k, N - integer, allocatable :: i_generator(:), subset(:) + integer :: i_generator, subset integer :: bsize ! Size of selection buffers logical :: sending PROVIDE global_selection_buffer global_selection_buffer_lock - 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 @@ -211,22 +209,17 @@ subroutine run_pt2_slave_large(thread,iproc,energy) 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 + done = task_id == 0 if (done) then n_tasks = n_tasks-1 endif if (n_tasks == 0) exit - do k=1,n_tasks - read (task(k),*) subset(k), i_generator(k), N - enddo + read (task,*) subset, i_generator, N if (b%N == 0) then ! Only first time bsize = min(N, (elec_alpha_num * (mo_num-elec_alpha_num))**2) @@ -238,15 +231,13 @@ subroutine run_pt2_slave_large(thread,iproc,energy) 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 pt2_alloc(pt2_data,N_states) + b%cur = 0 !double precision :: time2 !call wall_time(time2) - call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k))) + call select_connected(i_generator,energy,pt2_data,b,subset,pt2_F(i_generator)) !call wall_time(time1) !print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1)) - enddo call wall_time(time1) !print *, '-->', i_generator(1), time1-time0, n_tasks @@ -270,12 +261,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy) call push_pt2_results_async_send(zmq_socket_push, i_generator, pt2_data, b, task_id, n_tasks,sending) endif - do k=1,n_tasks - call pt2_dealloc(pt2_data(k)) - enddo -! ! 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 = 1 + call pt2_dealloc(pt2_data) end do call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending) @@ -322,8 +308,9 @@ subroutine push_pt2_results_async_send(zmq_socket_push, index, pt2_data, b, task integer, intent(in) :: n_tasks, index(n_tasks), task_id(n_tasks) type(selection_buffer), intent(inout) :: b logical, intent(inout) :: sending - integer :: rc + integer :: rc, i integer*8 :: rc8 + double precision, allocatable :: pt2_serialized(:,:) if (sending) then print *, irp_here, ': sending is true' @@ -351,12 +338,18 @@ subroutine push_pt2_results_async_send(zmq_socket_push, index, pt2_data, b, task endif - rc = f77_zmq_send( zmq_socket_push, pt2_data, pt2_type_size(N_states)*n_tasks, ZMQ_SNDMORE) + 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 /= pt2_type_size(N_states)*n_tasks) then + else if(rc /= size(pt2_serialized)*8) then stop 'push' endif @@ -468,6 +461,7 @@ subroutine pull_pt2_results(zmq_socket_pull, index, pt2_data, task_id, n_tasks, 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 @@ -485,14 +479,20 @@ subroutine pull_pt2_results(zmq_socket_pull, index, pt2_data, task_id, n_tasks, stop 'pull' endif - rc = f77_zmq_recv( zmq_socket_pull, pt2_data, pt2_type_size(N_states)*n_tasks, 0) + 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 /= pt2_type_size(N_states)*n_tasks) then + 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 diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index fe712c45..69a8a4c3 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -18,9 +18,7 @@ subroutine run_selection_slave(thread,iproc,energy) type(selection_buffer) :: buf, buf2 logical :: done, buffer_ready - double precision :: pt2(N_states) - double precision :: variance(N_states) - double precision :: norm(N_states) + 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 @@ -28,6 +26,7 @@ subroutine run_selection_slave(thread,iproc,energy) 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() @@ -42,9 +41,6 @@ subroutine run_selection_slave(thread,iproc,energy) buf%N = 0 buffer_ready = .False. ctask = 1 - pt2(:) = 0d0 - variance(:) = 0d0 - norm(:) = 0.d0 do integer, external :: get_task_from_taskserver @@ -69,7 +65,7 @@ subroutine run_selection_slave(thread,iproc,energy) stop '-1' end if end if - call select_connected(i_generator,energy,pt2,variance,norm,buf,subset,pt2_F(i_generator)) + call select_connected(i_generator,energy,pt2_data,buf,subset,pt2_F(i_generator)) endif integer, external :: task_done_to_taskserver @@ -88,12 +84,10 @@ subroutine run_selection_slave(thread,iproc,energy) if(ctask > 0) then call sort_selection_buffer(buf) ! call merge_selection_buffers(buf,buf2) -!print *, task_id(1), pt2(1), buf%cur, ctask - call push_selection_results(zmq_socket_push, pt2, variance, norm, buf, task_id(1), ctask) + 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 - pt2(:) = 0d0 - variance(:) = 0d0 - norm(:) = 0d0 buf%cur = 0 end if ctask = 0 @@ -106,11 +100,9 @@ subroutine run_selection_slave(thread,iproc,energy) if(ctask > 0) then call sort_selection_buffer(buf) ! call merge_selection_buffers(buf,buf2) - call push_selection_results(zmq_socket_push, pt2, variance, norm, buf, task_id(1), ctask) + call push_selection_results(zmq_socket_push, pt2_data, buf, task_id(1), ctask) + call pt2_dealloc(pt2_data) ! buf%mini = buf2%mini - pt2(:) = 0d0 - variance(:) = 0d0 - norm(:) = 0d0 buf%cur = 0 end if ctask = 0 @@ -129,18 +121,17 @@ subroutine run_selection_slave(thread,iproc,energy) end subroutine -subroutine push_selection_results(zmq_socket_push, pt2, variance, norm, b, task_id, ntask) +subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask) use f77_zmq use selection_types implicit none integer(ZMQ_PTR), intent(in) :: zmq_socket_push - double precision, intent(in) :: pt2(N_states) - double precision, intent(in) :: variance(N_states) - double precision, intent(in) :: norm(N_states) + type(pt2_type), intent(in) :: pt2_data type(selection_buffer), intent(inout) :: b integer, intent(in) :: ntask, 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 @@ -148,19 +139,19 @@ subroutine push_selection_results(zmq_socket_push, pt2, variance, norm, b, task_ endif - rc = f77_zmq_send( zmq_socket_push, pt2, 8*N_states, ZMQ_SNDMORE) - if(rc /= 8*N_states) then - print *, 'f77_zmq_send( zmq_socket_push, pt2, 8*N_states, ZMQ_SNDMORE)' - 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, variance, 8*N_states, ZMQ_SNDMORE) - if(rc /= 8*N_states) then - print *, 'f77_zmq_send( zmq_socket_push, variance, 8*N_states, ZMQ_SNDMORE)' - endif - - rc = f77_zmq_send( zmq_socket_push, norm, 8*N_states, ZMQ_SNDMORE) - if(rc /= 8*N_states) then - print *, 'f77_zmq_send( zmq_socket_push, norm, 8*N_states, ZMQ_SNDMORE)' + 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 if (b%cur > 0) then @@ -201,42 +192,36 @@ IRP_ENDIF end subroutine -subroutine pull_selection_results(zmq_socket_pull, pt2, variance, norm, val, det, N, task_id, ntask) +subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_id, ntask) use f77_zmq use selection_types implicit none integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - double precision, intent(inout) :: pt2(N_states) - double precision, intent(inout) :: variance(N_states) - double precision, intent(inout) :: norm(N_states) + 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, ntask, 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 - pt2(:) = 0.d0 - variance(:) = 0.d0 - norm(:) = 0.d0 - - rc = f77_zmq_recv( zmq_socket_pull, pt2, N_states*8, 0) - if(rc /= 8*N_states) then - print *, 'f77_zmq_recv( zmq_socket_pull, pt2, N_states*8, 0)' + 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 - rc = f77_zmq_recv( zmq_socket_pull, variance, N_states*8, 0) - if(rc /= 8*N_states) then - print *, 'f77_zmq_recv( zmq_socket_pull, variance, N_states*8, 0)' - endif - - rc = f77_zmq_recv( zmq_socket_pull, norm, N_states*8, 0) - if(rc /= 8*N_states) then - print *, 'f77_zmq_recv( zmq_socket_pull, norm, N_states*8, 0)' - endif + do i=1,n_tasks + call pt2_deserialize(pt2_data(i),N_states,pt2_serialized(1,i)) + enddo + deallocate(pt2_serialized) if (N>0) then rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, 0) diff --git a/src/cipsi/selection_types.f90 b/src/cipsi/selection_types.f90 index 52b84cf1..eef57aa5 100644 --- a/src/cipsi/selection_types.f90 +++ b/src/cipsi/selection_types.f90 @@ -24,7 +24,7 @@ module selection_types integer function pt2_type_size(N) implicit none integer, intent(in) :: N - pt2_type_size = 8*(8*n + 2*n*n) + pt2_type_size = (8*n + 2*n*n) end function end module diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 006e6578..789c7a26 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -129,12 +129,12 @@ subroutine ZMQ_selection(N_in, pt2_data) 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) - pt2_data % norm2(k) = pt2_data % norm2(k) * f(k) + pt2_data % pt2(k) = pt2_data % pt2(k) * f(k) + pt2_data % variance(k) = pt2_data % variance(k) * f(k) + pt2_data % norm2(k) = pt2_data % norm2(k) * f(k) - pt2_data % rpt2(k) = & - pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) + pt2_data % rpt2(k) = & + pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo call update_pt2_and_variance_weights(pt2_data, N_states) @@ -160,6 +160,7 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2_data) 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) @@ -190,15 +191,11 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2_data) pt2_data % pt2(:) = 0d0 pt2_data % variance(:) = 0.d0 pt2_data % norm2(:) = 0.d0 - pt2_mwen(:) = 0.d0 - variance_mwen(:) = 0.d0 - norm2_mwen(:) = 0.d0 + call pt2_alloc(pt2_data_tmp,N_states) do while (more == 1) - call pull_selection_results(zmq_socket_pull, pt2_mwen, variance_mwen, norm2_mwen, b2%val(1), b2%det(1,1,1), b2%cur, task_id, ntask) + call pull_selection_results(zmq_socket_pull, pt2_data_tmp, b2%val(1), b2%det(1,1,1), b2%cur, task_id, ntask) - pt2_data % pt2(:) += pt2_mwen(:) - pt2_data % variance(:) += variance_mwen(:) - pt2_data % norm2(:) += norm2_mwen(:) + 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 @@ -214,6 +211,7 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2_data) endif end do end do + call pt2_dealloc(pt2_data_tmp) call delete_selection_buffer(b2) From e9f1d7576a238fc8f82e4c81fa6b9468471e2419 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 Aug 2020 01:45:36 +0200 Subject: [PATCH 10/27] Removed error from type --- src/cipsi/cipsi.irp.f | 29 +++++++----- src/cipsi/pt2_stoch_routines.irp.f | 28 ++++++----- src/cipsi/pt2_type.irp.f | 73 +++++++---------------------- src/cipsi/run_pt2_slave.irp.f | 23 ++++----- src/cipsi/run_selection_slave.irp.f | 52 ++++++++++---------- src/cipsi/selection_types.f90 | 7 +-- src/cipsi/stochastic_cipsi.irp.f | 26 +++++----- src/fci/pt2.irp.f | 15 ++++-- src/iterations/print_summary.irp.f | 28 +++++------ 9 files changed, 125 insertions(+), 156 deletions(-) diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 44b87394..0f140240 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -6,7 +6,7 @@ subroutine run_cipsi ! stochastic PT2. END_DOC integer :: i,j,k - type(pt2_type) :: pt2_data + type(pt2_type) :: pt2_data, pt2_data_err double precision, allocatable :: zeros(:) integer :: to_select logical, external :: qp_stop @@ -25,6 +25,7 @@ subroutine run_cipsi allocate (zeros(N_states)) call pt2_alloc(pt2_data, N_states) + call pt2_alloc(pt2_data_err, N_states) double precision :: hf_energy_ref logical :: has @@ -79,16 +80,19 @@ subroutine run_cipsi to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) if (do_pt2) then - pt2_data % pt2 = 0.d0 - pt2_data % variance = 0.d0 - pt2_data % norm2 = 0.d0 + 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) threshold_generators_save = threshold_generators threshold_generators = 1.d0 SOFT_TOUCH threshold_generators - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error, 0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,pt2_data_err,relative_error, 0) ! Stochastic PT2 threshold_generators = threshold_generators_save SOFT_TOUCH threshold_generators else + call pt2_dealloc(pt2_data) + call pt2_alloc(pt2_data, N_states) call ZMQ_selection(to_select, pt2_data) endif @@ -98,7 +102,7 @@ subroutine run_cipsi call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep, & - pt2_data, N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) @@ -130,12 +134,13 @@ subroutine run_cipsi endif if (do_pt2) then - pt2_data % pt2(:) = 0.d0 - pt2_data % variance(:) = 0.d0 - pt2_data % norm2(:) = 0.d0 + 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) threshold_generators = 1d0 SOFT_TOUCH threshold_generators - call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, relative_error, 0) ! Stochastic PT2 SOFT_TOUCH threshold_generators endif print *, 'N_det = ', N_det @@ -145,9 +150,11 @@ subroutine run_cipsi call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) call print_summary(psi_energy_with_nucl_rep(1:N_states), & - pt2_data, N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call print_extrapolated_energy() endif + call pt2_dealloc(pt2_data) + call pt2_dealloc(pt2_data_err) end diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 4b781dd8..5b05a5a9 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -107,7 +107,7 @@ end function -subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in) +subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) use f77_zmq use selection_types @@ -117,7 +117,7 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in) 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 + type(pt2_type), intent(inout) :: pt2_data, pt2_data_err ! integer :: i, N @@ -298,13 +298,13 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in) i = omp_get_thread_num() if (i==0) then - call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, pt2_data, b, N) + 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 % norm2(pt2_stoch_istate)) !TODO : We should use here the correct formula for the error of X/Y - pt2_data % rpt2_err(pt2_stoch_istate) = & - pt2_data % pt2_err(pt2_stoch_istate)/(1.d0 + pt2_data % norm2(pt2_stoch_istate)) + pt2_data_err % rpt2(pt2_stoch_istate) = & + pt2_data_err % pt2(pt2_stoch_istate)/(1.d0 + pt2_data % norm2(pt2_stoch_istate)) else call pt2_slave_inproc(i) @@ -346,7 +346,7 @@ subroutine pt2_slave_inproc(i) end -subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) +subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_err, b, N_) use f77_zmq use selection_types use bitmasks @@ -355,7 +355,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) integer(ZMQ_PTR), intent(in) :: zmq_socket_pull double precision, intent(in) :: relative_error, E - type(pt2_type), intent(inout) :: pt2_data + type(pt2_type), intent(inout) :: pt2_data, pt2_data_err type(selection_buffer), intent(inout) :: b integer, intent(in) :: N_ @@ -414,9 +414,11 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) pt2_data % pt2(pt2_stoch_istate) = -huge(1.) - pt2_data % pt2_err(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 % norm2(pt2_stoch_istate) = 0.d0 + pt2_data_err % norm2(pt2_stoch_istate) = huge(1.) n = 1 t = 0 U = 0 @@ -479,8 +481,8 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) 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) % norm2(pt2_stoch_istate) / dble(c) @@ -498,22 +500,22 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_) 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 % pt2_err(pt2_stoch_istate) = eqt + 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 % variance_err(pt2_stoch_istate) = eqt + pt2_data_err % variance(pt2_stoch_istate) = eqt eqt = dabs((pt2_data_S2(t) % norm2(pt2_stoch_istate) / c) - (pt2_data_S(t) % norm2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability eqt = sqrt(eqt / (dble(c) - 1.5d0)) - pt2_data % norm2_err(pt2_stoch_istate) = eqt + pt2_data_err % norm2(pt2_stoch_istate) = eqt if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, F14.10, 2X, F14.10, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3, time-time0, '' if (stop_now .or. ( & - (do_exit .and. (dabs(pt2_data % pt2_err(pt2_stoch_istate)) / & + (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) diff --git a/src/cipsi/pt2_type.irp.f b/src/cipsi/pt2_type.irp.f index e6f31799..ef3fb882 100644 --- a/src/cipsi/pt2_type.irp.f +++ b/src/cipsi/pt2_type.irp.f @@ -6,27 +6,17 @@ subroutine pt2_alloc(pt2_data,N) integer :: k allocate(pt2_data % pt2(N) & - ,pt2_data % pt2_err(N) & ,pt2_data % variance(N) & - ,pt2_data % variance_err(N) & ,pt2_data % norm2(N) & - ,pt2_data % norm2_err(N) & ,pt2_data % rpt2(N) & - ,pt2_data % rpt2_err(N) & ,pt2_data % overlap(N,N) & - ,pt2_data % overlap_err(N,N) & ) pt2_data % pt2(:) = 0.d0 - pt2_data % pt2_err(:) = 0.d0 pt2_data % variance(:) = 0.d0 - pt2_data % variance_err(:) = 0.d0 pt2_data % norm2(:) = 0.d0 - pt2_data % norm2_err(:) = 0.d0 pt2_data % rpt2(:) = 0.d0 - pt2_data % rpt2_err(:) = 0.d0 pt2_data % overlap(:,:) = 0.d0 - pt2_data % overlap_err(:,:) = 0.d0 do k=1,N pt2_data % overlap(k,k) = 1.d0 @@ -38,15 +28,10 @@ subroutine pt2_dealloc(pt2_data) use selection_types type(pt2_type), intent(inout) :: pt2_data deallocate(pt2_data % pt2 & - ,pt2_data % pt2_err & ,pt2_data % variance & - ,pt2_data % variance_err& ,pt2_data % norm2 & - ,pt2_data % norm2_err & ,pt2_data % rpt2 & - ,pt2_data % rpt2_err & ,pt2_data % overlap & - ,pt2_data % overlap_err & ) end subroutine @@ -63,28 +48,18 @@ subroutine pt2_add(p1, w, p2) if (w == 1.d0) then p1 % pt2(:) = p1 % pt2(:) + p2 % pt2(:) - p1 % pt2_err(:) = p1 % pt2_err(:) + p2 % pt2_err(:) p1 % rpt2(:) = p1 % rpt2(:) + p2 % rpt2(:) - p1 % rpt2_err(:) = p1 % rpt2_err(:) + p2 % rpt2_err(:) p1 % variance(:) = p1 % variance(:) + p2 % variance(:) - p1 % variance_err(:) = p1 % variance_err(:) + p2 % variance_err(:) p1 % norm2(:) = p1 % norm2(:) + p2 % norm2(:) - p1 % norm2_err(:) = p1 % norm2_err(:) + p2 % norm2_err(:) p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:) - p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + p2 % overlap_err(:,:) else p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) - p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:) p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) - p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:) p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) - p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:) p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) - p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:) p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) - p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) endif @@ -104,28 +79,18 @@ subroutine pt2_add2(p1, w, p2) if (w == 1.d0) then p1 % pt2(:) = p1 % pt2(:) + p2 % pt2(:) * p2 % pt2(:) - p1 % pt2_err(:) = p1 % pt2_err(:) + p2 % pt2_err(:) * p2 % pt2_err(:) p1 % rpt2(:) = p1 % rpt2(:) + p2 % rpt2(:) * p2 % rpt2(:) - p1 % rpt2_err(:) = p1 % rpt2_err(:) + p2 % rpt2_err(:) * p2 % rpt2_err(:) p1 % variance(:) = p1 % variance(:) + p2 % variance(:) * p2 % variance(:) - p1 % variance_err(:) = p1 % variance_err(:) + p2 % variance_err(:) * p2 % variance_err(:) p1 % norm2(:) = p1 % norm2(:) + p2 % norm2(:) * p2 % norm2(:) - p1 % norm2_err(:) = p1 % norm2_err(:) + p2 % norm2_err(:) * p2 % norm2_err(:) p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:) * p2 % overlap(:,:) - p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + p2 % overlap_err(:,:) * p2 % overlap_err(:,:) else p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) * p2 % pt2(:) - p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:) * p2 % pt2_err(:) p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) * p2 % rpt2(:) - p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:) * p2 % rpt2_err(:) p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) * p2 % variance(:) - p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:) * p2 % variance_err(:) p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) * p2 % norm2(:) - p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:) * p2 % norm2_err(:) p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) * p2 % overlap(:,:) - p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) * p2 % overlap_err(:,:) endif @@ -142,18 +107,15 @@ subroutine pt2_serialize(pt2_data, n, x) integer :: i,k,n2 n2 = n*n - x(1:n) = pt2_data % pt2(1:n) - x(n+1:2*n) = pt2_data % pt2_err(1:n) - x(2*n+1:3*n) = pt2_data % rpt2(1:n) - x(3*n+1:4*n) = pt2_data % rpt2_err(1:n) - x(4*n+1:5*n) = pt2_data % variance(1:n) - x(5*n+1:6*n) = pt2_data % variance_err(1:n) - x(6*n+1:7*n) = pt2_data % norm2(1:n) - x(7*n+1:8*n) = pt2_data % norm2_err(1:n) - k=8*n - x(k+1:k+n2) = reshape(pt2_data % overlap(1:n,1:n), (/ n2 /)) - k=8*n+n2 - x(k+1:k+n2) = reshape(pt2_data % overlap_err(1:n,1:n), (/ n2 /)) + 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+n) = pt2_data % norm2(1:n) + k=k+n + x(k+1:k+n2) = reshape(pt2_data % overlap(1:n,1:n), (/ n2 /)) end @@ -168,16 +130,13 @@ subroutine pt2_deserialize(pt2_data, n, x) n2 = n*n pt2_data % pt2(1:n) = x(1:n) - pt2_data % pt2_err(1:n) = x(n+1:2*n) - pt2_data % rpt2(1:n) = x(2*n+1:3*n) - pt2_data % rpt2_err(1:n) = x(3*n+1:4*n) - pt2_data % variance(1:n) = x(4*n+1:5*n) - pt2_data % variance_err(1:n) = x(5*n+1:6*n) - pt2_data % norm2(1:n) = x(6*n+1:7*n) - pt2_data % norm2_err(1:n) = x(7*n+1:8*n) - k=8*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 % norm2(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 /)) - k=8*n+n2 - pt2_data % overlap_err(1:n,1:n) = reshape(x(k+1:k+n2), (/ n, n /)) end diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi/run_pt2_slave.irp.f index 3c72dac0..d3f4d45d 100644 --- a/src/cipsi/run_pt2_slave.irp.f +++ b/src/cipsi/run_pt2_slave.irp.f @@ -117,11 +117,11 @@ subroutine run_pt2_slave_small(thread,iproc,energy) 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 pt2_alloc(pt2_data(k),N_states) + b%cur = 0 !double precision :: time2 !call wall_time(time2) - call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k))) + call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k))) !call wall_time(time1) !print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1)) enddo @@ -157,6 +157,7 @@ subroutine run_pt2_slave_small(thread,iproc,energy) if (buffer_ready) then call delete_selection_buffer(b) endif + deallocate(pt2_data) end subroutine @@ -171,7 +172,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy) integer :: worker_id, ctask, ltask character*(512) :: task - integer :: task_id + integer :: task_id(1) integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket integer(ZMQ_PTR) :: zmq_to_qp_run_socket @@ -182,9 +183,9 @@ subroutine run_pt2_slave_large(thread,iproc,energy) type(selection_buffer) :: b logical :: done, buffer_ready - type(pt2_type) :: pt2_data + type(pt2_type) :: pt2_data(1) integer :: n_tasks, k, N - integer :: i_generator, subset + integer :: i_generator(1), subset integer :: bsize ! Size of selection buffers logical :: sending @@ -213,13 +214,13 @@ subroutine run_pt2_slave_large(thread,iproc,energy) if (get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task, n_tasks) == -1) then exit endif - done = task_id == 0 + done = task_id(1) == 0 if (done) then n_tasks = n_tasks-1 endif if (n_tasks == 0) exit - read (task,*) subset, i_generator, N + read (task,*) subset, i_generator(1), N if (b%N == 0) then ! Only first time bsize = min(N, (elec_alpha_num * (mo_num-elec_alpha_num))**2) @@ -231,11 +232,11 @@ subroutine run_pt2_slave_large(thread,iproc,energy) double precision :: time0, time1 call wall_time(time0) - call pt2_alloc(pt2_data,N_states) + call pt2_alloc(pt2_data(1),N_states) b%cur = 0 !double precision :: time2 !call wall_time(time2) - call select_connected(i_generator,energy,pt2_data,b,subset,pt2_F(i_generator)) + call select_connected(i_generator(1),energy,pt2_data(1),b,subset,pt2_F(i_generator(1))) !call wall_time(time1) !print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1)) call wall_time(time1) @@ -261,7 +262,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy) call push_pt2_results_async_send(zmq_socket_push, i_generator, pt2_data, b, task_id, n_tasks,sending) endif - call pt2_dealloc(pt2_data) + call pt2_dealloc(pt2_data(1)) end do call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending) diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index 69a8a4c3..650654be 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -101,11 +101,11 @@ subroutine run_selection_slave(thread,iproc,energy) 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) ! 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 @@ -121,7 +121,7 @@ subroutine run_selection_slave(thread,iproc,energy) end subroutine -subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask) +subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntasks) use f77_zmq use selection_types implicit none @@ -129,9 +129,9 @@ subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask) integer(ZMQ_PTR), intent(in) :: zmq_socket_push type(pt2_type), intent(in) :: pt2_data type(selection_buffer), intent(inout) :: b - integer, intent(in) :: ntask, task_id(*) + integer, intent(in) :: ntasks, task_id(*) integer :: rc - double precision, allocatable :: pt2_serialized(:,:) + double precision, allocatable :: pt2_serialized(:) rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE) if(rc /= 4) then @@ -139,13 +139,10 @@ subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask) 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 + 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) - deallocate(pt2_serialized) if (rc == -1) then print *, irp_here, ': error sending result' stop 3 @@ -153,6 +150,7 @@ subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask) else if(rc /= size(pt2_serialized)*8) then stop 'push' endif + deallocate(pt2_serialized) if (b%cur > 0) then @@ -168,14 +166,14 @@ subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask) endif - rc = f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE) + rc = f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE) if(rc /= 4) then - print *, 'f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE)' + print *, 'f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE)' endif - rc = f77_zmq_send( zmq_socket_push, task_id(1), ntask*4, 0) - if(rc /= 4*ntask) then - print *, 'f77_zmq_send( zmq_socket_push, task_id(1), ntask*4, 0)' + 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 @@ -192,7 +190,7 @@ IRP_ENDIF end subroutine -subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_id, ntask) +subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_id, ntasks) use f77_zmq use selection_types implicit none @@ -200,27 +198,25 @@ subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_i 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, ntask, task_id(*) + integer, intent(out) :: N, ntasks, task_id(*) integer :: rc, rn, i - double precision, allocatable :: pt2_serialized(:,:) + 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),n_tasks) ) - rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized)*n_tasks, 0) + allocate(pt2_serialized (pt2_type_size(N_states)) ) + rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized)*ntasks, 0) if (rc == -1) then - n_tasks = 1 + ntasks = 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 + call pt2_deserialize(pt2_data,N_states,pt2_serialized) deallocate(pt2_serialized) if (N>0) then @@ -235,14 +231,14 @@ subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_i endif endif - rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0) + rc = f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0) if(rc /= 4) then - print *, 'f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)' + print *, 'f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0)' endif - rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0) - if(rc /= 4*ntask) then - print *, 'f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0)' + 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 diff --git a/src/cipsi/selection_types.f90 b/src/cipsi/selection_types.f90 index eef57aa5..8df37653 100644 --- a/src/cipsi/selection_types.f90 +++ b/src/cipsi/selection_types.f90 @@ -8,15 +8,10 @@ module selection_types type pt2_type double precision, allocatable :: pt2(:) - double precision, allocatable :: pt2_err(:) double precision, allocatable :: rpt2(:) - double precision, allocatable :: rpt2_err(:) double precision, allocatable :: variance(:) - double precision, allocatable :: variance_err(:) double precision, allocatable :: norm2(:) - double precision, allocatable :: norm2_err(:) double precision, allocatable :: overlap(:,:) - double precision, allocatable :: overlap_err(:,:) endtype contains @@ -24,7 +19,7 @@ module selection_types integer function pt2_type_size(N) implicit none integer, intent(in) :: N - pt2_type_size = (8*n + 2*n*n) + pt2_type_size = (4*n + n*n) end function end module diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index e80bd856..953314b9 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -7,7 +7,7 @@ subroutine run_stochastic_cipsi integer :: i,j,k double precision, allocatable :: zeros(:) integer :: to_select - type(pt2_type) :: pt2_data + type(pt2_type) :: pt2_data, pt2_data_err logical, external :: qp_stop @@ -24,6 +24,7 @@ subroutine run_stochastic_cipsi allocate (zeros(N_states)) call pt2_alloc(pt2_data, N_states) + call pt2_alloc(pt2_data_err, N_states) double precision :: hf_energy_ref logical :: has @@ -79,10 +80,11 @@ subroutine run_stochastic_cipsi to_select = max(N_states_diag, to_select) - pt2_data % pt2 = 0.d0 - pt2_data % variance = 0.d0 - pt2_data % norm2 = 0.d0 - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,to_select) ! Stochastic PT2 and selection + 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(psi_energy_with_nucl_rep,pt2_data,pt2_data_err,relative_error,to_select) ! Stochastic PT2 and selection correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & (psi_energy_with_nucl_rep(1) + pt2_data % rpt2(1) - hf_energy_ref) @@ -90,7 +92,7 @@ subroutine run_stochastic_cipsi call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep, & - pt2_data, N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) @@ -121,17 +123,19 @@ subroutine run_stochastic_cipsi call save_energy(psi_energy_with_nucl_rep, zeros) endif - pt2_data % pt2(:) = 0.d0 - pt2_data % variance(:) = 0.d0 - pt2_data % norm2(:) = 0.d0 - call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2 + 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(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, relative_error, 0) ! Stochastic PT2 call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) call print_summary(psi_energy_with_nucl_rep, & - pt2_data , N_det, N_occ_pattern, N_states, psi_s2) + pt2_data , pt2_data_err, N_det, N_occ_pattern, N_states, psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call print_extrapolated_energy() endif call pt2_dealloc(pt2_data) + call pt2_dealloc(pt2_data_err) end diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f index c34705f5..eb11e28c 100644 --- a/src/fci/pt2.irp.f +++ b/src/fci/pt2.irp.f @@ -32,29 +32,34 @@ subroutine run integer :: i,j,k logical, external :: detEq - type(pt2_type) :: pt2_data + type(pt2_type) :: pt2_data, pt2_data_err integer :: degree integer :: n_det_before, to_select double precision :: threshold_davidson_in - double precision :: E_CI_before(N_states), relative_error + double precision :: relative_error + double precision, allocatable :: E_CI_before(:) + allocate ( E_CI_before(N_states)) call pt2_alloc(pt2_data, N_states) + call pt2_alloc(pt2_data_err, N_states) E_CI_before(:) = psi_energy(:) + nuclear_repulsion relative_error=PT2_relative_error if (do_pt2) then - call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, relative_error, 0) ! Stochastic PT2 else call ZMQ_selection(0, pt2_data) endif call print_summary(psi_energy_with_nucl_rep(1:N_states), & - pt2_data, N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2) - call save_energy(E_CI_before,pt2_data % pt2) + call save_energy(E_CI_before, pt2_data % pt2) call pt2_dealloc(pt2_data) + call pt2_dealloc(pt2_data_err) + deallocate(E_CI_before) end diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 32c07ce5..d2aa4282 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -1,4 +1,4 @@ -subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,n_st,s2_) +subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_occ_pattern_,n_st,s2_) use selection_types implicit none BEGIN_DOC @@ -7,7 +7,7 @@ subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,n_st,s2_) integer, intent(in) :: n_det_, n_occ_pattern_, n_st double precision, intent(in) :: e_(n_st), s2_(n_st) - type(pt2_type) , intent(in) :: pt2_data + type(pt2_type) , intent(in) :: pt2_data, pt2_data_err integer :: i, k integer :: N_states_p character*(9) :: pt2_string @@ -44,16 +44,16 @@ subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,n_st,s2_) write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0 endif write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))' - write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data % pt2_err(k), k=1,N_states_p) - write(*,fmt) '# rPT2'//pt2_string, (pt2_data % pt2(k)*f(k), pt2_data % pt2_err(k)*f(k), k=1,N_states_p) + write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p) + write(*,fmt) '# rPT2'//pt2_string, (pt2_data % pt2(k)*f(k), pt2_data_err % pt2(k)*f(k), k=1,N_states_p) write(*,'(A)') '#' - write(*,fmt) '# E+PT2 ', (e_(k)+pt2_data % pt2(k),pt2_data % pt2_err(k), k=1,N_states_p) - write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % pt2(k)*f(k),pt2_data % pt2_err(k)*f(k), k=1,N_states_p) + write(*,fmt) '# E+PT2 ', (e_(k)+pt2_data % pt2(k),pt2_data_err % pt2(k), k=1,N_states_p) + write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % pt2(k)*f(k),pt2_data_err % pt2(k)*f(k), k=1,N_states_p) if (N_states_p > 1) then write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), & - dsqrt(pt2_data % pt2_err(k)*pt2_data % pt2_err(k)+pt2_data % pt2_err(1)*pt2_data % pt2_err(1)), k=1,N_states_p) + dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p) write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, & - dsqrt(pt2_data % pt2_err(k)*pt2_data % pt2_err(k)+pt2_data % pt2_err(1)*pt2_data % pt2_err(1))*27.211396641308d0, k=1,N_states_p) + dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p) endif write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' write(*,fmt) @@ -70,12 +70,12 @@ subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,n_st,s2_) print*,'* State ',k print *, '< S^2 > = ', s2_(k) print *, 'E = ', e_(k) - print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data % variance_err(k) - print *, 'PT norm = ', dsqrt(pt2_data % norm2(k)), ' +/- ', 0.5d0*dsqrt(pt2_data % norm2(k)) * pt2_data % norm2_err(k) / pt2_data % norm2(k) - print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data % pt2_err(k) - print *, 'rPT2 = ', pt2_data % pt2(k)*f(k), ' +/- ', pt2_data % rpt2_err(k) - print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data % pt2_err(k) - print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % pt2(k)*f(k), ' +/- ', pt2_data % pt2_err(k)*f(k) + print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k) + print *, 'PT norm = ', dsqrt(pt2_data % norm2(k)), ' +/- ', 0.5d0*dsqrt(pt2_data % norm2(k)) * pt2_data_err % norm2(k) / pt2_data % norm2(k) + print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) + print *, 'rPT2 = ', pt2_data % pt2(k)*f(k), ' +/- ', pt2_data_err % rpt2(k) + print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) + print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % pt2(k)*f(k), ' +/- ', pt2_data_err % pt2(k)*f(k) print *, '' enddo From 614cf3b318ad5d4a9db5b07b890396bd3a9b5689 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 Aug 2020 11:37:00 +0200 Subject: [PATCH 11/27] Fixed NaN bug in DIIS --- src/scf_utils/roothaan_hall_scf.irp.f | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index d1236ce7..17458d9a 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -186,7 +186,7 @@ END_DOC implicit none - double precision,intent(in) :: Fock_matrix_DIIS(ao_num,ao_num,*),error_matrix_DIIS(ao_num,ao_num,*) + double precision,intent(in) :: Fock_matrix_DIIS(ao_num,ao_num,dim_DIIS),error_matrix_DIIS(ao_num,ao_num,dim_DIIS) integer,intent(in) :: iteration_SCF, size_Fock_matrix_AO double precision,intent(inout):: Fock_matrix_AO_(size_Fock_matrix_AO,ao_num) integer,intent(inout) :: dim_DIIS @@ -212,11 +212,12 @@ END_DOC ) ! Compute the matrices B and X + B_matrix_DIIS(:,:) = 0.d0 do j=1,dim_DIIS do i=1,dim_DIIS - j_DIIS = mod(iteration_SCF-j,max_dim_DIIS)+1 - i_DIIS = mod(iteration_SCF-i,max_dim_DIIS)+1 + j_DIIS = min(dim_DIIS,mod(iteration_SCF-j,max_dim_DIIS)+1) + i_DIIS = min(dim_DIIS,mod(iteration_SCF-i,max_dim_DIIS)+1) ! Compute product of two errors vectors @@ -229,7 +230,6 @@ END_DOC ! Compute Trace - B_matrix_DIIS(i,j) = 0.d0 do k=1,ao_num B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + scratch(k,k) enddo @@ -238,12 +238,11 @@ END_DOC ! Pad B matrix and build the X matrix + C_vector_DIIS(:) = 0.d0 do i=1,dim_DIIS B_matrix_DIIS(i,dim_DIIS+1) = -1.d0 B_matrix_DIIS(dim_DIIS+1,i) = -1.d0 - C_vector_DIIS(i) = 0.d0 enddo - B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1) = 0.d0 C_vector_DIIS(dim_DIIS+1) = -1.d0 deallocate(scratch) @@ -259,9 +258,10 @@ END_DOC allocate(AF(dim_DIIS+1,dim_DIIS+1)) allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+1)) ) allocate(scratch(lwork,1)) + scratch(:,1) = 0.d0 anorm = dlange('1', dim_DIIS+1, dim_DIIS+1, B_matrix_DIIS, & - size(B_matrix_DIIS,1), scratch) + size(B_matrix_DIIS,1), scratch(1,1)) AF(:,:) = B_matrix_DIIS(:,:) call dgetrf(dim_DIIS+1,dim_DIIS+1,AF,size(AF,1),ipiv,info) From 5cb15587a2f80a4bee44d303507e7b9161e4e5c7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 Aug 2020 15:14:58 +0200 Subject: [PATCH 12/27] Fix bug in zmq_pull --- src/cipsi/run_selection_slave.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index 650654be..c2ba2379 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -208,7 +208,7 @@ subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_i endif allocate(pt2_serialized (pt2_type_size(N_states)) ) - rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized)*ntasks, 0) + rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized), 0) if (rc == -1) then ntasks = 1 task_id(1) = 0 From 47e7e8869afd67057c492f889c6e32dc828492ad Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 Aug 2020 15:28:30 +0200 Subject: [PATCH 13/27] Initialize message in zmq_en_parallel --- src/zmq/utils.irp.f | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/zmq/utils.irp.f b/src/zmq/utils.irp.f index 503dedd2..07e3a88f 100644 --- a/src/zmq/utils.irp.f +++ b/src/zmq/utils.irp.f @@ -585,6 +585,7 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in) stop 'Wrong end of job' endif + message = repeat(' ',512) do i=360,1,-1 rc = f77_zmq_send(zmq_to_qp_run_socket, 'end_job '//trim(zmq_state),8+len(trim(zmq_state)),0) rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 512, 0) @@ -645,6 +646,7 @@ integer function connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) endif endif + message = repeat(' ',512) rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0) message = trim(message(1:rc)) if(message(1:5) == "error") then From 02a2695827b91f4cfe7dcf72adf2cb81042aa68e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 Aug 2020 22:39:40 +0200 Subject: [PATCH 14/27] Removed norm2 --- src/cipsi/cipsi.irp.f | 2 +- src/cipsi/pert_rdm_providers.irp.f | 11 ++++++++-- src/cipsi/pt2_stoch_routines.irp.f | 35 +++++++++++++++--------------- src/cipsi/pt2_type.irp.f | 14 ------------ src/cipsi/selection.irp.f | 3 --- src/cipsi/selection_types.f90 | 3 +-- src/cipsi/stochastic_cipsi.irp.f | 2 +- src/cipsi/zmq_selection.irp.f | 12 ++++++---- src/iterations/print_summary.irp.f | 19 ++++++---------- 9 files changed, 45 insertions(+), 56 deletions(-) diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 0f140240..34b16ff3 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -36,7 +36,7 @@ subroutine run_cipsi zeros = 0.d0 pt2_data % pt2 = -huge(1.e0) pt2_data % rpt2 = -huge(1.e0) - pt2_data % norm2 = 0.d0 + pt2_data % overlap(:,:) = 0.d0 pt2_data % variance = huge(1.e0) if (s2_eig) then diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f index 82ef1e63..caea57b2 100644 --- a/src/cipsi/pert_rdm_providers.irp.f +++ b/src/cipsi/pert_rdm_providers.irp.f @@ -47,7 +47,7 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo type(pt2_type), intent(inout) :: pt2_data type(selection_buffer), intent(inout) :: buf logical :: ok - integer :: s1, s2, p1, p2, ib, j, istate + integer :: s1, s2, p1, p2, ib, j, istate, jstate integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states) double precision, external :: diag_H_mat_elem_fock @@ -152,7 +152,14 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo print*,e_pert,coef,alpha_h_psi pt2_data % pt2(istate) += e_pert pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi - pt2_data % norm2(istate) = coef(istate) * coef(istate) + enddo + + do istate=1,N_states + alpha_h_psi = mat(istate, p1, p2) + e_pert = coef(istate) * alpha_h_psi + do jstate=1,N_states + pt2_data % overlap(jstate,jstate) = coef(istate) * coef(jstate) + enddo if (weight_selection /= 5) then ! Energy selection diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 5b05a5a9..a10cfebd 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -242,8 +242,8 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) 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)/8 & ! pt2_data_task - + N_det_generators*pt2_type_size(N_states)/8 & ! pt2_data_I + + 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 @@ -258,7 +258,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) 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, norm2 + + 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 @@ -300,11 +300,11 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) 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 % norm2(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 % norm2(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) @@ -377,7 +377,8 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ integer, allocatable :: task_id(:) integer, allocatable :: index(:) - double precision :: v, x, x2, x3, avg, avg2, avg3, eqt, E0, v0, n0 + 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(:) @@ -417,8 +418,8 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ 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 % norm2(pt2_stoch_istate) = 0.d0 - pt2_data_err % norm2(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 @@ -437,7 +438,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ n_tasks = 0 E0 = E v0 = 0.d0 - n0 = 0.d0 + n0(:) = 0.d0 more = 1 call wall_time(time0) time1 = time0 @@ -457,11 +458,11 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ t=t+1 E0 = 0.d0 v0 = 0.d0 - n0 = 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) % norm2(pt2_stoch_istate) + n0(:) += pt2_data_I(i) % overlap(:,pt2_stoch_istate) end do else exit @@ -485,7 +486,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ 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) % norm2(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 @@ -494,7 +495,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ endif pt2_data % pt2(pt2_stoch_istate) = avg pt2_data % variance(pt2_stoch_istate) = avg2 - pt2_data % norm2(pt2_stoch_istate) = avg3 + 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 @@ -506,14 +507,14 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ eqt = sqrt(eqt / (dble(c) - 1.5d0)) pt2_data_err % variance(pt2_stoch_istate) = eqt - eqt = dabs((pt2_data_S2(t) % norm2(pt2_stoch_istate) / c) - (pt2_data_S(t) % norm2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = sqrt(eqt / (dble(c) - 1.5d0)) - pt2_data_err % norm2(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 '(G10.3, 2X, F16.10, 2X, G10.3, 2X, F14.10, 2X, F14.10, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3, time-time0, '' + print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, G14.6, 2X, G14.6, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3(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 diff --git a/src/cipsi/pt2_type.irp.f b/src/cipsi/pt2_type.irp.f index ef3fb882..ee90d421 100644 --- a/src/cipsi/pt2_type.irp.f +++ b/src/cipsi/pt2_type.irp.f @@ -7,20 +7,15 @@ subroutine pt2_alloc(pt2_data,N) allocate(pt2_data % pt2(N) & ,pt2_data % variance(N) & - ,pt2_data % norm2(N) & ,pt2_data % rpt2(N) & ,pt2_data % overlap(N,N) & ) pt2_data % pt2(:) = 0.d0 pt2_data % variance(:) = 0.d0 - pt2_data % norm2(:) = 0.d0 pt2_data % rpt2(:) = 0.d0 pt2_data % overlap(:,:) = 0.d0 - do k=1,N - pt2_data % overlap(k,k) = 1.d0 - enddo end subroutine subroutine pt2_dealloc(pt2_data) @@ -29,7 +24,6 @@ subroutine pt2_dealloc(pt2_data) type(pt2_type), intent(inout) :: pt2_data deallocate(pt2_data % pt2 & ,pt2_data % variance & - ,pt2_data % norm2 & ,pt2_data % rpt2 & ,pt2_data % overlap & ) @@ -50,7 +44,6 @@ subroutine pt2_add(p1, w, p2) p1 % pt2(:) = p1 % pt2(:) + p2 % pt2(:) p1 % rpt2(:) = p1 % rpt2(:) + p2 % rpt2(:) p1 % variance(:) = p1 % variance(:) + p2 % variance(:) - p1 % norm2(:) = p1 % norm2(:) + p2 % norm2(:) p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:) else @@ -58,7 +51,6 @@ subroutine pt2_add(p1, w, p2) p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) - p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) endif @@ -81,7 +73,6 @@ subroutine pt2_add2(p1, w, p2) 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 % norm2(:) = p1 % norm2(:) + p2 % norm2(:) * p2 % norm2(:) p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:) * p2 % overlap(:,:) else @@ -89,7 +80,6 @@ subroutine pt2_add2(p1, w, p2) 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 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) * p2 % norm2(:) p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) * p2 % overlap(:,:) endif @@ -113,8 +103,6 @@ subroutine pt2_serialize(pt2_data, n, x) k=k+n x(k+1:k+n) = pt2_data % variance(1:n) k=k+n - x(k+1:k+n) = pt2_data % norm2(1:n) - k=k+n x(k+1:k+n2) = reshape(pt2_data % overlap(1:n,1:n), (/ n2 /)) end @@ -135,8 +123,6 @@ subroutine pt2_deserialize(pt2_data, n, x) k=k+n pt2_data % variance(1:n) = x(k+1:k+n) k=k+n - pt2_data % norm2(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/selection.irp.f b/src/cipsi/selection.irp.f index e68e4b79..e315a852 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -29,7 +29,6 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) type(pt2_type), intent(in) :: pt2_data double precision :: pt2(N_st) double precision :: variance(N_st) - double precision :: norm2(N_st) double precision :: avg, element, dt, x integer :: k @@ -39,7 +38,6 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) pt2(:) = pt2_data % pt2(:) variance(:) = pt2_data % variance(:) - norm2(:) = pt2_data % norm2(:) if (i_iter == 0) then allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax)) @@ -800,7 +798,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d e_pert = coef(istate) * alpha_h_psi pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi - pt2_data % norm2(istate) += coef(istate) * coef(istate) pt2_data % pt2(istate) += e_pert !!!DEBUG diff --git a/src/cipsi/selection_types.f90 b/src/cipsi/selection_types.f90 index 8df37653..58ce0e03 100644 --- a/src/cipsi/selection_types.f90 +++ b/src/cipsi/selection_types.f90 @@ -10,7 +10,6 @@ module selection_types double precision, allocatable :: pt2(:) double precision, allocatable :: rpt2(:) double precision, allocatable :: variance(:) - double precision, allocatable :: norm2(:) double precision, allocatable :: overlap(:,:) endtype @@ -19,7 +18,7 @@ module selection_types integer function pt2_type_size(N) implicit none integer, intent(in) :: N - pt2_type_size = (4*n + n*n) + pt2_type_size = (3*n + n*n) end function end module diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 953314b9..c529795e 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -35,7 +35,7 @@ subroutine run_stochastic_cipsi zeros = 0.d0 pt2_data % pt2 = -huge(1.e0) pt2_data % rpt2 = -huge(1.e0) - pt2_data % norm2 = 0.d0 + pt2_data % overlap= 0.d0 pt2_data % variance = huge(1.e0) if (s2_eig) then diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 789c7a26..e94fd422 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -7,7 +7,7 @@ subroutine ZMQ_selection(N_in, pt2_data) integer(ZMQ_PTR) :: zmq_to_qp_run_socket , zmq_socket_pull integer, intent(in) :: N_in type(selection_buffer) :: b - integer :: i, N + integer :: i, l, N integer, external :: omp_get_thread_num type(pt2_type), intent(inout) :: pt2_data @@ -131,10 +131,13 @@ subroutine ZMQ_selection(N_in, pt2_data) 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) - pt2_data % norm2(k) = pt2_data % norm2(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 % norm2(k)) + pt2_data % pt2(k)/(1.d0 + pt2_data % overlap(k,k)) enddo call update_pt2_and_variance_weights(pt2_data, N_states) @@ -182,6 +185,7 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2_data) 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) @@ -190,7 +194,7 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2_data) more = 1 pt2_data % pt2(:) = 0d0 pt2_data % variance(:) = 0.d0 - pt2_data % norm2(:) = 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) diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index d2aa4282..d04d8a93 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -12,7 +12,6 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_occ_pattern_,n_st,s2_ integer :: N_states_p character*(9) :: pt2_string character*(512) :: fmt - double precision :: f(n_st) if (do_pt2) then pt2_string = ' ' @@ -22,10 +21,6 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_occ_pattern_,n_st,s2_ N_states_p = min(N_det_,n_st) - do i=1,N_states_p - f(i) = 1.d0/(1.d0+pt2_data % norm2(i)) - enddo - print *, '' print '(A,I12)', 'Summary at N_det = ', N_det_ print '(A)', '-----------------------------------' @@ -45,10 +40,10 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_occ_pattern_,n_st,s2_ endif write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))' write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p) - write(*,fmt) '# rPT2'//pt2_string, (pt2_data % pt2(k)*f(k), pt2_data_err % pt2(k)*f(k), k=1,N_states_p) + write(*,fmt) '# rPT2'//pt2_string, (pt2_data % rpt2(k), pt2_data_err % rpt2(k), k=1,N_states_p) write(*,'(A)') '#' write(*,fmt) '# E+PT2 ', (e_(k)+pt2_data % pt2(k),pt2_data_err % pt2(k), k=1,N_states_p) - write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % pt2(k)*f(k),pt2_data_err % pt2(k)*f(k), k=1,N_states_p) + write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % rpt2(k),pt2_data_err % rpt2(k), k=1,N_states_p) if (N_states_p > 1) then write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), & dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p) @@ -71,11 +66,11 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_occ_pattern_,n_st,s2_ print *, '< S^2 > = ', s2_(k) print *, 'E = ', e_(k) print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k) - print *, 'PT norm = ', dsqrt(pt2_data % norm2(k)), ' +/- ', 0.5d0*dsqrt(pt2_data % norm2(k)) * pt2_data_err % norm2(k) / pt2_data % norm2(k) + print *, 'PT norm = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k)) print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) - print *, 'rPT2 = ', pt2_data % pt2(k)*f(k), ' +/- ', pt2_data_err % rpt2(k) + print *, 'rPT2 = ', pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k) print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) - print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % pt2(k)*f(k), ' +/- ', pt2_data_err % pt2(k)*f(k) + print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k) print *, '' enddo @@ -95,8 +90,8 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_occ_pattern_,n_st,s2_ print *, '-----' print*, 'Variational + renormalized perturbative Energy difference (au | eV)' do i=2, N_states_p - print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i)*f(i) - (e_(1) + pt2_data % pt2(1)*f(1))), & - (e_(i)+ pt2_data % pt2(i)*f(i) - (e_(1) + pt2_data % pt2(1)*f(1))) * 27.211396641308d0 + print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), & + (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0 enddo endif From 13abddef2d95041d2f33a34d1318a5a6670d55bb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 Aug 2020 23:04:34 +0200 Subject: [PATCH 15/27] Orthogonalized PT2 --- src/cipsi/energy.irp.f | 8 ++++++++ src/cipsi/pt2_stoch_routines.irp.f | 6 ++++++ src/cipsi/selection.irp.f | 10 +++++----- src/cipsi/zmq_selection.irp.f | 2 ++ src/scf_utils/roothaan_hall_scf.irp.f | 2 +- 5 files changed, 22 insertions(+), 6 deletions(-) diff --git a/src/cipsi/energy.irp.f b/src/cipsi/energy.irp.f index 37b29593..1d8c6bf5 100644 --- a/src/cipsi/energy.irp.f +++ b/src/cipsi/energy.irp.f @@ -37,3 +37,11 @@ BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ] 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/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index a10cfebd..f268430f 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -314,6 +314,12 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) print '(A)', '========== ================= =========== =============== =============== =================' + pt2_overlap(:,pt2_stoch_istate) = pt2_data % overlap(:,pt2_stoch_istate) +print *, 'Overlap' +print *, pt2_overlap(:,pt2_stoch_istate) +print *, '-------' + SOFT_TOUCH pt2_overlap + enddo FREE pt2_stoch_istate diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index e315a852..6813b5a1 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -787,11 +787,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d enddo ! Gram-Schmidt using input overlap matrix -! do istate=1,N_states -! do jstate=1,istate-1 -! coef(istate) = coef(istate) - pt2_overlap(jstate,istate)/(pt2_overlap(jstate,jstate)) * coef(jstate) -! enddo -! enddo + do istate=1,N_states + do jstate=1,istate-1 + coef(istate) = coef(istate) - pt2_overlap(istate,jstate)/(1.d0+pt2_overlap(jstate,jstate)) * coef(jstate) + enddo + enddo do istate=1,N_states alpha_h_psi = mat(istate, p1, p2) diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index e94fd422..9ce345cf 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -140,6 +140,8 @@ subroutine ZMQ_selection(N_in, pt2_data) pt2_data % pt2(k)/(1.d0 + pt2_data % overlap(k,k)) enddo + pt2_overlap(:,:) = pt2_data % overlap(:,:) + SOFT_TOUCH pt2_overlap call update_pt2_and_variance_weights(pt2_data, N_states) end subroutine diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 17458d9a..2b231c2b 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -186,10 +186,10 @@ END_DOC implicit none + integer,intent(inout) :: dim_DIIS double precision,intent(in) :: Fock_matrix_DIIS(ao_num,ao_num,dim_DIIS),error_matrix_DIIS(ao_num,ao_num,dim_DIIS) integer,intent(in) :: iteration_SCF, size_Fock_matrix_AO double precision,intent(inout):: Fock_matrix_AO_(size_Fock_matrix_AO,ao_num) - integer,intent(inout) :: dim_DIIS double precision,allocatable :: B_matrix_DIIS(:,:),X_vector_DIIS(:) double precision,allocatable :: C_vector_DIIS(:) From 7b6c0e13ebbd3b2696cc98f20de5e41b3d9a13bf Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Sep 2020 01:08:03 +0200 Subject: [PATCH 16/27] Changed normalization in PT2 --- src/cipsi/selection.irp.f | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 6813b5a1..9390a9bf 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -789,7 +789,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d ! Gram-Schmidt using input overlap matrix do istate=1,N_states do jstate=1,istate-1 - coef(istate) = coef(istate) - pt2_overlap(istate,jstate)/(1.d0+pt2_overlap(jstate,jstate)) * coef(jstate) + if ( (pt2_overlap(istate,jstate) == 0.d0).or.(pt2_overlap(jstate,jstate) == 0.d0) ) cycle + coef(istate) = coef(istate) - pt2_overlap(istate,jstate)/pt2_overlap(jstate,jstate) * coef(jstate) enddo enddo From 81628a6ae0f67f61d1411f75c39f715c28cecfc2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 2 Sep 2020 10:47:29 +0200 Subject: [PATCH 17/27] Printing overlap --- src/cipsi/pt2_stoch_routines.irp.f | 15 ++++++++++++--- src/cipsi/selection.irp.f | 16 ++++++++-------- 2 files changed, 20 insertions(+), 11 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index f268430f..f0bd8847 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -314,9 +314,18 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) print '(A)', '========== ================= =========== =============== =============== =================' - pt2_overlap(:,pt2_stoch_istate) = pt2_data % overlap(:,pt2_stoch_istate) -print *, 'Overlap' + do k=1,N_states + pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) + enddo +! ! The overlap is not exactly zero because of the guiding function. +! ! Remove the bias +! do k=1,pt2_stoch_istate-1 +! pt2_overlap(k,pt2_stoch_istate) -= pt2_data % overlap(k,pt2_stoch_istate) +! enddo +print *, 'Overlap before orthogonalization' print *, pt2_overlap(:,pt2_stoch_istate) +print *, 'Overlap after orthogonalization' +print *, pt2_overlap(pt2_stoch_istate,:) print *, '-------' SOFT_TOUCH pt2_overlap @@ -520,7 +529,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time - print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, G14.6, 2X, G14.6, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3(pt2_stoch_istate), time-time0, '' + print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, G14.6, 2X, G14.6, 2X, F10.4)', c, avg+E, eqt, avg2, avg3(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 diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 9390a9bf..071d9294 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -780,17 +780,17 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d endif enddo - do istate=1,N_states - do jstate=1,N_states - pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate) - enddo - enddo - ! Gram-Schmidt using input overlap matrix do istate=1,N_states do jstate=1,istate-1 - if ( (pt2_overlap(istate,jstate) == 0.d0).or.(pt2_overlap(jstate,jstate) == 0.d0) ) cycle - coef(istate) = coef(istate) - pt2_overlap(istate,jstate)/pt2_overlap(jstate,jstate) * coef(jstate) + if ( (pt2_overlap(jstate,istate) == 0.d0).or.(pt2_overlap(jstate,jstate) == 0.d0) ) cycle + coef(istate) = coef(istate) - pt2_overlap(jstate,istate)/pt2_overlap(jstate,jstate) * coef(jstate) + enddo + enddo + + do istate=1, N_states + do jstate=1,N_states + pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate) enddo enddo From d65c2cba57a1b794542432cef66d7a0619e17e90 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 2 Sep 2020 15:53:56 +0200 Subject: [PATCH 18/27] Fixed NaN bug in DIIS --- src/scf_utils/roothaan_hall_scf.irp.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index faf23a51..91c85f5e 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -208,8 +208,8 @@ END_DOC do j=1,dim_DIIS do i=1,dim_DIIS - j_DIIS = mod(iteration_SCF-j,max_dim_DIIS)+1 - i_DIIS = mod(iteration_SCF-i,max_dim_DIIS)+1 + j_DIIS = min(mod(iteration_SCF-j,max_dim_DIIS)+1,dim_DIIS) + i_DIIS = min(mod(iteration_SCF-i,max_dim_DIIS)+1,dim_DIIS) ! Compute product of two errors vectors From 202d00ca3a65f50ad50aa23a9a2f665001b79a50 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 2 Sep 2020 17:09:19 +0200 Subject: [PATCH 19/27] minor change --- src/scf_utils/roothaan_hall_scf.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 2b231c2b..dc44e262 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -214,9 +214,9 @@ END_DOC ! Compute the matrices B and X B_matrix_DIIS(:,:) = 0.d0 do j=1,dim_DIIS + j_DIIS = min(dim_DIIS,mod(iteration_SCF-j,max_dim_DIIS)+1) do i=1,dim_DIIS - j_DIIS = min(dim_DIIS,mod(iteration_SCF-j,max_dim_DIIS)+1) i_DIIS = min(dim_DIIS,mod(iteration_SCF-i,max_dim_DIIS)+1) ! Compute product of two errors vectors From 42b74b743f07739c31d84d5f2fe174ccfcb34d2f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Sep 2020 09:27:24 +0200 Subject: [PATCH 20/27] Stupid typo --- src/tools/molden.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f index 239baf18..417b25ad 100644 --- a/src/tools/molden.irp.f +++ b/src/tools/molden.irp.f @@ -17,7 +17,7 @@ program molden write(i_unit_output,'(A)') '[Molden Format]' - write(i_unit_output,'(A)') '[Atoms] Angs + write(i_unit_output,'(A)') '[Atoms] Angs' do i = 1, nucl_num write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') & trim(element_name(int(nucl_charge(i)))), & From c07232aae3ccbb640d29e444b70c558c287261f6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Sep 2020 11:09:25 +0200 Subject: [PATCH 21/27] Introduce in PT2 --- src/cipsi/selection.irp.f | 35 ++++++++++++++++++++++------------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 071d9294..22d737a1 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -674,7 +674,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d logical :: ok integer :: s1, s2, p1, p2, ib, j, istate, jstate integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert, delta_E, val, Hii, w, tmp, alpha_h_psi, coef(N_states) + double precision :: e_pert(N_states), coef(N_states), X(N_states) + double precision :: delta_E, val, Hii, w, tmp, alpha_h_psi double precision, external :: diag_H_mat_elem_fock double precision :: E_shift @@ -772,21 +773,26 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if (delta_E < 0.d0) then tmp = -tmp endif - e_pert = 0.5d0 * (tmp - delta_E) + e_pert(istate) = 0.5d0 * (tmp - delta_E) if (dabs(alpha_h_psi) > 1.d-4) then - coef(istate) = e_pert / alpha_h_psi + coef(istate) = e_pert(istate) / alpha_h_psi else coef(istate) = alpha_h_psi / delta_E endif + if (e_pert(istate) < 0.d0) then + X(istate) = -dsqrt(-e_pert(istate)) + else + X(istate) = dsqrt(e_pert(istate)) + endif enddo - ! Gram-Schmidt using input overlap matrix - do istate=1,N_states - do jstate=1,istate-1 - if ( (pt2_overlap(jstate,istate) == 0.d0).or.(pt2_overlap(jstate,jstate) == 0.d0) ) cycle - coef(istate) = coef(istate) - pt2_overlap(jstate,istate)/pt2_overlap(jstate,jstate) * coef(jstate) - enddo - enddo +! ! Gram-Schmidt using input overlap matrix +! do istate=1,N_states +! do jstate=1,istate-1 +! if ( (pt2_overlap(jstate,istate) == 0.d0).or.(pt2_overlap(jstate,jstate) == 0.d0) ) cycle +! coef(istate) = coef(istate) - pt2_overlap(jstate,istate)/pt2_overlap(jstate,jstate) * coef(jstate) +! enddo +! enddo do istate=1, N_states do jstate=1,N_states @@ -796,10 +802,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d do istate=1,N_states alpha_h_psi = mat(istate, p1, p2) - e_pert = coef(istate) * alpha_h_psi pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi - pt2_data % pt2(istate) += e_pert + pt2_data % pt2(istate) += e_pert(istate) !!!DEBUG ! delta_E = E0(istate) - Hii + E_shift @@ -831,7 +836,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d case default ! Energy selection - w = w + e_pert * selection_weight(istate) + w = w + e_pert(istate) * selection_weight(istate) + do jstate=1,N_states + if (istate == jstate) cycle + w = w - dabs(X(istate))*X(jstate) * sqrt(selection_weight(istate)*selection_weight(jstate)) + enddo end select end do From fd7825ea742163b80950c8495996e4b871c95c83 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Sep 2020 16:16:16 +0200 Subject: [PATCH 22/27] Changed default weight_selection to 1 --- src/cipsi/pt2_stoch_routines.irp.f | 15 +++++++++++---- src/cipsi/selection.irp.f | 23 ++++++++++++++++++----- src/determinants/EZFIO.cfg | 2 +- 3 files changed, 30 insertions(+), 10 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index f0bd8847..faced03a 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -287,9 +287,9 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) call omp_set_nested(.false.) - print '(A)', '========== ================= =========== =============== =============== =================' - print '(A)', ' Samples Energy Stat. Err Variance Norm^2 Seconds ' - print '(A)', '========== ================= =========== =============== =============== =================' + print '(A)', '========== ======================= ========= ========== ========== ========== ========== ' + print '(A)', ' Samples Energy Variance Norm^2 Seconds' + print '(A)', '========== ======================= ==================== ===================== ==========' PROVIDE global_selection_buffer @@ -529,7 +529,14 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time - print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, G14.6, 2X, G14.6, 2X, F10.4)', c, avg+E, eqt, avg2, avg3(pt2_stoch_istate), time-time0 + print '(I10, X, F11.7, X, G10.3, X, G10.3, X, G10.3, X, G10.3, 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 diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 22d737a1..e599737c 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -683,7 +683,12 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d double precision, allocatable :: values(:) integer, allocatable :: keys(:,:) integer :: nkeys - + double precision :: s_weight(N_states,N_states) + do jstate=1,N_states + do istate=1,N_states + s_weight(istate,jstate) = dsqrt(selection_weight(istate)*selection_weight(jstate)) + enddo + enddo if(sp == 3) then s1 = 1 @@ -829,17 +834,25 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d case(5) ! Variance selection - w = w - alpha_h_psi * alpha_h_psi * selection_weight(istate) + w = w - alpha_h_psi * alpha_h_psi * s_weight(istate,istate) + do jstate=1,N_states + if (istate == jstate) cycle + w = w + alpha_h_psi*mat(jstate,p1,p2) * s_weight(istate,jstate) + enddo case(6) - w = w - coef(istate) * coef(istate) * selection_weight(istate) + w = w - coef(istate) * coef(istate) * s_weight(istate,istate) + do jstate=1,N_states + if (istate == jstate) cycle + w = w + coef(istate)*coef(jstate) * s_weight(istate,jstate) + enddo case default ! Energy selection - w = w + e_pert(istate) * selection_weight(istate) + w = w + e_pert(istate) * s_weight(istate,istate) do jstate=1,N_states if (istate == jstate) cycle - w = w - dabs(X(istate))*X(jstate) * sqrt(selection_weight(istate)*selection_weight(jstate)) + w = w - dabs(X(istate))*X(jstate) * s_weight(istate,jstate) enddo end select diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index ef00080b..662c6fbb 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -44,7 +44,7 @@ default: 2 type: integer doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching, 6: CI coefficients 7: input state-average multiplied by variance and rPT2 matching 8: input state-average multiplied by rPT2 matching 9: input state-average multiplied by variance matching interface: ezfio,provider,ocaml -default: 2 +default: 1 [threshold_generators] type: Threshold From 8c75ad6cfab211da3519d9d7aad98676ac0a2b7f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Sep 2020 17:08:45 +0200 Subject: [PATCH 23/27] Prints --- src/cipsi/pt2_stoch_routines.irp.f | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index faced03a..72f00d70 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -287,9 +287,9 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) call omp_set_nested(.false.) - print '(A)', '========== ======================= ========= ========== ========== ========== ========== ' + print '(A)', '========== ====================== ===================== ===================== ===========' print '(A)', ' Samples Energy Variance Norm^2 Seconds' - print '(A)', '========== ======================= ==================== ===================== ==========' + print '(A)', '========== ====================== ===================== ===================== ===========' PROVIDE global_selection_buffer @@ -312,7 +312,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) !$OMP END PARALLEL call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') - print '(A)', '========== ================= =========== =============== =============== =================' + print '(A)', '========== ====================== ===================== ===================== ===========' do k=1,N_states pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) From 07785a6db19a93330d212c5e692797d7a9de2347 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Sep 2020 17:13:26 +0200 Subject: [PATCH 24/27] Print overlap in deterministic PT2 --- src/cipsi/pt2_stoch_routines.irp.f | 11 ++--------- src/cipsi/zmq_selection.irp.f | 6 ++++++ 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 72f00d70..4d6ea1e8 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -317,14 +317,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) do k=1,N_states pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) enddo -! ! The overlap is not exactly zero because of the guiding function. -! ! Remove the bias -! do k=1,pt2_stoch_istate-1 -! pt2_overlap(k,pt2_stoch_istate) -= pt2_data % overlap(k,pt2_stoch_istate) -! enddo -print *, 'Overlap before orthogonalization' -print *, pt2_overlap(:,pt2_stoch_istate) -print *, 'Overlap after orthogonalization' +print *, 'Overlap of perturbed states:' print *, pt2_overlap(pt2_stoch_istate,:) print *, '-------' SOFT_TOUCH pt2_overlap @@ -529,7 +522,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time - print '(I10, X, F11.7, X, G10.3, X, G10.3, X, G10.3, X, G10.3, X, G10.3, X, F10.4)', c, & + print '(I10, X, F11.7, X, G10.4, X, G10.3, X, G10.3, X, G10.3, 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), & diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 9ce345cf..448b409e 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -141,6 +141,12 @@ subroutine ZMQ_selection(N_in, pt2_data) 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) From 5fcdbe12dfa1ba81202ab001a2b222c2b22fe055 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Sep 2020 17:13:26 +0200 Subject: [PATCH 25/27] Print overlap in deterministic PT2 --- src/cipsi/pt2_stoch_routines.irp.f | 11 ++--------- src/cipsi/zmq_selection.irp.f | 6 ++++++ 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 72f00d70..7622e93f 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -317,14 +317,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) do k=1,N_states pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) enddo -! ! The overlap is not exactly zero because of the guiding function. -! ! Remove the bias -! do k=1,pt2_stoch_istate-1 -! pt2_overlap(k,pt2_stoch_istate) -= pt2_data % overlap(k,pt2_stoch_istate) -! enddo -print *, 'Overlap before orthogonalization' -print *, pt2_overlap(:,pt2_stoch_istate) -print *, 'Overlap after orthogonalization' +print *, 'Overlap of perturbed states:' print *, pt2_overlap(pt2_stoch_istate,:) print *, '-------' SOFT_TOUCH pt2_overlap @@ -529,7 +522,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time - print '(I10, X, F11.7, X, G10.3, X, G10.3, X, G10.3, X, G10.3, X, G10.3, X, F10.4)', c, & + print '(I10, X, F11.7, X, F10.6, X, G10.3, X, F10.6, X, G10.3, 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), & diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 9ce345cf..448b409e 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -141,6 +141,12 @@ subroutine ZMQ_selection(N_in, pt2_data) 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) From 0a5f3ac330cb9d75c641fc32ecf3b2c319558040 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Sep 2020 18:12:58 +0200 Subject: [PATCH 26/27] PT2 overlap --- src/cipsi/pt2_stoch_routines.irp.f | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 7622e93f..ab7e765a 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -317,14 +317,25 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) do k=1,N_states pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) enddo -print *, 'Overlap of perturbed states:' -print *, pt2_overlap(pt2_stoch_istate,:) -print *, '-------' 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 From f0454b5b344f011520e6f2c70ded66343fa9ded0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Sep 2020 23:14:08 +0200 Subject: [PATCH 27/27] Print format of PT2 --- src/cipsi/pt2_stoch_routines.irp.f | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 1d13c737..31f27e1d 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -287,9 +287,9 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) call omp_set_nested(.false.) - print '(A)', '========== ====================== ===================== ===================== ===========' + print '(A)', '========== ======================= ===================== ===================== ===========' print '(A)', ' Samples Energy Variance Norm^2 Seconds' - print '(A)', '========== ====================== ===================== ===================== ===========' + print '(A)', '========== ======================= ===================== ===================== ===========' PROVIDE global_selection_buffer @@ -312,7 +312,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) !$OMP END PARALLEL call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') - print '(A)', '========== ====================== ===================== ===================== ===========' + print '(A)', '========== ======================= ===================== ===================== ===========' do k=1,N_states pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) @@ -533,7 +533,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time - print '(I10, X, F11.7, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.4)', c, & + 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), &