diff --git a/.travis.yml b/.travis.yml index 5e032609..40c09bbc 100644 --- a/.travis.yml +++ b/.travis.yml @@ -26,7 +26,7 @@ python: script: - ./configure --production ./config/gfortran.cfg - - source ./quantum_package.rc ; qp_module.py install Full_CI Full_CI_ZMQ Hartree_Fock CAS_SD mrcepa0 All_singles + - source ./quantum_package.rc ; qp_module.py install Full_CI Full_CI_ZMQ Hartree_Fock CAS_SD_ZMQ mrcepa0 All_singles - source ./quantum_package.rc ; ninja - source ./quantum_package.rc ; cd ocaml ; make ; cd - - source ./quantum_package.rc ; cd tests ; ./run_tests.sh #-v diff --git a/plugins/CAS_SD_ZMQ/README.rst b/plugins/CAS_SD_ZMQ/README.rst new file mode 100644 index 00000000..45ba97e4 --- /dev/null +++ b/plugins/CAS_SD_ZMQ/README.rst @@ -0,0 +1,14 @@ +========== +CAS_SD_ZMQ +========== + +Selected CAS+SD module with Zero-MQ parallelization. + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f b/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f new file mode 100644 index 00000000..eb2d911f --- /dev/null +++ b/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f @@ -0,0 +1,234 @@ +program fci_zmq + implicit none + integer :: i,j,k + logical, external :: detEq + + double precision, allocatable :: pt2(:) + integer :: degree + + allocate (pt2(N_states)) + + pt2 = 1.d0 + diag_algorithm = "Lapack" + + if (N_det > N_det_max) then + call diagonalize_CI + call save_wavefunction + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = N_det_max + soft_touch N_det psi_det psi_coef + call diagonalize_CI + call save_wavefunction + 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+PT2 = ', CI_energy(k) + pt2(k) + print *, '-----' + enddo + endif + double precision :: E_CI_before(N_states) + + + integer :: n_det_before + print*,'Beginning the selection ...' + E_CI_before(1:N_states) = CI_energy(1:N_states) + + do while ( (N_det < N_det_max) .and. (maxval(abs(pt2(1:N_states))) > pt2_max) ) + n_det_before = N_det + call ZMQ_selection(max(256-N_det, N_det), pt2) + + PROVIDE psi_coef + PROVIDE psi_det + PROVIDE psi_det_sorted + + call diagonalize_CI + call save_wavefunction + + 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 *, '-----' + if(N_states.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_states + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_states + print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1)) + enddo + endif + E_CI_before(1:N_states) = CI_energy(1:N_states) + call ezfio_set_cas_sd_zmq_energy(CI_energy(1)) + enddo + + integer :: exc_max, degree_min + exc_max = 0 + print *, 'CAS determinants : ', N_det_cas + do i=1,min(N_det_cas,10) + do k=i,N_det_cas + call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int) + exc_max = max(exc_max,degree) + enddo + print *, psi_cas_coef(i,:) + call debug_det(psi_cas(1,1,i),N_int) + print *, '' + enddo + print *, 'Max excitation degree in the CAS :', exc_max + + if(do_pt2_end)then + print*,'Last iteration only to compute the PT2' + threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) + threshold_generators = max(threshold_generators,threshold_generators_pt2) + TOUCH threshold_selectors threshold_generators + E_CI_before(1:N_states) = CI_energy(1:N_states) + call ZMQ_selection(0, pt2) + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + do k=1,N_states + print *, 'State', k + print *, 'PT2 = ', pt2 + print *, 'E = ', E_CI_before + print *, 'E+PT2 = ', E_CI_before+pt2 + print *, '-----' + enddo + call ezfio_set_cas_sd_zmq_energy_pt2(E_CI_before+pt2) + endif + call save_wavefunction + call ezfio_set_cas_sd_zmq_energy(CI_energy(1)) + call ezfio_set_cas_sd_zmq_energy_pt2(E_CI_before+pt2) + +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= N_det_generators, 1, -step + i_generator_start = max(i-step+1,1) + i_generator_max = i + 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/e_corr_selectors.irp.f b/plugins/CAS_SD_ZMQ/e_corr_selectors.irp.f new file mode 100644 index 00000000..fec480f0 --- /dev/null +++ b/plugins/CAS_SD_ZMQ/e_corr_selectors.irp.f @@ -0,0 +1,79 @@ + +use bitmasks + BEGIN_PROVIDER [integer, exc_degree_per_selectors, (N_det_selectors)] +&BEGIN_PROVIDER [integer, double_index_selectors, (N_det_selectors)] +&BEGIN_PROVIDER [integer, n_double_selectors] + implicit none + BEGIN_DOC + ! degree of excitation respect to Hartree Fock for the wave function + ! + ! for the all the selectors determinants + ! + ! double_index_selectors = list of the index of the double excitations + ! + ! n_double_selectors = number of double excitations in the selectors determinants + END_DOC + integer :: i,degree + n_double_selectors = 0 + do i = 1, N_det_selectors + call get_excitation_degree(psi_selectors(1,1,i),ref_bitmask,degree,N_int) + exc_degree_per_selectors(i) = degree + if(degree==2)then + n_double_selectors += 1 + double_index_selectors(n_double_selectors) =i + endif + enddo +END_PROVIDER + + BEGIN_PROVIDER[double precision, coef_hf_selector] + &BEGIN_PROVIDER[double precision, inv_selectors_coef_hf] + &BEGIN_PROVIDER[double precision, inv_selectors_coef_hf_squared] + &BEGIN_PROVIDER[double precision, E_corr_per_selectors, (N_det_selectors)] + &BEGIN_PROVIDER[double precision, i_H_HF_per_selectors, (N_det_selectors)] + &BEGIN_PROVIDER[double precision, Delta_E_per_selector, (N_det_selectors)] + &BEGIN_PROVIDER[double precision, E_corr_double_only ] + &BEGIN_PROVIDER[double precision, E_corr_second_order ] + implicit none + BEGIN_DOC + ! energy of correlation per determinant respect to the Hartree Fock determinant + ! + ! for the all the double excitations in the selectors determinants + ! + ! E_corr_per_selectors(i) = * c(D_i)/c(HF) if |D_i> is a double excitation + ! + ! E_corr_per_selectors(i) = -1000.d0 if it is not a double excitation + ! + ! coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants + END_DOC + PROVIDE ref_bitmask_energy psi_selectors ref_bitmask N_int psi_selectors + integer :: i,degree + double precision :: hij,diag_H_mat_elem + E_corr_double_only = 0.d0 + E_corr_second_order = 0.d0 + do i = 1, N_det_selectors + if(exc_degree_per_selectors(i)==2)then + call i_H_j(ref_bitmask,psi_selectors(1,1,i),N_int,hij) + i_H_HF_per_selectors(i) = hij + E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij + E_corr_double_only += E_corr_per_selectors(i) +! E_corr_second_order += hij * hij /(ref_bitmask_energy - diag_H_mat_elem(psi_selectors(1,1,i),N_int)) + elseif(exc_degree_per_selectors(i) == 0)then + coef_hf_selector = psi_selectors_coef(i,1) + E_corr_per_selectors(i) = -1000.d0 + Delta_E_per_selector(i) = 0.d0 + else + E_corr_per_selectors(i) = -1000.d0 + endif + enddo + if (dabs(coef_hf_selector) > 1.d-8) then + inv_selectors_coef_hf = 1.d0/coef_hf_selector + inv_selectors_coef_hf_squared = inv_selectors_coef_hf * inv_selectors_coef_hf + else + inv_selectors_coef_hf = 0.d0 + inv_selectors_coef_hf_squared = 0.d0 + endif + do i = 1,n_double_selectors + E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf + enddo + E_corr_double_only = E_corr_double_only * inv_selectors_coef_hf + END_PROVIDER diff --git a/plugins/CAS_SD_ZMQ/energy.irp.f b/plugins/CAS_SD_ZMQ/energy.irp.f new file mode 100644 index 00000000..4999c176 --- /dev/null +++ b/plugins/CAS_SD_ZMQ/energy.irp.f @@ -0,0 +1,11 @@ +BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ] + implicit none + BEGIN_DOC + ! E0 in the denominator of the PT2 + END_DOC + pt2_E0_denominator(:) = CI_electronic_energy(:) +! pt2_E0_denominator(:) = HF_energy - nuclear_repulsion +! pt2_E0_denominator(:) = barycentric_electronic_energy(:) + call write_double(6,pt2_E0_denominator(1)+nuclear_repulsion, 'PT2 Energy denominator') +END_PROVIDER + diff --git a/plugins/CAS_SD_ZMQ/ezfio_interface.irp.f b/plugins/CAS_SD_ZMQ/ezfio_interface.irp.f new file mode 100644 index 00000000..8adab518 --- /dev/null +++ b/plugins/CAS_SD_ZMQ/ezfio_interface.irp.f @@ -0,0 +1,4 @@ +! DO NOT MODIFY BY HAND +! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py +! from file /home/scemama/quantum_package/src/CAS_SD_ZMQ/EZFIO.cfg + diff --git a/plugins/CAS_SD_ZMQ/run_selection_slave.irp.f b/plugins/CAS_SD_ZMQ/run_selection_slave.irp.f new file mode 100644 index 00000000..36550116 --- /dev/null +++ b/plugins/CAS_SD_ZMQ/run_selection_slave.irp.f @@ -0,0 +1,156 @@ + +subroutine run_selection_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, buf2 + logical :: done + double precision :: pt2(N_states) + + 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 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) + return + end if + buf%N = 0 + ctask = 1 + pt2 = 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_generator_start, i_generator_max, step, N + read (task,*) i_generator_start, i_generator_max, step, N + if(buf%N == 0) then + ! Only first time + call create_selection_buffer(N, N*2, buf) + call create_selection_buffer(N, N*3, buf2) + else + if(N /= buf%N) stop "N changed... wtf man??" + end if + !print *, "psi_selectors_coef ", psi_selectors_coef(N_det_selectors-5:N_det_selectors, 1) + !call debug_det(psi_selectors(1,1,N_det_selectors), N_int) + do i_generator=i_generator_start,i_generator_max,step + call select_connected(i_generator,energy,pt2,buf) + 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_selection_results(zmq_socket_push, pt2, buf, task_id(1), ctask) + do i=1,buf%cur + call add_to_selection_buffer(buf2, buf%det(1,1,i), buf%val(i)) + enddo + call sort_selection_buffer(buf2) + buf%mini = buf2%mini + pt2 = 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) +end subroutine + + +subroutine push_selection_results(zmq_socket_push, pt2, b, task_id, ntask) + use f77_zmq + use selection_types + implicit none + + integer(ZMQ_PTR), intent(in) :: zmq_socket_push + double precision, intent(in) :: pt2(N_states) + type(selection_buffer), intent(inout) :: b + integer, intent(in) :: ntask, task_id(*) + integer :: rc + + call sort_selection_buffer(b) + + rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE) + if(rc /= 4) stop "push" + rc = f77_zmq_send( zmq_socket_push, pt2, 8*N_states, ZMQ_SNDMORE) + if(rc /= 8*N_states) stop "push" + + rc = f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE) + if(rc /= 8*b%cur) stop "push" + + rc = f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE) + if(rc /= bit_kind*N_int*2*b%cur) 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(1), ntask*4, 0) + if(rc /= 4*ntask) stop "push" + +! Activate is zmq_socket_push is a REQ +! rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0) +end subroutine + + +subroutine pull_selection_results(zmq_socket_pull, pt2, val, det, N, task_id, ntask) + use f77_zmq + use selection_types + implicit none + integer(ZMQ_PTR), intent(in) :: zmq_socket_pull + double precision, intent(inout) :: pt2(N_states) + double precision, intent(out) :: val(*) + integer(bit_kind), intent(out) :: det(N_int, 2, *) + 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, pt2, N_states*8, 0) + if(rc /= 8*N_states) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, 0) + if(rc /= 8*N) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, 0) + if(rc /= bit_kind*N_int*2*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(1), ntask*4, 0) + if(rc /= 4*ntask) stop "pull" + +! Activate is zmq_socket_pull is a REP +! rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0) +end subroutine + + + diff --git a/plugins/CAS_SD_ZMQ/selection.irp.f b/plugins/CAS_SD_ZMQ/selection.irp.f new file mode 100644 index 00000000..6e7ba359 --- /dev/null +++ b/plugins/CAS_SD_ZMQ/selection.irp.f @@ -0,0 +1,1202 @@ +use bitmasks + + +double precision function integral8(i,j,k,l) + implicit none + + integer, intent(in) :: i,j,k,l + double precision, external :: get_mo_bielec_integral + integer :: ii + ii = l-mo_integrals_cache_min + ii = ior(ii, k-mo_integrals_cache_min) + ii = ior(ii, j-mo_integrals_cache_min) + ii = ior(ii, i-mo_integrals_cache_min) + if (iand(ii, -64) /= 0) then + integral8 = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + else + ii = l-mo_integrals_cache_min + ii = ior( ishft(ii,6), k-mo_integrals_cache_min) + ii = ior( ishft(ii,6), j-mo_integrals_cache_min) + ii = ior( ishft(ii,6), i-mo_integrals_cache_min) + integral8 = mo_integrals_cache(ii) + endif +end function + + +BEGIN_PROVIDER [ integer(1), psi_phasemask, (N_int*bit_kind_size, 2, N_det)] + use bitmasks + implicit none + + integer :: i + do i=1, N_det + call get_mask_phase(psi_selectors(1,1,i), psi_phasemask(1,1,i)) + end do +END_PROVIDER + + +subroutine assert(cond, msg) + character(*), intent(in) :: msg + logical, intent(in) :: cond + + if(.not. cond) then + print *, "assert fail: "//msg + stop + end if +end subroutine + + +subroutine get_mask_phase(det, phasemask) + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: det(N_int, 2) + integer(1), intent(out) :: phasemask(N_int*bit_kind_size, 2) + integer :: s, ni, i + logical :: change + + phasemask = 0_1 + do s=1,2 + change = .false. + do ni=1,N_int + do i=0,bit_kind_size-1 + if(BTEST(det(ni, s), i)) change = .not. change + if(change) phasemask((ni-1)*bit_kind_size + i + 1, s) = 1_1 + end do + end do + end do +end subroutine + + +subroutine select_connected(i_generator,E0,pt2,b) + use bitmasks + use selection_types + implicit none + integer, intent(in) :: i_generator + type(selection_buffer), intent(inout) :: b + double precision, intent(inout) :: pt2(N_states) + integer :: k,l + double precision, intent(in) :: E0(N_states) + + integer(bit_kind) :: hole_mask(N_int,2), particle_mask(N_int,2) + double precision :: fock_diag_tmp(2,mo_tot_num+1) + + call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int) + + do l=1,N_generators_bitmask + do k=1,N_int + hole_mask(k,1) = iand(generators_bitmask(k,1,s_hole,l), psi_det_generators(k,1,i_generator)) + hole_mask(k,2) = iand(generators_bitmask(k,2,s_hole,l), psi_det_generators(k,2,i_generator)) + particle_mask(k,1) = iand(generators_bitmask(k,1,s_part,l), not(psi_det_generators(k,1,i_generator)) ) + particle_mask(k,2) = iand(generators_bitmask(k,2,s_part,l), not(psi_det_generators(k,2,i_generator)) ) + + enddo + call select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b) + call select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b) + enddo +end subroutine + + +double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2) + use bitmasks + implicit none + + integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2) + integer, intent(in) :: s1, s2, h1, h2, p1, p2 + logical :: change + integer(1) :: np + double precision, parameter :: res(0:1) = (/1d0, -1d0/) + + np = phasemask(h1,s1) + phasemask(p1,s1) + phasemask(h2,s2) + phasemask(p2,s2) + if(p1 < h1) np = np + 1_1 + if(p2 < h2) np = np + 1_1 + + if(s1 == s2 .and. max(h1, p1) > min(h2, p2)) np = np + 1_1 + get_phase_bi = res(iand(np,1_1)) +end subroutine + + + +! Selection single +! ---------------- + +subroutine select_singles(i_gen,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf) + use bitmasks + use selection_types + implicit none + BEGIN_DOC +! Select determinants connected to i_det by H + END_DOC + integer, intent(in) :: i_gen + integer(bit_kind), intent(in) :: hole_mask(N_int,2), particle_mask(N_int,2) + double precision, intent(in) :: fock_diag_tmp(mo_tot_num) + double precision, intent(in) :: E0(N_states) + double precision, intent(inout) :: pt2(N_states) + type(selection_buffer), intent(inout) :: buf + + double precision :: vect(N_states, mo_tot_num) + logical :: bannedOrb(mo_tot_num) + integer :: i, j, k + integer :: h1,h2,s1,s2,i1,i2,ib,sp + integer(bit_kind) :: hole(N_int,2), particle(N_int,2), mask(N_int, 2) + logical :: fullMatch, ok + + + do k=1,N_int + hole (k,1) = iand(psi_det_generators(k,1,i_gen), hole_mask(k,1)) + hole (k,2) = iand(psi_det_generators(k,2,i_gen), hole_mask(k,2)) + particle(k,1) = iand(not(psi_det_generators(k,1,i_gen)), particle_mask(k,1)) + particle(k,2) = iand(not(psi_det_generators(k,2,i_gen)), particle_mask(k,2)) + enddo + + ! Create lists of holes and particles + ! ----------------------------------- + + integer :: N_holes(2), N_particles(2) + integer :: hole_list(N_int*bit_kind_size,2) + integer :: particle_list(N_int*bit_kind_size,2) + + call bitstring_to_list_ab(hole , hole_list , N_holes , N_int) + call bitstring_to_list_ab(particle, particle_list, N_particles, N_int) + + do sp=1,2 + do i=1, N_holes(sp) + h1 = hole_list(i,sp) + call apply_hole(psi_det_generators(1,1,i_gen), sp, h1, mask, ok, N_int) + bannedOrb = .true. + do j=1,N_particles(sp) + bannedOrb(particle_list(j, sp)) = .false. + end do + call spot_hasBeen(mask, sp, psi_selectors, i_gen, N_det, bannedOrb, fullMatch) + if(fullMatch) cycle + vect = 0d0 + call splash_p(mask, sp, psi_selectors(1,1,i_gen), psi_phasemask(1,1,i_gen), psi_selectors_coef_transp(1,i_gen), N_det_selectors - i_gen + 1, bannedOrb, vect) + call fill_buffer_single(i_gen, sp, h1, bannedOrb, fock_diag_tmp, E0, pt2, vect, buf) + end do + enddo +end subroutine + + +subroutine fill_buffer_single(i_generator, sp, h1, bannedOrb, fock_diag_tmp, E0, pt2, vect, buf) + use bitmasks + use selection_types + implicit none + + integer, intent(in) :: i_generator, sp, h1 + double precision, intent(in) :: vect(N_states, mo_tot_num) + logical, intent(in) :: bannedOrb(mo_tot_num) + double precision, intent(in) :: fock_diag_tmp(mo_tot_num) + double precision, intent(in) :: E0(N_states) + double precision, intent(inout) :: pt2(N_states) + type(selection_buffer), intent(inout) :: buf + logical :: ok + integer :: s1, s2, p1, p2, ib, istate + integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) + double precision :: e_pert, delta_E, val, Hii, max_e_pert, tmp + double precision, external :: diag_H_mat_elem_fock + + + call apply_hole(psi_det_generators(1,1,i_generator), sp, h1, mask, ok, N_int) + + do p1=1,mo_tot_num + if(bannedOrb(p1)) cycle + if(vect(1, p1) == 0d0) cycle + call apply_particle(mask, sp, p1, det, ok, N_int) + + + Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) + max_e_pert = 0d0 + + do istate=1,N_states + val = vect(istate, p1) + vect(istate, p1) + delta_E = E0(istate) - Hii + tmp = dsqrt(delta_E * delta_E + val * val) + if (delta_E < 0.d0) then + tmp = -tmp + endif + e_pert = 0.5d0 * ( tmp - delta_E) + pt2(istate) += e_pert + if(dabs(e_pert) > dabs(max_e_pert)) max_e_pert = e_pert + end do + + if(dabs(max_e_pert) > buf%mini) call add_to_selection_buffer(buf, det, max_e_pert) + end do +end subroutine + + +subroutine splash_p(mask, sp, det, phasemask, coefs, N_sel, bannedOrb, vect) + use bitmasks + implicit none + + integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int,2,N_sel) + integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2, N_sel) + double precision, intent(in) :: coefs(N_states, N_sel) + integer, intent(in) :: sp, N_sel + logical, intent(inout) :: bannedOrb(mo_tot_num) + double precision, intent(inout) :: vect(N_states, mo_tot_num) + + integer :: i, j, h(0:2,2), p(0:3,2), nt + integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2) + + do i=1,N_int + negMask(i,1) = not(mask(i,1)) + negMask(i,2) = not(mask(i,2)) + end do + + do i=1, N_sel + nt = 0 + do j=1,N_int + mobMask(j,1) = iand(negMask(j,1), det(j,1,i)) + mobMask(j,2) = iand(negMask(j,2), det(j,2,i)) + nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) + end do + + if(nt > 3) cycle + + do j=1,N_int + perMask(j,1) = iand(mask(j,1), not(det(j,1,i))) + perMask(j,2) = iand(mask(j,2), not(det(j,2,i))) + end do + + call bitstring_to_list(perMask(1,1), h(1,1), h(0,1), N_int) + call bitstring_to_list(perMask(1,2), h(1,2), h(0,2), N_int) + + call bitstring_to_list(mobMask(1,1), p(1,1), p(0,1), N_int) + call bitstring_to_list(mobMask(1,2), p(1,2), p(0,2), N_int) + + if(nt == 3) then + call get_m2(det(1,1,i), phasemask(1,1,i), bannedOrb, vect, mask, h, p, sp, coefs(1, i)) + else if(nt == 2) then + call get_m1(det(1,1,i), phasemask(1,1,i), bannedOrb, vect, mask, h, p, sp, coefs(1, i)) + else + call get_m0(det(1,1,i), phasemask(1,1,i), bannedOrb, vect, mask, h, p, sp, coefs(1, i)) + end if + end do +end subroutine + + +subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2) + integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2) + logical, intent(in) :: bannedOrb(mo_tot_num) + double precision, intent(in) :: coefs(N_states) + double precision, intent(inout) :: vect(N_states, mo_tot_num) + integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2) + integer :: i, j, h1, h2, p1, p2, sfix, hfix, pfix, hmob, pmob, puti + double precision :: hij + double precision, external :: get_phase_bi, integral8 + + integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) + integer, parameter :: turn2(2) = (/2,1/) + + if(h(0,sp) == 2) then + h1 = h(1, sp) + h2 = h(2, sp) + do i=1,3 + puti = p(i, sp) + if(bannedOrb(puti)) cycle + p1 = p(turn3_2(1,i), sp) + p2 = p(turn3_2(2,i), sp) + hij = integral8(p1, p2, h1, h2) - integral8(p2, p1, h1, h2) + hij *= get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2) + vect(:, puti) += hij * coefs + end do + else if(h(0,sp) == 1) then + sfix = turn2(sp) + hfix = h(1,sfix) + pfix = p(1,sfix) + hmob = h(1,sp) + do j=1,2 + puti = p(j, sp) + if(bannedOrb(puti)) cycle + pmob = p(turn2(j), sp) + hij = integral8(pfix, pmob, hfix, hmob) + hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix) + vect(:, puti) += hij * coefs + end do + else + puti = p(1,sp) + if(.not. bannedOrb(puti)) then + sfix = turn2(sp) + p1 = p(1,sfix) + p2 = p(2,sfix) + h1 = h(1,sfix) + h2 = h(2,sfix) + hij = (integral8(p1,p2,h1,h2) - integral8(p2,p1,h1,h2)) + hij *= get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2) + vect(:, puti) += hij * coefs + end if + end if +end subroutine + + + +subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2) + integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2) + logical, intent(in) :: bannedOrb(mo_tot_num) + double precision, intent(in) :: coefs(N_states) + double precision, intent(inout) :: vect(N_states, mo_tot_num) + integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2) + integer :: i, hole, p1, p2, sh + logical :: ok, lbanned(mo_tot_num) + integer(bit_kind) :: det(N_int, 2) + double precision :: hij + double precision, external :: get_phase_bi, integral8 + + lbanned = bannedOrb + sh = 1 + if(h(0,2) == 1) sh = 2 + hole = h(1, sh) + lbanned(p(1,sp)) = .true. + if(p(0,sp) == 2) lbanned(p(2,sp)) = .true. + !print *, "SPm1", sp, sh + + p1 = p(1, sp) + + if(sp == sh) then + p2 = p(2, sp) + lbanned(p2) = .true. + + do i=1,hole-1 + if(lbanned(i)) cycle + hij = (integral8(p1, p2, i, hole) - integral8(p2, p1, i, hole)) + hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2) + vect(:,i) += hij * coefs + end do + do i=hole+1,mo_tot_num + if(lbanned(i)) cycle + hij = (integral8(p1, p2, hole, i) - integral8(p2, p1, hole, i)) + hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2) + vect(:,i) += hij * coefs + end do + + call apply_particle(mask, sp, p2, det, ok, N_int) + call i_h_j(gen, det, N_int, hij) + vect(:, p2) += hij * coefs + else + p2 = p(1, sh) + do i=1,mo_tot_num + if(lbanned(i)) cycle + hij = integral8(p1, p2, i, hole) + hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2) + vect(:,i) += hij * coefs + end do + end if + + call apply_particle(mask, sp, p1, det, ok, N_int) + call i_h_j(gen, det, N_int, hij) + vect(:, p1) += hij * coefs +end subroutine + + +subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2) + integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2) + logical, intent(in) :: bannedOrb(mo_tot_num) + double precision, intent(in) :: coefs(N_states) + double precision, intent(inout) :: vect(N_states, mo_tot_num) + integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2) + integer :: i + logical :: ok, lbanned(mo_tot_num) + integer(bit_kind) :: det(N_int, 2) + double precision :: hij + + lbanned = bannedOrb + lbanned(p(1,sp)) = .true. + do i=1,mo_tot_num + if(lbanned(i)) cycle + call apply_particle(mask, sp, i, det, ok, N_int) + call i_h_j(gen, det, N_int, hij) + vect(:, i) += hij * coefs + end do +end subroutine + + +subroutine spot_hasBeen(mask, sp, det, i_gen, N, banned, fullMatch) + use bitmasks + implicit none + + integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N) + integer, intent(in) :: i_gen, N, sp + logical, intent(inout) :: banned(mo_tot_num) + logical, intent(out) :: fullMatch + + + integer :: i, j, na, nb, list(3), nt + integer(bit_kind) :: myMask(N_int, 2), negMask(N_int, 2) + + fullMatch = .false. + + do i=1,N_int + negMask(i,1) = not(mask(i,1)) + negMask(i,2) = not(mask(i,2)) + end do + + do i=1, N + nt = 0 + + do j=1, N_int + myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1)) + myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2)) + nt += popcnt(myMask(j, 1)) + popcnt(myMask(j, 2)) + end do + + if(nt > 3) cycle + + if(nt <= 2 .and. i < i_gen) then + fullMatch = .true. + return + end if + + call bitstring_to_list(myMask(1,sp), list(1), na, N_int) + + if(nt == 3 .and. i < i_gen) then + do j=1,na + banned(list(j)) = .true. + end do + else if(nt == 1 .and. na == 1) then + banned(list(1)) = .true. + end if + end do +end subroutine + + + + +! Selection double +! ---------------- + +subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf) + use bitmasks + use selection_types + implicit none + + integer, intent(in) :: i_generator + integer(bit_kind), intent(in) :: hole_mask(N_int,2), particle_mask(N_int,2) + double precision, intent(in) :: fock_diag_tmp(mo_tot_num) + double precision, intent(in) :: E0(N_states) + double precision, intent(inout) :: pt2(N_states) + type(selection_buffer), intent(inout) :: buf + + double precision :: mat(N_states, mo_tot_num, mo_tot_num) + integer :: h1,h2,s1,s2,s3,i1,i2,ib,sp,k,i,j,nt,ii + integer(bit_kind) :: hole(N_int,2), particle(N_int,2), mask(N_int, 2), pmask(N_int, 2) + logical :: fullMatch, ok + + integer(bit_kind) :: mobMask(N_int, 2), negMask(N_int, 2) + integer,allocatable :: preinteresting(:), prefullinteresting(:), interesting(:), fullinteresting(:) + integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) + + allocate(minilist(N_int, 2, N_det_selectors), fullminilist(N_int, 2, N_det)) + allocate(preinteresting(0:N_det_selectors), prefullinteresting(0:N_det), interesting(0:N_det_selectors), fullinteresting(0:N_det)) + + do k=1,N_int + hole (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1)) + hole (k,2) = iand(psi_det_generators(k,2,i_generator), hole_mask(k,2)) + particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), particle_mask(k,1)) + particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), particle_mask(k,2)) + enddo + + integer :: N_holes(2), N_particles(2) + integer :: hole_list(N_int*bit_kind_size,2) + integer :: particle_list(N_int*bit_kind_size,2) + + call bitstring_to_list_ab(hole , hole_list , N_holes , N_int) + call bitstring_to_list_ab(particle, particle_list, N_particles, N_int) + + + preinteresting(0) = 0 + prefullinteresting(0) = 0 + + do i=1,N_int + negMask(i,1) = not(psi_det_generators(i,1,i_generator)) + negMask(i,2) = not(psi_det_generators(i,2,i_generator)) + end do + + do i=1,N_det + nt = 0 + do j=1,N_int + mobMask(j,1) = iand(negMask(j,1), psi_selectors(j,1,i)) + mobMask(j,2) = iand(negMask(j,2), psi_selectors(j,2,i)) + nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) + end do + + if(nt <= 4) then + if(i <= N_det_selectors) then + preinteresting(0) += 1 + preinteresting(preinteresting(0)) = i + else if(nt <= 2) then + prefullinteresting(0) += 1 + prefullinteresting(prefullinteresting(0)) = i + end if + end if + end do + + + do s1=1,2 + do i1=N_holes(s1),1,-1 ! Generate low excitations first + h1 = hole_list(i1,s1) + call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int) + + do i=1,N_int + negMask(i,1) = not(pmask(i,1)) + negMask(i,2) = not(pmask(i,2)) + end do + + interesting(0) = 0 + fullinteresting(0) = 0 + + do ii=1,preinteresting(0) + i = preinteresting(ii) + nt = 0 + do j=1,N_int + mobMask(j,1) = iand(negMask(j,1), psi_selectors(j,1,i)) + mobMask(j,2) = iand(negMask(j,2), psi_selectors(j,2,i)) + nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) + end do + + if(nt <= 4) then + interesting(0) += 1 + interesting(interesting(0)) = i + minilist(:,:,interesting(0)) = psi_selectors(:,:,i) + if(nt <= 2) then + fullinteresting(0) += 1 + fullinteresting(fullinteresting(0)) = i + fullminilist(:,:,fullinteresting(0)) = psi_selectors(:,:,i) + end if + end if + end do + + do ii=1,prefullinteresting(0) + i = prefullinteresting(ii) + nt = 0 + do j=1,N_int + mobMask(j,1) = iand(negMask(j,1), psi_selectors(j,1,i)) + mobMask(j,2) = iand(negMask(j,2), psi_selectors(j,2,i)) + nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) + end do + + if(nt <= 2) then + fullinteresting(0) += 1 + fullinteresting(fullinteresting(0)) = i + fullminilist(:,:,fullinteresting(0)) = psi_selectors(:,:,i) + end if + end do + + do s2=s1,2 + sp = s1 + if(s1 /= s2) sp = 3 + + ib = 1 + if(s1 == s2) ib = i1+1 + do i2=N_holes(s2),ib,-1 ! Generate low excitations first + + h2 = hole_list(i2,s2) + call apply_hole(pmask, s2,h2, mask, ok, N_int) + + logical :: banned(mo_tot_num, mo_tot_num,2) + logical :: bannedOrb(mo_tot_num, 2) + + banned = .false. + + call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting) + + if(fullMatch) cycle + + bannedOrb(1:mo_tot_num, 1:2) = .true. + do s3=1,2 + do i=1,N_particles(s3) + bannedOrb(particle_list(i,s3), s3) = .false. + enddo + enddo + + mat = 0d0 + call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) + call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf) + enddo + enddo + enddo + enddo +end subroutine + + +subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf) + use bitmasks + use selection_types + implicit none + + integer, intent(in) :: i_generator, sp, h1, h2 + double precision, intent(in) :: mat(N_states, mo_tot_num, mo_tot_num) + logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num) + double precision, intent(in) :: fock_diag_tmp(mo_tot_num) + double precision, intent(in) :: E0(N_states) + double precision, intent(inout) :: pt2(N_states) + type(selection_buffer), intent(inout) :: buf + logical :: ok + integer :: s1, s2, p1, p2, ib, j, istate + integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) + double precision :: e_pert, delta_E, val, Hii, max_e_pert,tmp + double precision, external :: diag_H_mat_elem_fock + + logical, external :: detEq + + + if(sp == 3) then + s1 = 1 + s2 = 2 + else + s1 = sp + s2 = sp + end if + + call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) + + do p1=1,mo_tot_num + if(bannedOrb(p1, s1)) cycle + ib = 1 + if(sp /= 3) ib = p1+1 + do p2=ib,mo_tot_num + if(bannedOrb(p2, s2)) cycle + if(banned(p1,p2)) cycle + if(mat(1, p1, p2) == 0d0) cycle + call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) + + + Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) + max_e_pert = 0d0 + + do istate=1,N_states + delta_E = E0(istate) - Hii + val = mat(istate, p1, p2) + mat(istate, p1, p2) + tmp = dsqrt(delta_E * delta_E + val * val) + if (delta_E < 0.d0) then + tmp = -tmp + endif + e_pert = 0.5d0 * ( tmp - delta_E) + pt2(istate) = pt2(istate) + e_pert + max_e_pert = min(e_pert,max_e_pert) + end do + + if(dabs(max_e_pert) > buf%mini) then + call add_to_selection_buffer(buf, det, max_e_pert) + end if + end do + end do +end subroutine + + +subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting) + use bitmasks + implicit none + + integer, intent(in) :: interesting(0:N_sel) + + integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel) + integer, intent(in) :: sp, i_gen, N_sel + logical, intent(inout) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2) + double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) + + integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt + integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2) +! logical :: bandon +! +! bandon = .false. + mat = 0d0 + + do i=1,N_int + negMask(i,1) = not(mask(i,1)) + negMask(i,2) = not(mask(i,2)) + end do + + do i=1, N_sel ! interesting(0) + !i = interesting(ii) + + nt = 0 + do j=1,N_int + mobMask(j,1) = iand(negMask(j,1), det(j,1,i)) + mobMask(j,2) = iand(negMask(j,2), det(j,2,i)) + nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) + end do + + if(nt > 4) cycle + + do j=1,N_int + perMask(j,1) = iand(mask(j,1), not(det(j,1,i))) + perMask(j,2) = iand(mask(j,2), not(det(j,2,i))) + end do + + call bitstring_to_list(perMask(1,1), h(1,1), h(0,1), N_int) + call bitstring_to_list(perMask(1,2), h(1,2), h(0,2), N_int) + + call bitstring_to_list(mobMask(1,1), p(1,1), p(0,1), N_int) + call bitstring_to_list(mobMask(1,2), p(1,2), p(0,2), N_int) + + if(interesting(i) < i_gen) then + if(nt == 4) call past_d2(banned, p, sp) + if(nt == 3) call past_d1(bannedOrb, p) + else + if(interesting(i) == i_gen) then +! bandon = .true. + if(sp == 3) then + banned(:,:,2) = transpose(banned(:,:,1)) + else + do k=1,mo_tot_num + do l=k+1,mo_tot_num + banned(l,k,1) = banned(k,l,1) + end do + end do + end if + end if + if(nt == 4) then + call get_d2(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) + else if(nt == 3) then + call get_d1(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) + else + call get_d0(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) + end if + end if + end do +end subroutine + + +subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2) + integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2) + logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2) + double precision, intent(in) :: coefs(N_states) + double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) + integer, intent(in) :: h(0:2,2), p(0:4,2), sp + + double precision, external :: get_phase_bi, integral8 + + integer :: i, j, tip, ma, mi, puti, putj + integer :: h1, h2, p1, p2, i1, i2 + double precision :: hij, phase + + integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0, 3,4, 0,0, 0,0, 2,4, 1,4, 0,0, 2,3, 1,3, 1,2 /), (/2,3,4/)) + integer, parameter :: turn2(2) = (/2, 1/) + integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) + + integer :: bant + bant = 1 + + tip = p(0,1) * p(0,2) + + ma = sp + if(p(0,1) > p(0,2)) ma = 1 + if(p(0,1) < p(0,2)) ma = 2 + mi = mod(ma, 2) + 1 + + if(sp == 3) then + if(ma == 2) bant = 2 + + if(tip == 3) then + puti = p(1, mi) + do i = 1, 3 + putj = p(i, ma) + if(banned(putj,puti,bant)) cycle + i1 = turn3(1,i) + i2 = turn3(2,i) + p1 = p(i1, ma) + p2 = p(i2, ma) + h1 = h(1, ma) + h2 = h(2, ma) + + hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) + if(ma == 1) then + mat(:, putj, puti) += coefs * hij + else + mat(:, puti, putj) += coefs * hij + end if + end do + else + do i = 1,2 + do j = 1,2 + puti = p(i, 1) + putj = p(j, 2) + + if(banned(puti,putj,bant)) cycle + p1 = p(turn2(i), 1) + p2 = p(turn2(j), 2) + h1 = h(1,1) + h2 = h(1,2) + + hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) + mat(:, puti, putj) += coefs * hij + end do + end do + end if + + else + if(tip == 0) then + h1 = h(1, ma) + h2 = h(2, ma) + do i=1,3 + puti = p(i, ma) + do j=i+1,4 + putj = p(j, ma) + if(banned(puti,putj,1)) cycle + + i1 = turn2d(1, i, j) + i2 = turn2d(2, i, j) + p1 = p(i1, ma) + p2 = p(i2, ma) + hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) + mat(:, puti, putj) += coefs * hij + end do + end do + else if(tip == 3) then + h1 = h(1, mi) + h2 = h(1, ma) + p1 = p(1, mi) + do i=1,3 + puti = p(turn3(1,i), ma) + putj = p(turn3(2,i), ma) + if(banned(puti,putj,1)) cycle + p2 = p(i, ma) + + hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2) + mat(:, min(puti, putj), max(puti, putj)) += coefs * hij + end do + else ! tip == 4 + puti = p(1, sp) + putj = p(2, sp) + if(.not. banned(puti,putj,1)) then + p1 = p(1, mi) + p2 = p(2, mi) + h1 = h(1, mi) + h2 = h(2, mi) + hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2) + mat(:, puti, putj) += coefs * hij + end if + end if + end if +end subroutine + + +subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2) + integer(1),intent(in) :: phasemask(N_int*bit_kind_size, 2) + logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2) + integer(bit_kind) :: det(N_int, 2) + double precision, intent(in) :: coefs(N_states) + double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) + double precision :: hij, tmp_row(N_states, mo_tot_num), tmp_row2(N_states, mo_tot_num) + double precision, external :: get_phase_bi, integral8 + + logical :: lbanned(mo_tot_num, 2), ok + integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j, hfix, pfix, h1, h2, p1, p2, ib + + integer, intent(in) :: h(0:2,2), p(0:4,2), sp + + integer, parameter :: turn2(2) = (/2,1/) + integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) + + integer :: bant + + + lbanned = bannedOrb + + do i=1, p(0,1) + lbanned(p(i,1), 1) = .true. + end do + do i=1, p(0,2) + lbanned(p(i,2), 2) = .true. + end do + + ma = 1 + if(p(0,2) >= 2) ma = 2 + mi = turn2(ma) + + bant = 1 + + if(sp == 3) then + !move MA + if(ma == 2) bant = 2 + puti = p(1,mi) + hfix = h(1,ma) + p1 = p(1,ma) + p2 = p(2,ma) + if(.not. bannedOrb(puti, mi)) then + tmp_row = 0d0 + do putj=1, hfix-1 + if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle + hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) + tmp_row(1:N_states,putj) += hij * coefs(1:N_states) + end do + do putj=hfix+1, mo_tot_num + if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle + hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) + tmp_row(1:N_states,putj) += hij * coefs(1:N_states) + end do + + if(ma == 1) then + mat(1:N_states,1:mo_tot_num,puti) += tmp_row(1:N_states,1:mo_tot_num) + else + mat(1:N_states,puti,1:mo_tot_num) += tmp_row(1:N_states,1:mo_tot_num) + end if + end if + + !MOVE MI + pfix = p(1,mi) + tmp_row = 0d0 + tmp_row2 = 0d0 + do puti=1,mo_tot_num + if(lbanned(puti,mi)) cycle + !p1 fixed + putj = p1 + if(.not. banned(putj,puti,bant)) then + hij = integral8(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix) + tmp_row(:,puti) += hij * coefs + end if + + putj = p2 + if(.not. banned(putj,puti,bant)) then + hij = integral8(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix) + tmp_row2(:,puti) += hij * coefs + end if + end do + + if(mi == 1) then + mat(:,:,p1) += tmp_row(:,:) + mat(:,:,p2) += tmp_row2(:,:) + else + mat(:,p1,:) += tmp_row(:,:) + mat(:,p2,:) += tmp_row2(:,:) + end if + else + if(p(0,ma) == 3) then + do i=1,3 + hfix = h(1,ma) + puti = p(i, ma) + p1 = p(turn3(1,i), ma) + p2 = p(turn3(2,i), ma) + tmp_row = 0d0 + do putj=1,hfix-1 + if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle + hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2) + tmp_row(:,putj) += hij * coefs + end do + do putj=hfix+1,mo_tot_num + if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle + hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2) + tmp_row(:,putj) += hij * coefs + end do + + mat(:, :puti-1, puti) += tmp_row(:,:puti-1) + mat(:, puti, puti:) += tmp_row(:,puti:) + end do + else + hfix = h(1,mi) + pfix = p(1,mi) + p1 = p(1,ma) + p2 = p(2,ma) + tmp_row = 0d0 + tmp_row2 = 0d0 + do puti=1,mo_tot_num + if(lbanned(puti,ma)) cycle + putj = p2 + if(.not. banned(puti,putj,1)) then + hij = integral8(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1) + tmp_row(:,puti) += hij * coefs + end if + + putj = p1 + if(.not. banned(puti,putj,1)) then + hij = integral8(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2) + tmp_row2(:,puti) += hij * coefs + end if + end do + mat(:,:p2-1,p2) += tmp_row(:,:p2-1) + mat(:,p2,p2:) += tmp_row(:,p2:) + mat(:,:p1-1,p1) += tmp_row2(:,:p1-1) + mat(:,p1,p1:) += tmp_row2(:,p1:) + end if + end if + + !! MONO + if(sp == 3) then + s1 = 1 + s2 = 2 + else + s1 = sp + s2 = sp + end if + + do i1=1,p(0,s1) + ib = 1 + if(s1 == s2) ib = i1+1 + do i2=ib,p(0,s2) + p1 = p(i1,s1) + p2 = p(i2,s2) + if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle + call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) + call i_h_j(gen, det, N_int, hij) + mat(:, p1, p2) += coefs * hij + end do + end do +end subroutine + + + + +subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2) + integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2) + logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2) + integer(bit_kind) :: det(N_int, 2) + double precision, intent(in) :: coefs(N_states) + double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) + integer, intent(in) :: h(0:2,2), p(0:4,2), sp + + integer :: i, j, s, h1, h2, p1, p2, puti, putj + double precision :: hij, phase + double precision, external :: get_phase_bi, integral8 + logical :: ok + + integer :: bant + bant = 1 + + + if(sp == 3) then ! AB + h1 = p(1,1) + h2 = p(1,2) + do p1=1, mo_tot_num + if(bannedOrb(p1, 1)) cycle + do p2=1, mo_tot_num + if(bannedOrb(p2,2)) cycle + if(banned(p1, p2, bant)) cycle ! rentable? + if(p1 == h1 .or. p2 == h2) then + call apply_particles(mask, 1,p1,2,p2, det, ok, N_int) + call i_h_j(gen, det, N_int, hij) + else + hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) + phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) + end if + mat(:, p1, p2) += coefs(:) * hij + end do + end do + else ! AA BB + p1 = p(1,sp) + p2 = p(2,sp) + do puti=1, mo_tot_num + if(bannedOrb(puti, sp)) cycle + do putj=puti+1, mo_tot_num + if(bannedOrb(putj, sp)) cycle + if(banned(puti, putj, bant)) cycle ! rentable? + if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then + call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int) + call i_h_j(gen, det, N_int, hij) + else + hij = (integral8(p1, p2, puti, putj) - integral8(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2) + end if + mat(:, puti, putj) += coefs(:) * hij + end do + end do + end if +end subroutine + + +subroutine past_d1(bannedOrb, p) + use bitmasks + implicit none + + logical, intent(inout) :: bannedOrb(mo_tot_num, 2) + integer, intent(in) :: p(0:4, 2) + integer :: i,s + + do s = 1, 2 + do i = 1, p(0, s) + bannedOrb(p(i, s), s) = .true. + end do + end do +end subroutine + + +subroutine past_d2(banned, p, sp) + use bitmasks + implicit none + + logical, intent(inout) :: banned(mo_tot_num, mo_tot_num) + integer, intent(in) :: p(0:4, 2), sp + integer :: i,j + + if(sp == 3) then + do i=1,p(0,1) + do j=1,p(0,2) + banned(p(i,1), p(j,2)) = .true. + end do + end do + else + do i=1,p(0, sp) + do j=1,i-1 + banned(p(j,sp), p(i,sp)) = .true. + banned(p(i,sp), p(j,sp)) = .true. + end do + end do + end if +end subroutine + + + +subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting) + use bitmasks + implicit none + + integer, intent(in) :: interesting(0:N) + integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N) + integer, intent(in) :: i_gen, N + logical, intent(inout) :: banned(mo_tot_num, mo_tot_num) + logical, intent(out) :: fullMatch + + + integer :: i, j, na, nb, list(3) + integer(bit_kind) :: myMask(N_int, 2), negMask(N_int, 2) + + fullMatch = .false. + + do i=1,N_int + negMask(i,1) = not(mask(i,1)) + negMask(i,2) = not(mask(i,2)) + end do + + genl : do i=1, N + do j=1, N_int + if(iand(det(j,1,i), mask(j,1)) /= mask(j, 1)) cycle genl + if(iand(det(j,2,i), mask(j,2)) /= mask(j, 2)) cycle genl + end do + + if(interesting(i) < i_gen) then + fullMatch = .true. + return + end if + + do j=1, N_int + myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1)) + myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2)) + end do + + call bitstring_to_list(myMask(1,1), list(1), na, N_int) + call bitstring_to_list(myMask(1,2), list(na+1), nb, N_int) + banned(list(1), list(2)) = .true. + end do genl +end subroutine + diff --git a/plugins/CAS_SD_ZMQ/selection_buffer.irp.f b/plugins/CAS_SD_ZMQ/selection_buffer.irp.f new file mode 100644 index 00000000..2bcb11d3 --- /dev/null +++ b/plugins/CAS_SD_ZMQ/selection_buffer.irp.f @@ -0,0 +1,70 @@ + +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 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(dabs(val) >= b%mini) then + b%cur += 1 + b%det(:,:,b%cur) = det(:,:) + b%val(b%cur) = val + if(b%cur == size(b%val)) then + call sort_selection_buffer(b) + end if + end if +end subroutine + + +subroutine sort_selection_buffer(b) + use selection_types + implicit none + + type(selection_buffer), intent(inout) :: b + double precision, allocatable :: vals(:), absval(:) + integer, allocatable :: iorder(:) + integer(bit_kind), allocatable :: detmp(:,:,:) + integer :: i, nmwen + logical, external :: detEq + nmwen = min(b%N, b%cur) + + + allocate(iorder(b%cur), detmp(N_int, 2, nmwen), absval(b%cur), vals(nmwen)) + absval = -dabs(b%val(:b%cur)) + do i=1,b%cur + iorder(i) = i + end do + call dsort(absval, iorder, b%cur) + + do i=1, nmwen + detmp(:,:,i) = b%det(:,:,iorder(i)) + vals(i) = b%val(iorder(i)) + end do + b%det(:,:,:nmwen) = detmp(:,:,:) + b%det(:,:,nmwen+1:) = 0_bit_kind + b%val(:nmwen) = vals(:) + b%val(nmwen+1:) = 0d0 + b%mini = max(b%mini,dabs(b%val(b%N))) + b%cur = nmwen +end subroutine + diff --git a/plugins/CAS_SD_ZMQ/selection_types.f90 b/plugins/CAS_SD_ZMQ/selection_types.f90 new file mode 100644 index 00000000..9506629c --- /dev/null +++ b/plugins/CAS_SD_ZMQ/selection_types.f90 @@ -0,0 +1,9 @@ +module selection_types + type selection_buffer + integer :: N, cur + integer(8), allocatable :: det(:,:,:) + double precision, allocatable :: val(:) + double precision :: mini + endtype +end module + diff --git a/plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES b/plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES index cb6ff46e..7ff203d4 100644 --- a/plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES +++ b/plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_full ZMQ Full_CI +Perturbation Selectors_full Generators_full ZMQ diff --git a/plugins/Full_CI_ZMQ/fci_zmq.irp.f b/plugins/Full_CI_ZMQ/fci_zmq.irp.f index e03db458..8b9488d2 100644 --- a/plugins/Full_CI_ZMQ/fci_zmq.irp.f +++ b/plugins/Full_CI_ZMQ/fci_zmq.irp.f @@ -79,7 +79,7 @@ program fci_zmq enddo endif E_CI_before(1:N_states) = CI_energy(1:N_states) - call ezfio_set_full_ci_energy(CI_energy) + call ezfio_set_full_ci_zmq_energy(CI_energy) enddo if(do_pt2_end)then @@ -99,7 +99,7 @@ program fci_zmq print *, 'E+PT2 = ', E_CI_before+pt2 print *, '-----' enddo - call ezfio_set_full_ci_energy_pt2(E_CI_before+pt2) + call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before+pt2) endif call save_wavefunction end diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index d8a18437..68a7a050 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -350,8 +350,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] integral = ao_bielec_integral(1,1,1,1) real :: map_mb - print*, 'read_ao_integrals',read_ao_integrals - print*, 'disk_access_ao_integrals',disk_access_ao_integrals + PROVIDE read_ao_integrals disk_access_ao_integrals if (read_ao_integrals) then print*,'Reading the AO integrals' call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) diff --git a/tests/bats/cassd.bats b/tests/bats/cassd.bats index 07d79f1a..44b44ee6 100644 --- a/tests/bats/cassd.bats +++ b/tests/bats/cassd.bats @@ -3,15 +3,15 @@ source $QP_ROOT/tests/bats/common.bats.sh @test "CAS_SD H2O cc-pVDZ" { - test_exe cas_sd_selected || skip + test_exe cassd_zmq || skip INPUT=h2o.ezfio qp_edit -c $INPUT ezfio set_file $INPUT ezfio set perturbation do_pt2_end False - ezfio set determinants n_det_max 1000 + ezfio set determinants n_det_max 2000 qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]" - qp_run cas_sd_selected $INPUT - energy="$(ezfio get cas_sd energy)" + qp_run cassd_zmq $INPUT + energy="$(ezfio get cas_sd_zmq energy)" eq $energy -76.2221842108163 1.E-5 } diff --git a/tests/bats/fci.bats b/tests/bats/fci.bats index 174c8f61..79ff91ab 100644 --- a/tests/bats/fci.bats +++ b/tests/bats/fci.bats @@ -20,7 +20,7 @@ function run_FCI() { function run_FCI_ZMQ() { thresh=5.e-5 - test_exe full_ci || skip + test_exe fci_zmq || skip qp_edit -c $1 ezfio set_file $1 ezfio set perturbation do_pt2_end True @@ -28,9 +28,9 @@ function run_FCI_ZMQ() { ezfio set davidson threshold_davidson 1.e-10 qp_run fci_zmq $1 - energy="$(ezfio get full_ci energy)" + energy="$(ezfio get full_ci_zmq energy)" eq $energy $3 $thresh - energy_pt2="$(ezfio get full_ci energy_pt2)" + energy_pt2="$(ezfio get full_ci_zmq energy_pt2)" eq $energy_pt2 $4 $thresh } diff --git a/tests/bats/pseudo.bats b/tests/bats/pseudo.bats index 8cccf229..a20b0842 100644 --- a/tests/bats/pseudo.bats +++ b/tests/bats/pseudo.bats @@ -23,7 +23,7 @@ function run_HF() { function run_FCI_ZMQ() { thresh=5.e-5 - test_exe full_ci || skip + test_exe fci_zmq|| skip qp_edit -c $1 ezfio set_file $1 ezfio set perturbation do_pt2_end True @@ -31,9 +31,9 @@ function run_FCI_ZMQ() { ezfio set davidson threshold_davidson 1.e-10 qp_run fci_zmq $1 - energy="$(ezfio get full_ci energy)" + energy="$(ezfio get full_ci_zmq energy)" eq $energy $3 $thresh - energy_pt2="$(ezfio get full_ci energy_pt2)" + energy_pt2="$(ezfio get full_ci_zmq energy_pt2)" eq $energy_pt2 $4 $thresh }