BEGIN_PROVIDER [ double precision, integral8, (mo_tot_num, mo_tot_num, mo_tot_num, mo_tot_num) ] integral8 = 0d0 integer :: h1, h2 do h1=1, mo_tot_num do h2=1, mo_tot_num call get_mo_bielec_integrals_ij(h1, h2 ,mo_tot_num,integral8(1,1,h1,h2),mo_integrals_map) end do end do END_PROVIDER subroutine select_connected(i_generator,E0,pt2,b) !use f77_zmq 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)) ) 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_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 select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf) !use f77_zmq use bitmasks use selection_types implicit none BEGIN_DOC ! Select determinants connected to i_det by H END_DOC integer, intent(in) :: i_generator double precision, intent(in) :: fock_diag_tmp(mo_tot_num) double precision, intent(inout) :: pt2(N_states) 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 msg_size = bit_kind*N_int*2 ! Apply hole and particle masks ! ----------------------------- integer(bit_kind) :: hole(N_int,2), particle(N_int,2) 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 ! 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) ! Create excited determinants ! --------------------------- integer :: ispin, other_spin integer(bit_kind) :: exc_det(N_int,2), ion_det(N_int,2) 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) ion_det(k,1) = psi_det_generators(k,1,i_generator) ion_det(k,2) = psi_det_generators(k,2,i_generator) enddo 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 i=1, N_holes(ispin) ion_det(:,:) = psi_det_generators(:,:,i_generator) integer :: i_hole i_hole = hole_list(i,ispin) ! Apply the hole 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(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_transp(:,idx_microlist(j)) enddo if(ptr_microlist(mo_tot_num * 2 + 1) == 1) then cycle endif do j=1,N_particles(ispin) exc_det(:,:) = ion_det(:,:) integer :: i_particle i_particle = particle_list(j,ispin) 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(exc_det(k_particle,ispin),j_particle) logical, external :: is_in_wavefunction logical :: nok if (.not. is_in_wavefunction(exc_det,N_int)) then double precision :: i_H_psi_value(N_states), i_H_psi_value2(N_states) i_H_psi_value = 0d0 i_H_psi_value2 = 0d0 integer :: sporb nok = .false. 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 (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 call add_to_selection_buffer(buf, exc_det, e_pertm) endif ! Reset exc_det exc_det(k_particle,ispin) = psi_det_generators(k_particle,ispin,i_generator) enddo ! j ! Reset ion_det ion_det(k_hole,ispin) = psi_det_generators(k_hole,ispin,i_generator) enddo ! i enddo ! ispin end subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf) !use f77_zmq use bitmasks use selection_types implicit none BEGIN_DOC ! Select determinants connected to i_det by H END_DOC integer, intent(in) :: i_generator double precision, intent(in) :: fock_diag_tmp(mo_tot_num) double precision, intent(inout) :: pt2(N_states) 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 logical :: isinwf(mo_tot_num*2, mo_tot_num*2) double precision :: d0s(mo_tot_num, mo_tot_num, N_states) integer :: i,j,k,l,j1,j2,i1,i2,ib,jb integer :: msg_size msg_size = bit_kind*N_int*2 ! Apply hole and particle masks ! ----------------------------- integer(bit_kind) :: hole(N_int,2), particle(N_int,2) 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 ! 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) ! 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) ion_det(k,1) = psi_det_generators(k,1,i_generator) ion_det(k,2) = psi_det_generators(k,2,i_generator) enddo do ispin1=1,2 do ispin2=1,ispin1 integer :: i_hole1, i_hole2, j_hole, k_hole do i1=N_holes(ispin1),1,-1 ! Generate low excitations first if(ispin1 == ispin2) then ib = i1+1 else ib = 1 endif do i2=N_holes(ispin2),ib,-1 ! Generate low excitations first ion_det(:,:) = psi_det_generators(:,:,i_generator) 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(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_transp(:,idx_microlist(j)) enddo do j=1, ptr_tmicrolist(mo_tot_num * 2 + 1) - 1 psi_coef_tmicrolist(j,:) = psi_selectors_coef_transp(:,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)) < 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)) < i_generator) cycle if(isinwf(p1, p2)) cycle exc_det = ion_det if(N_microlist(p1) < N_microlist(p2)) then sporb = p1 lorb = p2 else 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) logical, external :: is_in_wavefunction logical :: nok ! 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 ! do j=1, N_states ! i_H_psi_value(j) = d0s(mod(p1-1, mo_tot_num)+1, mod(p2-1, mo_tot_num)+1, j) ! end do 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)) ! if(idx_microlist(c1) < idx_microlist(c2)) then ! c1 += 1 ! else if(idx_microlist(c1) > idx_microlist(c2)) then ! c2 += 1 ! else ! call i_H_j(exc_det,microlist(1,1,c1),N_int,hij) ! do j = 1, N_states ! i_H_psi_value2(j) = i_H_psi_value2(j) + psi_coef_microlist(c1,j)*hij ! end do ! c1 += 1 ! c2 += 1 ! endif ! end do ! else 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 if(.true.) then ! DET DRIVEN 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 end if 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 e_pertm += dabs(e_pert(k)) ! 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 enddo enddo enddo enddo end subroutine create_futur_ptr(ptr_microlist, idx_microlist, ptr_futur_microlist, N_futur_microlist, i_generator) 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 ptr_futur_microlist(i) = j N_futur_microlist(i) = ptr_microlist(i+1) - j exit end if end do end do end subroutine subroutine create_microlist_single(minilist, i_cur, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint) 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 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)) 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 !!!!!!!!!!!!!!!!!!!!! DESACTIVADO N_microlist(:) = 0 !!!! PAST LINKED TO EVERYBODY! ptr_microlist(:) = 1 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 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) end do 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 microlist(k,1,cur_microlist(0)) = minilist(k,1,i) microlist(k,2,cur_microlist(0)) = minilist(k,2,i) enddo cur_microlist(0) = cur_microlist(0) + 1 else if(n_element(1) + n_element(2) == 3) then do s = 1, 2 do j=1,n_element(s) nt = list(j,s) + mo_tot_num * (s-1) idx_microlist(cur_microlist(nt)) = i do k=1,Nint microlist(k,1,cur_microlist(nt)) = minilist(k,1,i) microlist(k,2,cur_microlist(nt)) = minilist(k,2,i) enddo cur_microlist(nt) = cur_microlist(nt) + 1 end do end do end if end do 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) logical,intent(inout) :: isinwf(mo_tot_num*2, mo_tot_num*2) 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 end subroutine subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, & tmicrolist, idx_tmicrolist, N_tmicrolist, ptr_tmicrolist, & isinwf, d0s, Nint) use bitmasks 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 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), tmp_det(Nint, 2) integer :: mo_tot_num_2, pwen(4), pweni 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) logical :: localbanned(mo_tot_num*2), banned(mo_tot_num*2), banned_pair(mo_tot_num*2, mo_tot_num*2), ok banned = .false. banned_pair = .false. 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 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)) 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 ptr_microlist = 1 N_tmicrolist = 0 ptr_tmicrolist = 1 return else N_microlist(0) = N_microlist(0) + 1 endif else do s=1,2 do j=1,n_element(s,idx(0)) k = list(j,s,idx(0)) + mo_tot_num * (s-1) if(nt == 4) N_microlist(k) = N_microlist(k) + 1 if(nt == 3) then N_tmicrolist(k) = N_tmicrolist(k) + 1 if(idx(i) < i_cur) banned(nt) = .true. end if end do end do endif end do ptr_microlist(0) = 1 ptr_tmicrolist(0) = 1 do i=1,mo_tot_num_2+1 ptr_microlist(i) = ptr_microlist(i-1) + N_microlist(i-1) ptr_tmicrolist(i) = ptr_tmicrolist(i-1) + N_tmicrolist(i-1) end do do i=0,mo_tot_num_2+1 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) == 2) cycle pweni = 0 do s = 1, 2 do j=1,n_element(s,i) nt = list(j,s,i) + mo_tot_num * (s-1) pweni += 1 pwen(pweni) = nt if(n_element(1,i) + n_element(2,i) == 4) then idx_microlist(cur_microlist(nt)) = idx(i) do k=1,Nint microlist(k,1,cur_microlist(nt)) = minilist(k,1,idx(i)) microlist(k,2,cur_microlist(nt)) = minilist(k,2,idx(i)) enddo cur_microlist(nt) = cur_microlist(nt) + 1 else idx_tmicrolist(cur_tmicrolist(nt)) = idx(i) do k=1,Nint tmicrolist(k,1,cur_tmicrolist(nt)) = minilist(k,1,idx(i)) tmicrolist(k,2,cur_tmicrolist(nt)) = minilist(k,2,idx(i)) enddo cur_tmicrolist(nt) = cur_tmicrolist(nt) + 1 endif end do end do if(idx(i) < i_cur .and. pweni == 4) then do j=1,4 do k=j+1,4 banned_pair(pwen(j), pwen(k)) = .true. banned_pair(pwen(k), pwen(j)) = .true. end do end do end if end do do i=1, idx(0) if(n_element(1, i) + n_element(2, i) <= 2) then idx_microlist(cur_microlist(0)) = idx(i) do k=1,Nint microlist(k,1,cur_microlist(0)) = minilist(k,1,idx(i)) 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)), banned, banned_pair, d0s, 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, psi_selectors_coef_transp(1,idx(i))) ! do j=1, N_states ! do nt2=1, mo_tot_num ! do nt=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 if(.false. .and. n_element(1, i) + n_element(2, i) == 3) then ! INTEGRAL DRIVEN ! -459.6399263191298 pweni = 0 do s = 1, 2 do j=1,n_element(s,i) nt = list(j,s,i) + mo_tot_num * (s-1) pweni += 1 pwen(pweni) = nt end do end do call get_d1(minilist(1,1,idx(i)), banned, banned_pair, d0s, key_mask, pwen, psi_selectors_coef_transp(1,idx(i))) ! do k=1, N_states ! do nt2=1, mo_tot_num ! do nt=1, mo_tot_num ! d0s(nt,nt2,k) = d0s(nt,nt2,k) + (integ(nt,nt2) * psi_selectors_coef(idx(i), k)) !!! SUPPOSE MINILIST = SELECTORS !!!! ! end do ! end do ! end do end if end do end subroutine subroutine get_d1(gen, banned, banned_pair, mat, mask, pwen, coefs) use bitmasks implicit none integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2) logical, intent(in) :: banned(mo_tot_num*2), banned_pair(mo_tot_num*2, mo_tot_num*2) integer(bit_kind) :: deth(N_int, 2), det(N_int, 2), i8 double precision, intent(in) :: coefs(N_states) double precision, intent(inout) :: mat(mo_tot_num, mo_tot_num, N_states) double precision :: hij, phase, inv, inv2 integer, intent(in) :: pwen(3) integer :: s(3), p(3), i, j, k, h1, h2, ns(2), sm, mwen, a1, a2, pwens(2), sp, st integer :: sfix, pfix integer :: exc(0:2, 2, 2) logical :: lbanned(mo_tot_num*2) logical :: ok, mono, ab integer :: tmp_array(4) lbanned = banned !mat = 0d0 pwens = 0 ns = 0 do sp=1,2 do i=1, N_int ns(sp) += popcnt(gen(i, sp)) - popcnt(mask(i, sp)) i8 = iand(not(gen(i, sp)), mask(i, sp)) if(i8 /= 0_8) then sfix = sp pfix = 1+trailz(i8) + bit_kind*8*(i-1) end if end do end do do i=1,3 s(i) = 1+(pwen(i)-1)/mo_tot_num p(i) = 1+mod(pwen(i)-1, mo_tot_num) pwens(s(i)) += 1 end do do i=1,3 if(s(i) == 1 .and. ns(1) == 0) cycle if(s(i) == 2 .and. ns(2) == 0) cycle if(lbanned(pwen(i))) cycle ab = pwens(s(i)) == 2 if(ns(1) == 1) sm = mod(s(i), 2) + 1 if(ns(1) == 2) sm = 1 if(ns(2) == 2) sm = 2 lbanned(pwen(i)) = .true. if(ab) then if(s(mod(i,3)+1) == 2) then a1 = mod(i, 3) + 1 a2 = mod(i+1, 3) + 1 else a2 = mod(i,3)+1 a1 = mod(i+1,3)+1 end if exc(0, :, 1) = 1 exc(0, :, 2) = 1 exc(1, 1, 1) = p(a2) exc(1, 1, 2) = p(a1) exc(1, 2, sfix) = pfix call apply_particle(mask, 0, 0 ,s(i), p(i), deth, ok, N_int) do j=1,mo_tot_num mwen = j + (sm-1)*mo_tot_num if(lbanned(mwen)) cycle tmp_array = (/0,0,sm,j/) call apply_particle(deth, 0, 0 ,s(i), p(i), det, ok, N_int) if(.not. ok) cycle mono = mwen == pwen(a1) .or. mwen == pwen(a2) if(mono) then call i_h_j(gen, det, N_int, hij) else exc(1, 2, sm) = j call get_double_excitation_phase(gen, det, exc, phase, N_int) if(sfix == 1) hij = integral8(j, pfix, p(a1), p(a2)) * phase if(sfix == 2) hij = integral8(pfix, j, p(a1), p(a2)) * phase end if if(ns(1) == 1) then do st=1, N_states if(sm == 2) mat(j, p(i), st) = mat(j, p(i), st) + hij * coefs(st) if(sm == 1) mat(p(i), j, st) = mat(p(i), j, st) + hij * coefs(st) end do else do st=1, N_states mat(j, p(i), st) += hij * coefs(st) mat(p(i), j, st) += hij * coefs(st) end do end if end do else !! AA / BB a1 = mod(i,3)+1 a2 = mod(i+1,3)+1 h1 = p(a1) h2 = p(a2) inv = 1d0 if(h1 > h2) inv = -1d0 if(pwens(s(i)) == 1) sp = mod(s(i), 2)+1 if(pwens(s(i)) == 3) sp = s(i) exc(0, :, sp) = 2 exc(0, :, mod(sp, 2)+1) = 0 exc(1, 1, sp) = min(h1, h2) exc(2, 1, sp) = max(h1, h2) call apply_particle(mask, 0, 0 ,s(i), p(i) , deth, ok, N_int) do j=1,mo_tot_num if(j == pfix) inv = -inv mwen = j + (sm-1)*mo_tot_num if(lbanned(mwen)) cycle call apply_particle(deth, 0, 0 ,s(i), p(i) , det, ok, N_int) if(.not. ok) cycle mono = mwen == pwen(a1) .or. mwen == pwen(a2) if(mono) then call i_h_j(gen, det, N_int, hij) else exc(1, 2, sfix) = min(j,pfix) exc(2, 2, sp) = max(j,pfix) call get_double_excitation_phase(gen, det, exc, phase, N_int) hij = (integral8(j, pfix, h1, h2) - integral8(pfix,j, h1, h2))*phase*inv end if if(ns(1) == 1) then do st=1, N_states if(sm == 2) mat(j, p(i), st) = mat(j, p(i), st) + hij * coefs(st) if(sm == 1) mat(p(i), j, st) = mat(p(i), j, st) + hij * coefs(st) end do else do st=1, N_states mat(j, p(i), st) += hij * coefs(st) mat(p(i), j, st) += hij * coefs(st) end do end if end do end if end do end subroutine subroutine get_d0(gen, banned, banned_pair, mat, mask, s1, s2, h1, h2, coefs) use bitmasks implicit none double precision, intent(inout) :: mat(mo_tot_num, mo_tot_num, N_states) logical, intent(in) :: banned(mo_tot_num*2), banned_pair(mo_tot_num*2, mo_tot_num*2) double precision :: mat_mwen(mo_tot_num, mo_tot_num) double precision, intent(in) :: coefs(N_states) 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, hij integer :: p1, p2, hmi, hma, ns1, ns2, st logical, external :: detEq integer :: exc(0:2, 2, 2), exc2(0:2,2,2) integer :: tmp_array(4) exc = 0 ! mat_mwen = integral8(:,:,h1,h2) !call get_mo_bielec_integrals_ij(h1, h2 ,mo_tot_num,mat_mwen,mo_integrals_map) ns1 = mo_tot_num*(s1-1) ns2 = mo_tot_num*(s2-1) !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 p2=1,mo_tot_num if(banned(p2 + ns2)) cycle do p1=1,mo_tot_num if(banned(p1 + ns1)) cycle if(p1 == p2) cycle if(banned_pair(p1 + ns1, p2 + ns2)) cycle !tmp_array = (/s1,p1,s2,p2/) 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, hij) do st=1, N_states mat(p1, p2, st) += hij * coefs(st) end do 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) do st=1, N_states mat(p1, p2, st) += coefs(st) * inv * (integral8(p1, p2, h1, h2) - integral8(p2, p1, h1, h2)) * phase end do 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 p2=1, mo_tot_num if(banned(p2 + ns2)) cycle do p1=1, mo_tot_num if(banned(p1 + ns1)) cycle if(banned_pair(p1 + ns1, p2 + ns2)) cycle !tmp_array = (/s1,p1,s2,p2/) 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, hij) do st=1, N_states mat(p1, p2, st) += hij * coefs(st) end do else exc(1, 2, s1) = p1 exc(1, 2, s2) = p2 call get_double_excitation_phase(gen, det2, exc, phase, N_int) do st=1, N_states mat(p1, p2, st) += coefs(st) * integral8(p1, p2, h1, h2) * phase end do !mat(p1, p2) = integral8(p1, p2, h1, h2) * phase end if end do end do end if 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=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))) end do if(s <= 4) then ok = .true. return end if end do end subroutine 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 do ni=1,Nint s += popcnt(xor(det(ni,1), list(ni,1,i))) + popcnt(xor(det(ni,2), list(ni,2,i))) end do if(s <= 4) then ok = .true. return end if end do end subroutine