LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 03:23:29 +01:00
Code Releases Activity

Cleaned rPT2

Browse Source
This commit is contained in:
Anthony Scemama 2020-08-28 16:05:53 +02:00
parent 7bde6f7451
commit cf2d78fced
5 changed files with 45 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
src/cipsi/cipsi.irp.f
View File

@ -7,7 +7,7 @@ subroutine run_cipsi
END_DOC END_DOC
integer :: i,j,k integer :: i,j,k
type(pt2_type) :: pt2_data type(pt2_type) :: pt2_data
double precision, allocatable :: rpt2(:), zeros(:) double precision, allocatable :: zeros(:)
integer :: to_select integer :: to_select
logical, external :: qp_stop logical, external :: qp_stop
@ -23,8 +23,8 @@ subroutine run_cipsi
rss = memory_of_double(N_states)*4.d0 rss = memory_of_double(N_states)*4.d0
call check_mem(rss,irp_here) call check_mem(rss,irp_here)
allocate (zeros(N_states), rpt2(N_states)) allocate (zeros(N_states))
allocate (pt2_data % pt2(N_states), pt2_data % norm2(N_states), pt2_data % variance(N_states)) call pt2_alloc(pt2_data, N_states)
double precision :: hf_energy_ref double precision :: hf_energy_ref
logical :: has logical :: has
@ -33,8 +33,8 @@ subroutine run_cipsi
relative_error=PT2_relative_error relative_error=PT2_relative_error
zeros = 0.d0 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 % norm2 = 0.d0
pt2_data % variance = huge(1.e0) pt2_data % variance = huge(1.e0)
@ -92,21 +92,17 @@ subroutine run_cipsi
call ZMQ_selection(to_select, pt2_data) call ZMQ_selection(to_select, pt2_data)
endif 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) / & 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) correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
call write_double(6,correlation_energy_ratio, 'Correlation ratio') call write_double(6,correlation_energy_ratio, 'Correlation ratio')
call print_summary(psi_energy_with_nucl_rep, & call print_summary(psi_energy_with_nucl_rep, &
pt2_data, 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) 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() call print_extrapolated_energy()
N_iter += 1 N_iter += 1
@ -147,15 +143,10 @@ subroutine run_cipsi
print *, 'N_states = ', N_states print *, 'N_states = ', N_states
print*, 'correlation_ratio = ', correlation_energy_ratio print*, 'correlation_ratio = ', correlation_energy_ratio
call save_energy(psi_energy_with_nucl_rep, pt2_data % 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 print_summary(psi_energy_with_nucl_rep(1:N_states), & call print_summary(psi_energy_with_nucl_rep(1:N_states), &
pt2_data, 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 save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call print_extrapolated_energy() call print_extrapolated_energy()
endif endif

38
src/cipsi/pt2_stoch_routines.irp.f
View File

@ -292,16 +292,19 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in)
print '(A)', '========== ================= =========== =============== =============== =================' print '(A)', '========== ================= =========== =============== =============== ================='
PROVIDE global_selection_buffer PROVIDE global_selection_buffer
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc_target+1) & !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc_target+1) &
!$OMP PRIVATE(i) !$OMP PRIVATE(i)
i = omp_get_thread_num() i = omp_get_thread_num()
if (i==0) then 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) call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, pt2_data, b, N)
pt2_data % pt2(pt2_stoch_istate) = w(pt2_stoch_istate,1) pt2_data % rpt2(pt2_stoch_istate) = &
pt2_data % pt2_err(pt2_stoch_istate) = w(pt2_stoch_istate,2) pt2_data % pt2(pt2_stoch_istate)/(1.d0 + pt2_data % norm2(pt2_stoch_istate))
pt2_data % variance(pt2_stoch_istate) = w(pt2_stoch_istate,3)
pt2_data % norm2(pt2_stoch_istate) = w(pt2_stoch_istate,4) !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 else
call pt2_slave_inproc(i) call pt2_slave_inproc(i)
@ -343,7 +346,7 @@ subroutine pt2_slave_inproc(i)
end 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 f77_zmq
use selection_types use selection_types
use bitmasks 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 integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
double precision, intent(in) :: relative_error, E double precision, intent(in) :: relative_error, E
double precision, intent(out) :: pt2(N_states), error(N_states) type(pt2_type), intent(inout) :: pt2_data
double precision, intent(out) :: variance(N_states), norm2(N_states)
type(selection_buffer), intent(inout) :: b type(selection_buffer), intent(inout) :: b
integer, intent(in) :: N_ 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) call create_selection_buffer(N_, N_*2, b2)
pt2(:) = -huge(1.) pt2_data % pt2(pt2_stoch_istate) = -huge(1.)
error(:) = huge(1.) pt2_data % pt2_err(pt2_stoch_istate) = huge(1.)
variance(:) = huge(1.) pt2_data % variance(pt2_stoch_istate) = huge(1.)
norm2(:) = 0.d0 pt2_data % norm2(pt2_stoch_istate) = 0.d0
S(:) = 0d0 S(:) = 0d0
S2(:) = 0d0 S2(:) = 0d0
T2(:) = 0d0 T2(:) = 0d0
@ -486,21 +488,21 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error, varianc
if (qp_stop()) then if (qp_stop()) then
stop_now = .True. stop_now = .True.
endif endif
pt2(pt2_stoch_istate) = avg pt2_data % pt2(pt2_stoch_istate) = avg
variance(pt2_stoch_istate) = avg2 pt2_data % variance(pt2_stoch_istate) = avg2
norm2(pt2_stoch_istate) = avg3 pt2_data % norm2(pt2_stoch_istate) = avg3
call wall_time(time) call wall_time(time)
! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969) ! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969)
if(c > 2) then if(c > 2) then
eqt = dabs((S2(t) / c) - (S(t)/c)**2) ! dabs for numerical stability eqt = dabs((S2(t) / c) - (S(t)/c)**2) ! dabs for numerical stability
eqt = sqrt(eqt / (dble(c) - 1.5d0)) 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 if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then
time1 = time 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, F14.10, 2X, F14.10, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3, time-time0, ''
if (stop_now .or. ( & if (stop_now .or. ( &
(do_exit .and. (dabs(error(pt2_stoch_istate)) / & (do_exit .and. (dabs(pt2_data % pt2_err(pt2_stoch_istate)) / &
(1.d-20 + dabs(pt2(pt2_stoch_istate)) ) <= relative_error))) ) then (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then
if (zmq_abort(zmq_to_qp_run_socket) == -1) then if (zmq_abort(zmq_to_qp_run_socket) == -1) then
call sleep(10) call sleep(10)
if (zmq_abort(zmq_to_qp_run_socket) == -1) then if (zmq_abort(zmq_to_qp_run_socket) == -1) then

7
src/cipsi/selection.irp.f
View File

@ -31,7 +31,7 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st)
double precision :: variance(N_st) double precision :: variance(N_st)
double precision :: norm2(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 :: k
integer, save :: i_iter=0 integer, save :: i_iter=0
integer, parameter :: i_itermax = 1 integer, parameter :: i_itermax = 1
@ -54,11 +54,6 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st)
dt = 2.0d0 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 avg = sum(pt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero
do k=1,N_st do k=1,N_st
element = exp(dt*(pt2(k)/avg -1.d0)) element = exp(dt*(pt2(k)/avg -1.d0))

24
src/cipsi/stochastic_cipsi.irp.f
View File

@ -5,7 +5,7 @@ subroutine run_stochastic_cipsi
! Selected Full Configuration Interaction with Stochastic selection and PT2. ! Selected Full Configuration Interaction with Stochastic selection and PT2.
END_DOC END_DOC
integer :: i,j,k integer :: i,j,k
double precision, allocatable :: rpt2(:), zeros(:) double precision, allocatable :: zeros(:)
integer :: to_select integer :: to_select
type(pt2_type) :: pt2_data type(pt2_type) :: pt2_data
logical, external :: qp_stop logical, external :: qp_stop
@ -22,7 +22,7 @@ subroutine run_stochastic_cipsi
rss = memory_of_double(N_states)*4.d0 rss = memory_of_double(N_states)*4.d0
call check_mem(rss,irp_here) call check_mem(rss,irp_here)
allocate (zeros(N_states), rpt2(N_states)) allocate (zeros(N_states))
call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data, N_states)
double precision :: hf_energy_ref double precision :: hf_energy_ref
@ -32,8 +32,8 @@ subroutine run_stochastic_cipsi
relative_error=PT2_relative_error relative_error=PT2_relative_error
zeros = 0.d0 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 % norm2 = 0.d0
pt2_data % variance = huge(1.e0) pt2_data % variance = huge(1.e0)
@ -84,21 +84,17 @@ subroutine run_stochastic_cipsi
pt2_data % norm2 = 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 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) / & 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) correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
call write_double(6,correlation_energy_ratio, 'Correlation ratio') call write_double(6,correlation_energy_ratio, 'Correlation ratio')
call print_summary(psi_energy_with_nucl_rep, & call print_summary(psi_energy_with_nucl_rep, &
pt2_data, 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) 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() call print_extrapolated_energy()
N_iter += 1 N_iter += 1
@ -130,14 +126,10 @@ subroutine run_stochastic_cipsi
pt2_data % norm2(:) = 0.d0 pt2_data % norm2(:) = 0.d0
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, relative_error, 0) ! Stochastic PT2
do k=1,N_states call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
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, & call print_summary(psi_energy_with_nucl_rep, &
pt2_data , 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 save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call print_extrapolated_energy() call print_extrapolated_energy()
endif endif
call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data)

5
src/cipsi/zmq_selection.irp.f
View File

@ -128,13 +128,14 @@ subroutine ZMQ_selection(N_in, pt2_data)
endif endif
call delete_selection_buffer(b) 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 % pt2(k) = pt2_data % pt2(k) * f(k)
pt2_data % variance(k) = pt2_data % variance(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 % norm2(k) = pt2_data % norm2(k) * f(k)
pt2_data % rpt2(k) = &
pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k))
enddo enddo
endif
call update_pt2_and_variance_weights(pt2_data, N_states) call update_pt2_and_variance_weights(pt2_data, N_states)