missing files.......
This commit is contained in:
Yann Garniron
parent
b5b333904f
commit
2f2840ebb5
|
|
@ -0,0 +1,23 @@
|
|||
BEGIN_PROVIDER [ logical, initialize_mrcc_E0_denominator ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! If true, initialize mrcc_E0_denominator
|
||||
END_DOC
|
||||
initialize_mrcc_E0_denominator = .True.
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, mrcc_E0_denominator, (N_states) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! E0 in the denominator of the mrcc
|
||||
END_DOC
|
||||
if (initialize_mrcc_E0_denominator) then
|
||||
mrcc_E0_denominator(1:N_states) = psi_energy(1:N_states)
|
||||
! mrcc_E0_denominator(1:N_states) = HF_energy - nuclear_repulsion
|
||||
! mrcc_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
|
||||
call write_double(6,mrcc_E0_denominator(1)+nuclear_repulsion, 'mrcc Energy denominator')
|
||||
else
|
||||
mrcc_E0_denominator = -huge(1.d0)
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
|
|
@ -0,0 +1,80 @@
|
|||
! DO NOT MODIFY BY HAND
|
||||
! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
|
||||
! from file /home/garniron/quantum_package/src/mrcepa0/EZFIO.cfg
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ logical, perturbative_triples ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Compute perturbative contribution of the Triples
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
|
||||
call ezfio_has_mrcepa0_perturbative_triples(has)
|
||||
if (has) then
|
||||
call ezfio_get_mrcepa0_perturbative_triples(perturbative_triples)
|
||||
else
|
||||
print *, 'mrcepa0/perturbative_triples not found in EZFIO file'
|
||||
stop 1
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, thresh_dressed_ci ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Threshold on the convergence of the dressed CI energy
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
|
||||
call ezfio_has_mrcepa0_thresh_dressed_ci(has)
|
||||
if (has) then
|
||||
call ezfio_get_mrcepa0_thresh_dressed_ci(thresh_dressed_ci)
|
||||
else
|
||||
print *, 'mrcepa0/thresh_dressed_ci not found in EZFIO file'
|
||||
stop 1
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, n_it_max_dressed_ci ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Maximum number of dressed CI iterations
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
|
||||
call ezfio_has_mrcepa0_n_it_max_dressed_ci(has)
|
||||
if (has) then
|
||||
call ezfio_get_mrcepa0_n_it_max_dressed_ci(n_it_max_dressed_ci)
|
||||
else
|
||||
print *, 'mrcepa0/n_it_max_dressed_ci not found in EZFIO file'
|
||||
stop 1
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, lambda_type ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! lambda type
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
|
||||
call ezfio_has_mrcepa0_lambda_type(has)
|
||||
if (has) then
|
||||
call ezfio_get_mrcepa0_lambda_type(lambda_type)
|
||||
else
|
||||
print *, 'mrcepa0/lambda_type not found in EZFIO file'
|
||||
stop 1
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
|
@ -0,0 +1,76 @@
|
|||
program mrcc_slave
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Helper program to compute the mrcc in distributed mode.
|
||||
END_DOC
|
||||
|
||||
read_wf = .False.
|
||||
distributed_davidson = .False.
|
||||
SOFT_TOUCH read_wf distributed_davidson
|
||||
call provide_everything
|
||||
call switch_qp_run_to_master
|
||||
call run_wf
|
||||
end
|
||||
|
||||
subroutine provide_everything
|
||||
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
|
||||
end
|
||||
|
||||
subroutine run_wf
|
||||
use f77_zmq
|
||||
implicit none
|
||||
|
||||
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
double precision :: energy(N_states_diag)
|
||||
character*(64) :: states(1)
|
||||
integer :: rc, i
|
||||
|
||||
call provide_everything
|
||||
|
||||
zmq_context = f77_zmq_ctx_new ()
|
||||
states(1) = 'mrcc'
|
||||
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
|
||||
do
|
||||
|
||||
call wait_for_states(states,zmq_state,1)
|
||||
|
||||
if(trim(zmq_state) == 'Stopped') then
|
||||
|
||||
exit
|
||||
|
||||
else if (trim(zmq_state) == 'mrcc') then
|
||||
|
||||
! Selection
|
||||
! ---------
|
||||
|
||||
print *, 'mrcc'
|
||||
call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
|
||||
|
||||
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
|
||||
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
|
||||
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
|
||||
PROVIDE psi_bilinear_matrix_transp_order
|
||||
|
||||
!$OMP PARALLEL PRIVATE(i)
|
||||
i = omp_get_thread_num()
|
||||
call mrcc_slave_tcp(i, energy)
|
||||
!$OMP END PARALLEL
|
||||
print *, 'mrcc done'
|
||||
|
||||
endif
|
||||
|
||||
end do
|
||||
end
|
||||
|
||||
subroutine mrcc_slave_tcp(i,energy)
|
||||
implicit none
|
||||
double precision, intent(in) :: energy(N_states_diag)
|
||||
integer, intent(in) :: i
|
||||
logical :: lstop
|
||||
lstop = .False.
|
||||
call run_mrcc_slave(0,i,energy,lstop)
|
||||
end
|
||||
|
||||
|
|
@ -0,0 +1,38 @@
|
|||
program mrcc_stoch
|
||||
implicit none
|
||||
read_wf = .True.
|
||||
SOFT_TOUCH read_wf
|
||||
PROVIDE mo_bielec_integrals_in_map
|
||||
call run
|
||||
end
|
||||
|
||||
subroutine run
|
||||
implicit none
|
||||
integer :: i,j,k
|
||||
logical, external :: detEq
|
||||
|
||||
double precision, allocatable :: mrcc(:)
|
||||
integer :: degree
|
||||
integer :: n_det_before, to_select
|
||||
double precision :: threshold_davidson_in
|
||||
|
||||
double precision :: E_CI_before, relative_error
|
||||
|
||||
allocate (mrcc(N_states))
|
||||
mrcc = 0.d0
|
||||
|
||||
E_CI_before = mrcc_E0_denominator(1) + nuclear_repulsion
|
||||
threshold_selectors = 1.d0
|
||||
threshold_generators = 1d0
|
||||
relative_error = 1.d-3
|
||||
call ZMQ_mrcc(E_CI_before, mrcc, relative_error)
|
||||
!print *, 'Final step'
|
||||
!print *, 'N_det = ', N_det
|
||||
print *, 'mrcc = ', mrcc
|
||||
!print *, 'E = ', E_CI_before
|
||||
!print *, 'E+mrcc = ', E_CI_before+mrcc
|
||||
!print *, '-----'
|
||||
!call ezfio_set_full_ci_zmq_energy_mrcc(E_CI_before+mrcc(1))
|
||||
end
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,530 @@
|
|||
BEGIN_PROVIDER [ integer, fragment_first ]
|
||||
implicit none
|
||||
fragment_first = first_det_of_teeth(1)
|
||||
END_PROVIDER
|
||||
|
||||
subroutine ZMQ_mrcc(E, mrcc,relative_error)
|
||||
use f77_zmq
|
||||
use selection_types
|
||||
|
||||
implicit none
|
||||
|
||||
character(len=64000) :: task
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_to_qp_run_socket2
|
||||
type(selection_buffer) :: b
|
||||
integer, external :: omp_get_thread_num
|
||||
double precision, intent(in) :: relative_error, E
|
||||
double precision, intent(out) :: mrcc(N_states)
|
||||
|
||||
|
||||
double precision, allocatable :: mrcc_detail(:,:), comb(:)
|
||||
logical, allocatable :: computed(:)
|
||||
integer, allocatable :: tbc(:)
|
||||
integer :: i, j, k, Ncomb, generator_per_task, i_generator_end
|
||||
integer, external :: mrcc_find
|
||||
|
||||
double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
|
||||
double precision, external :: omp_get_wtime
|
||||
double precision :: time
|
||||
|
||||
!if (N_det < 10) then
|
||||
! !call ZMQ_selection(0, mrcc)
|
||||
! return
|
||||
!else
|
||||
|
||||
allocate(mrcc_detail(N_states, N_det_generators+1), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
|
||||
sumabove = 0d0
|
||||
sum2above = 0d0
|
||||
Nabove = 0d0
|
||||
|
||||
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral mrcc_weight psi_selectors
|
||||
|
||||
computed = .false.
|
||||
|
||||
tbc(0) = first_det_of_comb - 1
|
||||
do i=1, tbc(0)
|
||||
tbc(i) = i
|
||||
computed(i) = .true.
|
||||
end do
|
||||
|
||||
mrcc_detail = 0d0
|
||||
generator_per_task = 1
|
||||
print *, '========== ================= ================= ================='
|
||||
print *, ' Samples Energy Stat. Error Seconds '
|
||||
print *, '========== ================= ================= ================='
|
||||
|
||||
call new_parallel_job(zmq_to_qp_run_socket,'mrcc')
|
||||
call zmq_put_psi(zmq_to_qp_run_socket,1,mrcc_e0_denominator,size(mrcc_e0_denominator))
|
||||
call create_selection_buffer(1, 1*2, b)
|
||||
|
||||
Ncomb=size(comb)
|
||||
call get_carlo_workbatch(computed, comb, Ncomb, tbc)
|
||||
|
||||
do i=1,comb_teeth
|
||||
print *, "TOOTH", first_det_of_teeth(i+1) - first_det_of_teeth(i)
|
||||
end do
|
||||
!stop
|
||||
|
||||
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
|
||||
integer :: ipos
|
||||
ipos=1
|
||||
|
||||
do i=1,tbc(0)
|
||||
if(tbc(i) > fragment_first) then
|
||||
write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') 0, tbc(i)
|
||||
ipos += 20
|
||||
if (ipos > 63980) then
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos)))
|
||||
ipos=1
|
||||
endif
|
||||
else
|
||||
do j=1,fragment_count
|
||||
write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') j, tbc(i)
|
||||
ipos += 20
|
||||
if (ipos > 63980) then
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos)))
|
||||
ipos=1
|
||||
endif
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
if (ipos > 1) then
|
||||
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos)))
|
||||
endif
|
||||
|
||||
call zmq_set_running(zmq_to_qp_run_socket)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc+1) &
|
||||
!$OMP PRIVATE(i)
|
||||
i = omp_get_thread_num()
|
||||
if (i==0) then
|
||||
call mrcc_collector(E, b, tbc, comb, Ncomb, computed, mrcc_detail, sumabove, sum2above, Nabove, relative_error, mrcc)
|
||||
else
|
||||
call mrcc_slave_inproc(i)
|
||||
endif
|
||||
!$OMP END PARALLEL
|
||||
call delete_selection_buffer(b)
|
||||
call end_parallel_job(zmq_to_qp_run_socket, 'mrcc')
|
||||
|
||||
print *, '========== ================= ================= ================='
|
||||
|
||||
!endif
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine do_carlo(tbc, Ncomb, comb, mrcc_detail, computed, sumabove, sum2above, Nabove)
|
||||
integer, intent(in) :: tbc(0:size_tbc), Ncomb
|
||||
logical, intent(in) :: computed(N_det_generators)
|
||||
double precision, intent(in) :: comb(Ncomb), mrcc_detail(N_states, N_det_generators)
|
||||
double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
|
||||
integer :: i, dets(comb_teeth)
|
||||
double precision :: myVal, myVal2
|
||||
|
||||
mainLoop : do i=1,Ncomb
|
||||
call get_comb(comb(i), dets, comb_teeth)
|
||||
do j=1,comb_teeth
|
||||
if(.not.(computed(dets(j)))) then
|
||||
exit mainLoop
|
||||
end if
|
||||
end do
|
||||
|
||||
myVal = 0d0
|
||||
myVal2 = 0d0
|
||||
do j=comb_teeth,1,-1
|
||||
myVal += mrcc_detail(1, dets(j)) * mrcc_weight_inv(dets(j)) * comb_step
|
||||
sumabove(j) += myVal
|
||||
sum2above(j) += myVal*myVal
|
||||
Nabove(j) += 1
|
||||
end do
|
||||
end do mainLoop
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine mrcc_slave_inproc(i)
|
||||
implicit none
|
||||
integer, intent(in) :: i
|
||||
|
||||
call run_mrcc_slave(1,i,mrcc_e0_denominator)
|
||||
end
|
||||
|
||||
subroutine mrcc_collector(E, b, tbc, comb, Ncomb, computed, mrcc_detail, sumabove, sum2above, Nabove, relative_error, mrcc)
|
||||
use f77_zmq
|
||||
use selection_types
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
|
||||
integer, intent(in) :: Ncomb
|
||||
double precision, intent(inout) :: mrcc_detail(N_states, N_det_generators)
|
||||
double precision, intent(in) :: comb(Ncomb), relative_error, E
|
||||
logical, intent(inout) :: computed(N_det_generators)
|
||||
integer, intent(in) :: tbc(0:size_tbc)
|
||||
double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
|
||||
double precision, intent(out) :: mrcc(N_states)
|
||||
|
||||
|
||||
type(selection_buffer), intent(inout) :: b
|
||||
double precision, allocatable :: mrcc_mwen(:,:)
|
||||
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
|
||||
integer(ZMQ_PTR), external :: new_zmq_pull_socket
|
||||
integer(ZMQ_PTR) :: zmq_socket_pull
|
||||
|
||||
integer :: msg_size, rc, more
|
||||
integer :: acc, i, j, robin, N, ntask
|
||||
double precision, allocatable :: val(:)
|
||||
integer(bit_kind), allocatable :: det(:,:,:)
|
||||
integer, allocatable :: task_id(:)
|
||||
integer :: Nindex
|
||||
integer, allocatable :: index(:)
|
||||
double precision, save :: time0 = -1.d0
|
||||
double precision :: time, timeLast, Nabove_old
|
||||
double precision, external :: omp_get_wtime
|
||||
integer :: tooth, firstTBDcomb, orgTBDcomb
|
||||
integer, allocatable :: parts_to_get(:)
|
||||
logical, allocatable :: actually_computed(:)
|
||||
integer :: ncomputed
|
||||
|
||||
ncomputed = 0
|
||||
|
||||
character*(512) :: task
|
||||
Nabove_old = -1.d0
|
||||
|
||||
allocate(actually_computed(N_det_generators), parts_to_get(N_det_generators), &
|
||||
mrcc_mwen(N_states, N_det_generators) )
|
||||
mrcc_mwen(1:N_states, 1:N_det_generators) =0.d0
|
||||
do i=1,N_det_generators
|
||||
actually_computed(i) = computed(i)
|
||||
enddo
|
||||
|
||||
parts_to_get(:) = 1
|
||||
if(fragment_first > 0) then
|
||||
do i=1,fragment_first
|
||||
parts_to_get(i) = fragment_count
|
||||
enddo
|
||||
endif
|
||||
|
||||
do i=1,tbc(0)
|
||||
actually_computed(tbc(i)) = .false.
|
||||
end do
|
||||
|
||||
orgTBDcomb = int(Nabove(1))
|
||||
firstTBDcomb = 1
|
||||
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
zmq_socket_pull = new_zmq_pull_socket()
|
||||
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det_generators), index(1))
|
||||
more = 1
|
||||
if (time0 < 0.d0) then
|
||||
call wall_time(time0)
|
||||
endif
|
||||
timeLast = time0
|
||||
|
||||
call get_first_tooth(actually_computed, tooth)
|
||||
Nabove_old = Nabove(tooth)
|
||||
|
||||
pullLoop : do while (more == 1)
|
||||
|
||||
call pull_mrcc_results(zmq_socket_pull, Nindex, index, mrcc_mwen, task_id, ntask)
|
||||
do i=1,Nindex
|
||||
mrcc_detail(1:N_states, index(i)) += mrcc_mwen(1:N_states,i)
|
||||
parts_to_get(index(i)) -= 1
|
||||
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 ??"
|
||||
end if
|
||||
if(parts_to_get(index(i)) == 0) then
|
||||
actually_computed(index(i)) = .true.
|
||||
ncomputed += 1
|
||||
end if
|
||||
end do
|
||||
|
||||
do i=1, ntask
|
||||
if(task_id(i) == 0) then
|
||||
print *, "Error in collector"
|
||||
endif
|
||||
call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more)
|
||||
end do
|
||||
|
||||
time = omp_get_wtime()
|
||||
|
||||
if(time - timeLast > 1d0 .or. more /= 1) then
|
||||
timeLast = time
|
||||
do i=1, first_det_of_teeth(1)-1
|
||||
if(.not.(actually_computed(i))) then
|
||||
cycle pullLoop
|
||||
end if
|
||||
end do
|
||||
|
||||
double precision :: E0, avg, eqt, prop
|
||||
call do_carlo(tbc, Ncomb+1-firstTBDcomb, comb(firstTBDcomb), mrcc_detail, actually_computed, sumabove, sum2above, Nabove)
|
||||
firstTBDcomb = int(Nabove(1)) - orgTBDcomb + 1
|
||||
if(Nabove(1) < 5d0) cycle
|
||||
call get_first_tooth(actually_computed, tooth)
|
||||
|
||||
E0 = sum(mrcc_detail(1,:first_det_of_teeth(tooth)-1))
|
||||
if (tooth <= comb_teeth) then
|
||||
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - mrcc_cweight(first_det_of_teeth(tooth)-1))
|
||||
prop = prop * mrcc_weight_inv(first_det_of_teeth(tooth))
|
||||
E0 += mrcc_detail(1,first_det_of_teeth(tooth)) * prop
|
||||
avg = E0 + (sumabove(tooth) / Nabove(tooth))
|
||||
eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
|
||||
else
|
||||
eqt = 0.d0
|
||||
endif
|
||||
call wall_time(time)
|
||||
if (dabs(eqt/avg) < relative_error) then
|
||||
! Termination
|
||||
mrcc(1) = avg
|
||||
print '(G10.3, 2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, ''
|
||||
print *, avg
|
||||
call zmq_abort(zmq_to_qp_run_socket)
|
||||
else
|
||||
if (Nabove(tooth) > Nabove_old) then
|
||||
print '(G10.3, 2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, ''
|
||||
print *, avg
|
||||
Nabove_old = Nabove(tooth)
|
||||
endif
|
||||
endif
|
||||
end if
|
||||
end do pullLoop
|
||||
|
||||
E0 = sum(mrcc_detail(1,:first_det_of_teeth(tooth)-1))
|
||||
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - mrcc_cweight(first_det_of_teeth(tooth)-1))
|
||||
prop = prop * mrcc_weight_inv(first_det_of_teeth(tooth))
|
||||
E0 += mrcc_detail(1,first_det_of_teeth(tooth)) * prop
|
||||
mrcc(1) = E0 + (sumabove(tooth) / Nabove(tooth))
|
||||
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
call end_zmq_pull_socket(zmq_socket_pull)
|
||||
call sort_selection_buffer(b)
|
||||
end subroutine
|
||||
|
||||
integer function mrcc_find(v, w, sze, imin, imax)
|
||||
implicit none
|
||||
integer, intent(in) :: sze, imin, imax
|
||||
double precision, intent(in) :: v, w(sze)
|
||||
integer :: i,l,h
|
||||
integer, parameter :: block=64
|
||||
|
||||
l = imin
|
||||
h = imax-1
|
||||
|
||||
do while(h-l >= block)
|
||||
i = ishft(h+l,-1)
|
||||
if(w(i+1) > v) then
|
||||
h = i-1
|
||||
else
|
||||
l = i+1
|
||||
end if
|
||||
end do
|
||||
!DIR$ LOOP COUNT (64)
|
||||
do mrcc_find=l,h
|
||||
if(w(mrcc_find) >= v) then
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
end function
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, comb_teeth ]
|
||||
implicit none
|
||||
comb_teeth = 100
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
subroutine get_first_tooth(computed, first_teeth)
|
||||
implicit none
|
||||
logical, intent(in) :: computed(N_det_generators)
|
||||
integer, intent(out) :: first_teeth
|
||||
integer :: i, first_det
|
||||
|
||||
first_det = 1
|
||||
first_teeth = 1
|
||||
do i=first_det_of_comb, N_det_generators
|
||||
if(.not.(computed(i))) then
|
||||
first_det = i
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
do i=comb_teeth, 1, -1
|
||||
if(first_det_of_teeth(i) < first_det) then
|
||||
first_teeth = i
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, size_tbc ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Size of the tbc array
|
||||
END_DOC
|
||||
size_tbc = (comb_teeth+1)*N_det_generators + fragment_count*fragment_first
|
||||
END_PROVIDER
|
||||
|
||||
subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
|
||||
implicit none
|
||||
integer, intent(inout) :: Ncomb
|
||||
double precision, intent(out) :: comb(Ncomb)
|
||||
integer, intent(inout) :: tbc(0:size_tbc)
|
||||
logical, intent(inout) :: computed(N_det_generators)
|
||||
integer :: i, j, last_full, dets(comb_teeth)
|
||||
integer :: icount, n
|
||||
integer :: k, l
|
||||
l=first_det_of_comb
|
||||
call RANDOM_NUMBER(comb)
|
||||
do i=1,size(comb)
|
||||
comb(i) = comb(i) * comb_step
|
||||
!DIR$ FORCEINLINE
|
||||
call add_comb(comb(i), computed, tbc, size_tbc, comb_teeth)
|
||||
Ncomb = i
|
||||
if (tbc(0) == N_det_generators) return
|
||||
do while (computed(l))
|
||||
l=l+1
|
||||
enddo
|
||||
k=tbc(0)+1
|
||||
tbc(k) = l
|
||||
computed(l) = .True.
|
||||
tbc(0) = k
|
||||
enddo
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
|
||||
subroutine get_comb(stato, dets, ct)
|
||||
implicit none
|
||||
integer, intent(in) :: ct
|
||||
double precision, intent(in) :: stato
|
||||
integer, intent(out) :: dets(ct)
|
||||
double precision :: curs
|
||||
integer :: j
|
||||
integer, external :: mrcc_find
|
||||
|
||||
curs = 1d0 - stato
|
||||
do j = comb_teeth, 1, -1
|
||||
!DIR$ FORCEINLINE
|
||||
dets(j) = mrcc_find(curs, mrcc_cweight,size(mrcc_cweight), first_det_of_teeth(j), first_det_of_teeth(j+1))
|
||||
curs -= comb_step
|
||||
end do
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine add_comb(comb, computed, tbc, stbc, ct)
|
||||
implicit none
|
||||
integer, intent(in) :: stbc, ct
|
||||
double precision, intent(in) :: comb
|
||||
logical, intent(inout) :: computed(N_det_generators)
|
||||
integer, intent(inout) :: tbc(0:stbc)
|
||||
integer :: i, k, l, dets(ct)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call get_comb(comb, dets, ct)
|
||||
|
||||
k=tbc(0)+1
|
||||
do i = 1, ct
|
||||
l = dets(i)
|
||||
if(.not.(computed(l))) then
|
||||
tbc(k) = l
|
||||
k = k+1
|
||||
computed(l) = .true.
|
||||
end if
|
||||
end do
|
||||
tbc(0) = k-1
|
||||
end subroutine
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, mrcc_weight, (N_det_generators) ]
|
||||
&BEGIN_PROVIDER [ double precision, mrcc_cweight, (N_det_generators) ]
|
||||
&BEGIN_PROVIDER [ double precision, mrcc_cweight_cache, (N_det_generators) ]
|
||||
&BEGIN_PROVIDER [ double precision, comb_step ]
|
||||
&BEGIN_PROVIDER [ integer, first_det_of_teeth, (comb_teeth+1) ]
|
||||
&BEGIN_PROVIDER [ integer, first_det_of_comb ]
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: norm_left, stato
|
||||
integer, external :: mrcc_find
|
||||
|
||||
mrcc_weight(1) = psi_coef_generators(1,1)**2
|
||||
mrcc_cweight(1) = psi_coef_generators(1,1)**2
|
||||
|
||||
do i=1,N_det_generators
|
||||
mrcc_weight(i) = psi_coef_generators(i,1)**2
|
||||
enddo
|
||||
|
||||
! Important to loop backwards for numerical precision
|
||||
mrcc_cweight(N_det_generators) = mrcc_weight(N_det_generators)
|
||||
do i=N_det_generators-1,1,-1
|
||||
mrcc_cweight(i) = mrcc_weight(i) + mrcc_cweight(i+1)
|
||||
end do
|
||||
|
||||
do i=1,N_det_generators
|
||||
mrcc_weight(i) = mrcc_weight(i) / mrcc_cweight(1)
|
||||
mrcc_cweight(i) = mrcc_cweight(i) / mrcc_cweight(1)
|
||||
enddo
|
||||
|
||||
do i=1,N_det_generators-1
|
||||
mrcc_cweight(i) = 1.d0 - mrcc_cweight(i+1)
|
||||
end do
|
||||
mrcc_cweight(N_det_generators) = 1.d0
|
||||
|
||||
norm_left = 1d0
|
||||
|
||||
comb_step = 1d0/dfloat(comb_teeth)
|
||||
first_det_of_comb = 1
|
||||
do i=1,N_det_generators
|
||||
if(mrcc_weight(i)/norm_left < .5d0*comb_step) then
|
||||
first_det_of_comb = i
|
||||
exit
|
||||
end if
|
||||
norm_left -= mrcc_weight(i)
|
||||
end do
|
||||
first_det_of_comb = max(2,first_det_of_comb)
|
||||
call write_int(6, first_det_of_comb-1, 'Size of deterministic set')
|
||||
|
||||
comb_step = (1d0 - mrcc_cweight(first_det_of_comb-1)) * comb_step
|
||||
|
||||
stato = 1d0 - comb_step
|
||||
iloc = N_det_generators
|
||||
do i=comb_teeth, 1, -1
|
||||
integer :: iloc
|
||||
iloc = mrcc_find(stato, mrcc_cweight, N_det_generators, 1, iloc)
|
||||
first_det_of_teeth(i) = iloc
|
||||
stato -= comb_step
|
||||
end do
|
||||
first_det_of_teeth(comb_teeth+1) = N_det_generators + 1
|
||||
first_det_of_teeth(1) = first_det_of_comb
|
||||
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, mrcc_weight_inv, (N_det_generators) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Inverse of mrcc_weight array
|
||||
END_DOC
|
||||
integer :: i
|
||||
do i=1,N_det_generators
|
||||
mrcc_weight_inv(i) = 1.d0/mrcc_weight(i)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,215 @@
|
|||
|
||||
BEGIN_PROVIDER [ integer, fragment_count ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of fragments for the deterministic part
|
||||
END_DOC
|
||||
fragment_count = 1 ! (elec_alpha_num-n_core_orb)**2
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine run_mrcc_slave(thread,iproc,energy)
|
||||
use f77_zmq
|
||||
use selection_types
|
||||
implicit none
|
||||
|
||||
double precision, intent(in) :: energy(N_states_diag)
|
||||
integer, intent(in) :: thread, iproc
|
||||
integer :: rc, i
|
||||
|
||||
integer :: worker_id, task_id(1), ctask, ltask
|
||||
character*(512) :: task
|
||||
|
||||
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
|
||||
integer(ZMQ_PTR), external :: new_zmq_push_socket
|
||||
integer(ZMQ_PTR) :: zmq_socket_push
|
||||
|
||||
!type(selection_buffer) :: buf
|
||||
logical :: done
|
||||
|
||||
double precision,allocatable :: mrcc_detail(:,:)
|
||||
integer :: index
|
||||
integer :: Nindex
|
||||
|
||||
integer(bit_kind),allocatable :: abuf(:,:,:)
|
||||
integer(bit_kind) :: mask(N_int,2), omask(N_int,2)
|
||||
|
||||
double precision,allocatable :: delta_ij_loc(:,:,:)
|
||||
double precision,allocatable :: delta_ii_loc(:,:)
|
||||
double precision,allocatable :: delta_ij_s2_loc(:,:,:)
|
||||
double precision,allocatable :: delta_ii_s2_loc(:,:)
|
||||
integer :: h,p,n
|
||||
logical :: ok
|
||||
double precision :: contrib(N_states)
|
||||
|
||||
allocate(delta_ij_loc(N_states,N_det_non_ref,N_det_ref) &
|
||||
,delta_ii_loc(N_states,N_det_ref) &
|
||||
,delta_ij_s2_loc(N_states,N_det_non_ref,N_det_ref) &
|
||||
,delta_ii_s2_loc(N_states, N_det_ref))
|
||||
|
||||
|
||||
allocate(abuf(N_int, 2, N_det_non_ref))
|
||||
|
||||
allocate(mrcc_detail(N_states, N_det_generators))
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
zmq_socket_push = new_zmq_push_socket(thread)
|
||||
call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
|
||||
if(worker_id == -1) then
|
||||
print *, "WORKER -1"
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
call end_zmq_push_socket(zmq_socket_push,thread)
|
||||
return
|
||||
end if
|
||||
!buf%N = 0
|
||||
ctask = 1
|
||||
Nindex=1
|
||||
mrcc_detail = 0d0
|
||||
do
|
||||
call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task)
|
||||
|
||||
done = task_id(ctask) == 0
|
||||
if (done) then
|
||||
ctask = ctask - 1
|
||||
else
|
||||
integer :: i_generator, i_i_generator, subset
|
||||
read (task,*) subset, index
|
||||
|
||||
! if(buf%N == 0) then
|
||||
! ! Only first time
|
||||
! call create_selection_buffer(1, 2, buf)
|
||||
! end if
|
||||
do i_i_generator=1, Nindex
|
||||
contrib = 0d0
|
||||
i_generator = index
|
||||
!call select_connected(i_generator,energy,mrcc_detail(1, i_i_generator),buf,subset)
|
||||
|
||||
!!!!!!!!!!!!!!!!!!!!!!
|
||||
!if(mod(gen, 1000) == 0) print *, "mrcc ", gen, "/", N_det_generators
|
||||
do h=1, hh_shortcut(0)
|
||||
call apply_hole_local(psi_det_generators(1,1,i_generator), hh_exists(1, h), mask, ok, N_int)
|
||||
if(.not. ok) cycle
|
||||
omask = 0_bit_kind
|
||||
if(hh_exists(1, h) /= 0) omask = mask
|
||||
n = 1
|
||||
do p=hh_shortcut(h), hh_shortcut(h+1)-1
|
||||
call apply_particle_local(mask, pp_exists(1, p), abuf(1,1,n), ok, N_int)
|
||||
if(ok) n = n + 1
|
||||
if(n > N_det_non_ref) stop "Buffer too small in MRCC..."
|
||||
end do
|
||||
n = n - 1
|
||||
|
||||
if(n /= 0) then
|
||||
call mrcc_part_dress(delta_ij_loc, delta_ii_loc, delta_ij_s2_loc, delta_ii_s2_loc, &
|
||||
i_generator,n,abuf,N_int,omask,1d0,contrib)
|
||||
endif
|
||||
end do
|
||||
!!!!!!!!!!!!!!!!!!!!!!
|
||||
!print *, "contrib", i_generator, contrib
|
||||
mrcc_detail(:, i_i_generator) = contrib
|
||||
enddo
|
||||
endif
|
||||
|
||||
if(done .or. (ctask == size(task_id)) ) then
|
||||
! if(buf%N == 0 .and. ctask > 0) stop "uninitialized selection_buffer"
|
||||
do i=1, ctask
|
||||
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i))
|
||||
end do
|
||||
if(ctask > 0) then
|
||||
call push_mrcc_results(zmq_socket_push, Nindex, index, mrcc_detail, task_id(1), ctask)
|
||||
mrcc_detail(:,:Nindex) = 0d0
|
||||
! buf%cur = 0
|
||||
end if
|
||||
ctask = 0
|
||||
end if
|
||||
|
||||
if(done) exit
|
||||
ctask = ctask + 1
|
||||
end do
|
||||
call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
call end_zmq_push_socket(zmq_socket_push,thread)
|
||||
! call delete_selection_buffer(buf)
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine push_mrcc_results(zmq_socket_push, N, index, mrcc_detail, task_id, ntask)
|
||||
use f77_zmq
|
||||
use selection_types
|
||||
implicit none
|
||||
|
||||
integer(ZMQ_PTR), intent(in) :: zmq_socket_push
|
||||
double precision, intent(in) :: mrcc_detail(N_states, N_det_generators)
|
||||
integer, intent(in) :: ntask, N, index, task_id(*)
|
||||
integer :: rc
|
||||
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, N, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4) stop "push"
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, index, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4*N) stop "push"
|
||||
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, mrcc_detail, 8*N_states*N, ZMQ_SNDMORE)
|
||||
if(rc /= 8*N_states*N) stop "push"
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4) stop "push"
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, task_id, ntask*4, 0)
|
||||
if(rc /= 4*ntask) stop "push"
|
||||
|
||||
! Activate is zmq_socket_push is a REQ
|
||||
IRP_IF ZMQ_PUSH
|
||||
IRP_ELSE
|
||||
character*(2) :: ok
|
||||
rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0)
|
||||
IRP_ENDIF
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine pull_mrcc_results(zmq_socket_pull, N, index, mrcc_detail, task_id, ntask)
|
||||
use f77_zmq
|
||||
use selection_types
|
||||
implicit none
|
||||
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
|
||||
double precision, intent(inout) :: mrcc_detail(N_states, N_det_generators)
|
||||
integer, intent(out) :: index
|
||||
integer, intent(out) :: N, ntask, task_id(*)
|
||||
integer :: rc, rn, i
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, N, 4, 0)
|
||||
if(rc /= 4) stop "pull"
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, index, 4, 0)
|
||||
if(rc /= 4*N) stop "pull"
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, mrcc_detail, N_states*8*N, 0)
|
||||
if(rc /= 8*N_states*N) stop "pull"
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)
|
||||
if(rc /= 4) stop "pull"
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, task_id, ntask*4, 0)
|
||||
if(rc /= 4*ntask) stop "pull"
|
||||
|
||||
! Activate is zmq_socket_pull is a REP
|
||||
IRP_IF ZMQ_PUSH
|
||||
IRP_ELSE
|
||||
rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0)
|
||||
IRP_ENDIF
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, mrcc_workload, (N_det_generators) ]
|
||||
integer :: i
|
||||
do i=1,N_det_generators
|
||||
mrcc_workload(i) = dfloat(N_det_generators - i + 1)**2
|
||||
end do
|
||||
mrcc_workload = mrcc_workload / sum(mrcc_workload)
|
||||
END_PROVIDER
|
||||
|
||||
|
|
@ -0,0 +1,133 @@
|
|||
|
||||
subroutine create_selection_buffer(N, siz, res)
|
||||
use selection_types
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: N, siz
|
||||
type(selection_buffer), intent(out) :: res
|
||||
|
||||
allocate(res%det(N_int, 2, siz), res%val(siz))
|
||||
|
||||
res%val = 0d0
|
||||
res%det = 0_8
|
||||
res%N = N
|
||||
res%mini = 0d0
|
||||
res%cur = 0
|
||||
end subroutine
|
||||
|
||||
subroutine delete_selection_buffer(b)
|
||||
use selection_types
|
||||
implicit none
|
||||
type(selection_buffer), intent(inout) :: b
|
||||
if (associated(b%det)) then
|
||||
deallocate(b%det)
|
||||
endif
|
||||
if (associated(b%val)) then
|
||||
deallocate(b%val)
|
||||
endif
|
||||
end
|
||||
|
||||
|
||||
subroutine add_to_selection_buffer(b, det, val)
|
||||
use selection_types
|
||||
implicit none
|
||||
|
||||
type(selection_buffer), intent(inout) :: b
|
||||
integer(bit_kind), intent(in) :: det(N_int, 2)
|
||||
double precision, intent(in) :: val
|
||||
integer :: i
|
||||
|
||||
if(val <= b%mini) then
|
||||
b%cur += 1
|
||||
b%det(1:N_int,1:2,b%cur) = det(1:N_int,1:2)
|
||||
b%val(b%cur) = val
|
||||
if(b%cur == size(b%val)) then
|
||||
call sort_selection_buffer(b)
|
||||
end if
|
||||
end if
|
||||
end subroutine
|
||||
|
||||
subroutine merge_selection_buffers(b1, b2)
|
||||
use selection_types
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Merges the selection buffers b1 and b2 into b2
|
||||
END_DOC
|
||||
type(selection_buffer), intent(inout) :: b1
|
||||
type(selection_buffer), intent(inout) :: b2
|
||||
integer(bit_kind), pointer :: detmp(:,:,:)
|
||||
double precision, pointer :: val(:)
|
||||
integer :: i, i1, i2, k, nmwen
|
||||
if (b1%cur == 0) return
|
||||
do while (b1%val(b1%cur) > b2%mini)
|
||||
b1%cur = b1%cur-1
|
||||
if (b1%cur == 0) then
|
||||
return
|
||||
endif
|
||||
enddo
|
||||
nmwen = min(b1%N, b1%cur+b2%cur)
|
||||
allocate( val(size(b1%val)), detmp(N_int, 2, size(b1%det,3)) )
|
||||
i1=1
|
||||
i2=1
|
||||
do i=1,nmwen
|
||||
if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then
|
||||
exit
|
||||
else if (i1 > b1%cur) then
|
||||
val(i) = b2%val(i2)
|
||||
detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2)
|
||||
detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2)
|
||||
i2=i2+1
|
||||
else if (i2 > b2%cur) then
|
||||
val(i) = b1%val(i1)
|
||||
detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1)
|
||||
detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1)
|
||||
i1=i1+1
|
||||
else
|
||||
if (b1%val(i1) <= b2%val(i2)) then
|
||||
val(i) = b1%val(i1)
|
||||
detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1)
|
||||
detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1)
|
||||
i1=i1+1
|
||||
else
|
||||
val(i) = b2%val(i2)
|
||||
detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2)
|
||||
detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2)
|
||||
i2=i2+1
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
deallocate(b2%det, b2%val)
|
||||
b2%det => detmp
|
||||
b2%val => val
|
||||
b2%mini = min(b2%mini,b2%val(b2%N))
|
||||
b2%cur = nmwen
|
||||
end
|
||||
|
||||
|
||||
subroutine sort_selection_buffer(b)
|
||||
use selection_types
|
||||
implicit none
|
||||
|
||||
type(selection_buffer), intent(inout) :: b
|
||||
integer, allocatable :: iorder(:)
|
||||
integer(bit_kind), pointer :: detmp(:,:,:)
|
||||
integer :: i, nmwen
|
||||
logical, external :: detEq
|
||||
if (b%cur == 0) return
|
||||
nmwen = min(b%N, b%cur)
|
||||
|
||||
allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)))
|
||||
do i=1,b%cur
|
||||
iorder(i) = i
|
||||
end do
|
||||
call dsort(b%val, iorder, b%cur)
|
||||
do i=1, nmwen
|
||||
detmp(1:N_int,1,i) = b%det(1:N_int,1,iorder(i))
|
||||
detmp(1:N_int,2,i) = b%det(1:N_int,2,iorder(i))
|
||||
end do
|
||||
deallocate(b%det,iorder)
|
||||
b%det => detmp
|
||||
b%mini = min(b%mini,b%val(b%N))
|
||||
b%cur = nmwen
|
||||
end subroutine
|
||||
|
||||
Loading…
Reference in New Issue