10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00

Merge branch 'master' of github.com:scemama/quantum_package

This commit is contained in:
Anthony Scemama 2017-02-01 11:29:38 +01:00
commit f56246b456
12 changed files with 340 additions and 178 deletions

View File

@ -102,10 +102,10 @@ program fci_zmq
threshold_selectors = 1.d0 threshold_selectors = 1.d0
threshold_generators = 1d0 threshold_generators = 1d0
E_CI_before(1:N_states) = CI_energy(1:N_states) E_CI_before(1:N_states) = CI_energy(1:N_states)
!call ZMQ_selection(0, pt2)! pour non-stochastic
double precision :: relative_error double precision :: relative_error
relative_error=1.d-3 relative_error=1.d-3
call ZMQ_pt2(pt2,relative_error) !call ZMQ_pt2(pt2,relative_error)
call ZMQ_selection(0, pt2)! pour non-stochastic
print *, 'Final step' print *, 'Final step'
print *, 'N_det = ', N_det print *, 'N_det = ', N_det
print *, 'N_states = ', N_states print *, 'N_states = ', N_states

View File

@ -1,7 +1,5 @@
BEGIN_PROVIDER [ integer, fragment_first ]
BEGIN_PROVIDER [ integer, fragment_count ] implicit none
&BEGIN_PROVIDER [ integer, fragment_first ]
fragment_count = (elec_alpha_num-n_core_orb)**2
fragment_first = first_det_of_teeth(1) fragment_first = first_det_of_teeth(1)
END_PROVIDER END_PROVIDER
@ -12,7 +10,7 @@ subroutine ZMQ_pt2(pt2,relative_error)
implicit none implicit none
character*(512) :: task character*(512) :: task
integer(ZMQ_PTR) :: zmq_to_qp_run_socket 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 double precision, intent(in) :: relative_error
@ -29,16 +27,17 @@ subroutine ZMQ_pt2(pt2,relative_error)
double precision, external :: omp_get_wtime double precision, external :: omp_get_wtime
double precision :: time0, time double precision :: time0, time
allocate(pt2_detail(N_states, N_det_generators), comb(100000), computed(N_det_generators), tbc(0:N_det_generators)) allocate(pt2_detail(N_states, N_det_generators), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
sumabove = 0d0 sumabove = 0d0
sum2above = 0d0 sum2above = 0d0
Nabove = 0d0 Nabove = 0d0
provide nproc provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral
call random_seed() !call random_seed()
computed = .false. computed = .false.
tbc(0) = first_det_of_comb - 1 tbc(0) = first_det_of_comb - 1
do i=1, tbc(0) do i=1, tbc(0)
tbc(i) = i tbc(i) = i
@ -46,41 +45,65 @@ subroutine ZMQ_pt2(pt2,relative_error)
end do end do
pt2_detail = 0d0 pt2_detail = 0d0
time0 = omp_get_wtime() time0 = omp_get_wtime()
print *, "grep - time - avg - err - n_combs" print *, "grep - time - avg - err - n_combs"
generator_per_task = 1
do while(.true.) do while(.true.)
call write_time(6)
call new_parallel_job(zmq_to_qp_run_socket,"pt2") call new_parallel_job(zmq_to_qp_run_socket,"pt2")
call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator)) call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
call zmq_set_running(zmq_to_qp_run_socket)
call create_selection_buffer(1, 1*2, b) call create_selection_buffer(1, 1*2, b)
! TODO PARAMETER : 1.d-2 Ncomb=size(comb)
call get_carlo_workbatch(1d-2, computed, comb, Ncomb, tbc) call get_carlo_workbatch(computed, comb, Ncomb, tbc)
generator_per_task = 1
do i=1,tbc(0) call write_time(6)
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2, relative_error) NUM_THREADS(nproc+1) &
!$OMP PRIVATE(i,zmq_to_qp_run_socket2,i_generator_end,task,j)
zmq_to_qp_run_socket2 = new_zmq_to_qp_run_socket()
!$OMP DO SCHEDULE(static,1)
do i=1,min(2000,tbc(0))
i_generator_end = min(i+generator_per_task-1, tbc(0)) i_generator_end = min(i+generator_per_task-1, tbc(0))
if(tbc(i) > fragment_first) then if(tbc(i) > fragment_first) then
integer :: zero write(task,*) (i_generator_end-i+1), 0, tbc(i:i_generator_end)
zero = 0 call add_task_to_taskserver(zmq_to_qp_run_socket2,task)
write(task,*) (i_generator_end-i+1), zero, tbc(i:i_generator_end)
call add_task_to_taskserver(zmq_to_qp_run_socket,task)
else else
do j=1,fragment_count do j=1,fragment_count
write(task,*) (i_generator_end-i+1), j, tbc(i:i_generator_end) write(task,*) (i_generator_end-i+1), j, tbc(i:i_generator_end)
call add_task_to_taskserver(zmq_to_qp_run_socket,task) call add_task_to_taskserver(zmq_to_qp_run_socket2,task)
end do end do
end if end if
end do end do
!$OMP END DO NOWAIT
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2, relative_error) PRIVATE(i) NUM_THREADS(nproc+1)
i = omp_get_thread_num() i = omp_get_thread_num()
if (i==0) then if (i==0) then
call zmq_set_running(zmq_to_qp_run_socket)
call pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, pt2) call pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, pt2)
else if (i==1) then
do i=2001,tbc(0)
i_generator_end = min(i+generator_per_task-1, tbc(0))
if(tbc(i) > fragment_first) then
write(task,*) (i_generator_end-i+1), 0, tbc(i:i_generator_end)
call add_task_to_taskserver(zmq_to_qp_run_socket2,task)
else
do j=1,fragment_count
write(task,*) (i_generator_end-i+1), j, tbc(i:i_generator_end)
call add_task_to_taskserver(zmq_to_qp_run_socket2,task)
end do
end if
end do
call pt2_slave_inproc(1)
else else
call pt2_slave_inproc(i) call pt2_slave_inproc(i)
endif endif
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket2)
!$OMP END PARALLEL !$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, 'pt2') call end_parallel_job(zmq_to_qp_run_socket, 'pt2')
tbc(0) = 0 tbc(0) = 0
@ -93,7 +116,7 @@ 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)
integer, intent(in) :: tbc(0:N_det_generators), Ncomb integer, intent(in) :: tbc(0:size_tbc), Ncomb
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)
@ -103,7 +126,7 @@ subroutine do_carlo(tbc, Ncomb, comb, pt2_detail, computed, sumabove, sum2above,
mainLoop : do i=1,Ncomb mainLoop : do i=1,Ncomb
call get_comb(comb(i), dets) call get_comb(comb(i), dets)
do j=1,comb_teeth do j=1,comb_teeth
if(not(computed(dets(j)))) then if(.not.(computed(dets(j)))) then
exit mainLoop exit mainLoop
end if end if
end do end do
@ -111,7 +134,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)) / weight(dets(j)) * comb_step myVal += pt2_detail(1, dets(j)) * pt2_weight_inv(dets(j)) * comb_step
sumabove(j) += myVal sumabove(j) += myVal
sum2above(j) += myVal**2 sum2above(j) += myVal**2
Nabove(j) += 1 Nabove(j) += 1
@ -138,7 +161,7 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
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) double precision, intent(in) :: comb(Ncomb)
logical, intent(inout) :: computed(N_det_generators) logical, intent(inout) :: computed(N_det_generators)
integer, intent(in) :: tbc(0:N_det_generators) integer, intent(in) :: tbc(0:size_tbc)
double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth), relative_error double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth), relative_error
double precision, intent(out) :: pt2(N_states) double precision, intent(out) :: pt2(N_states)
@ -166,7 +189,7 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
logical, allocatable :: actually_computed(:) logical, allocatable :: actually_computed(:)
allocate(actually_computed(N_det_generators), parts_to_get(N_det_generators)) allocate(actually_computed(N_det_generators), parts_to_get(N_det_generators))
actually_computed = computed actually_computed(:) = computed(:)
parts_to_get(:) = 1 parts_to_get(:) = 1
if(fragment_first > 0) parts_to_get(1:fragment_first) = fragment_count if(fragment_first > 0) parts_to_get(1:fragment_first) = fragment_count
@ -194,6 +217,9 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
pt2_detail(:, index(i)) += pt2_mwen(:,i) pt2_detail(:, index(i)) += pt2_mwen(:,i)
parts_to_get(index(i)) -= 1 parts_to_get(index(i)) -= 1
if(parts_to_get(index(i)) < 0) then if(parts_to_get(index(i)) < 0) then
print *, i, index(i), parts_to_get(index(i)), Nindex
print *, "PARTS ??"
print *, parts_to_get
stop "PARTS ??" stop "PARTS ??"
end if end if
if(parts_to_get(index(i)) == 0) actually_computed(index(i)) = .true. if(parts_to_get(index(i)) == 0) actually_computed(index(i)) = .true.
@ -208,10 +234,10 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
time = omp_get_wtime() time = omp_get_wtime()
if(time - timeLast > 30.0 .or. more /= 1) then if(time - timeLast > 1d1 .or. more /= 1) then
timeLast = time timeLast = time
do i=1, first_det_of_teeth(1)-1 do i=1, first_det_of_teeth(1)-1
if(not(actually_computed(i))) then if(.not.(actually_computed(i))) then
print *, "PT2 : deterministic part not finished" print *, "PT2 : deterministic part not finished"
cycle pullLoop cycle pullLoop
end if end if
@ -220,7 +246,7 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
double precision :: E0, avg, eqt, prop double precision :: E0, avg, eqt, prop
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)
firstTBDcomb = Nabove(1) - orgTBDcomb + 1 firstTBDcomb = Nabove(1) - orgTBDcomb + 1
if(Nabove(1) < 2.0) cycle if(Nabove(1) < 2d0) cycle
call get_first_tooth(actually_computed, tooth) call get_first_tooth(actually_computed, tooth)
done = 0 done = 0
@ -229,8 +255,8 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
end do end do
E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1)) E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1))
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - 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 / weight(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(1,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))
@ -251,16 +277,17 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
end subroutine end subroutine
integer function pt2_find(v, w) integer function pt2_find(v, w, sze)
implicit none implicit none
double precision :: v, w(N_det) integer, intent(in) :: sze
double precision, intent(in) :: v, w(sze)
integer :: i,l,h integer :: i,l,h
l = 0 l = 0
h = N_det-1 h = N_det-1
do while(h >= l) do while(h >= l)
i = (h+l)/2 i = ishft(h+l,-1)
if(w(i+1) > v) then if(w(i+1) > v) then
h = i-1 h = i-1
else else
@ -287,7 +314,7 @@ subroutine get_first_tooth(computed, first_teeth)
first_det = 1 first_det = 1
first_teeth = 1 first_teeth = 1
do i=first_det_of_comb, N_det_generators do i=first_det_of_comb, N_det_generators
if(not(computed(i))) then if(.not.(computed(i))) then
first_det = i first_det = i
exit exit
end if end if
@ -310,10 +337,10 @@ subroutine get_last_full_tooth(computed, last_tooth)
integer :: i, j, missing integer :: i, j, missing
last_tooth = 0 last_tooth = 0
combLoop : do i=comb_teeth-1, 1, -1 combLoop : do i=comb_teeth, 1, -1
missing = 1+ (first_det_of_teeth(i+1)-first_det_of_teeth(i))/100 missing = 1+ ishft(first_det_of_teeth(i+1)-first_det_of_teeth(i),-6) ! /64
do j=first_det_of_teeth(i), first_det_of_teeth(i+1)-1 do j=first_det_of_teeth(i), first_det_of_teeth(i+1)-1
if(not(computed(j))) then if(.not.computed(j)) then
missing -= 1 missing -= 1
if(missing < 0) cycle combLoop if(missing < 0) cycle combLoop
end if end if
@ -324,106 +351,131 @@ subroutine get_last_full_tooth(computed, last_tooth)
end subroutine end subroutine
BEGIN_PROVIDER [ integer, size_tbc ]
size_tbc = N_det_generators + fragment_count*fragment_first
END_PROVIDER
subroutine get_carlo_workbatch(maxWorkload, computed, comb, Ncomb, tbc) subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
implicit none implicit none
double precision, intent(in) :: maxWorkload double precision, intent(out) :: comb(Ncomb)
double precision, intent(out) :: comb(N_det_generators) integer, intent(inout) :: tbc(0:size_tbc)
integer, intent(inout) :: tbc(0:N_det_generators) integer, intent(inout) :: Ncomb
integer, intent(out) :: Ncomb
logical, intent(inout) :: computed(N_det_generators) logical, intent(inout) :: computed(N_det_generators)
integer :: i, j, last_full, dets(comb_teeth) integer :: i, j, last_full, dets(comb_teeth), tbc_save
double precision :: myWorkload
myWorkload = 0d0 call RANDOM_NUMBER(comb)
do j=1,size(comb),100
do i=1,size(comb) do i=j,min(size(comb),j+99)
call RANDOM_NUMBER(comb(i))
comb(i) = comb(i) * comb_step comb(i) = comb(i) * comb_step
call add_comb(comb(i), computed, tbc, myWorkload) tbc_save = tbc(0)
!DIR$ FORCEINLINE
call add_comb(comb(i), computed, tbc, size_tbc, comb_teeth)
if (tbc(0) < size(tbc)) then
Ncomb = i Ncomb = i
else
tbc(0) = tbc_save
return
endif
end do
call get_filling_teeth(computed, tbc)
enddo
end subroutine
subroutine get_filling_teeth(computed, tbc)
implicit none
integer, intent(inout) :: tbc(0:size_tbc)
logical, intent(inout) :: computed(N_det_generators)
integer :: i, j, k, last_full, dets(comb_teeth)
call get_last_full_tooth(computed, last_full) call get_last_full_tooth(computed, last_full)
if(Ncomb >= 30 .and. last_full /= 0) then if(last_full /= 0) then
if (tbc(0) > size(tbc) - first_det_of_teeth(last_full+1) -2) then
return
endif
k = tbc(0)+1
do j=1,first_det_of_teeth(last_full+1)-1 do j=1,first_det_of_teeth(last_full+1)-1
if(not(computed(j))) then if(.not.(computed(j))) then
tbc(0) += 1 tbc(k) = j
tbc(tbc(0)) = j k=k+1
computed(j) = .true. computed(j) = .true.
myWorkload += comb_workload(j) ! print *, "filled ", j, "to reach tooth", last_full, "ending at", first_det_of_teeth(last_full+1)
print *, "filled ", j, "to reach tooth", last_full, "ending at", first_det_of_teeth(last_full+1)
end if end if
end do end do
tbc(0) = k-1
end if end if
if(myWorkload > maxWorkload .and. i >= 100) exit
end do
end subroutine end subroutine
subroutine reorder_tbc(tbc) subroutine reorder_tbc(tbc)
implicit none implicit none
integer, intent(inout) :: tbc(0:N_det_generators) integer, intent(inout) :: tbc(0:size_tbc)
logical, allocatable :: ltbc(:) logical, allocatable :: ltbc(:)
integer :: i, ci integer :: i, ci
allocate(ltbc(N_det_generators)) allocate(ltbc(size_tbc))
ltbc = .false. ltbc(:) = .false.
do i=1,tbc(0) do i=1,tbc(0)
ltbc(tbc(i)) = .true. ltbc(tbc(i)) = .true.
end do end do
ci = 0 ci = 0
do i=1,N_det_generators do i=1,size_tbc
if(ltbc(i)) then if(ltbc(i)) then
ci += 1 ci = ci+1
tbc(ci) = i tbc(ci) = i
end if end if
end do end do
end subroutine end subroutine
subroutine get_comb(stato, dets) subroutine get_comb(stato, dets, ct)
implicit none implicit none
integer, intent(in) :: ct
double precision, intent(in) :: stato double precision, intent(in) :: stato
integer, intent(out) :: dets(comb_teeth) integer, intent(out) :: dets(ct)
double precision :: curs double precision :: curs
integer :: j integer :: j
integer, external :: pt2_find integer, external :: pt2_find
curs = 1d0 - stato curs = 1d0 - stato
do j = comb_teeth, 1, -1 do j = comb_teeth, 1, -1
dets(j) = pt2_find(curs, cweight) !DIR$ FORCEINLINE
dets(j) = pt2_find(curs, pt2_cweight,size(pt2_cweight))
curs -= comb_step curs -= comb_step
end do end do
end subroutine end subroutine
subroutine add_comb(comb, computed, tbc, workload) subroutine add_comb(comb, computed, tbc, stbc, ct)
implicit none implicit none
integer, intent(in) :: stbc, ct
double precision, intent(in) :: comb double precision, intent(in) :: comb
logical, intent(inout) :: computed(N_det_generators) logical, intent(inout) :: computed(N_det_generators)
double precision, intent(inout) :: workload integer, intent(inout) :: tbc(0:stbc)
integer, intent(inout) :: tbc(0:N_det_generators) integer :: i, k, l, dets(ct)
integer :: i, dets(comb_teeth)
call get_comb(comb, dets) !DIR$ FORCEINLINE
call get_comb(comb, dets, ct)
do i = 1, comb_teeth k=tbc(0)+1
if(not(computed(dets(i)))) then do i = 1, ct
tbc(0) += 1 l = dets(i)
tbc(tbc(0)) = dets(i) if(.not.(computed(l))) then
workload += comb_workload(dets(i)) tbc(k) = l
computed(dets(i)) = .true. k = k+1
computed(l) = .true.
end if end if
end do end do
tbc(0) = k-1
end subroutine end subroutine
BEGIN_PROVIDER [ double precision, weight, (N_det_generators) ] BEGIN_PROVIDER [ double precision, pt2_weight, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, cweight, (N_det_generators) ] &BEGIN_PROVIDER [ double precision, pt2_cweight, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, comb_workload, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, comb_step ] &BEGIN_PROVIDER [ double precision, comb_step ]
&BEGIN_PROVIDER [ integer, first_det_of_teeth, (comb_teeth+1) ] &BEGIN_PROVIDER [ integer, first_det_of_teeth, (comb_teeth+1) ]
&BEGIN_PROVIDER [ integer, first_det_of_comb ] &BEGIN_PROVIDER [ integer, first_det_of_comb ]
@ -432,39 +484,53 @@ end subroutine
double precision :: norm_left, stato double precision :: norm_left, stato
integer, external :: pt2_find integer, external :: pt2_find
weight(1) = psi_coef_generators(1,1)**2 pt2_weight(1) = psi_coef_generators(1,1)**2
cweight(1) = psi_coef_generators(1,1)**2 pt2_cweight(1) = psi_coef_generators(1,1)**2
do i=2,N_det_generators do i=2,N_det_generators
weight(i) = psi_coef_generators(i,1)**2 pt2_weight(i) = psi_coef_generators(i,1)**2
cweight(i) = cweight(i-1) + psi_coef_generators(i,1)**2 pt2_cweight(i) = pt2_cweight(i-1) + psi_coef_generators(i,1)**2
end do end do
weight = weight / cweight(N_det_generators) pt2_weight = pt2_weight / pt2_cweight(N_det_generators)
cweight = cweight / cweight(N_det_generators) pt2_cweight = pt2_cweight / pt2_cweight(N_det_generators)
comb_workload = 1d0 / dfloat(N_det_generators)
norm_left = 1d0 norm_left = 1d0
comb_step = 1d0/dfloat(comb_teeth) comb_step = 1d0/dfloat(comb_teeth)
do i=1,N_det_generators do i=1,N_det_generators
if(weight(i)/norm_left < comb_step/2d0) then if(pt2_weight(i)/norm_left < comb_step*.5d0) then
first_det_of_comb = i first_det_of_comb = i
exit exit
end if end if
norm_left -= weight(i) norm_left -= pt2_weight(i)
end do end do
comb_step = 1d0 / dfloat(comb_teeth) * (1d0 - cweight(first_det_of_comb-1)) comb_step = (1d0 - pt2_cweight(first_det_of_comb-1)) * comb_step
stato = 1d0 - comb_step! + 1d-5 stato = 1d0 - comb_step! + 1d-5
do i=comb_teeth, 1, -1 do i=comb_teeth, 1, -1
first_det_of_teeth(i) = pt2_find(stato, cweight) first_det_of_teeth(i) = pt2_find(stato, pt2_cweight, N_det_generators)
stato -= comb_step stato -= comb_step
end do end do
first_det_of_teeth(comb_teeth+1) = N_det_generators + 1 first_det_of_teeth(comb_teeth+1) = N_det_generators + 1
first_det_of_teeth(1) = first_det_of_comb first_det_of_teeth(1) = first_det_of_comb
if(first_det_of_teeth(1) /= first_det_of_comb) stop "comb provider" if(first_det_of_teeth(1) /= first_det_of_comb) then
print *, 'Error in ', irp_here
stop "comb provider"
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, pt2_weight_inv, (N_det_generators) ]
implicit none
BEGIN_DOC
! Inverse of pt2_weight array
END_DOC
integer :: i
do i=1,N_det_generators
pt2_weight_inv(i) = 1.d0/pt2_weight(i)
enddo
END_PROVIDER END_PROVIDER

