From cf2d78fced92a660ff145b33d0a79b5f70318976 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Aug 2020 16:05:53 +0200 Subject: [PATCH] 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)