PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-26 15:12:19 +02:00
Code Releases Activity

Removed error from type

Browse Source
This commit is contained in:
Anthony Scemama 2020-08-31 01:45:36 +02:00
parent 93fc49000c
commit e9f1d7576a
9 changed files with 125 additions and 156 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -6,7 +6,7 @@ subroutine run_cipsi
! stochastic PT2.
END_DOC
integer :: i,j,k
type(pt2_type) :: pt2_data
type(pt2_type) :: pt2_data, pt2_data_err
double precision, allocatable :: zeros(:)
integer :: to_select
logical, external :: qp_stop
@ -25,6 +25,7 @@ subroutine run_cipsi
allocate (zeros(N_states))
call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states)
double precision :: hf_energy_ref
logical :: has
@ -79,16 +80,19 @@ subroutine run_cipsi
to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
to_select = max(N_states_diag, to_select)
if (do_pt2) then
pt2_data % pt2 = 0.d0
pt2_data % variance = 0.d0
pt2_data % norm2 = 0.d0
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states)
threshold_generators_save = threshold_generators
threshold_generators = 1.d0
SOFT_TOUCH threshold_generators
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error, 0) ! Stochastic PT2
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,pt2_data_err,relative_error, 0) ! Stochastic PT2
threshold_generators = threshold_generators_save
SOFT_TOUCH threshold_generators
else
call pt2_dealloc(pt2_data)
call pt2_alloc(pt2_data, N_states)
call ZMQ_selection(to_select, pt2_data)
endif
@ -98,7 +102,7 @@ subroutine run_cipsi
call write_double(6,correlation_energy_ratio, 'Correlation ratio')
call print_summary(psi_energy_with_nucl_rep, &
pt2_data, N_det,N_occ_pattern,N_states,psi_s2)
pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2)
call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
@ -130,12 +134,13 @@ subroutine run_cipsi
endif
if (do_pt2) then
pt2_data % pt2(:) = 0.d0
pt2_data % variance(:) = 0.d0
pt2_data % norm2(:) = 0.d0
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states)
threshold_generators = 1d0
SOFT_TOUCH threshold_generators
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, relative_error, 0) ! Stochastic PT2
SOFT_TOUCH threshold_generators
endif
print *, 'N_det = ', N_det
@ -145,9 +150,11 @@ subroutine run_cipsi
call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
call print_summary(psi_energy_with_nucl_rep(1:N_states), &
pt2_data, N_det,N_occ_pattern,N_states,psi_s2)
pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call print_extrapolated_energy()
endif
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
end

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