View File

@ -124,7 +124,7 @@ subroutine push_pt2_results(zmq_socket_push, N, index, pt2_detail, task_id, ntas
if(rc /= 4*ntask) stop "push" if(rc /= 4*ntask) stop "push"
! Activate is zmq_socket_push is a REQ ! Activate is zmq_socket_push is a REQ
! rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0) rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
end subroutine end subroutine
@ -154,7 +154,7 @@ subroutine pull_pt2_results(zmq_socket_pull, N, index, pt2_detail, task_id, ntas
if(rc /= 4*ntask) stop "pull" if(rc /= 4*ntask) stop "pull"
! Activate is zmq_socket_pull is a REP ! Activate is zmq_socket_pull is a REP
! rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0) rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
end subroutine end subroutine

View File

@ -115,7 +115,7 @@ subroutine push_selection_results(zmq_socket_push, pt2, b, task_id, ntask)
if(rc /= 4*ntask) stop "push" if(rc /= 4*ntask) stop "push"
! Activate is zmq_socket_push is a REQ ! Activate is zmq_socket_push is a REQ
! rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0) rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
end subroutine end subroutine
@ -149,7 +149,7 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, val, det, N, task_id, nt
if(rc /= 4*ntask) stop "pull" if(rc /= 4*ntask) stop "pull"
! Activate is zmq_socket_pull is a REP ! Activate is zmq_socket_pull is a REP
! rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0) rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
end subroutine end subroutine

