pick_occupied_orbital_hubbard Subroutine

    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