@ -107,7 +107,7 @@ end function
subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in)
subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
use f77_zmq
use selection_types
@ -117,7 +117,7 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in)
integer, intent(in) :: N_in
! integer, intent(inout) :: N_in
double precision, intent(in) :: relative_error, E(N_states)
type(pt2_type), intent(inout) :: pt2_data
type(pt2_type), intent(inout) :: pt2_data, pt2_data_err
!
integer :: i, N
@ -298,13 +298,13 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in)
i = omp_get_thread_num()
if (i==0) then
call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, pt2_data, b, N)
call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, pt2_data, pt2_data_err, b, N)
pt2_data % rpt2(pt2_stoch_istate) = &
pt2_data % pt2(pt2_stoch_istate)/(1.d0 + pt2_data % norm2(pt2_stoch_istate))
!TODO : We should use here the correct formula for the error of X/Y
pt2_data % rpt2_err(pt2_stoch_istate) = &
pt2_data % pt2_err(pt2_stoch_istate)/(1.d0 + pt2_data % norm2(pt2_stoch_istate))
pt2_data_err % rpt2(pt2_stoch_istate) = &
pt2_data_err % pt2(pt2_stoch_istate)/(1.d0 + pt2_data % norm2(pt2_stoch_istate))
else
call pt2_slave_inproc(i)
@ -346,7 +346,7 @@ subroutine pt2_slave_inproc(i)
end
subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_err, b, N_)
use f77_zmq
use selection_types
use bitmasks
@ -355,7 +355,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
double precision, intent(in) :: relative_error, E
type(pt2_type), intent(inout) :: pt2_data
type(pt2_type), intent(inout) :: pt2_data, pt2_data_err
type(selection_buffer), intent(inout) :: b
integer, intent(in) :: N_
@ -414,9 +414,11 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
pt2_data % pt2(pt2_stoch_istate) = -huge(1.)
pt2_data % pt2_err(pt2_stoch_istate) = huge(1.)
pt2_data_err % pt2(pt2_stoch_istate) = huge(1.)
pt2_data % variance(pt2_stoch_istate) = huge(1.)
pt2_data_err % variance(pt2_stoch_istate) = huge(1.)
pt2_data % norm2(pt2_stoch_istate) = 0.d0
pt2_data_err % norm2(pt2_stoch_istate) = huge(1.)
n = 1
t = 0
U = 0
@ -479,8 +481,8 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
call pt2_add ( pt2_data_S(p), 1.d0, pt2_data_teeth )
call pt2_add2( pt2_data_S2(p), 1.d0, pt2_data_teeth )
enddo
call pt2_dealloc(pt2_data_teeth)
avg = E0 + pt2_data_S(t) % pt2(pt2_stoch_istate) / dble(c)
avg2 = v0 + pt2_data_S(t) % variance(pt2_stoch_istate) / dble(c)
avg3 = n0 + pt2_data_S(t) % norm2(pt2_stoch_istate) / dble(c)
@ -498,22 +500,22 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
if(c > 2) then
eqt = dabs((pt2_data_S2(t) % pt2(pt2_stoch_istate) / c) - (pt2_data_S(t) % pt2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability
eqt = sqrt(eqt / (dble(c) - 1.5d0))
pt2_data % pt2_err(pt2_stoch_istate) = eqt
pt2_data_err % pt2(pt2_stoch_istate) = eqt
eqt = dabs((pt2_data_S2(t) % variance(pt2_stoch_istate) / c) - (pt2_data_S(t) % variance(pt2_stoch_istate)/c)**2) ! dabs for numerical stability
eqt = sqrt(eqt / (dble(c) - 1.5d0))
pt2_data % variance_err(pt2_stoch_istate) = eqt
pt2_data_err % variance(pt2_stoch_istate) = eqt
eqt = dabs((pt2_data_S2(t) % norm2(pt2_stoch_istate) / c) - (pt2_data_S(t) % norm2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability
eqt = sqrt(eqt / (dble(c) - 1.5d0))
pt2_data % norm2_err(pt2_stoch_istate) = eqt
pt2_data_err % norm2(pt2_stoch_istate) = eqt
if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then
time1 = time
print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, F14.10, 2X, F14.10, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3, time-time0, ''
if (stop_now .or. ( &
(do_exit .and. (dabs(pt2_data % pt2_err(pt2_stoch_istate)) / &
(do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) / &
(1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then
if (zmq_abort(zmq_to_qp_run_socket) == -1) then
call sleep(10)

73
src/cipsi/pt2_type.irp.f
View File

@ -6,27 +6,17 @@ subroutine pt2_alloc(pt2_data,N)
integer :: k
allocate(pt2_data % pt2(N) &
,pt2_data % pt2_err(N) &
,pt2_data % variance(N) &
,pt2_data % variance_err(N) &
,pt2_data % norm2(N) &
,pt2_data % norm2_err(N) &
,pt2_data % rpt2(N) &
,pt2_data % rpt2_err(N) &
,pt2_data % overlap(N,N) &
,pt2_data % overlap_err(N,N) &
)
pt2_data % pt2(:) = 0.d0
pt2_data % pt2_err(:) = 0.d0
pt2_data % variance(:) = 0.d0
pt2_data % variance_err(:) = 0.d0
pt2_data % norm2(:) = 0.d0
pt2_data % norm2_err(:) = 0.d0
pt2_data % rpt2(:) = 0.d0
pt2_data % rpt2_err(:) = 0.d0
pt2_data % overlap(:,:) = 0.d0
pt2_data % overlap_err(:,:) = 0.d0
do k=1,N
pt2_data % overlap(k,k) = 1.d0
@ -38,15 +28,10 @@ subroutine pt2_dealloc(pt2_data)
use selection_types
type(pt2_type), intent(inout) :: pt2_data
deallocate(pt2_data % pt2 &
,pt2_data % pt2_err &
,pt2_data % variance &
,pt2_data % variance_err&
,pt2_data % norm2 &
,pt2_data % norm2_err &
,pt2_data % rpt2 &
,pt2_data % rpt2_err &
,pt2_data % overlap &
,pt2_data % overlap_err &
)
end subroutine
@ -63,28 +48,18 @@ subroutine pt2_add(p1, w, p2)
if (w == 1.d0) then
p1 % pt2(:) = p1 % pt2(:) + p2 % pt2(:)
p1 % pt2_err(:) = p1 % pt2_err(:) + p2 % pt2_err(:)
p1 % rpt2(:) = p1 % rpt2(:) + p2 % rpt2(:)
p1 % rpt2_err(:) = p1 % rpt2_err(:) + p2 % rpt2_err(:)
p1 % variance(:) = p1 % variance(:) + p2 % variance(:)
p1 % variance_err(:) = p1 % variance_err(:) + p2 % variance_err(:)
p1 % norm2(:) = p1 % norm2(:) + p2 % norm2(:)
p1 % norm2_err(:) = p1 % norm2_err(:) + p2 % norm2_err(:)
p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:)
p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + p2 % overlap_err(:,:)
else
p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:)
p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:)
p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:)
p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:)
p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:)
p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:)
p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:)
p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:)
p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:)
p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:)
endif
@ -104,28 +79,18 @@ subroutine pt2_add2(p1, w, p2)
if (w == 1.d0) then
p1 % pt2(:) = p1 % pt2(:) + p2 % pt2(:) * p2 % pt2(:)
p1 % pt2_err(:) = p1 % pt2_err(:) + p2 % pt2_err(:) * p2 % pt2_err(:)
p1 % rpt2(:) = p1 % rpt2(:) + p2 % rpt2(:) * p2 % rpt2(:)
p1 % rpt2_err(:) = p1 % rpt2_err(:) + p2 % rpt2_err(:) * p2 % rpt2_err(:)
p1 % variance(:) = p1 % variance(:) + p2 % variance(:) * p2 % variance(:)
p1 % variance_err(:) = p1 % variance_err(:) + p2 % variance_err(:) * p2 % variance_err(:)
p1 % norm2(:) = p1 % norm2(:) + p2 % norm2(:) * p2 % norm2(:)
p1 % norm2_err(:) = p1 % norm2_err(:) + p2 % norm2_err(:) * p2 % norm2_err(:)
p1 % overlap(:,:) = p1 % overlap(:,:) + p2 % overlap(:,:) * p2 % overlap(:,:)
p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + p2 % overlap_err(:,:) * p2 % overlap_err(:,:)
else
p1 % pt2(:) = p1 % pt2(:) + w * p2 % pt2(:) * p2 % pt2(:)
p1 % pt2_err(:) = p1 % pt2_err(:) + w * p2 % pt2_err(:) * p2 % pt2_err(:)
p1 % rpt2(:) = p1 % rpt2(:) + w * p2 % rpt2(:) * p2 % rpt2(:)
p1 % rpt2_err(:) = p1 % rpt2_err(:) + w * p2 % rpt2_err(:) * p2 % rpt2_err(:)
p1 % variance(:) = p1 % variance(:) + w * p2 % variance(:) * p2 % variance(:)
p1 % variance_err(:) = p1 % variance_err(:) + w * p2 % variance_err(:) * p2 % variance_err(:)
p1 % norm2(:) = p1 % norm2(:) + w * p2 % norm2(:) * p2 % norm2(:)
p1 % norm2_err(:) = p1 % norm2_err(:) + w * p2 % norm2_err(:) * p2 % norm2_err(:)
p1 % overlap(:,:) = p1 % overlap(:,:) + w * p2 % overlap(:,:) * p2 % overlap(:,:)
p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) * p2 % overlap_err(:,:)
endif
@ -142,18 +107,15 @@ subroutine pt2_serialize(pt2_data, n, x)
integer :: i,k,n2
n2 = n*n
x(1:n) = pt2_data % pt2(1:n)
x(n+1:2*n) = pt2_data % pt2_err(1:n)
x(2*n+1:3*n) = pt2_data % rpt2(1:n)
x(3*n+1:4*n) = pt2_data % rpt2_err(1:n)
x(4*n+1:5*n) = pt2_data % variance(1:n)
x(5*n+1:6*n) = pt2_data % variance_err(1:n)
x(6*n+1:7*n) = pt2_data % norm2(1:n)
x(7*n+1:8*n) = pt2_data % norm2_err(1:n)
k=8*n
x(k+1:k+n2) = reshape(pt2_data % overlap(1:n,1:n), (/ n2 /))
k=8*n+n2
x(k+1:k+n2) = reshape(pt2_data % overlap_err(1:n,1:n), (/ n2 /))
x(1:n) = pt2_data % pt2(1:n)
k=n
x(k+1:k+n) = pt2_data % rpt2(1:n)
k=k+n
x(k+1:k+n) = pt2_data % variance(1:n)
k=k+n
x(k+1:k+n) = pt2_data % norm2(1:n)
k=k+n
x(k+1:k+n2) = reshape(pt2_data % overlap(1:n,1:n), (/ n2 /))
end
@ -168,16 +130,13 @@ subroutine pt2_deserialize(pt2_data, n, x)
n2 = n*n
pt2_data % pt2(1:n) = x(1:n)
pt2_data % pt2_err(1:n) = x(n+1:2*n)
pt2_data % rpt2(1:n) = x(2*n+1:3*n)
pt2_data % rpt2_err(1:n) = x(3*n+1:4*n)
pt2_data % variance(1:n) = x(4*n+1:5*n)
pt2_data % variance_err(1:n) = x(5*n+1:6*n)
pt2_data % norm2(1:n) = x(6*n+1:7*n)
pt2_data % norm2_err(1:n) = x(7*n+1:8*n)
k=8*n
k=n
pt2_data % rpt2(1:n) = x(k+1:k+n)
k=k+n
pt2_data % variance(1:n) = x(k+1:k+n)
k=k+n
pt2_data % norm2(1:n) = x(k+1:k+n)
k=k+n
pt2_data % overlap(1:n,1:n) = reshape(x(k+1:k+n2), (/ n, n /))
k=8*n+n2
pt2_data % overlap_err(1:n,1:n) = reshape(x(k+1:k+n2), (/ n, n /))
end