View File

@ -1,5 +1,10 @@
use bitmasks use bitmasks
BEGIN_PROVIDER [ integer, fragment_count ]
implicit none
fragment_count = (elec_alpha_num-n_core_orb)**2
END_PROVIDER
double precision function integral8(i,j,k,l) double precision function integral8(i,j,k,l)
implicit none implicit none
@ -356,20 +361,6 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
integer :: nb_count integer :: nb_count
do s1=1,2 do s1=1,2
do i1=N_holes(s1),1,-1 ! Generate low excitations first do i1=N_holes(s1),1,-1 ! Generate low excitations first
! will_compute = (subset == 0)
! nb_count = 0
! if (s1==1) then
! nb_count = N_holes(1)-i1 + N_holes(2)
! else
! nb_count = N_holes(2)-i1
! endif
! maskInd = 12345
! fragment_count = 400
! subset = 3
! nb_count = 100
! if( nb_count >= (fragment_count - mod(maskInd+1, fragment_count) + subset-1) ) then
! will_compute = .true.
! end if
h1 = hole_list(i1,s1) h1 = hole_list(i1,s1)
call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int) call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int)

View File

@ -77,6 +77,7 @@ END_PROVIDER
norm_psi_ref(j) += psi_ref_coef(i,j) * psi_ref_coef(i,j) norm_psi_ref(j) += psi_ref_coef(i,j) * psi_ref_coef(i,j)
enddo enddo
inv_norm_psi_ref(j) = 1.d0/(dsqrt(norm_psi_Ref(j))) inv_norm_psi_ref(j) = 1.d0/(dsqrt(norm_psi_Ref(j)))
print *, inv_norm_psi_ref(j)
enddo enddo
END_PROVIDER END_PROVIDER

