From a117ac1a6bee5f9c7c329b09d408ce48f6135b2b Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Mon, 25 Jul 2016 14:10:28 +0200 Subject: [PATCH] ready for integral driven version - dirty copypastas --- plugins/Full_CI_ZMQ/fci_zmq.irp.f | 25 +- plugins/Full_CI_ZMQ/selection.irp.f | 920 ++++++++++++++++------------ 2 files changed, 531 insertions(+), 414 deletions(-) diff --git a/plugins/Full_CI_ZMQ/fci_zmq.irp.f b/plugins/Full_CI_ZMQ/fci_zmq.irp.f index e221533b..b3a89b74 100644 --- a/plugins/Full_CI_ZMQ/fci_zmq.irp.f +++ b/plugins/Full_CI_ZMQ/fci_zmq.irp.f @@ -7,8 +7,8 @@ program fci_zmq double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) integer :: N_st, degree - integer :: it, mit(0:5) - mit = (/1, 268, 1517, 10018, 45096, 100000/) + integer :: it, mit(0:6) + mit = (/1, 246, 1600, 17528, 112067, 519459, 2685970/) it = 0 N_st = N_states allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) @@ -41,18 +41,18 @@ program fci_zmq E_CI_before = CI_energy do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) n_det_before = N_det -! call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) + ! call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) it += 1 - if(it > 5) stop - call ZMQ_selection(mit(it) - mit(it-1), pt2) !(max(N_det*3, 1000-N_det)) + if(it > 6) stop + call ZMQ_selection(mit(it) - mit(it-1), pt2) ! max(1000-N_det, N_det), pt2) - do i=1, N_det + !do i=1, N_det !if(popcnt(psi_det(1,1,i)) + popcnt(psi_det(2,1,i)) /= 23) stop "ZZ1" -2099.2504682049275 !if(popcnt(psi_det(1,2,i)) + popcnt(psi_det(2,2,i)) /= 23) stop "ZZ2" - do k=1,i-1 - if(detEq(psi_det(1,1,i), psi_det(1,1,k), N_int)) stop "ATRRGRZER" - end do - end do + ! do k=1,i-1 + ! if(detEq(psi_det(1,1,i), psi_det(1,1,k), N_int)) stop "ATRRGRZER" + ! end do + !end do PROVIDE psi_coef PROVIDE psi_det PROVIDE psi_det_sorted @@ -128,11 +128,10 @@ subroutine ZMQ_selection(N, pt2) integer :: i integer, external :: omp_get_thread_num double precision, intent(out) :: pt2(N_states) - call flip_generators() + !call flip_generators() call new_parallel_job(zmq_to_qp_run_socket,'selection') call create_selection_buffer(N, N*2, b) - do i= N_det_generators, 1, -1 write(task,*) i, N call add_task_to_taskserver(zmq_to_qp_run_socket,task) @@ -149,7 +148,7 @@ subroutine ZMQ_selection(N, pt2) endif !$OMP END PARALLEL call end_parallel_job(zmq_to_qp_run_socket, 'selection') - call flip_generators() + !call flip_generators() call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN call copy_H_apply_buffer_to_wf() end subroutine diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index 5e4eef36..d0c7cbba 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/plugins/Full_CI_ZMQ/selection.irp.f @@ -4,31 +4,32 @@ subroutine selection_slave(thread,iproc) use f77_zmq use selection_types implicit none - + integer, intent(in) :: thread, iproc integer :: rc, i integer :: worker_id, task_id(100), 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 - + 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 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) - + buf%N = 0 ctask = 1 pt2 = 0d0 - do + + do call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task) done = task_id(ctask) == 0 if (.not. done) then @@ -41,30 +42,30 @@ subroutine selection_slave(thread,iproc) end if call select_connected(i_generator,ci_electronic_energy,pt2,buf) !! ci_electronic_energy ?? end if - + if(done) ctask = ctask - 1 - - if(done .or. ctask == size(task_id)) then - if(ctask > 0 .and. buf%N /= 0) then + + 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) pt2 = 0d0 buf%cur = 0 end if - - do i=1, ctask - call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i)) - end do + + 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 @@ -72,27 +73,32 @@ 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 + + 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) - - rc = f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE) - - rc = f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE) - - rc = f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE) - - rc = f77_zmq_send( zmq_socket_push, task_id(1), ntask*4, 0) - + 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" + end subroutine @@ -100,24 +106,30 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, val, det, N, task_id, nt use f77_zmq use selection_types implicit none - integer(ZMQ_PTR), intent(in) :: zmq_socket_pull + 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, ZMQ_SNDMORE) - - rc = f77_zmq_recv( zmq_socket_pull, pt2, N_states*8, ZMQ_SNDMORE) - - rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, ZMQ_SNDMORE) - - rc = f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, ZMQ_SNDMORE) - - rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, ZMQ_SNDMORE) - - rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0) + + rc = f77_zmq_recv( zmq_socket_pull, N, 4, ZMQ_SNDMORE) + if(rc /= 4) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, pt2, N_states*8, ZMQ_SNDMORE) + if(rc /= 8*N_states) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, ZMQ_SNDMORE) + if(rc /= 8*N) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, ZMQ_SNDMORE) + if(rc /= bit_kind*N_int*2*N) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, ZMQ_SNDMORE) + if(rc /= 4) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0) + if(rc /= 4*ntask) stop "pull" end subroutine @@ -134,37 +146,37 @@ subroutine select_connected(i_generator,E0,pt2,b) 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)) + 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)) ) - + hole_mask(k,1) = ior(generators_bitmask(k,1,s_hole,l), generators_bitmask(k,1,s_part,l)) hole_mask(k,2) = ior(generators_bitmask(k,2,s_hole,l), generators_bitmask(k,2,s_part,l)) particle_mask(k,:) = hole_mask(k,:) enddo - + call select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b) call select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b) enddo end - + 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 @@ -176,12 +188,12 @@ 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(:,:) @@ -196,7 +208,7 @@ 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(:) @@ -204,15 +216,15 @@ subroutine sort_selection_buffer(b) 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)) @@ -231,47 +243,53 @@ subroutine selection_collector(b, pt2) 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(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(:) - - zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() - zmq_socket_pull = new_zmq_pull_socket() + 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)) - pt2 = 0d0 - more = 1 - do while (more == 1) + 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 - call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more) - endif + if(task_id(i) == 0) stop "collector" + 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 - + subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf) use f77_zmq @@ -287,7 +305,7 @@ subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p integer(bit_kind), intent(in) :: hole_mask(N_int,2), particle_mask(N_int,2) double precision, intent(in) :: E0(N_states) type(selection_buffer), intent(inout) :: buf - + integer :: i,j,k,l integer :: msg_size @@ -313,7 +331,7 @@ subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p call bitstring_to_list_ab(hole , hole_list , N_holes , N_int) call bitstring_to_list_ab(particle, particle_list, N_particles, N_int) - + ! Create excited determinants ! --------------------------- @@ -326,26 +344,18 @@ subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p ion_det(k,1) = psi_det_generators(k,1,i_generator) ion_det(k,2) = psi_det_generators(k,2,i_generator) enddo - - - ! Create the mini wave function where = - ! -------------------------------------------------------------- -! integer(bit_kind) :: psi_det_connected(N_int,2,psi_selectors_size) -! double precision :: psi_coef_connected(psi_selectors_size,N_states) - + + integer :: ptr_microlist(0:mo_tot_num * 2 + 1), N_microlist(0:mo_tot_num * 2) integer, allocatable :: idx_microlist(:) integer(bit_kind), allocatable :: microlist(:, :, :) double precision, allocatable :: psi_coef_microlist(:,:) - + allocate(microlist(N_int, 2, N_det_selectors * 3), psi_coef_microlist(psi_selectors_size * 3, N_states), idx_microlist(N_det_selectors * 3)) - + do ispin=1,2 -! do k=1,N_int -! ion_det(k,ispin) = psi_det_generators(k,ispin,i_generator) -! enddo - + do i=1, N_holes(ispin) ion_det(:,:) = psi_det_generators(:,:,i_generator) @@ -356,81 +366,64 @@ subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p integer :: j_hole, k_hole k_hole = ishft(i_hole-1,-bit_kind_shift)+1 ! N_int j_hole = i_hole-ishft(k_hole-1,bit_kind_shift)-1 ! bit index -! ion_det(k_hole,ispin) = ibclr(psi_det_generators(k_hole,ispin,i_generator),j_hole) ion_det(k_hole,ispin) = ibclr(ion_det(k_hole,ispin),j_hole) - + call create_microlist_single(psi_selectors, i_generator, N_det_selectors, ion_det, microlist, idx_microlist, N_microlist, ptr_microlist, N_int) - + do j=1, ptr_microlist(mo_tot_num * 2 + 1) - 1 psi_coef_microlist(j,:) = psi_selectors_coef(idx_microlist(j),:) enddo - + if(ptr_microlist(mo_tot_num * 2 + 1) == 1) then cycle endif - ! Create particles - ! ---------------- do j=1,N_particles(ispin) -! exc_det(k_hole,ispin) = ion_det(k_hole,ispin) exc_det(:,:) = ion_det(:,:) - + integer :: i_particle i_particle = particle_list(j,ispin) - ! Apply the particle integer :: j_particle, k_particle k_particle = ishft(i_particle-1,-bit_kind_shift)+1 ! N_int j_particle = i_particle-ishft(k_particle-1,bit_kind_shift)-1 ! bit index -! exc_det(k_particle,ispin) = ibset(ion_det(k_particle,ispin),j_particle) exc_det(k_particle,ispin) = ibset(exc_det(k_particle,ispin),j_particle) - - ! TODO - + + logical, external :: is_in_wavefunction logical :: nok - ! TODO : Check connected to ref if (.not. is_in_wavefunction(exc_det,N_int)) then - ! Compute perturbative contribution and select determinant double precision :: i_H_psi_value(N_states), i_H_psi_value2(N_states) - double precision :: i_H_full(N_states) i_H_psi_value = 0d0 i_H_psi_value2 = 0d0 - i_H_full = 0d0 integer :: sporb - -! call i_H_psi(exc_det,psi_selectors,psi_selectors_coef,N_int,N_det_selectors,psi_selectors_size,N_states,i_H_full) -! + + nok = .false. - sporb = i_particle + (ispin - 1) * mo_tot_num -! ! ! subroutine check_past(det, list, idx, N, cur, ok, Nint) - if(N_microlist(sporb) > 0) call check_past(exc_det, microlist(1,1,ptr_microlist(sporb)), idx_microlist(ptr_microlist(sporb)), N_microlist(sporb), i_generator,nok, N_int) + sporb = i_particle + (ispin - 1) * mo_tot_num + + if(N_microlist(sporb) > 0) call check_past(exc_det, microlist(1,1,ptr_microlist(sporb)), idx_microlist(ptr_microlist(sporb)), N_microlist(sporb), i_generator, nok, N_int) if(nok) cycle -! + if(N_microlist(0) > 0) call i_H_psi(exc_det,microlist,psi_coef_microlist,N_int,N_microlist(0),psi_selectors_size*3,N_states,i_H_psi_value) if(N_microlist(sporb) > 0) call i_H_psi(exc_det,microlist(1,1,ptr_microlist(sporb)),psi_coef_microlist(ptr_microlist(sporb), 1),N_int,N_microlist(sporb),psi_selectors_size*3,N_states,i_H_psi_value2) - - i_H_psi_value(:) = i_H_psi_value(:) + i_H_psi_value2(:) double precision :: Hii, diag_H_mat_elem_fock Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),exc_det,fock_diag_tmp,N_int) - + double precision :: delta_E, e_pert(N_states), e_pertm e_pert(:) = 0d0 e_pertm = 0d0 - + do k=1,N_states -! if(dabs(1d0 - i_H_psi_value(k)/i_H_full(k)) > 1d-6) then -! stop "PAS BON, PAS BOOOOON!! (single)" -! endif if (i_H_psi_value(k) == 0.d0) cycle delta_E = E0(k) - Hii if (delta_E < 0.d0) then - e_pert(k) = 0.5d0 * (-dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E) + e_pert(k) = 0.5d0 * (-dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E) else - e_pert(k) = 0.5d0 * ( dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E) + e_pert(k) = 0.5d0 * ( dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E) endif if(dabs(e_pert(k)) > dabs(e_pertm)) e_pertm = e_pert(k) pt2(k) += e_pert(k) @@ -466,11 +459,7 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p type(selection_buffer), intent(inout) :: buf logical :: isinwf(mo_tot_num*2, mo_tot_num*2) double precision :: d0s(mo_tot_num, mo_tot_num, N_states) - d0s = 0d0 -! double precision, save :: d0 = 0d0 -! double precision, save :: d1 = 0d0 -! double precision, save :: d2 = 0d0 - + integer :: i,j,k,l,j1,j2,i1,i2,ib,jb integer :: msg_size @@ -496,31 +485,31 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p call bitstring_to_list_ab(hole , hole_list , N_holes , N_int) call bitstring_to_list_ab(particle, particle_list, N_particles, N_int) - + ! Create excited determinants ! --------------------------- integer :: ispin1, ispin2, other_spin integer(bit_kind) :: exc_det(N_int,2), ion_det(N_int,2) - - + + integer :: ptr_microlist(0:mo_tot_num * 2 + 1), N_microlist(0:mo_tot_num * 2) double precision, allocatable :: psi_coef_microlist(:,:) - + integer :: ptr_tmicrolist(0:mo_tot_num * 2 + 1), N_tmicrolist(0:mo_tot_num * 2) double precision, allocatable :: psi_coef_tmicrolist(:,:) - + integer, allocatable :: idx_tmicrolist(:), idx_microlist(:) integer(bit_kind), allocatable :: microlist(:,:,:), tmicrolist(:,:,:) - + integer :: ptr_futur_microlist(0:mo_tot_num * 2 + 1), ptr_futur_tmicrolist(0:mo_tot_num * 2 + 1) integer :: N_futur_microlist(0:mo_tot_num * 2), N_futur_tmicrolist(0:mo_tot_num * 2) logical :: pastlink - + allocate(idx_tmicrolist(N_det_selectors * 3), idx_microlist(N_det_selectors * 4)) allocate(microlist(N_int, 2, N_det_selectors * 4), tmicrolist(N_int, 2, N_det_selectors * 3)) allocate(psi_coef_tmicrolist(psi_selectors_size * 3, N_states), psi_coef_microlist(psi_selectors_size * 4, N_states)) - + do k=1,N_int exc_det(k,1) = psi_det_generators(k,1,i_generator) exc_det(k,2) = psi_det_generators(k,2,i_generator) @@ -536,62 +525,62 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p if(ispin1 == ispin2) ib = i1+1 do i2=ib, N_holes(ispin2) ion_det(:,:) = psi_det_generators(:,:,i_generator) -! call set_hole(ion_det, hole_list(i1,ispin1), ispin1, hole_list(i1,ispin1), ispin1, Nint) + i_hole1 = hole_list(i1,ispin1) k_hole = ishft(i_hole1-1,-bit_kind_shift)+1 ! N_int j_hole = i_hole1-ishft(k_hole-1,bit_kind_shift)-1 ! bit index ion_det(k_hole,ispin1) = ibclr(ion_det(k_hole,ispin1),j_hole) - + i_hole2 = hole_list(i2,ispin2) k_hole = ishft(i_hole2-1,-bit_kind_shift)+1 ! N_int j_hole = i_hole2-ishft(k_hole-1,bit_kind_shift)-1 ! bit index ion_det(k_hole,ispin2) = ibclr(ion_det(k_hole,ispin2),j_hole) - + call create_microlist_double(psi_selectors, i_generator, N_det_selectors, ion_det, & microlist, idx_microlist, N_microlist, ptr_microlist, & tmicrolist, idx_tmicrolist, N_tmicrolist, ptr_tmicrolist, & isinwf, d0s, N_int) - if(N_microlist(0) > 0 .and. idx_microlist(1) > i_generator) stop "wtf..." + if(ptr_microlist(mo_tot_num * 2 + 1) == 1 .and. ptr_tmicrolist(mo_tot_num * 2 + 1) == 1) cycle - + call finish_isinwf(ion_det, psi_det_sorted(1,1,N_det_selectors+1), N_det - N_det_selectors, isinwf) - - + + call create_futur_ptr(ptr_microlist, idx_microlist, ptr_futur_microlist, N_futur_microlist, i_generator) call create_futur_ptr(ptr_tmicrolist, idx_tmicrolist, ptr_futur_tmicrolist, N_futur_tmicrolist, i_generator) - - + + do j=1, ptr_microlist(mo_tot_num * 2 + 1) - 1 psi_coef_microlist(j,:) = psi_selectors_coef(idx_microlist(j),:) enddo do j=1, ptr_tmicrolist(mo_tot_num * 2 + 1) - 1 psi_coef_tmicrolist(j,:) = psi_selectors_coef(idx_tmicrolist(j),:) enddo - - + + ! Create particles ! ---------------- integer :: i_particle1, i_particle2, k_particle, j_particle integer :: p1, p2, sporb, lorb - + do j1=1,N_particles(ispin1) i_particle1 = particle_list(j1, ispin1) p1 = i_particle1 + (ispin1 - 1) * mo_tot_num - if(N_tmicrolist(p1) > 0 .and. idx_tmicrolist(ptr_tmicrolist(p1+1)-1) > i_generator) cycle + if(N_tmicrolist(p1) > 0 .and. idx_tmicrolist(ptr_tmicrolist(p1)) < i_generator) cycle jb = 1 if(ispin1 == ispin2) jb = j1+1 do j2=jb,N_particles(ispin2) - + i_particle2 = particle_list(j2, ispin2) - - - + + + p2 = i_particle2 + (ispin2 - 1) * mo_tot_num - if(N_tmicrolist(p2) > 0 .and. idx_tmicrolist(ptr_tmicrolist(p2+1)-1) > i_generator) cycle + if(N_tmicrolist(p2) > 0 .and. idx_tmicrolist(ptr_tmicrolist(p2)) < i_generator) cycle if(isinwf(p1, p2)) cycle exc_det = ion_det - - + + if(N_microlist(p1) < N_microlist(p2)) then sporb = p1 lorb = p2 @@ -599,75 +588,40 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p sporb = p2 lorb = p1 endif - - + + ! Apply the particle k_particle = ishft(i_particle2-1,-bit_kind_shift)+1 ! N_int j_particle = i_particle2-ishft(k_particle-1,bit_kind_shift)-1 ! bit index exc_det(k_particle,ispin2) = ibset(exc_det(k_particle,ispin2),j_particle) - + ! Apply the particle k_particle = ishft(i_particle1-1,-bit_kind_shift)+1 ! N_int j_particle = i_particle1-ishft(k_particle-1,bit_kind_shift)-1 ! bit index exc_det(k_particle,ispin1) = ibset(exc_det(k_particle,ispin1),j_particle) - - -! if(.false. .or. (is_in_wavefunction(exc_det,N_int) .and. .not. isinwf(p1,p2))) then -! print *, p1, p2 -! call debug_det(ion_det, N_int) -! call debug_det(exc_det, N_int) -! ! do i=1,mo_tot_num*2 -! ! print *, isinwf(:, i) -! ! end do -! print *, isinwf(p1, p2) -! stop "isw" -! end if - - ! TODO - + logical, external :: is_in_wavefunction logical :: nok - ! TODO : Check connected to ref -! if (.not. is_in_wavefunction(exc_det,N_int)) then + ! Compute perturbative contribution and select determinant + double precision :: i_H_psi_value(N_states), i_H_psi_value2(N_states) + i_H_psi_value = 0d0 + i_H_psi_value2 = 0d0 + + nok = .false. + call check_past_s(exc_det, microlist(1,1,ptr_microlist(sporb)), N_microlist(sporb) - N_futur_microlist(sporb), nok, N_int) + if(nok) cycle + !DET DRIVEN + if(N_futur_microlist(0) > 0) then + call i_H_psi(exc_det,microlist(1,1,ptr_futur_microlist(0)),psi_coef_microlist(ptr_futur_microlist(0), 1),N_int,N_futur_microlist(0),psi_selectors_size*4,N_states,i_H_psi_value) + end if + !INTEGRAL DRIVEN +! i_H_psi_value = d0s(mod(p1-1, mo_tot_num)+1, mod(p2-1, mo_tot_num)+1, :) + - - ! Compute perturbative contribution and select determinant - double precision :: i_H_psi_value(N_states), i_H_psi_value2(N_states) - double precision :: i_H_full(N_states) - i_H_psi_value = 0d0 - i_H_psi_value2 = 0d0 - i_H_full = 0d0 - -! call i_H_psi(exc_det,psi_selectors,psi_selectors_coef,N_int,N_det_selectors,psi_selectors_size,N_states,i_H_full) - - - -! call check_past(exc_det, microlist, idx_microlist, N_microlist(0), i_generator, nok, N_int) -! if(nok) cycle - - nok = .false. - - call check_futur(exc_det, microlist(1,1,ptr_futur_microlist(sporb)), N_futur_microlist(sporb), nok, N_int) - - if(nok) cycle - - if(N_microlist(0)-N_futur_microlist(0) > 0) then - call i_H_psi(exc_det,microlist(1,1,ptr_microlist(0)),psi_coef_microlist(ptr_microlist(0), 1),N_int,N_microlist(0)-N_futur_microlist(0),psi_selectors_size*4,N_states,i_H_psi_value) -! if(i_H_psi_value(1) /= d0s(p1, p2, 1) .and. d0s(p1, p2, 1) /= 0d0) then -! print *, d0s(p1, p2, 1), i_H_psi_value(1) -! print *, d0s(:3, :3, 1) -! stop "SKSL" -! end if - end if -! if(N_futur_microlist(sporb) > 0) call i_H_psi(exc_det,microlist(1,1,ptr_futur_microlist(sporb)),psi_coef_microlist(ptr_futur_microlist(sporb), 1),N_int,N_futur_microlist(sporb),psi_selectors_size*4,N_states,i_H_psi_value2) -! !$OMP ATOMIC -! d0 += dabs(i_H_psi_value(1)) -! d2 += N_futur_microlist(sporb) - - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - if(N_microlist(sporb)-N_futur_microlist(sporb) > 0) then -! ! if(dfloat(N_futur_microlist(lorb)) / dfloat(N_futur_microlist(sporb)) < 2d0) then + if(N_futur_microlist(sporb) > 0) then + !!! COMPUTE INTERSECTION + !!!!!!!!!!!!! +! if(dfloat(N_futur_microlist(lorb)) / dfloat(N_futur_microlist(sporb)) < 2d0) then ! c1 = ptr_futur_microlist(p1) ! c2 = ptr_futur_microlist(p2) ! do while(c1 < ptr_microlist(p1+1) .and. c2 < ptr_microlist(p2+1)) @@ -685,112 +639,72 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p ! endif ! end do ! else - call i_H_psi(exc_det,microlist(1,1,ptr_microlist(sporb)),psi_coef_microlist(ptr_microlist(sporb), 1),N_int,N_microlist(sporb)-N_futur_microlist(sporb),psi_selectors_size*4,N_states,i_H_psi_value2) - i_H_psi_value = i_H_psi_value + i_H_psi_value2 -! !$OMP ATOMIC -! d2 += dabs(i_H_psi_value2(1)) -! d2 += N_futur_microlist(sporb) -! ! end if - end if - - -! enddo 2099.3283623955813 - - - !!!!!!!!!!!!!!!!!! -! if(N_microlist(0) > 0) call i_H_psi(exc_det,microlist,psi_coef_microlist(ptr_microlist(0), 1),N_int,N_microlist(0),psi_selectors_size*4,N_states,i_H_psi_value) -! if(N_microlist(sporb) > 0) call i_H_psi(exc_det,microlist(1,1,ptr_microlist(sporb)),psi_coef_microlist(ptr_microlist(sporb), 1),N_int,N_microlist(sporb),psi_selectors_size*4,N_states,i_H_psi_value2) - -! i_H_psi_value = i_H_psi_value + i_H_psi_value2 - - integer :: c1, c2 - double precision :: hij - c1 = ptr_tmicrolist(p1) - c2 = ptr_tmicrolist(p2) - do while(.true.) - if(c1 >= ptr_futur_tmicrolist(p1) .or. c2 >= ptr_futur_tmicrolist(p2)) then - if(ptr_futur_tmicrolist(p1) /= c1) then - call i_H_psi(exc_det,tmicrolist(1,1,c1),psi_coef_tmicrolist(c1, 1),N_int, ptr_futur_tmicrolist(p1)-c1 ,psi_selectors_size*3,N_states,i_H_psi_value2) - i_H_psi_value = i_H_psi_value + i_H_psi_value2 -! ! !$OMP ATOMIC -! d1 += dabs(i_H_psi_value2(1)) - end if - - if(ptr_futur_tmicrolist(p2) /= c2) then - call i_H_psi(exc_det,tmicrolist(1,1,c2),psi_coef_tmicrolist(c2, 1),N_int, ptr_futur_tmicrolist(p2)-c2 ,psi_selectors_size*3,N_states,i_H_psi_value2) - i_H_psi_value = i_H_psi_value + i_H_psi_value2 -! !$OMP ATOMIC -! d1 += dabs(i_H_psi_value2(1)) - endif - - exit - endif - - if(idx_tmicrolist(c1) < idx_tmicrolist(c2)) then - call i_H_j(exc_det,tmicrolist(1,1,c1),N_int,hij) - do j = 1, N_states - i_H_psi_value(j) = i_H_psi_value(j) + psi_coef_tmicrolist(c1,j)*hij - enddo -! !$OMP ATOMIC -! d1 += dabs(psi_coef_tmicrolist(c1,1)*hij) - c1 += 1 - else - call i_H_j(exc_det,tmicrolist(1,1,c2),N_int,hij) - do j = 1, N_states - i_H_psi_value(j) = i_H_psi_value(j) + psi_coef_tmicrolist(c2,j)*hij - enddo - if(idx_tmicrolist(c1) == idx_tmicrolist(c2)) c1 = c1 + 1 -! !$OMP ATOMIC -! d1 += dabs(psi_coef_tmicrolist(c2,1)*hij) - c2 += 1 - end if - enddo - - double precision :: Hii, diag_H_mat_elem_fock - Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),exc_det,fock_diag_tmp,N_int) - double precision :: delta_E, e_pert(N_states), e_pertm - e_pert(:) = 0d0 - e_pertm = 0d0 - - do k=1,N_states -! if(dabs(1d0 - i_H_psi_value(k)/i_H_full(k)) > 1d-6) then -! print *, i_H_psi_value(k), i_H_full(k), i_H_psi_value(k)/i_H_full(k) -! stop "PAS BON, PAS BOOON (double)" -! -! endif - - if (i_H_psi_value(k) == 0.d0) cycle - delta_E = E0(k) - Hii - if (delta_E < 0.d0) then - e_pert(k) = 0.5d0 * (-dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E) - else - e_pert(k) = 0.5d0 * ( dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E) - endif - if(dabs(e_pert(k)) > dabs(e_pertm)) e_pertm = e_pert(k) - pt2(k) += e_pert(k) - enddo - if(dabs(e_pertm) > dabs(buf%mini)) then - call add_to_selection_buffer(buf, exc_det, e_pertm) - end if -! endif ! iwf + call i_H_psi(exc_det,microlist(1,1,ptr_futur_microlist(sporb)),psi_coef_microlist(ptr_futur_microlist(sporb), 1),N_int,N_futur_microlist(sporb),psi_selectors_size*4,N_states,i_H_psi_value2) + i_H_psi_value = i_H_psi_value + i_H_psi_value2 + end if - - - - - - - ! Reset exc_det -! exc_det(k_particle,ispin) = psi_det_generators(k_particle,ispin,i_generator) - enddo ! j + + integer :: c1, c2 + double precision :: hij + c1 = ptr_futur_tmicrolist(p1) + c2 = ptr_futur_tmicrolist(p2) + do while(.true.) + if(c1 >= ptr_tmicrolist(p1+1) .or. c2 >= ptr_tmicrolist(p2+1)) then + if(ptr_tmicrolist(p1+1) /= c1) then + call i_H_psi(exc_det,tmicrolist(1,1,c1),psi_coef_tmicrolist(c1, 1),N_int, ptr_tmicrolist(p1+1)-c1 ,psi_selectors_size*3,N_states,i_H_psi_value2) + i_H_psi_value = i_H_psi_value + i_H_psi_value2 + end if + + if(ptr_tmicrolist(p2+1) /= c2) then + call i_H_psi(exc_det,tmicrolist(1,1,c2),psi_coef_tmicrolist(c2, 1),N_int, ptr_tmicrolist(p2+1)-c2 ,psi_selectors_size*3,N_states,i_H_psi_value2) + i_H_psi_value = i_H_psi_value + i_H_psi_value2 + endif + + exit + endif + + if(idx_tmicrolist(c1) < idx_tmicrolist(c2)) then + call i_H_j(exc_det,tmicrolist(1,1,c1),N_int,hij) + do j = 1, N_states + i_H_psi_value(j) = i_H_psi_value(j) + psi_coef_tmicrolist(c1,j)*hij + enddo + c1 += 1 + else + call i_H_j(exc_det,tmicrolist(1,1,c2),N_int,hij) + do j = 1, N_states + i_H_psi_value(j) = i_H_psi_value(j) + psi_coef_tmicrolist(c2,j)*hij + enddo + if(idx_tmicrolist(c1) == idx_tmicrolist(c2)) c1 = c1 + 1 + c2 += 1 + end if + enddo + + double precision :: Hii, diag_H_mat_elem_fock + Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),exc_det,fock_diag_tmp,N_int) + double precision :: delta_E, e_pert(N_states), e_pertm + e_pert(:) = 0d0 + e_pertm = 0d0 + + do k=1,N_states + if (i_H_psi_value(k) == 0.d0) cycle + delta_E = E0(k) - Hii + if (delta_E < 0.d0) then + e_pert(k) = 0.5d0 * (-dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E) + else + e_pert(k) = 0.5d0 * ( dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E) + endif + if(dabs(e_pert(k)) > dabs(e_pertm)) e_pertm = e_pert(k) + pt2(k) += e_pert(k) + enddo + if(dabs(e_pertm) > dabs(buf%mini)) then + call add_to_selection_buffer(buf, exc_det, e_pertm) + end if + enddo enddo - ! Reset ion_det -! ion_det(k_hole,ispin) = psi_det_generators(k_hole,ispin,i_generator) - enddo ! i enddo - enddo ! ispin + enddo + enddo enddo - !print *, "D ::: ", d0/1000000, d1/1000000, d2/1000000 end @@ -799,12 +713,12 @@ subroutine create_futur_ptr(ptr_microlist, idx_microlist, ptr_futur_microlist, N integer, intent(in) :: ptr_microlist(0:mo_tot_num * 2 + 1), idx_microlist(*), i_generator integer, intent(out) :: ptr_futur_microlist(0:mo_tot_num * 2 + 1), N_futur_microlist(0:mo_tot_num * 2) integer :: i, j - + N_futur_microlist = 0 do i=0,mo_tot_num*2 ptr_futur_microlist(i) = ptr_microlist(i+1) do j=ptr_microlist(i), ptr_microlist(i+1) - 1 - if(idx_microlist(j) > i_generator) then + if(idx_microlist(j) >= i_generator) then ptr_futur_microlist(i) = j N_futur_microlist(i) = ptr_microlist(i+1) - j exit @@ -818,26 +732,26 @@ subroutine create_microlist_single(minilist, i_cur, N_minilist, key_mask, microl use bitmasks integer, intent(in) :: Nint, i_cur, N_minilist integer(bit_kind), intent(in) :: minilist(Nint,2,N_minilist), key_mask(Nint,2) - + integer, intent(out) :: N_microlist(0:mo_tot_num*2), ptr_microlist(0:mo_tot_num*2+1), idx_microlist(N_minilist*4) integer(bit_kind), intent(out) :: microlist(Nint,2,N_minilist*4) - + integer :: i,j,k,s,nt,n_element(2) integer :: list(Nint*bit_kind_size,2), cur_microlist(0:mo_tot_num*2+1) integer(bit_kind) :: key_mask_neg(Nint,2), mobileMask(Nint,2) integer :: mo_tot_num_2 mo_tot_num_2 = mo_tot_num+mo_tot_num - + do i=1,Nint key_mask_neg(i,1) = not(key_mask(i,1)) key_mask_neg(i,2) = not(key_mask(i,2)) end do - + do i=0,mo_tot_num_2 N_microlist(i) = 0 enddo - + do i=1, N_minilist nt = 0 do j=1,Nint @@ -845,30 +759,30 @@ subroutine create_microlist_single(minilist, i_cur, N_minilist, key_mask, microl mobileMask(j,2) = iand(key_mask_neg(j,2), minilist(j,2,i)) nt += popcnt(mobileMask(j, 1)) + popcnt(mobileMask(j, 2)) end do - + if(nt > 3) then !! TOO MANY DIFFERENCES continue else if(nt < 3) then - if(i > i_cur) then + if(i < i_cur) then !!!!!!!!!!!!!!!!!!!!! DESACTIVADO N_microlist(:) = 0 !!!! PAST LINKED TO EVERYBODY! ptr_microlist(:) = 1 - return + return else !! FUTUR LINKED TO EVERYBODY N_microlist(0) = N_microlist(0) + 1 endif else call bitstring_to_list(mobileMask(1,1), list(1,1), n_element(1), Nint) call bitstring_to_list(mobileMask(1,2), list(1,2), n_element(2), Nint) - + do s=1,2 do j=1,n_element(s) - nt = list(j,s) + mo_tot_num * (s-1) - N_microlist(nt) = N_microlist(nt) + 1 + nt = list(j,s) + mo_tot_num * (s-1) + N_microlist(nt) = N_microlist(nt) + 1 end do end do endif end do - + ptr_microlist(0) = 1 do i=1,mo_tot_num_2+1 ptr_microlist(i) = ptr_microlist(i-1) + N_microlist(i-1) @@ -877,18 +791,18 @@ subroutine create_microlist_single(minilist, i_cur, N_minilist, key_mask, microl do i=0,mo_tot_num_2+1 cur_microlist(i) = ptr_microlist(i) end do - - + + do i=1, N_minilist do j=1,Nint mobileMask(j,1) = iand(key_mask_neg(j,1), minilist(j,1,i)) mobileMask(j,2) = iand(key_mask_neg(j,2), minilist(j,2,i)) end do - + call bitstring_to_list(mobileMask(1,1), list(1,1), n_element(1), Nint) call bitstring_to_list(mobileMask(1,2), list(1,2), n_element(2), Nint) - + if(n_element(1) + n_element(2) < 3) then idx_microlist(cur_microlist(0)) = i do k=1,Nint @@ -916,7 +830,7 @@ end subroutine subroutine finish_isinwf(key_mask, keys, N_keys, isinwf) use bitmasks implicit none - + integer(bit_kind), intent(in) :: key_mask(N_int, 2), keys(N_int, 2, N_keys) integer(bit_kind) :: key_mask_neg(N_int, 2) integer(bit_kind) :: mobileMask(N_int, 2) @@ -924,31 +838,31 @@ subroutine finish_isinwf(key_mask, keys, N_keys, isinwf) integer, intent(in) :: N_keys integer :: i,j,nt,nt2,list(2,2), n_element(2) logical, external :: detEq - + do i=1,N_int key_mask_neg(i,1) = not(key_mask(i,1)) key_mask_neg(i,2) = not(key_mask(i,2)) end do - + do i=1, N_keys nt = 0 - + do j=1,N_int mobileMask(j,1) = iand(key_mask_neg(j,1), keys(j,1,i)) mobileMask(j,2) = iand(key_mask_neg(j,2), keys(j,2,i)) nt += popcnt(mobileMask(j, 1)) + popcnt(mobileMask(j, 2)) end do - + if(nt /= 2) cycle - + call bitstring_to_list(mobileMask(1,1), list(1,1), n_element(1), N_int) call bitstring_to_list(mobileMask(1,2), list(1,2), n_element(2), N_int) - + if(n_element(1) >= 1) nt = list(1,1) if(n_element(1) == 2) nt2 = list(2,1) if(n_element(2) == 2) nt = list(2,2) + mo_tot_num if(n_element(2) >= 1) nt2 = list(1,2) + mo_tot_num - + isinwf(nt, nt2) = .true. isinwf(nt2, nt) = .true. end do @@ -962,37 +876,41 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl implicit none integer, intent(in) :: Nint, i_cur, N_minilist integer(bit_kind), intent(in) :: minilist(Nint,2,N_minilist), key_mask(Nint,2) - + integer, intent(out) :: N_microlist(0:mo_tot_num*2), ptr_microlist(0:mo_tot_num*2+1), idx_microlist(N_minilist*4) integer(bit_kind), intent(out) :: microlist(Nint,2,N_minilist*4) - + integer, intent(out) :: N_tmicrolist(0:mo_tot_num*2), ptr_tmicrolist(0:mo_tot_num*2+1), idx_tmicrolist(N_minilist*4) integer(bit_kind), intent(out) :: tmicrolist(Nint,2,N_minilist*4) - - - integer :: i,j,k,s,nt,nt2,n_element(2,N_minilist), idx(0:N_minilist) - integer :: list(4,2,N_minilist), cur_microlist(0:mo_tot_num*2+1), cur_tmicrolist(0:mo_tot_num*2+1) + + + integer :: i,j,k,s,nt,nt2 + integer, allocatable :: n_element(:,:), idx(:), list(:,:,:) + integer :: cur_microlist(0:mo_tot_num*2+1), cur_tmicrolist(0:mo_tot_num*2+1) integer(bit_kind) :: key_mask_neg(Nint,2), mobileMask(Nint,2) integer :: mo_tot_num_2 logical,intent(out) :: isinwf(mo_tot_num*2, mo_tot_num*2) double precision, intent(out) :: d0s(mo_tot_num, mo_tot_num, N_states) double precision :: integ(mo_tot_num, mo_tot_num) + + allocate(list(4,2,N_minilist), n_element(2,N_minilist), idx(0:N_minilist)) + isinwf = .false. integ = 0d0 d0s = 0d0 mo_tot_num_2 = mo_tot_num+mo_tot_num - + idx(0) = 0 do i=1,Nint key_mask_neg(i,1) = not(key_mask(i,1)) key_mask_neg(i,2) = not(key_mask(i,2)) end do - + do i=0,mo_tot_num_2 N_microlist(i) = 0 N_tmicrolist(i) = 0 enddo - + do i=1, N_minilist nt = 0 do j=1,Nint @@ -1000,44 +918,24 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl mobileMask(j,2) = iand(key_mask_neg(j,2), minilist(j,2,i)) nt += popcnt(mobileMask(j, 1)) + popcnt(mobileMask(j, 2)) end do - + if(nt > 4) cycle !! TOO MANY DIFFERENCES - idx(0) += 1 idx(idx(0)) = i - + call bitstring_to_list(mobileMask(1,1), list(1,1,idx(0)), n_element(1, idx(0)), Nint) call bitstring_to_list(mobileMask(1,2), list(1,2,idx(0)), n_element(2, idx(0)), Nint) - - - if(nt <= 2) then - if(i > i_cur) then - N_microlist = 0 + + + if(nt == 2) then + if(i < i_cur) then + N_microlist(:) = 0 ptr_microlist = 1 - N_tmicrolist = 0 + N_tmicrolist = 0 ptr_tmicrolist = 1 - return + return else - !n_element(:, idx(0)) = (/2, 0/) N_microlist(0) = N_microlist(0) + 1 - - if(n_element(1,idx(0)) >= 1) nt = list(1,1,idx(0)) - if(n_element(1,idx(0)) == 2) nt2 = list(2,1,idx(0)) - if(n_element(2,idx(0)) == 2) nt = list(2,2,idx(0)) + mo_tot_num - if(n_element(2,idx(0)) >= 1) nt2 = list(1,2,idx(0)) + mo_tot_num - - isinwf(nt, nt2) = .true. - isinwf(nt2, nt) = .true. - double precision, external :: get_mo_bielec_integral - nt = mod(nt, mo_tot_num) - nt2 = mod(nt2, mo_tot_num) - ! call get_mo_bielec_integrals_ij(nt, nt2 ,mo_tot_num,integ,mo_integrals_map) - -! do j=1, N_states -! call i_h_j -! d0s(:,:,j) += integ(:,:) * psi_selectors_coef(i,j) !!!!!!!!!!!!!!!!!!! MOOOOOCHE !!!!! suppose que minilist = psi_selectors ..... -! end do -! print *, "TO ", integ(mod(nt, mo_tot_num), mod(nt2, mo_tot_num)) endif else do s=1,2 @@ -1049,7 +947,7 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl end do endif end do - + ptr_microlist(0) = 1 ptr_tmicrolist(0) = 1 do i=1,mo_tot_num_2+1 @@ -1061,10 +959,9 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl cur_microlist(i) = ptr_microlist(i) cur_tmicrolist(i) = ptr_tmicrolist(i) end do - - + + do i=1, idx(0) - if(n_element(1, i) + n_element(2, i) > 4) stop "wired" if(n_element(1, i) + n_element(2, i) <= 2) then idx_microlist(cur_microlist(0)) = idx(i) do k=1,Nint @@ -1072,11 +969,31 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl microlist(k,2,cur_microlist(0)) = minilist(k,2,idx(i)) enddo cur_microlist(0) = cur_microlist(0) + 1 + + if(n_element(1,i) >= 1) nt = list(1,1,i) + if(n_element(1,i) == 2) nt2 = list(2,1,i) + if(n_element(2,i) == 2) nt = list(2,2,i) + mo_tot_num + if(n_element(2,i) >= 1) nt2 = list(1,2,i) + mo_tot_num + + isinwf(nt, nt2) = .true. + isinwf(nt2, nt) = .true. + !!!! INTEGRAL DRIVEN + !!!!!!!!!!!!!!!!!!!! +! call get_d0(minilist(1,1,idx(i)), integ, key_mask, 1+(nt2-1)/mo_tot_num, 1+(nt-1)/mo_tot_num, & +! mod(nt2-1, mo_tot_num)+1, mod(nt-1, mo_tot_num)+1) +! +! do j=1, N_states +! do nt=1, mo_tot_num +! do nt2=1, mo_tot_num +! d0s(nt,nt2,j) = d0s(nt,nt2,j) + (integ(nt,nt2) * psi_selectors_coef(idx(i), j)) !!! SUPPOSE MINILIST = SELECTORS !!!! +! end do +! end do +! end do else do s = 1, 2 do j=1,n_element(s,i) nt = list(j,s,i) + mo_tot_num * (s-1) - + if(n_element(1,i) + n_element(2,i) == 4) then idx_microlist(cur_microlist(nt)) = idx(i) do k=1,Nint @@ -1102,15 +1019,15 @@ end subroutine subroutine check_past(det, list, idx, N, cur, ok, Nint) implicit none use bitmasks - + integer(bit_kind), intent(in) :: det(Nint, 2), list(Nint, 2, N) integer, intent(in) :: Nint, idx(N), N, cur logical, intent(out) :: ok integer :: i,s,ni - + ok = .false. - do i=N,1,-1 - if(idx(i) <= cur) exit + do i=1,N + if(idx(i) >= cur) exit s = 0 do ni=1,Nint s += popcnt(xor(det(ni,1), list(ni,1,i))) + popcnt(xor(det(ni,2), list(ni,2,i))) @@ -1123,15 +1040,15 @@ subroutine check_past(det, list, idx, N, cur, ok, Nint) end subroutine -subroutine check_futur(det, list, N, ok, Nint) +subroutine check_past_s(det, list, N, ok, Nint) implicit none use bitmasks - + integer(bit_kind), intent(in) :: det(Nint, 2), list(Nint, 2, N) integer, intent(in) :: Nint, N logical, intent(out) :: ok integer :: i,s,ni - + ok = .false. do i=1,N s = 0 @@ -1145,3 +1062,204 @@ subroutine check_futur(det, list, N, ok, Nint) end do end subroutine + +subroutine get_d0(gen, mat, mask, s1, s2, h1, h2) + use bitmasks + implicit none + + double precision, intent(out) :: mat(mo_tot_num, mo_tot_num) + double precision :: mat_mwen(mo_tot_num, mo_tot_num) + integer, intent(in) :: h1, h2, s1, s2 + integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2) + integer(bit_kind) :: det1(N_int, 2), det2(N_int, 2) + logical :: ok, mono + double precision :: phase, phase2, inv + integer :: p1, p2, hmi, hma + logical, external :: detEq + integer :: exc(0:2, 2, 2), exc2(0:2,2,2) + + exc = 0 + call get_mo_bielec_integrals_ij(h1, h2 ,mo_tot_num,mat_mwen,mo_integrals_map) + mat = 0d0 + if(s1 == s2) then + hmi = min(h1, h2) + hma = max(h1, h2) + inv = 1d0 + if(h1 > h2) inv = -1d0 + exc(0, :, s1) = 2 + exc(1, 1, s1) = hmi + exc(2, 1, s1) = hma + do p1=1,mo_tot_num + do p2=1,mo_tot_num + if(p1 == p2) cycle + call apply_particle(mask, (/s1,p1,s2,p2/), det2, ok, N_int) + if(.not. ok) cycle + mono = (hmi == p1 .or. hma == p2 .or. hmi == p2 .or. hma == p1) + if(mono) then + call i_h_j(gen, det2, N_int, mat(p1, p2)) + !mat(p1, p2) = phase + else + exc(1, 2, s1) = min(p1, p2) + exc(2, 2, s2) = max(p2, p1) + call get_double_excitation_phase(gen, det2, exc, phase, N_int) + mat(p1, p2) = inv * (mat_mwen(p1, p2) - mat_mwen(p2, p1)) * phase + end if + end do + end do + else + exc(0, :, 1) = 1 + exc(0, :, 2) = 1 + if(s1 /= 2) stop "alpha beta inversified" + exc(1, 1, 1) = h2 + exc(1, 1, 2) = h1 + + do p1=1, mo_tot_num + do p2=1, mo_tot_num + call apply_particle(mask, (/s1,p1,s2,p2/), det2, ok, N_int) + if(.not. ok) cycle + mono = (h1 == p1 .or. h2 == p2) + if(mono) then + call i_h_j(gen, det2, N_int, phase) + mat(p1, p2) = phase + else + exc(1, 2, s1) = p1 + exc(1, 2, s2) = p2 + call get_double_excitation_phase(gen, det2, exc, phase, N_int) + mat(p1, p2) = mat_mwen(p1, p2) * phase + end if + end do + end do + end if +end subroutine + + +subroutine apply_particle(det, exc, res, ok, Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer, intent(in) :: exc(4) + integer :: s1, s2, p1, p2 + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: ii, pos + + ok = .false. + s1 = exc(1) + p1 = exc(2) + s2 = exc(3) + p2 = exc(4) + res = det + + if(p1 /= 0) then + ii = (p1-1)/bit_kind_size + 1 + pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1) + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s1) = ibset(res(ii, s1), pos) + end if + + ii = (p2-1)/bit_kind_size + 1 + pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s2) = ibset(res(ii, s2), pos) + + ok = .true. +end subroutine + + +subroutine get_double_excitation_phase(det1,det2,exc,phase,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Returns the two excitation operators between two doubly excited determinants and the phase + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det1(Nint,2) + integer(bit_kind), intent(in) :: det2(Nint,2) + integer, intent(in) :: exc(0:2,2,2) + double precision, intent(out) :: phase + integer :: tz + integer :: l, ispin, idx_hole, idx_particle, ishift + integer :: nperm + integer :: i,j,k,m,n + integer :: high, low + integer :: a,b,c,d + integer(bit_kind) :: hole, particle, tmp + double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /) + + ASSERT (Nint > 0) + nperm = 0 + do ispin = 1,2 + select case (exc(0,1,ispin)) + case(0) + cycle + + case(1) + low = min(exc(1,1,ispin), exc(1,2,ispin)) + high = max(exc(1,1,ispin), exc(1,2,ispin)) + + ASSERT (low > 0) + j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint) + n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size) + ASSERT (high > 0) + k = ishft(high-1,-bit_kind_shift)+1 + m = iand(high-1,bit_kind_size-1)+1 + + if (j==k) then + nperm = nperm + popcnt(iand(det1(j,ispin), & + iand( ibset(0_bit_kind,m-1)-1_bit_kind, & + ibclr(-1_bit_kind,n)+1_bit_kind ) )) + else + nperm = nperm + popcnt(iand(det1(k,ispin), & + ibset(0_bit_kind,m-1)-1_bit_kind)) + if (n < bit_kind_size) then + nperm = nperm + popcnt(iand(det1(j,ispin), ibclr(-1_bit_kind,n) +1_bit_kind)) + endif + do i=j+1,k-1 + nperm = nperm + popcnt(det1(i,ispin)) + end do + endif + + case (2) + + do i=1,2 + low = min(exc(i,1,ispin), exc(i,2,ispin)) + high = max(exc(i,1,ispin), exc(i,2,ispin)) + + ASSERT (low > 0) + j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint) + n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size) + ASSERT (high > 0) + k = ishft(high-1,-bit_kind_shift)+1 + m = iand(high-1,bit_kind_size-1)+1 + + if (j==k) then + nperm = nperm + popcnt(iand(det1(j,ispin), & + iand( ibset(0_bit_kind,m-1)-1_bit_kind, & + ibclr(-1_bit_kind,n)+1_bit_kind ) )) + else + nperm = nperm + popcnt(iand(det1(k,ispin), & + ibset(0_bit_kind,m-1)-1_bit_kind)) + if (n < bit_kind_size) then + nperm = nperm + popcnt(iand(det1(j,ispin), ibclr(-1_bit_kind,n) +1_bit_kind)) + endif + do l=j+1,k-1 + nperm = nperm + popcnt(det1(l,ispin)) + end do + endif + + enddo + + a = min(exc(1,1,ispin), exc(1,2,ispin)) + b = max(exc(1,1,ispin), exc(1,2,ispin)) + c = min(exc(2,1,ispin), exc(2,2,ispin)) + d = max(exc(2,1,ispin), exc(2,2,ispin)) + if (c>a .and. cb) then + nperm = nperm + 1 + endif + exit + end select + + enddo + phase = phase_dble(iand(nperm,1)) +end