LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-01-03 10:05:57 +01:00
Code Issues Releases Wiki Activity

Multi-state PT2 OK

Browse Source
This commit is contained in:
Anthony Scemama 2017-11-20 15:19:00 +01:00
parent 31557a1180
commit 14eb5b4a42
3 changed files with 80 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
plugins/Full_CI_ZMQ/fci_zmq.irp.f
View File

@ -72,7 +72,7 @@ program fci_zmq
if (do_pt2) then if (do_pt2) then
pt2_string = ' ' pt2_string = ' '
pt2 = 0.d0 pt2 = 0.d0
if (N_states == 1) then ! if (N_states == 1) then
threshold_selectors = 1.d0 threshold_selectors = 1.d0
threshold_generators = 1d0 threshold_generators = 1d0
SOFT_TOUCH threshold_selectors threshold_generators SOFT_TOUCH threshold_selectors threshold_generators
@ -80,12 +80,12 @@ program fci_zmq
threshold_selectors = threshold_selectors_save threshold_selectors = threshold_selectors_save
threshold_generators = threshold_generators_save threshold_generators = threshold_generators_save
SOFT_TOUCH threshold_selectors threshold_generators SOFT_TOUCH threshold_selectors threshold_generators
else ! else
threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) ! threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
threshold_generators = max(threshold_generators,threshold_generators_pt2) ! threshold_generators = max(threshold_generators,threshold_generators_pt2)
SOFT_TOUCH threshold_selectors threshold_generators ! SOFT_TOUCH threshold_selectors threshold_generators
call ZMQ_selection(0, pt2) ! Deterministic PT2 ! call ZMQ_selection(0, pt2) ! Deterministic PT2
endif ! endif
else else
pt2_string = '(approx)' pt2_string = '(approx)'
endif endif
@ -104,24 +104,26 @@ program fci_zmq
print*,'State ',k print*,'State ',k
print *, 'PT2 = ', pt2(k) print *, 'PT2 = ', pt2(k)
print *, 'E = ', CI_energy(k) print *, 'E = ', CI_energy(k)
if (N_states==1) then ! if (N_states==1) then
print *, 'E+PT2'//pt2_string//' = ', CI_energy(k)+pt2(k), ' +/- ', error print *, 'E+PT2'//pt2_string//' = ', CI_energy(k)+pt2(k), ' +/- ', error
else ! else
print *, 'E+PT2'//pt2_string//' = ', CI_energy(k)+pt2(k) ! print *, 'E+PT2'//pt2_string//' = ', CI_energy(k)+pt2(k)
endif ! endif
enddo enddo
print *, '-----' print *, '-----'
if(N_states.gt.1)then if(N_states.gt.1)then
print*,'Variational Energy difference' print*,'Variational Energy difference (au | eV)'
do i = 2, N_states do i = 2, N_states
print*,'Delta E = ',CI_energy(i) - CI_energy(1) print*,'Delta E = ', (CI_energy(i) - CI_energy(1)), &
(CI_energy(i) - CI_energy(1)) * 27.2107362681d0
enddo enddo
endif endif
if(N_states.gt.1)then if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference' print*,'Variational + perturbative Energy difference (au | eV)'
do i = 2, N_states do i = 2, N_states
print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1)) print*,'Delta E = ', (CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1))), &
(CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1))) * 27.2107362681d0
enddo enddo
endif endif
@ -157,7 +159,7 @@ program fci_zmq
if (do_pt2) then if (do_pt2) then
pt2 = 0.d0 pt2 = 0.d0
if (N_states == 1) then ! if (N_states == 1) then
threshold_selectors = 1.d0 threshold_selectors = 1.d0
threshold_generators = 1d0 threshold_generators = 1d0
SOFT_TOUCH threshold_selectors threshold_generators SOFT_TOUCH threshold_selectors threshold_generators
@ -165,12 +167,12 @@ program fci_zmq
threshold_selectors = threshold_selectors_save threshold_selectors = threshold_selectors_save
threshold_generators = threshold_generators_save threshold_generators = threshold_generators_save
SOFT_TOUCH threshold_selectors threshold_generators SOFT_TOUCH threshold_selectors threshold_generators
else ! else
threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) ! threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
threshold_generators = max(threshold_generators,threshold_generators_pt2) ! threshold_generators = max(threshold_generators,threshold_generators_pt2)
SOFT_TOUCH threshold_selectors threshold_generators ! SOFT_TOUCH threshold_selectors threshold_generators
call ZMQ_selection(0, pt2) ! Deterministic PT2 ! call ZMQ_selection(0, pt2) ! Deterministic PT2
endif ! endif
call ezfio_set_full_ci_zmq_energy(CI_energy(1)) call ezfio_set_full_ci_zmq_energy(CI_energy(1))
call ezfio_set_full_ci_zmq_energy_pt2(CI_energy(1)+pt2(1)) call ezfio_set_full_ci_zmq_energy_pt2(CI_energy(1)+pt2(1))
call dump_fci_iterations_value(N_det,CI_energy(1),pt2(1)) ! This call automatically appends data call dump_fci_iterations_value(N_det,CI_energy(1),pt2(1)) ! This call automatically appends data
@ -183,11 +185,11 @@ program fci_zmq
print*,'State ',k print*,'State ',k
print *, 'PT2 = ', pt2(k) print *, 'PT2 = ', pt2(k)
print *, 'E = ', CI_energy(k) print *, 'E = ', CI_energy(k)
if (N_states==1) then ! if (N_states==1) then
print *, 'E+PT2'//pt2_string//' = ', CI_energy(k)+pt2(k), ' +/- ', error print *, 'E+PT2'//pt2_string//' = ', CI_energy(k)+pt2(k), ' +/- ', error
else ! else
print *, 'E+PT2'//pt2_string//' = ', CI_energy(k)+pt2(k) ! print *, 'E+PT2'//pt2_string//' = ', CI_energy(k)+pt2(k)
endif ! endif
enddo enddo
print *, '-----' print *, '-----'

