# pick_occupied_orbital_hubbard Subroutine

## public subroutine pick_occupied_orbital_hubbard(ilutI, part_type, pgen, orb, calc_pgen)

### Arguments

Type IntentOptional Attributes Name
integer(kind=n_int), intent(in) :: ilutI(0:niftot)
integer, intent(in) :: part_type
real(kind=dp), intent(out) :: pgen
integer, intent(out) :: orb
logical, intent(in), optional :: calc_pgen

## Source Code

    subroutine pick_occupied_orbital_hubbard(ilutI, part_type, pgen, orb, calc_pgen)
! routine to pick one possible orbital from the occupied manifold
! thats the easiest of all implementations actually..
integer, intent(in) :: part_type
integer(n_int), intent(in) :: ilutI(0:niftot)
real(dp), intent(out) :: pgen
integer, intent(out) :: orb
logical, intent(in), optional :: calc_pgen
character(*), parameter :: this_routine = "pick_occupied_orbital_hubbard"
integer :: n_valid, j, occ_orbs(nel), ind, i

integer :: orb_a
n_valid = 0
j = 1
occ_orbs = -1
orb = -1

! i also have to include the whole symmetry shabang in the
! picker here or?? wtf
do i = 1, nel
orb_a = projedet(i, part_type_to_run(part_type))

! what am i testing here, actually
! i want to loop over the reference det and check if it is
! no in nI! that can stay..
! or i could just pass also ilut into this function and then
! check if ref(i) is occupied..
! and i also should include k-point symmetry here or??
! why didn't i do that and why does it work anyway..
! nah.. in the hubbard this just works fine..
if (IsNotOcc(ilutI, orb_a)) then
n_valid = n_valid + 1
occ_orbs(j) = orb_a
j = j + 1
end if
end do

if (n_valid == 0) then
orb = 0
pgen = 0.0_dp
return
end if

ASSERT(n_valid > 0)

pgen = 1.0_dp / real(n_valid, dp)
if (present(calc_pgen)) then
if (calc_pgen) return
end if

ind = 1 + int(n_valid * genrand_real2_dSFMT())
orb = occ_orbs(ind)

ASSERT(orb > 0)
ASSERT(orb <= nbasis)

end subroutine pick_occupied_orbital_hubbard