From 9d2c209c057ecddd1fd05ade1f84877f729ac7d8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 14 Feb 2017 15:44:39 +0100 Subject: [PATCH] target_pt2_ratio_cassd.irp.f --- plugins/CAS_SD_ZMQ/cassd_zmq.irp.f | 121 ------------------ plugins/CAS_SD_ZMQ/selection.irp.f | 118 +++++++++++++++++ .../CAS_SD_ZMQ/target_pt2_ratio_cassd.irp.f | 109 ++++++++++++++++ 3 files changed, 227 insertions(+), 121 deletions(-) create mode 100644 plugins/CAS_SD_ZMQ/target_pt2_ratio_cassd.irp.f diff --git a/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f b/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f index 881f74c3..5b364400 100644 --- a/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f +++ b/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f @@ -132,124 +132,3 @@ program fci_zmq call ezfio_set_cas_sd_zmq_energy_pt2(E_CI_before(1)+pt2(1)) end - - - - -subroutine ZMQ_selection(N_in, pt2) - use f77_zmq - use selection_types - - implicit none - - character*(512) :: task - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - integer, intent(in) :: N_in - type(selection_buffer) :: b - integer :: i, N - integer, external :: omp_get_thread_num - double precision, intent(out) :: pt2(N_states) - - - if (.True.) then - PROVIDE pt2_e0_denominator - N = max(N_in,1) - provide nproc - call new_parallel_job(zmq_to_qp_run_socket,"selection") - call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator)) - call zmq_set_running(zmq_to_qp_run_socket) - call create_selection_buffer(N, N*2, b) - endif - - integer :: i_generator, i_generator_start, i_generator_max, step -! step = int(max(1.,10*elec_num/mo_tot_num) - - step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num )) - step = max(1,step) - do i= 1, N_det_generators,step - i_generator_start = i - i_generator_max = min(i+step-1,N_det_generators) - write(task,*) i_generator_start, i_generator_max, 1, N - call add_task_to_taskserver(zmq_to_qp_run_socket,task) - end do - - !$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1) - i = omp_get_thread_num() - if (i==0) then - call selection_collector(b, pt2) - else - call selection_slave_inproc(i) - endif - !$OMP END PARALLEL - call end_parallel_job(zmq_to_qp_run_socket, 'selection') - if (N_in > 0) then - call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN - call copy_H_apply_buffer_to_wf() - if (s2_eig) then - call make_s2_eigenfunction - endif - endif -end subroutine - - -subroutine selection_slave_inproc(i) - implicit none - integer, intent(in) :: i - - call run_selection_slave(1,i,pt2_e0_denominator) -end - -subroutine selection_collector(b, pt2) - use f77_zmq - use selection_types - use bitmasks - implicit none - - - type(selection_buffer), intent(inout) :: b - double precision, intent(out) :: pt2(N_states) - double precision :: pt2_mwen(N_states) - 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 :: done - real :: time, time0 - 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)) - done = 0 - more = 1 - pt2(:) = 0d0 - call CPU_TIME(time0) - do while (more == 1) - call pull_selection_results(zmq_socket_pull, pt2_mwen, val(1), det(1,1,1), N, task_id, ntask) - pt2 += pt2_mwen - do i=1, N - call add_to_selection_buffer(b, det(1,1,i), val(i)) - 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 - done += ntask - call CPU_TIME(time) -! print *, "DONE" , done, time - time0 - end do - - - 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 - diff --git a/plugins/CAS_SD_ZMQ/selection.irp.f b/plugins/CAS_SD_ZMQ/selection.irp.f index db8ebbf0..70230e9e 100644 --- a/plugins/CAS_SD_ZMQ/selection.irp.f +++ b/plugins/CAS_SD_ZMQ/selection.irp.f @@ -1205,3 +1205,121 @@ subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting) end do genl end subroutine + +subroutine ZMQ_selection(N_in, pt2) + use f77_zmq + use selection_types + + implicit none + + character*(512) :: task + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + integer, intent(in) :: N_in + type(selection_buffer) :: b + integer :: i, N + integer, external :: omp_get_thread_num + double precision, intent(out) :: pt2(N_states) + + + if (.True.) then + PROVIDE pt2_e0_denominator + N = max(N_in,1) + provide nproc + call new_parallel_job(zmq_to_qp_run_socket,"selection") + call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator)) + call zmq_set_running(zmq_to_qp_run_socket) + call create_selection_buffer(N, N*2, b) + endif + + integer :: i_generator, i_generator_start, i_generator_max, step +! step = int(max(1.,10*elec_num/mo_tot_num) + + step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num )) + step = max(1,step) + do i= 1, N_det_generators,step + i_generator_start = i + i_generator_max = min(i+step-1,N_det_generators) + write(task,*) i_generator_start, i_generator_max, 1, N + call add_task_to_taskserver(zmq_to_qp_run_socket,task) + end do + + !$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call selection_collector(b, pt2) + else + call selection_slave_inproc(i) + endif + !$OMP END PARALLEL + call end_parallel_job(zmq_to_qp_run_socket, 'selection') + if (N_in > 0) then + call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN + call copy_H_apply_buffer_to_wf() + if (s2_eig) then + call make_s2_eigenfunction + endif + endif +end subroutine + + +subroutine selection_slave_inproc(i) + implicit none + integer, intent(in) :: i + + call run_selection_slave(1,i,pt2_e0_denominator) +end + +subroutine selection_collector(b, pt2) + use f77_zmq + use selection_types + use bitmasks + implicit none + + + type(selection_buffer), intent(inout) :: b + double precision, intent(out) :: pt2(N_states) + double precision :: pt2_mwen(N_states) + 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 :: done + real :: time, time0 + 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)) + done = 0 + more = 1 + pt2(:) = 0d0 + call CPU_TIME(time0) + do while (more == 1) + call pull_selection_results(zmq_socket_pull, pt2_mwen, val(1), det(1,1,1), N, task_id, ntask) + pt2 += pt2_mwen + do i=1, N + call add_to_selection_buffer(b, det(1,1,i), val(i)) + 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 + done += ntask + call CPU_TIME(time) +! print *, "DONE" , done, time - time0 + end do + + + 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 + diff --git a/plugins/CAS_SD_ZMQ/target_pt2_ratio_cassd.irp.f b/plugins/CAS_SD_ZMQ/target_pt2_ratio_cassd.irp.f new file mode 100644 index 00000000..cf934a46 --- /dev/null +++ b/plugins/CAS_SD_ZMQ/target_pt2_ratio_cassd.irp.f @@ -0,0 +1,109 @@ +program fci_zmq + implicit none + integer :: i,j,k + logical, external :: detEq + + double precision, allocatable :: pt2(:) + integer :: Nmin, Nmax + integer :: n_det_before, to_select + double precision :: threshold_davidson_in, ratio, E_ref + + double precision, allocatable :: psi_coef_ref(:,:) + integer(bit_kind), allocatable :: psi_det_ref(:,:,:) + + + allocate (pt2(N_states)) + + pt2 = 1.d0 + threshold_davidson_in = threshold_davidson + threshold_davidson = threshold_davidson_in * 100.d0 + SOFT_TOUCH threshold_davidson + + ! Stopping criterion is the PT2max + + double precision :: E_CI_before(N_states) + do while (dabs(pt2(1)) > pt2_max) + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + do k=1, N_states + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', CI_energy(k) + print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k) + enddo + print *, '-----' + E_CI_before(1:N_states) = CI_energy(1:N_states) + call ezfio_set_cas_sd_zmq_energy(CI_energy(1)) + + n_det_before = N_det + to_select = N_det + to_select = max(64-to_select, to_select) + call ZMQ_selection(to_select, pt2) + + PROVIDE psi_coef + PROVIDE psi_det + PROVIDE psi_det_sorted + + call diagonalize_CI + call save_wavefunction + call ezfio_set_cas_sd_zmq_energy(CI_energy(1)) + enddo + + threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) + threshold_generators = max(threshold_generators,threshold_generators_pt2) + threshold_davidson = threshold_davidson_in + TOUCH threshold_selectors threshold_generators threshold_davidson + call diagonalize_CI + call ZMQ_selection(0, pt2) + + E_ref = CI_energy(1) + pt2(1) + print *, 'Est FCI = ', E_ref + + Nmax = N_det + Nmin = 2 + allocate (psi_coef_ref(size(psi_coef_sorted,1),size(psi_coef_sorted,2))) + allocate (psi_det_ref(N_int,2,size(psi_det_sorted,3))) + psi_coef_ref = psi_coef_sorted + psi_det_ref = psi_det_sorted + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + TOUCH psi_coef psi_det + do while (Nmax-Nmin > 1) + psi_coef = psi_coef_ref + psi_det = psi_det_ref + TOUCH psi_det psi_coef + call diagonalize_CI + ratio = (CI_energy(1) - HF_energy) / (E_ref - HF_energy) + if (ratio < var_pt2_ratio) then + Nmin = N_det + else + Nmax = N_det + psi_coef_ref = psi_coef + psi_det_ref = psi_det + TOUCH psi_det psi_coef + endif + N_det = Nmin + (Nmax-Nmin)/2 + print *, '-----' + print *, 'Det min, Det max: ', Nmin, Nmax + print *, 'Ratio : ', ratio, ' ~ ', var_pt2_ratio + print *, 'N_det = ', N_det + print *, 'E = ', CI_energy(1) + call save_wavefunction + enddo + call ZMQ_selection(0, pt2) + print *, '------' + print *, 'HF_energy = ', HF_energy + print *, 'Est FCI = ', E_ref + print *, 'E = ', CI_energy(1) + print *, 'PT2 = ', pt2(1) + print *, 'E+PT2 = ', CI_energy(1)+pt2(1) + + E_CI_before(1:N_states) = CI_energy(1:N_states) + call save_wavefunction + call ezfio_set_cas_sd_zmq_energy(CI_energy(1)) + call ezfio_set_cas_sd_zmq_energy_pt2(E_CI_before(1)+pt2(1)) +end + + + +