pick_occupied_orbital Subroutine

    subroutine pick_occupied_orbital(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
        ! routine to pick an orbital of the occupied manifold of the
        ! reference determinant uniformly
        ! to be compatible with the rest of the 4ind-weighted-2
        ! excitation generators i have to be carefull with the cum_lists
        ! and stuff..
        character(*), parameter :: this_routine = "pick_occupied_orbital"
        integer :: occ_orbs(nel), n_valid, j, ind, i, orb_a

        ! soo what do i need?
        ! i have to check if any of the possible orbitals for nI is occupied
        ! in the reference determinant!

        ! better idea:
        n_valid = 0
        j = 1
        occ_orbs = -1
        orb = -1
        ! loop over ref det
        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 (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
                    ! there is some weird shenanigan in the gen_a_orb_cum_list
                    ! if the spins are anti-parallel.. why?
                    ! this "only" has to do with the weighting of the
                    ! matrix element.. so it does not affect me here i guess
                    ! so here all the orbitals are alowed..
                    ! UPDATE: nope this also implies that (a) is always a
                    ! beta orbital for anti-parallel spin excitations
                    ! i do not know why exactly, but somebody decided to do
                    ! it this way.. so just to be sure, also do it like that
                    ! in the back-spawn method
                    if (is_beta(orb_a)) then
                        n_valid = n_valid + 1
                        occ_orbs(j) = orb_a
                        j = j + 1
                    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)

        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