# pick_ab_orbitals_hubbard Subroutine

## public subroutine pick_ab_orbitals_hubbard(nI, ilutI, src, orbs, p_orb)

### Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer(kind=n_int), intent(in) :: ilutI(0:niftot)
integer, intent(in) :: src(2)
integer, intent(out) :: orbs(2)
real(kind=dp), intent(out) :: p_orb

## Source Code

    subroutine pick_ab_orbitals_hubbard(nI, ilutI, src, orbs, p_orb)
! depending on the already picked electrons (ij) pick an orbital
! (a) and the connected orbital (b)
integer, intent(in) :: nI(nel), src(2)
integer(n_int), intent(in) :: ilutI(0:niftot)
integer, intent(out) :: orbs(2)
real(dp), intent(out) :: p_orb
#ifdef DEBUG_
real(dp) :: test
integer :: ex(2, 2)
#endif
real(dp) :: cum_arr(nbasis)
real(dp) :: cum_sum
integer :: orb_list(nbasis, 2)
integer :: ind

! without transcorrelation factor this is uniform, but with a
! transcorrelation factor the matrix element might change and so also
! the pgen should change.
call create_ab_list_hubbard(nI, ilutI, src, orb_list, cum_arr, cum_sum)

if (cum_sum < EPS) then
orbs(1) = ABORT_EXCITATION
return
end if

! this stuff is also written so often i should finally make a routine
! out of that
call pick_from_cum_list(cum_arr, cum_sum, ind, p_orb)

orbs = orb_list(ind, :)

#ifdef DEBUG_
! check that the other way of picking the orbital has the same
! probability..
call create_ab_list_hubbard(nI, ilutI, src, orb_list, cum_arr, cum_sum, &
orbs(2), test)

if (abs(test - p_orb) > 1.e-8) then
print *, "pgen assumption wrong:!"
print *, "p_orb: ", p_orb
print *, "test: ", test
print *, "orbs: ", orbs
end if

!todo: also call the calc_pgen_k_space_hubbard here and check
! pgens
ex(1, :) = src
ex(2, :) = orbs

#endif

! do i have to recalc. the pgen the other way around? yes!
! effectively reuse the above functionality
! i am pretty sure i just have to find the position in the
! list.. OR: since in the hubbard it is just twice the
! probability or? i am pretty sure yes.. but for all of them..
! so in the end it shouldnt matter again..
p_orb = 2.0_dp * p_orb

end subroutine pick_ab_orbitals_hubbard