LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-04 10:25:58 +02:00
Code Issues Releases Wiki Activity
ad69f39f99
quantum_package/plugins/dress_zmq/dress_stoch_routines.irp.f
Yann Garniron ad69f39f99 dress_zmq re-implemented
2018-08-29 20:54:58 +02:00

526 lines
15 KiB
Fortran
Raw Blame History

BEGIN_PROVIDER [ integer, dress_stoch_istate ]
implicit none
dress_stoch_istate = 1
END_PROVIDER
BEGIN_PROVIDER [ integer, pt2_N_teeth ]
&BEGIN_PROVIDER [ integer, pt2_minDetInFirstTeeth ]
&BEGIN_PROVIDER [ integer, pt2_n_tasks_max ]
&BEGIN_PROVIDER [ integer, pt2_F, (N_det_generators) ]
implicit none
pt2_F(:) = 1
pt2_F(:N_det_generators/100 + 1) = 1
pt2_n_tasks_max = N_det_generators/100 + 1
if(N_det_generators < 256) then
pt2_minDetInFirstTeeth = 1
pt2_N_teeth = 1
else
pt2_minDetInFirstTeeth = 5
pt2_N_teeth = 16
end if
END_PROVIDER
BEGIN_PROVIDER [ integer, dress_N_cp_max ]
dress_N_cp_max = 100
END_PROVIDER
BEGIN_PROVIDER[ integer, dress_M_m, (dress_N_cp_max)]
implicit none
integer :: i
do i=1,dress_N_cp_max-1
dress_M_m(i) = N_det_generators * i / (dress_N_cp_max+1)
end do
dress_M_m(1) = 1
dress_M_m(dress_N_cp_max) = N_det_generators+1
END_PROVIDER
BEGIN_PROVIDER[ integer, pt2_J, (N_det_generators)]
&BEGIN_PROVIDER[ double precision, pt2_u, (N_det_generators)]
&BEGIN_PROVIDER[ integer, dress_R, (0:N_det_generators)]
&BEGIN_PROVIDER[ integer, dress_R1, (0:N_det_generators)]
&BEGIN_PROVIDER[ double precision, dress_M_mi, (dress_N_cp_max, N_det_generators+1)]
&BEGIN_PROVIDER [ integer, dress_P, (N_det_generators) ]
&BEGIN_PROVIDER [ integer, dress_T, (N_det_generators) ]
&BEGIN_PROVIDER [ integer, dress_N_cp ]
implicit none
integer :: N_c, N_j, U, t, i, m
double precision :: v
double precision, allocatable :: tilde_M(:)
logical, allocatable :: d(:)
integer, external :: dress_find_sample
allocate(d(N_det_generators), tilde_M(N_det_generators))
dress_M_mi = 0d0
tilde_M = 0d0
dress_R(:) = 0
dress_R1(:) = 0
N_c = 0
N_j = pt2_n_0(1)
d(:) = .false.
do i=1,N_j
d(i) = .true.
pt2_J(i) = i
end do
call random_seed(put=(/3211,64,6566,321,65,321,654,65,321,6321,654,65,321,621,654,65,321,65,654,65,321,65/))
call RANDOM_NUMBER(pt2_u)
call RANDOM_NUMBER(pt2_u)
U = 0
m = 1
do while(N_j < N_det_generators)
!ADD_COMB
N_c += 1
do t=0, pt2_N_teeth-1
v = pt2_u_0 + pt2_W_T * (dble(t) + pt2_u(N_c))
i = dress_find_sample(v, pt2_cW)
tilde_M(i) += 1d0
if(.not. d(i)) then
N_j += 1
pt2_J(N_j) = i
d(i) = .true.
end if
end do
!FILL_TOOTH
do while(U < N_det_generators)
U += 1
if(.not. d(U)) then
N_j += 1
pt2_J(N_j) = U
d(U) = .true.
exit;
end if
end do
if(N_c == dress_M_m(m)) then
dress_R1(m) = N_j
dress_R(N_j) = N_c
dress_M_mi(m, :N_det_generators) = tilde_M(:)
m += 1
end if
enddo
dress_N_cp = m-1
dress_R1(dress_N_cp) = N_j
!!!!!!!!!!!!!!
do m=1,dress_N_cp
do i=dress_R1(m-1)+1, dress_R1(m)
dress_P(pt2_J(i)) = m
end do
end do
do i=1, pt2_n_0(1)
dress_T(i) = 0
end do
do t=2,pt2_N_teeth+1
do i=pt2_n_0(t-1)+1, pt2_n_0(t)
dress_T(i) = t-1
end do
end do
!!!!!!!!!!!!!
END_PROVIDER
! BEGIN_PROVIDER [double precision, dress_e, (N_det_generators, dress_N_cp)]
!&BEGIN_PROVIDER [integer, dress_dot_t, (0:dress_N_cp)]
!&BEGIN_PROVIDER [integer, dress_dot_n_0, (0:dress_N_cp)]
! implicit none
! dress_e(:,:) = 1d0
! dress_dot_t(:) = 0
! dress_dot_n_0(:) = 0
!
! integer :: U, m, t, i
!
! U = pt2_n_0(1)+1
!
! do m=1,dress_N_cp
! do while(dress_M_mi(m, U) /= 0d0)
! U = U+1
! end do
! dress_dot_t(m) = pt2_N_teeth + 1
! dress_dot_n_0(m) = N_det_generators
!!
! do t = 2, pt2_N_teeth+1
! if(U <= pt2_n_0(t)) then
! dress_dot_t(m) = t-1
! dress_dot_n_0(m) = pt2_n_0(t-1)
! exit
! end if
! end do
! do t=dress_dot_t(m), pt2_N_teeth
! do i=pt2_n_0(t)+1, pt2_n_0(t+1)
! dress_e(i,m) = pt2_W_T * dress_M_mi(m,i) / pt2_w(i)
! end do
! end do
! end do
! do m=dress_N_cp, 2, -1
! dress_e(:,m) -= dress_e(:,m-1)
! end do
!END_PROVIDER
subroutine ZMQ_dress(E, dress, delta_out, delta_s2_out, relative_error)
use f77_zmq
implicit none
character(len=64000) :: task
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
integer, external :: omp_get_thread_num
double precision, intent(in) :: E(N_states), relative_error
double precision, intent(out) :: dress(N_states)
double precision, intent(out) :: delta_out(N_states, N_det)
double precision, intent(out) :: delta_s2_out(N_states, N_det)
double precision, allocatable :: delta(:,:)
double precision, allocatable :: delta_s2(:,:)
integer :: i, j, k, Ncp
integer, external :: add_task_to_taskserver
double precision :: state_average_weight_save(N_states)
task(:) = CHAR(0)
allocate(delta(N_states,N_det), delta_s2(N_states, N_det))
state_average_weight_save(:) = state_average_weight(:)
do dress_stoch_istate=1,N_states
SOFT_TOUCH dress_stoch_istate
state_average_weight(:) = 0.d0
state_average_weight(dress_stoch_istate) = 1.d0
TOUCH state_average_weight
!provide psi_coef_generators
provide nproc mo_bielec_integrals_in_map mo_mono_elec_integral psi_selectors
!print *, dress_e0_denominator
print *, '========== ================= ================= ================='
print *, ' Samples Energy Stat. Error Seconds '
print *, '========== ================= ================= ================='
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull, 'dress')
integer, external :: zmq_put_psi
integer, external :: zmq_put_N_det_generators
integer, external :: zmq_put_N_det_selectors
integer, external :: zmq_put_dvector
integer, external :: zmq_set_running
if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then
stop 'Unable to put psi on ZMQ server'
endif
if (zmq_put_N_det_generators(zmq_to_qp_run_socket, 1) == -1) then
stop 'Unable to put N_det_generators on ZMQ server'
endif
if (zmq_put_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) then
stop 'Unable to put N_det_selectors on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',dress_e0_denominator,size(dress_e0_denominator)) == -1) then
stop 'Unable to put energy on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,"state_average_weight",state_average_weight,N_states) == -1) then
stop 'Unable to put state_average_weight on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,"dress_stoch_istate",real(dress_stoch_istate,8),1) == -1) then
stop 'Unable to put dress_stoch_istate on ZMQ server'
endif
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
do i=1,N_det_generators
do j=1,pt2_F(i) !!!!!!!!!!!!
write(task(1:20),'(I9,1X,I9''|'')') j, pt2_J(i)
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:20))) == -1) then
stop 'Unable to add task to task server'
endif
end do
end do
if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
print *, irp_here, ': Failed in zmq_set_running'
endif
!!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc) &
! !$OMP PRIVATE(i)
!i = omp_get_thread_num()
!if (i==0) then
call dress_collector(zmq_socket_pull,E, relative_error, delta, delta_s2, dress,&
dress_stoch_istate)
!else
! call dress_slave_inproc(i)
!endif
!!$OMP END PARALLEL
delta_out(dress_stoch_istate,1:N_det) = delta(dress_stoch_istate,1:N_det)
delta_s2_out(dress_stoch_istate,1:N_det) = delta_s2(dress_stoch_istate,1:N_det)
call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'dress')
print *, '========== ================= ================= ================='
enddo
FREE dress_stoch_istate
state_average_weight(:) = state_average_weight_save(:)
TOUCH state_average_weight
deallocate(delta,delta_s2)
end subroutine
subroutine dress_slave_inproc(i)
implicit none
integer, intent(in) :: i
call run_dress_slave(1,i,dress_e0_denominator)
end
BEGIN_PROVIDER [double precision, dress_e, (N_det_generators, dress_N_cp)]
&BEGIN_PROVIDER [integer, dress_dot_t, (0:dress_N_cp)]
&BEGIN_PROVIDER [integer, dress_dot_n_0, (0:dress_N_cp)]
&BEGIN_PROVIDER [integer, dress_dot_F, (dress_N_cp)]
implicit none
logical, allocatable :: d(:)
integer :: U, m, t, i
allocate(d(N_det_generators+1))
dress_e(:,:) = 1d0
dress_dot_t(:) = 0
dress_dot_n_0(:) = 0
dress_dot_F = 0
d(:) = .false.
U=0
do m=1,dress_N_cp
do i=dress_R1(m-1)+1,dress_R1(m)
dress_dot_F(m) += pt2_F(pt2_J(i))
d(pt2_J(i)) = .true.
end do
do while(d(U+1))
U += 1
end do
dress_dot_t(m) = pt2_N_teeth + 1
dress_dot_n_0(m) = N_det_generators
do t = 2, pt2_N_teeth+1
if(U < pt2_n_0(t)) then
dress_dot_t(m) = t-1
dress_dot_n_0(m) = pt2_n_0(t-1)
exit
end if
end do
do t=dress_dot_t(m), pt2_N_teeth
do i=pt2_n_0(t)+1, pt2_n_0(t+1)
dress_e(i,m) = pt2_W_T * dress_M_mi(m,i) / pt2_w(i)
end do
end do
end do
do m=dress_N_cp, 2, -1
dress_e(:,m) -= dress_e(:,m-1)
end do
END_PROVIDER
subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, dress, istate)
use f77_zmq
use bitmasks
implicit none
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
integer, intent(in) :: istate
double precision, intent(in) :: relative_error, E(N_states)
double precision, intent(out) :: dress(N_states)
double precision, allocatable :: cp(:,:,:,:)
double precision, intent(out) :: delta(N_states, N_det)
double precision, intent(out) :: delta_s2(N_states, N_det)
double precision, allocatable :: breve_delta_m(:,:,:), S(:), S2(:)
double precision, allocatable :: edI(:,:), edI_task(:,:)
integer, allocatable :: edI_index(:)
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 :: more
integer :: i, c, j, k, f, t, m, p, m_task
integer :: task_id, n_tasks
double precision :: E0, error, x, v, time, time0
double precision :: avg, eqt
double precision, external :: omp_get_wtime
integer, allocatable :: dot_f(:)
integer, external :: zmq_delete_tasks, dress_find_sample
delta = 0d0
delta_s2 = 0d0
allocate(cp(N_states, N_det, dress_N_cp, 2), edI(N_states, N_det))
allocate(edI_task(N_states, N_det), edI_index(N_det))
allocate(breve_delta_m(N_states, N_det, 2))
allocate(dot_f(dress_N_cp))
allocate(S(pt2_N_teeth+1), S2(pt2_N_teeth+1))
edI = -100000d0
cp = 0d0
dot_f(:) = dress_dot_F(:)
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
more = 1
m = 1
c = 0
S(:) = 0d0
S2(:) = 0d0
time0 = omp_get_wtime()
do while (m <= dress_N_cp)
if(dot_f(m) == 0) then
E0 = 0
do i=dress_dot_n_0(m),1,-1
E0 += edI(dress_stoch_istate, i)
end do
do while(c < dress_M_m(m))
c = c+1
x = 0d0
do p=pt2_N_teeth, 1, -1
v = pt2_u_0 + pt2_W_T * (pt2_u(c) + dble(p-1))
i = dress_find_sample(v, pt2_cW)
x += edI(dress_stoch_istate, i) * pt2_W_T / pt2_w(i)
S(p) += x
S2(p) += x**2
end do
end do
t = dress_dot_t(m)
avg = S(t) / dble(c)
eqt = (S2(t) / c) - (S(t)/c)**2
eqt = sqrt(eqt / dble(c-1))
error = eqt
time = omp_get_wtime()
print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', c, avg+E0+E(dress_stoch_istate), eqt, time-time0, ''
!end do
m += 1
else
task_id = 0
do
call pull_dress_results(zmq_socket_pull, m_task, f, edI_task, edI_index, breve_delta_m, task_id, n_tasks)
if(task_id == 0) exit
if (zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,n_tasks,more) == -1) then
stop 'Unable to delete tasks'
endif
end do
do i=1,n_tasks
edI(:, edI_index(i)) = edI_task(:, i) !!!!!!!!!!!!!!! += !!!!!
end do
cp(:,:,m_task,1) += breve_delta_m(:,:,1)
cp(:,:,m_task,2) += breve_delta_m(:,:,2)
dot_f(m_task) -= f
end if
end do
delta(:,:) = cp(:,:,m-1,1)
delta_s2(:,:) = cp(:,:,m-1,2)
dress(istate) = E(istate)+E0+avg
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
end subroutine
integer function dress_find_sample(v, w)
implicit none
double precision, intent(in) :: v, w(0:N_det_generators)
integer :: i,l,h
integer, parameter :: block=64
l = 0
h = N_det_generators
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 dress_find_sample=l,h
if(w(dress_find_sample) >= v) then
exit
end if
end do
end function
BEGIN_PROVIDER [ double precision, pt2_w, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, pt2_cW, (0:N_det_generators) ]
&BEGIN_PROVIDER [ double precision, pt2_W_T ]
&BEGIN_PROVIDER [ double precision, pt2_u_0 ]
&BEGIN_PROVIDER [ integer, pt2_n_0, (pt2_N_teeth+1) ]
implicit none
integer :: i, t
double precision, allocatable :: tilde_w(:), tilde_cW(:)
double precision :: r, tooth_width
integer, external :: dress_find_sample
allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators))
tilde_cW(0) = 0d0
do i=1,N_det_generators
tilde_w(i) = psi_coef_generators(i,dress_stoch_istate)**2
tilde_cW(i) = tilde_cW(i-1) + tilde_w(i)
enddo
pt2_n_0(1) = 0
do
pt2_u_0 = tilde_cW(pt2_n_0(1))
r = tilde_cW(pt2_n_0(1) + pt2_minDetInFirstTeeth)
pt2_W_T = (1d0 - pt2_u_0) / dble(pt2_N_teeth)
if(pt2_W_T >= r - pt2_u_0) then
exit
end if
pt2_n_0(1) += 1
if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then
stop "teeth building failed"
end if
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do t=2, pt2_N_teeth
r = pt2_u_0 + pt2_W_T * dble(t-1)
pt2_n_0(t) = dress_find_sample(r, tilde_cW)
end do
pt2_n_0(pt2_N_teeth+1) = N_det_generators
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
pt2_w(:pt2_n_0(1)) = tilde_w(:pt2_n_0(1))
do t=1, pt2_N_teeth
tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t))
do i=pt2_n_0(t)+1, pt2_n_0(t+1)
pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width
end do
end do
pt2_cW(0) = 0d0
do i=1,N_det_generators
pt2_cW(i) = pt2_cW(i-1) + pt2_w(i)
end do
pt2_n_0(pt2_N_teeth+1) = N_det_generators
END_PROVIDER
Reference in New Issue View Git Blame Copy Permalink