PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00
Code Releases Activity

Factorized pt2 data

Browse Source
This commit is contained in:
Anthony Scemama 2020-08-28 15:39:01 +02:00
parent e0f17d516b
commit 7bde6f7451
8 changed files with 148 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -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

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

@ -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)

74
src/cipsi/pt2_type.irp.f Normal file
View File

@ -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

9
src/cipsi/selection_types.f90
View File

@ -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

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

@ -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

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

@ -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

18
src/fci/pt2.irp.f
View File

@ -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

44
src/iterations/print_summary.irp.f
View File

@ -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