pick_hole_from_active_space Function

    function pick_hole_from_active_space(this, ilutI, nI, iGAS, ms, r, pgen) result(tgt)
        ! pick a hole of ilutI with spin ms from the active space with index srcGASInd
        ! the random number is to be supplied as r
        class(GAS_disc_ExcGenerator_t), intent(in) :: this
        integer(n_int), intent(in) :: ilutI(0:NIfD)
        integer, intent(in) :: nI(nel)
        integer, intent(in) :: ms, iGAS
        real(dp), intent(in) :: r
        real(dp), intent(out) :: pgen
        integer :: tgt
        integer :: nEmpty, nOrb

        ! this sum only converts an array of size 1 to a scalar
        nEmpty = this%GAS_size(iGAS) - sum(popCnt(iand(ilutI(0:NIfD), this%GAS_spin_orbs(0:NIfD, iGAS, ms))))

        ! if there are no empyty orbitals in this gas (can happen when allowing for
        ! spin-exchange), the excitation is invalid
        if (nEmpty == 0) then
            tgt = 0
            pGen = 1.0_dp
            return
        end if
        ! adjust pgen
        pgen = 1.0_dp / real(nEmpty, dp)
        ! We can safely assume there is always an empty orb in each active space

        ! index of the target orb
        nOrb = int(r * nEmpty) + 1

        tgt = this%GAS_spin_orb_list(map_to_unoccupied(nOrb, iGAS, ms), iGAS, ms)

    contains

        !> Return the i-th unoccupied orbital in the iGAS space with spin ms.
        function map_to_unoccupied(i, iGAS, ms) result(new_orb_idx)
            integer, intent(in) :: i, iGAS, ms
            integer :: new_orb_idx

            integer :: j

            new_orb_idx = i
            do j = 1, nEl
                if (nI(j) > this%GAS_spin_orb_list(new_orb_idx, iGAS, ms)) return
                if (this%GAS_table(nI(j)) == iGAS .and. get_spin(nI(j)) == ms) then
                    ! if yes, we skip this orbital, increase by 1
                    new_orb_idx = new_orb_idx + 1
                end if
            end do
        end function
    end function pick_hole_from_active_space