23
src/cipsi/run_pt2_slave.irp.f
View File

@ -117,11 +117,11 @@ subroutine run_pt2_slave_small(thread,iproc,energy)
double precision :: time0, time1
call wall_time(time0)
do k=1,n_tasks
call pt2_alloc(pt2_data(k),N_states)
b%cur = 0
call pt2_alloc(pt2_data(k),N_states)
b%cur = 0
!double precision :: time2
!call wall_time(time2)
call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k)))
call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k)))
!call wall_time(time1)
!print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1))
enddo
@ -157,6 +157,7 @@ subroutine run_pt2_slave_small(thread,iproc,energy)
if (buffer_ready) then
call delete_selection_buffer(b)
endif
deallocate(pt2_data)
end subroutine
@ -171,7 +172,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
integer :: worker_id, ctask, ltask
character*(512) :: task
integer :: task_id
integer :: task_id(1)
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
@ -182,9 +183,9 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
type(selection_buffer) :: b
logical :: done, buffer_ready
type(pt2_type) :: pt2_data
type(pt2_type) :: pt2_data(1)
integer :: n_tasks, k, N
integer :: i_generator, subset
integer :: i_generator(1), subset
integer :: bsize ! Size of selection buffers
logical :: sending
@ -213,13 +214,13 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
if (get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task, n_tasks) == -1) then
exit
endif
done = task_id == 0
done = task_id(1) == 0
if (done) then
n_tasks = n_tasks-1
endif
if (n_tasks == 0) exit
read (task,*) subset, i_generator, N
read (task,*) subset, i_generator(1), N
if (b%N == 0) then
! Only first time
bsize = min(N, (elec_alpha_num * (mo_num-elec_alpha_num))**2)
@ -231,11 +232,11 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
double precision :: time0, time1
call wall_time(time0)
call pt2_alloc(pt2_data,N_states)
call pt2_alloc(pt2_data(1),N_states)
b%cur = 0
!double precision :: time2
!call wall_time(time2)
call select_connected(i_generator,energy,pt2_data,b,subset,pt2_F(i_generator))
call select_connected(i_generator(1),energy,pt2_data(1),b,subset,pt2_F(i_generator(1)))
!call wall_time(time1)
!print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1))
call wall_time(time1)
@ -261,7 +262,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
call push_pt2_results_async_send(zmq_socket_push, i_generator, pt2_data, b, task_id, n_tasks,sending)
endif
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data(1))
end do
call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending)

