pick_occupied_orbital_ueg Subroutine

    subroutine pick_occupied_orbital_ueg(ilutI, src, ispn, part_type, cpt, &
                                         cum_sum, orb, calc_pgen)
        integer, intent(in) :: src(2), ispn, part_type
        integer(n_int), intent(in) :: ilutI(0:niftot)
        real(dp), intent(out) :: cpt, cum_sum
        integer, intent(out) :: orb
        logical, intent(in), optional :: calc_pgen
        character(*), parameter :: this_routine = "pick_occupied_orbital_ueg"

        integer :: occ_orbs(nel), n_valid, j, ind, i, orb_a, orb_b
        ! the only difference in the ueg orbital picker is that it is not
        ! restricted to pick a beta-orbital first in the case of
        ! a opposite spin excitation

        ! better idea:
        n_valid = 0
        j = 1
        occ_orbs = -1
        orb = -1

        ! loop over ref det
        ! in this routine i also have to check if the momentum is
        ! allowed! why didn't i do that before??
        do i = 1, nel
            orb_a = projedet(i, part_type_to_run(part_type))
            ! check if ref-det electron is NOT in nI
            if (IsNotOcc(ilutI, orb_a)) then
                ! check the symmetry here.. or atleast the spin..
                ! if we are parallel i have to ensure the orbital has the
                ! same spin
                if (is_allowed_ueg_k_vector(src(1), src(2), orb_a)) then
                    ! i also have to check if orb b is not occupied ofc..
                    ! what the fuck? why didn't i do that before
                    orb_b = get_orb_from_kpoints(src(1), src(2), orb_a)

                    if (IsNotOcc(ilutI, orb_b) .and. (orb_a /= orb_b)) then
                        if (ispn /= 2) then
                            if (is_beta(orb_a) .eqv. is_beta(src(1))) then
                                ! this is a valid orbital i guess..
                                n_valid = n_valid + 1
                                occ_orbs(j) = orb_a
                                j = j + 1
                            end if
                        else
                            ! in the ueg case we can pick the first orbital freely
                            ! for a ispn == 2 excitation..
                            n_valid = n_valid + 1
                            occ_orbs(j) = orb_a
                            j = j + 1
                        end if
                    end if
                end if
            end if
        end do

        ! so now we have a list of the possible orbitals in occ_orbs
        ! this has to be atleast 2, or otherwise we won't find a second
        ! orbital.. well no! since the second orbital can be picked from
        ! all the orbitals!
        if (n_valid == 0) then
            orb = 0
            cpt = 0.0_dp
            ! can i set cum_sum to 0 here, or does this invoke divisions by zero?
            cum_sum = 1.0_dp
            return
        end if

        ASSERT(n_valid > 0)
        ! and now the cum_sums and pgens..
        cpt = 1.0_dp / real(n_valid, dp)
        cum_sum = 1.0_dp

        if (present(calc_pgen)) then
            if (calc_pgen) return
        end if

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

        orb = occ_orbs(ind)

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

    end subroutine pick_occupied_orbital_ueg