View File

@ -4,44 +4,95 @@
implicit none implicit none
integer :: i,j,k,l integer :: i,j,k,l
integer, allocatable :: idx(:) integer, allocatable :: idx(:)
integer, allocatable :: holes_part(:,:)
double precision, allocatable :: e_corr(:,:) double precision, allocatable :: e_corr(:,:)
double precision, allocatable :: accu(:) double precision, allocatable :: accu(:)
double precision, allocatable :: ihpsi_current(:) double precision, allocatable :: ihpsi_current(:)
double precision, allocatable :: H_jj(:),H_jj_total(:),S2_jj(:) double precision, allocatable :: H_jj(:),H_jj_total(:),S2_jj(:)
integer :: number_of_particles, number_of_holes, n_h,n_p
allocate(e_corr(N_det_non_ref,N_states),ihpsi_current(N_states),accu(N_states),H_jj(N_det_non_ref),idx(0:N_det_non_ref)) allocate(e_corr(N_det_non_ref,N_states),ihpsi_current(N_states),accu(N_states),H_jj(N_det_non_ref),idx(0:N_det_non_ref))
allocate(H_jj_total(N_det),S2_jj(N_det)) allocate(H_jj_total(N_det),S2_jj(N_det))
allocate(holes_part(N_det,2))
accu = 0.d0 accu = 0.d0
do i = 1, N_det_non_ref do i = 1, N_det_non_ref
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef_interm_norm, N_int, N_det_ref,& holes_part(i,1) = number_of_holes(psi_non_ref(1,1,i))
holes_part(i,2) = number_of_particles(psi_non_ref(1,1,i))
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
size(psi_ref_coef_interm_norm,1), N_states,ihpsi_current) size(psi_ref_coef_interm_norm,1), N_states,ihpsi_current)
do j = 1, N_states do j = 1, N_states
e_corr(i,j) = psi_non_ref_coef_interm_norm(i,j) * ihpsi_current(j) e_corr(i,j) = psi_non_ref_coef(i,j) * ihpsi_current(j) * inv_norm_psi_ref(j)
accu(j) += e_corr(i,j) accu(j) += e_corr(i,j)
enddo enddo
enddo enddo
print *, 'accu = ',accu
double precision :: hjj,diag_h_mat_elem double precision :: hjj,diag_h_mat_elem
do i = 1, N_det_non_ref do i = 1, N_det_non_ref
call filter_not_connected(psi_non_ref,psi_non_ref(1,1,i),N_int,N_det_non_ref,idx)
H_jj(i) = 0.d0 H_jj(i) = 0.d0
n_h = holes_part(i,1)
n_p = holes_part(i,2)
integer :: degree
! do j = 1, N_det_non_ref
! call get_excitation_degree(psi_non_ref(1,1,i),psi_non_ref(1,1,j),degree,N_int)
! if(degree .gt. 2)then
! if(n_h + holes_part(j,1) .gt. 2 .or. n_p + holes_part(j,2) .gt. 2 ) then
! H_jj(i) += e_corr(j,1)
! endif
! endif
! enddo
call filter_not_connected(psi_non_ref,psi_non_ref(1,1,i),N_int,N_det_non_ref,idx)
do j = 1, idx(0) do j = 1, idx(0)
if(n_h + holes_part(idx(j),1) .gt. 2 .or. n_p + holes_part(idx(j),2) .gt. 2 ) then
H_jj(i) += e_corr(idx(j),1) H_jj(i) += e_corr(idx(j),1)
endif
enddo enddo
enddo enddo
do i=1,N_Det do i=1,N_Det
H_jj_total(i) = diag_h_mat_elem(psi_det(1,1,i),N_int) H_jj_total(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
call get_s2(psi_det(1,1,i),psi_det(1,1,i),N_int,S2_jj(i)) call get_s2(psi_det(1,1,i),psi_det(1,1,i),N_int,S2_jj(i))
enddo enddo
do i=1, N_det_non_ref do i = 1, N_det_non_ref
H_jj_total(idx_non_ref(i)) += H_jj(i) H_jj_total(idx_non_ref(i)) += H_jj(i)
enddo enddo
print *, 'coef'
call davidson_diag_hjj_sjj(psi_det,CI_eigenvectors_sc2_no_amp,H_jj_total,S2_jj,CI_electronic_energy_sc2_no_amp,size(CI_eigenvectors_sc2_no_amp,1),N_Det,N_states,N_states_diag,N_int,6) call davidson_diag_hjj_sjj(psi_det,CI_eigenvectors_sc2_no_amp,H_jj_total,S2_jj,CI_electronic_energy_sc2_no_amp,size(CI_eigenvectors_sc2_no_amp,1),N_Det,N_states,N_states_diag,N_int,6)
do i = 1, N_det
hjj = diag_h_mat_elem(psi_det(1,1,i),N_int)
! if(hjj<-210.d0)then
! call debug_det(psi_det(1,1,i),N_int)
! print *, CI_eigenvectors_sc2_no_amp((i),1),hjj, H_jj_total(i)
! endif
enddo
print *, 'ref',N_det_ref
do i =1, N_det_ref
call debug_det(psi_det(1,1,idx_ref(i)),N_int)
print *, CI_eigenvectors_sc2_no_amp(idx_ref(i),1), H_jj_total(idx_ref(i))
enddo
print *, 'non ref',N_det_non_ref
do i=1, N_det_non_ref
hjj = diag_h_mat_elem(psi_non_ref(1,1,i),N_int)
! print *, CI_eigenvectors_sc2_no_amp(idx_non_ref(i),1),H_jj_total(idx_non_ref(i)), H_jj(i)
! if(dabs(CI_eigenvectors_sc2_no_amp(idx_non_ref(i),1)).gt.1.d-1)then
! if(hjj<-210.d0)then
! call debug_det(psi_det(1,1,idx_non_ref(i)),N_int)
! write(*,'(10(F16.10,X))') CI_eigenvectors_sc2_no_amp(idx_non_ref(i),1),hjj, H_jj(i),H_jj_total(idx_non_ref(i))
! endif
enddo
! do i = 1, N_det
! print *, CI_eigenvectors_sc2_no_amp(i,1)
! enddo
do i=1,N_states_diag do i=1,N_states_diag
CI_eigenvectors_s2_sc2_no_amp(i) = S2_jj(i) CI_eigenvectors_s2_sc2_no_amp(i) = S2_jj(i)
enddo enddo
deallocate(e_corr,ihpsi_current,accu,H_jj,idx,H_jj_total,s2_jj) deallocate(e_corr,ihpsi_current,accu,H_jj,idx,H_jj_total,s2_jj,holes_part)
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_energy_sc2_no_amp, (N_states_diag) ] BEGIN_PROVIDER [ double precision, CI_energy_sc2_no_amp, (N_states_diag) ]