52
src/cipsi/run_selection_slave.irp.f
View File

@ -101,11 +101,11 @@ subroutine run_selection_slave(thread,iproc,energy)
call sort_selection_buffer(buf)
! call merge_selection_buffers(buf,buf2)
call push_selection_results(zmq_socket_push, pt2_data, buf, task_id(1), ctask)
call pt2_dealloc(pt2_data)
! buf%mini = buf2%mini
buf%cur = 0
end if
ctask = 0
call pt2_dealloc(pt2_data)
integer, external :: disconnect_from_taskserver
if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then
@ -121,7 +121,7 @@ subroutine run_selection_slave(thread,iproc,energy)
end subroutine
subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask)
subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntasks)
use f77_zmq
use selection_types
implicit none
@ -129,9 +129,9 @@ subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask)
integer(ZMQ_PTR), intent(in) :: zmq_socket_push
type(pt2_type), intent(in) :: pt2_data
type(selection_buffer), intent(inout) :: b
integer, intent(in) :: ntask, task_id(*)
integer, intent(in) :: ntasks, task_id(*)
integer :: rc
double precision, allocatable :: pt2_serialized(:,:)
double precision, allocatable :: pt2_serialized(:)
rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE)
if(rc /= 4) then
@ -139,13 +139,10 @@ subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask)
endif
allocate(pt2_serialized (pt2_type_size(N_states),n_tasks) )
do i=1,n_tasks
call pt2_serialize(pt2_data(i),N_states,pt2_serialized(1,i))
enddo
allocate(pt2_serialized (pt2_type_size(N_states)) )
call pt2_serialize(pt2_data,N_states,pt2_serialized)
rc = f77_zmq_send( zmq_socket_push, pt2_serialized, size(pt2_serialized)*8, ZMQ_SNDMORE)
deallocate(pt2_serialized)
if (rc == -1) then
print *, irp_here, ': error sending result'
stop 3
@ -153,6 +150,7 @@ subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask)
else if(rc /= size(pt2_serialized)*8) then
stop 'push'
endif
deallocate(pt2_serialized)
if (b%cur > 0) then
@ -168,14 +166,14 @@ subroutine push_selection_results(zmq_socket_push, pt2_data, b, task_id, ntask)
endif
rc = f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE)
rc = f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE)
if(rc /= 4) then
print *, 'f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE)'
print *, 'f77_zmq_send( zmq_socket_push, ntasks, 4, ZMQ_SNDMORE)'
endif
rc = f77_zmq_send( zmq_socket_push, task_id(1), ntask*4, 0)
if(rc /= 4*ntask) then
print *, 'f77_zmq_send( zmq_socket_push, task_id(1), ntask*4, 0)'
rc = f77_zmq_send( zmq_socket_push, task_id(1), ntasks*4, 0)
if(rc /= 4*ntasks) then
print *, 'f77_zmq_send( zmq_socket_push, task_id(1), ntasks*4, 0)'
endif
! Activate is zmq_socket_push is a REQ
@ -192,7 +190,7 @@ IRP_ENDIF
end subroutine
subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_id, ntask)
subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_id, ntasks)
use f77_zmq
use selection_types
implicit none
@ -200,27 +198,25 @@ subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_i
type(pt2_type), intent(inout) :: pt2_data
double precision, intent(out) :: val(*)
integer(bit_kind), intent(out) :: det(N_int, 2, *)
integer, intent(out) :: N, ntask, task_id(*)
integer, intent(out) :: N, ntasks, task_id(*)
integer :: rc, rn, i
double precision, allocatable :: pt2_serialized(:,:)
double precision, allocatable :: pt2_serialized(:)
rc = f77_zmq_recv( zmq_socket_pull, N, 4, 0)
if(rc /= 4) then
print *, 'f77_zmq_recv( zmq_socket_pull, N, 4, 0)'
endif
allocate(pt2_serialized (pt2_type_size(N_states),n_tasks) )
rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized)*n_tasks, 0)
allocate(pt2_serialized (pt2_type_size(N_states)) )
rc = f77_zmq_recv( zmq_socket_pull, pt2_serialized, 8*size(pt2_serialized)*ntasks, 0)
if (rc == -1) then
n_tasks = 1
ntasks = 1
task_id(1) = 0
else if(rc /= 8*size(pt2_serialized)) then
stop 'pull'
endif
do i=1,n_tasks
call pt2_deserialize(pt2_data(i),N_states,pt2_serialized(1,i))
enddo
call pt2_deserialize(pt2_data,N_states,pt2_serialized)
deallocate(pt2_serialized)
if (N>0) then
@ -235,14 +231,14 @@ subroutine pull_selection_results(zmq_socket_pull, pt2_data, val, det, N, task_i
endif
endif
rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)
rc = f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0)
if(rc /= 4) then
print *, 'f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)'
print *, 'f77_zmq_recv( zmq_socket_pull, ntasks, 4, 0)'
endif
rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0)
if(rc /= 4*ntask) then
print *, 'f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0)'
rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntasks*4, 0)
if(rc /= 4*ntasks) then
print *, 'f77_zmq_recv( zmq_socket_pull, task_id(1), ntasks*4, 0)'
endif
! Activate is zmq_socket_pull is a REP

