subroutine pick_orbitals_guga_heisenberg(ilut, nI, csf_i, excitInfo, orb_pgen)
! i "just" need to implement a custom orbital picker for the
! spin-free Heisenberg exchange
integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
integer, intent(in) :: nI(nel)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(out) :: excitInfo
real(dp), intent(out) :: orb_pgen
integer :: elec, src, id, ind, tgt, start, ende
real(dp) :: p_elec, cum_sum, r, p_orb
integer, allocatable :: neighbors(:)
real(dp), allocatable :: cum_arr(:)
! Exists for the function pointer interface
unused_var(ilut)
! here i want to only pick nearest neighbor electrons, where a
! spin recoupling is possible
! i need to pick one orbital at random: since it should be
! general for the t-J model too, pick an occupied orbital
elec = 1 + int(genrand_real2_dsfmt() * nel)
p_elec = 1.0_dp / real(nel, dp)
src = nI(elec)
! spatial orbital:
id = gtID(src)
neighbors = lat%get_neighbors(lat%get_site_index(id))
call gen_guga_heisenberg_cum_list(csf_i, id, cum_arr)
cum_sum = cum_arr(size(neighbors))
if (cum_sum < EPS) then
excitInfo%valid = .false.
orb_pgen = 0.0_dp
return
end if
r = genrand_real2_dsfmt() * cum_sum
ind = binary_search_first_ge(cum_arr, r)
tgt = neighbors(ind)
if (ind == 1) then
p_orb = cum_arr(1) / cum_sum
else
p_orb = (cum_arr(ind) - cum_arr(ind - 1)) / cum_sum
end if
! and now we have to to the GUGA excitation:
! were I just realised, everything will get recalculated anyway..
! so i have to adapt the pgen recalculation.. to fit this
! scheme above..
start = min(id, tgt)
ende = max(id, tgt)
excitInfo = assign_excitInfo_values_double(excit_type%fullstart_stop_mixed, &
gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
start, ende, ende, start, start, start, ende, ende, 0, 2, 1.0_dp, 1.0_dp)
orb_pgen = p_elec * p_orb
end subroutine pick_orbitals_guga_heisenberg