View File

@ -1,9 +1,14 @@
program pouet program pouet
provide ao_bielec_integrals_in_map
call bla
end
subroutine bla
implicit none implicit none
integer :: i integer :: i
do i = 1, 10 do i = 1, 10
call diagonalize_CI_sc2_no_amp call diagonalize_CI_sc2_no_amp
TOUCH psi_coef TOUCH psi_coef
enddo enddo
print *, "E+PT2 = ", ci_energy_sc2_no_amp(:)
end end

View File

@ -28,32 +28,32 @@ subroutine routine
if(degree == 0)then if(degree == 0)then
print*,'Reference determinant ' print*,'Reference determinant '
else else
call i_H_j(psi_det(1,1,i),psi_det(1,1,1),N_int,hij) call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij)
call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int) call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*,'phase = ',phase print*,'phase = ',phase
if(degree == 1)then ! if(degree == 1)then
print*,'s1',s1 ! print*,'s1',s1
print*,'h1,p1 = ',h1,p1 ! print*,'h1,p1 = ',h1,p1
if(s1 == 1)then ! if(s1 == 1)then
norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1)) ! norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1))
else ! else
norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1)) ! norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
endif ! endif
print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,list_act(1),list_act(1),p1,mo_integrals_map) ! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,list_act(1),list_act(1),p1,mo_integrals_map)
double precision :: hmono,hdouble ! double precision :: hmono,hdouble
call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble) ! call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble)
print*,'hmono = ',hmono ! print*,'hmono = ',hmono
print*,'hdouble = ',hdouble ! print*,'hdouble = ',hdouble
print*,'hmono+hdouble = ',hmono+hdouble ! print*,'hmono+hdouble = ',hmono+hdouble
print*,'hij = ',hij ! print*,'hij = ',hij
else ! else
print*,'s1',s1 ! print*,'s1',s1
print*,'h1,p1 = ',h1,p1 ! print*,'h1,p1 = ',h1,p1
print*,'s2',s2 ! print*,'s2',s2
print*,'h2,p2 = ',h2,p2 ! print*,'h2,p2 = ',h2,p2
print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map) ! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
endif ! endif
print*,'<Ref| H |D_I> = ',hij print*,'<Ref| H |D_I> = ',hij
endif endif