7
src/cipsi/selection_types.f90
View File

@ -8,15 +8,10 @@ module selection_types
type pt2_type
double precision, allocatable :: pt2(:)
double precision, allocatable :: pt2_err(:)
double precision, allocatable :: rpt2(:)
double precision, allocatable :: rpt2_err(:)
double precision, allocatable :: variance(:)
double precision, allocatable :: variance_err(:)
double precision, allocatable :: norm2(:)
double precision, allocatable :: norm2_err(:)
double precision, allocatable :: overlap(:,:)
double precision, allocatable :: overlap_err(:,:)
endtype
contains
@ -24,7 +19,7 @@ module selection_types
integer function pt2_type_size(N)
implicit none
integer, intent(in) :: N
pt2_type_size = (8*n + 2*n*n)
pt2_type_size = (4*n + n*n)
end function
end module

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

@ -7,7 +7,7 @@ subroutine run_stochastic_cipsi
integer :: i,j,k
double precision, allocatable :: zeros(:)
integer :: to_select
type(pt2_type) :: pt2_data
type(pt2_type) :: pt2_data, pt2_data_err
logical, external :: qp_stop
@ -24,6 +24,7 @@ subroutine run_stochastic_cipsi
allocate (zeros(N_states))
call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states)
double precision :: hf_energy_ref
logical :: has
@ -79,10 +80,11 @@ subroutine run_stochastic_cipsi
to_select = max(N_states_diag, to_select)
pt2_data % pt2 = 0.d0
pt2_data % variance = 0.d0
pt2_data % norm2 = 0.d0
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,to_select) ! Stochastic PT2 and selection
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states)
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,pt2_data_err,relative_error,to_select) ! Stochastic PT2 and selection
correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / &
(psi_energy_with_nucl_rep(1) + pt2_data % rpt2(1) - hf_energy_ref)
@ -90,7 +92,7 @@ subroutine run_stochastic_cipsi
call write_double(6,correlation_energy_ratio, 'Correlation ratio')
call print_summary(psi_energy_with_nucl_rep, &
pt2_data, N_det,N_occ_pattern,N_states,psi_s2)
pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2)
call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
@ -121,17 +123,19 @@ subroutine run_stochastic_cipsi
call save_energy(psi_energy_with_nucl_rep, zeros)
endif
pt2_data % pt2(:) = 0.d0
pt2_data % variance(:) = 0.d0
pt2_data % norm2(:) = 0.d0
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states)
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, relative_error, 0) ! Stochastic PT2
call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
call print_summary(psi_energy_with_nucl_rep, &
pt2_data , N_det, N_occ_pattern, N_states, psi_s2)
pt2_data , pt2_data_err, N_det, N_occ_pattern, N_states, psi_s2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call print_extrapolated_energy()
endif
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
end

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