40
plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
View File

@ -13,7 +13,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, eqt)
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_to_qp_run_socket2 integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_to_qp_run_socket2
type(selection_buffer) :: b type(selection_buffer) :: b
integer, external :: omp_get_thread_num integer, external :: omp_get_thread_num
double precision, intent(in) :: relative_error, absolute_error, E double precision, intent(in) :: relative_error, absolute_error, E(N_states)
double precision, intent(out) :: pt2(N_states),eqt double precision, intent(out) :: pt2(N_states),eqt
@ -26,13 +26,20 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, eqt)
double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth) double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
double precision, external :: omp_get_wtime double precision, external :: omp_get_wtime
double precision :: time double precision :: time
double precision :: w(N_states)
integer :: istate
if (N_det < 10) then if (N_det < 10) then
call ZMQ_selection(0, pt2) call ZMQ_selection(0, pt2)
return return
else else
allocate(pt2_detail(N_states, N_det_generators+1), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc)) do istate=1,N_states
w(:) = 0.d0
w(istate) = 1.d0
call update_psi_average_norm_contrib(w)
allocate(pt2_detail(N_states,N_det_generators+1), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
sumabove = 0d0 sumabove = 0d0
sum2above = 0d0 sum2above = 0d0
Nabove = 0d0 Nabove = 0d0
@ -93,7 +100,8 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, eqt)
!$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(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, absolute_error, pt2,eqt) call pt2_collector(E(istate), b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, absolute_error,w,eqt,istate)
pt2(istate) = w(istate)
else else
call pt2_slave_inproc(i) call pt2_slave_inproc(i)
endif endif
@ -104,15 +112,16 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, eqt)
print *, '========== ================= ================= =================' print *, '========== ================= ================= ================='
deallocate(pt2_detail, comb, computed, tbc) deallocate(pt2_detail, comb, computed, tbc)
enddo
endif endif
end subroutine end subroutine
subroutine do_carlo(tbc, Ncomb, comb, pt2_detail, computed, sumabove, sum2above, Nabove) subroutine do_carlo(tbc, Ncomb, comb, pt2_detail, computed, sumabove, sum2above, Nabove,istate)
integer, intent(in) :: tbc(0:size_tbc), Ncomb integer, intent(in) :: tbc(0:size_tbc), Ncomb, istate
logical, intent(in) :: computed(N_det_generators) logical, intent(in) :: computed(N_det_generators)
double precision, intent(in) :: comb(Ncomb), pt2_detail(N_states, N_det_generators) double precision, intent(in) :: comb(Ncomb), pt2_detail(N_states,N_det_generators)
double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth) double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
integer :: i, dets(comb_teeth) integer :: i, dets(comb_teeth)
double precision :: myVal, myVal2 double precision :: myVal, myVal2
@ -128,7 +137,7 @@ subroutine do_carlo(tbc, Ncomb, comb, pt2_detail, computed, sumabove, sum2above,
myVal = 0d0 myVal = 0d0
myVal2 = 0d0 myVal2 = 0d0
do j=comb_teeth,1,-1 do j=comb_teeth,1,-1
myVal += pt2_detail(1, dets(j)) * pt2_weight_inv(dets(j)) * comb_step myVal += pt2_detail(istate,dets(j)) * pt2_weight_inv(dets(j)) * comb_step
sumabove(j) += myVal sumabove(j) += myVal
sum2above(j) += myVal*myVal sum2above(j) += myVal*myVal
Nabove(j) += 1 Nabove(j) += 1
@ -144,7 +153,7 @@ subroutine pt2_slave_inproc(i)
call run_pt2_slave(1,i,pt2_e0_denominator) call run_pt2_slave(1,i,pt2_e0_denominator)
end end
subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, absolute_error, pt2,eqt) subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, absolute_error, pt2,eqt,istate)
use f77_zmq use f77_zmq
use selection_types use selection_types
use bitmasks use bitmasks
@ -152,6 +161,7 @@ subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove,
integer, intent(in) :: Ncomb integer, intent(in) :: Ncomb
integer, intent(in) :: istate
double precision, intent(inout) :: pt2_detail(N_states, N_det_generators) double precision, intent(inout) :: pt2_detail(N_states, N_det_generators)
double precision, intent(in) :: comb(Ncomb), relative_error, absolute_error, E double precision, intent(in) :: comb(Ncomb), relative_error, absolute_error, E
logical, intent(inout) :: computed(N_det_generators) logical, intent(inout) :: computed(N_det_generators)
@ -254,16 +264,16 @@ subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove,
call zmq_abort(zmq_to_qp_run_socket) call zmq_abort(zmq_to_qp_run_socket)
exit pullLoop exit pullLoop
endif endif
call do_carlo(tbc, Ncomb+1-firstTBDcomb, comb(firstTBDcomb), pt2_detail, actually_computed, sumabove, sum2above, Nabove) call do_carlo(tbc, Ncomb+1-firstTBDcomb, comb(firstTBDcomb), pt2_detail, actually_computed, sumabove, sum2above, Nabove, istate)
firstTBDcomb = int(Nabove(1)) - orgTBDcomb + 1 firstTBDcomb = int(Nabove(1)) - orgTBDcomb + 1
if(Nabove(1) < 5d0) cycle if(Nabove(1) < 5d0) cycle
call get_first_tooth(actually_computed, tooth) call get_first_tooth(actually_computed, tooth)
E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1)) E0 = sum(pt2_detail(istate,:first_det_of_teeth(tooth)-1))
if (tooth <= comb_teeth) then if (tooth <= comb_teeth) then
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1)) prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1))
prop = prop * pt2_weight_inv(first_det_of_teeth(tooth)) prop = prop * pt2_weight_inv(first_det_of_teeth(tooth))
E0 += pt2_detail(1,first_det_of_teeth(tooth)) * prop E0 += pt2_detail(istate,first_det_of_teeth(tooth)) * prop
avg = E0 + (sumabove(tooth) / Nabove(tooth)) avg = E0 + (sumabove(tooth) / Nabove(tooth))
eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2)) eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
else else
@ -272,7 +282,7 @@ subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove,
call wall_time(time) call wall_time(time)
if ( (dabs(eqt/avg) < relative_error) .or. (dabs(eqt) < absolute_error) ) then if ( (dabs(eqt/avg) < relative_error) .or. (dabs(eqt) < absolute_error) ) then
! Termination ! Termination
pt2(1) = avg pt2(istate) = avg
print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, '' print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, ''
call zmq_abort(zmq_to_qp_run_socket) call zmq_abort(zmq_to_qp_run_socket)
else else
@ -284,11 +294,11 @@ subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove,
end if end if
end do pullLoop end do pullLoop
E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1)) E0 = sum(pt2_detail(istate,:first_det_of_teeth(tooth)-1))
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1)) prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1))
prop = prop * pt2_weight_inv(first_det_of_teeth(tooth)) prop = prop * pt2_weight_inv(first_det_of_teeth(tooth))
E0 += pt2_detail(1,first_det_of_teeth(tooth)) * prop E0 += pt2_detail(istate,first_det_of_teeth(tooth)) * prop
pt2(1) = E0 + (sumabove(tooth) / Nabove(tooth)) pt2(istate) = E0 + (sumabove(tooth) / Nabove(tooth))
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_pull_socket(zmq_socket_pull) call end_zmq_pull_socket(zmq_socket_pull)

27
src/Determinants/determinants.irp.f
View File

@ -179,6 +179,33 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
END_PROVIDER END_PROVIDER
subroutine update_psi_average_norm_contrib(w)
implicit none
BEGIN_DOC
! Compute psi_average_norm_contrib for different state average weights w(:)
END_DOC
double precision, intent(in) :: w(N_states)
double precision :: w0(N_states), f
w0(:) = w(:)/sum(w(:))
integer :: i,j,k
do i=1,N_det
psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*w(1)
enddo
do k=2,N_states
do i=1,N_det
psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + &
psi_coef(i,k)*psi_coef(i,k)*w(k)
enddo
enddo
f = 1.d0/sum(psi_average_norm_contrib(1:N_det))
do i=1,N_det
psi_average_norm_contrib(i) = psi_average_norm_contrib(i)*f
enddo
SOFT_TOUCH psi_average_norm_contrib
end subroutine
BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ] BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ]
implicit none implicit none