View File

@ -57,12 +57,12 @@ subroutine push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value,
endif endif
! Activate is zmq_socket_push is a REQ ! Activate is zmq_socket_push is a REQ
! integer :: idummy integer :: idummy
! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0) rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
! if (rc /= 4) then if (rc /= 4) then
! print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)' print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
! stop 'error' stop 'error'
! endif endif
end end
@ -187,11 +187,11 @@ subroutine ao_bielec_integrals_in_map_collector
rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
! Activate if zmq_socket_pull is a REP ! Activate if zmq_socket_pull is a REP
! rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0) rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
! if (rc /= 4) then if (rc /= 4) then
! print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...' print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
! stop 'error' stop 'error'
! endif endif
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)

View File

@ -73,10 +73,11 @@ subroutine map_load_from_disk(filename,map)
implicit none implicit none
character*(*), intent(in) :: filename character*(*), intent(in) :: filename
type(map_type), intent(inout) :: map type(map_type), intent(inout) :: map
double precision :: x
type(c_ptr) :: c_pointer(3) type(c_ptr) :: c_pointer(3)
integer :: fd(3) integer :: fd(3)
integer*8 :: i,k integer*8 :: i,k, l
integer :: n_elements integer :: n_elements, j
@ -95,7 +96,9 @@ subroutine map_load_from_disk(filename,map)
call mmap(trim(filename)//'_consolidated_value', (/ map % n_elements /), integral_kind, fd(3), .True., c_pointer(3)) call mmap(trim(filename)//'_consolidated_value', (/ map % n_elements /), integral_kind, fd(3), .True., c_pointer(3))
call c_f_pointer(c_pointer(3),map % consolidated_value, (/ map % n_elements /)) call c_f_pointer(c_pointer(3),map % consolidated_value, (/ map % n_elements /))
l = 0_8
k = 1_8 k = 1_8
x = 0.d0
do i=0_8, map % map_size do i=0_8, map % map_size
deallocate(map % map(i) % value) deallocate(map % map(i) % value)
deallocate(map % map(i) % key) deallocate(map % map(i) % key)
@ -106,9 +109,15 @@ subroutine map_load_from_disk(filename,map)
k = map % consolidated_idx (i+2) k = map % consolidated_idx (i+2)
map % map(i) % map_size = n_elements map % map(i) % map_size = n_elements
map % map(i) % n_elements = n_elements map % map(i) % n_elements = n_elements
! Load memory from disk
do j=1,n_elements
x = x + map % map(i) % value(j)
l = iand(l,map % map(i) % key(j))
enddo enddo
enddo
map % sorted = x>0 .or. l == 0_8
map % n_elements = k-1 map % n_elements = k-1
map % sorted = .True. map % sorted = map % sorted .or. .True.
map % consolidated = .True. map % consolidated = .True.
end end

View File

@ -235,8 +235,8 @@ function new_zmq_pull_socket()
if (zmq_context == 0_ZMQ_PTR) then if (zmq_context == 0_ZMQ_PTR) then
stop 'zmq_context is uninitialized' stop 'zmq_context is uninitialized'
endif endif
new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL) ! new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL)
! new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_REP) new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_REP)
call omp_unset_lock(zmq_lock) call omp_unset_lock(zmq_lock)
if (new_zmq_pull_socket == 0_ZMQ_PTR) then if (new_zmq_pull_socket == 0_ZMQ_PTR) then
stop 'Unable to create zmq pull socket' stop 'Unable to create zmq pull socket'
@ -312,8 +312,8 @@ function new_zmq_push_socket(thread)
if (zmq_context == 0_ZMQ_PTR) then if (zmq_context == 0_ZMQ_PTR) then
stop 'zmq_context is uninitialized' stop 'zmq_context is uninitialized'
endif endif
new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_PUSH) ! new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_PUSH)
! new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_REQ) new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_REQ)
call omp_unset_lock(zmq_lock) call omp_unset_lock(zmq_lock)
if (new_zmq_push_socket == 0_ZMQ_PTR) then if (new_zmq_push_socket == 0_ZMQ_PTR) then
stop 'Unable to create zmq push socket' stop 'Unable to create zmq push socket'
@ -696,8 +696,7 @@ subroutine add_task_to_taskserver(zmq_to_qp_run_socket,task)
endif endif
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0) rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
message = trim(message(1:rc)) if (message(1:rc) /= 'ok') then
if (trim(message) /= 'ok') then
print *, trim(task) print *, trim(task)
print *, 'Unable to add the next task' print *, 'Unable to add the next task'
stop -1 stop -1
@ -705,6 +704,47 @@ subroutine add_task_to_taskserver(zmq_to_qp_run_socket,task)
end end
subroutine add_task_to_taskserver_send(zmq_to_qp_run_socket,task)
use f77_zmq
implicit none
BEGIN_DOC
! Get a task from the task server
END_DOC
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
character*(*), intent(in) :: task
integer :: rc, sze
character*(512) :: message
write(message,*) 'add_task '//trim(zmq_state)//' '//trim(task)
sze = len(trim(message))
rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
if (rc /= sze) then
print *, rc, sze
print *, irp_here,': f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)'
stop 'error'
endif
end
subroutine add_task_to_taskserver_recv(zmq_to_qp_run_socket)
use f77_zmq
implicit none
BEGIN_DOC
! Get a task from the task server
END_DOC
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
integer :: rc, sze
character*(512) :: message
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
if (message(1:rc) /= 'ok') then
print *, 'Unable to add the next task'
stop -1
endif
end
subroutine task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_id) subroutine task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_id)
use f77_zmq use f77_zmq
implicit none implicit none
@ -726,8 +766,7 @@ subroutine task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_id)
endif endif
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0) rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
message = trim(message(1:rc)) if (trim(message(1:rc)) /= 'ok') then
if (trim(message) /= 'ok') then
print *, 'Unable to send task_done message' print *, 'Unable to send task_done message'
stop -1 stop -1
endif endif