@ -32,29 +32,34 @@ subroutine run
integer :: i,j,k
logical, external :: detEq
type(pt2_type) :: pt2_data
type(pt2_type) :: pt2_data, pt2_data_err
integer :: degree
integer :: n_det_before, to_select
double precision :: threshold_davidson_in
double precision :: E_CI_before(N_states), relative_error
double precision :: relative_error
double precision, allocatable :: E_CI_before(:)
allocate ( E_CI_before(N_states))
call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states)
E_CI_before(:) = psi_energy(:) + nuclear_repulsion
relative_error=PT2_relative_error
if (do_pt2) then
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, 0) ! Stochastic PT2
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, relative_error, 0) ! Stochastic PT2
else
call ZMQ_selection(0, pt2_data)
endif
call print_summary(psi_energy_with_nucl_rep(1:N_states), &
pt2_data, N_det,N_occ_pattern,N_states,psi_s2)
pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2)
call save_energy(E_CI_before,pt2_data % pt2)
call save_energy(E_CI_before, pt2_data % pt2)
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
deallocate(E_CI_before)
end

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

@ -1,4 +1,4 @@
subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,n_st,s2_)
subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_occ_pattern_,n_st,s2_)
use selection_types
implicit none
BEGIN_DOC
@ -7,7 +7,7 @@ subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,n_st,s2_)
integer, intent(in) :: n_det_, n_occ_pattern_, n_st
double precision, intent(in) :: e_(n_st), s2_(n_st)
type(pt2_type) , intent(in) :: pt2_data
type(pt2_type) , intent(in) :: pt2_data, pt2_data_err
integer :: i, k
integer :: N_states_p
character*(9) :: pt2_string
@ -44,16 +44,16 @@ subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,n_st,s2_)
write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0
endif
write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))'
write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data % pt2_err(k), k=1,N_states_p)
write(*,fmt) '# rPT2'//pt2_string, (pt2_data % pt2(k)*f(k), pt2_data % pt2_err(k)*f(k), k=1,N_states_p)
write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p)
write(*,fmt) '# rPT2'//pt2_string, (pt2_data % pt2(k)*f(k), pt2_data_err % pt2(k)*f(k), k=1,N_states_p)
write(*,'(A)') '#'
write(*,fmt) '# E+PT2 ', (e_(k)+pt2_data % pt2(k),pt2_data % pt2_err(k), k=1,N_states_p)
write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % pt2(k)*f(k),pt2_data % pt2_err(k)*f(k), k=1,N_states_p)
write(*,fmt) '# E+PT2 ', (e_(k)+pt2_data % pt2(k),pt2_data_err % pt2(k), k=1,N_states_p)
write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % pt2(k)*f(k),pt2_data_err % pt2(k)*f(k), k=1,N_states_p)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), &
dsqrt(pt2_data % pt2_err(k)*pt2_data % pt2_err(k)+pt2_data % pt2_err(1)*pt2_data % pt2_err(1)), k=1,N_states_p)
dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p)
write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, &
dsqrt(pt2_data % pt2_err(k)*pt2_data % pt2_err(k)+pt2_data % pt2_err(1)*pt2_data % pt2_err(1))*27.211396641308d0, k=1,N_states_p)
dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p)
endif
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
@ -70,12 +70,12 @@ subroutine print_summary(e_,pt2_data,n_det_,n_occ_pattern_,n_st,s2_)
print*,'* State ',k
print *, '< S^2 > = ', s2_(k)
print *, 'E = ', e_(k)
print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data % variance_err(k)
print *, 'PT norm = ', dsqrt(pt2_data % norm2(k)), ' +/- ', 0.5d0*dsqrt(pt2_data % norm2(k)) * pt2_data % norm2_err(k) / pt2_data % norm2(k)
print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data % pt2_err(k)
print *, 'rPT2 = ', pt2_data % pt2(k)*f(k), ' +/- ', pt2_data % rpt2_err(k)
print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data % pt2_err(k)
print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % pt2(k)*f(k), ' +/- ', pt2_data % pt2_err(k)*f(k)
print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k)
print *, 'PT norm = ', dsqrt(pt2_data % norm2(k)), ' +/- ', 0.5d0*dsqrt(pt2_data % norm2(k)) * pt2_data_err % norm2(k) / pt2_data % norm2(k)
print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)
print *, 'rPT2 = ', pt2_data % pt2(k)*f(k), ' +/- ', pt2_data_err % rpt2(k)
print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)
print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % pt2(k)*f(k), ' +/- ', pt2_data_err % pt2(k)*f(k)
print *, ''
enddo