diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 44b87394..0f140240 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -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 diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 4b781dd8..5b05a5a9 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -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) diff --git a/src/cipsi/pt2_type.irp.f b/src/cipsi/pt2_type.irp.f index e6f31799..ef3fb882 100644 --- a/src/cipsi/pt2_type.irp.f +++ b/src/cipsi/pt2_type.irp.f @@ -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 diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi/run_pt2_slave.irp.f index 3c72dac0..d3f4d45d 100644 --- a/src/cipsi/run_pt2_slave.irp.f +++ b/src/cipsi/run_pt2_slave.irp.f @@ -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) diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index 69a8a4c3..650654be 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -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 diff --git a/src/cipsi/selection_types.f90 b/src/cipsi/selection_types.f90 index eef57aa5..8df37653 100644 --- a/src/cipsi/selection_types.f90 +++ b/src/cipsi/selection_types.f90 @@ -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 diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index e80bd856..953314b9 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -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 diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f index c34705f5..eb11e28c 100644 --- a/src/fci/pt2.irp.f +++ b/src/fci/pt2.irp.f @@ -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 diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 32c07ce5..d2aa4282 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -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