pick_three_orbs_sym Subroutine

    subroutine pick_three_orbs_sym(nI, src, tgt, pgen, ms)
        ! picks three random unoccupied orbitals, given the occupied orbitals
        integer, intent(in) :: nI(nel), ms, src(3)
        integer, intent(out) :: tgt(3)
        real(dp), intent(inout) :: pgen

        integer :: i, msCur, msOrb
        integer :: cc_occ(ScratchSize), cc_unocc(ScratchSize), tgt_sym
        ! pool to draw each orbital from
        integer, allocatable :: pool(:)
        ! size of the pools
        integer :: pool_sizes(3)
        real(dp) :: pgen_pick

        msCur = ms
        pgen_pick = pgen
        ! First, pick two arbitrary orbitals
        do i = 0, 1
            ! get the ms of this orb
            ! we take ms = 1 until the total leftover ms is negative
            ! This ensures the first two electrons have the same spin
            if (msCur > 0) then
                msOrb = 1
            else
                ! then we start taking ms = -1
                msOrb = -1
            end if
            call create_full_pool(nI, tgt, msOrb, pool, i)
            call get_orb_from_pool(tgt, pool, i, pgen_pick)
            ! the remaining ms
            msCur = msCur - msOrb
            ! keep the size of the pool in memory
            pool_sizes(i + 1) = size(pool)
        end do

        ! Get the number of orbitals per symmetry/occupation
        call construct_class_counts(nI, cc_occ, cc_unocc)
        ! Now, pick the last one according to symmetry
        tgt_sym = get_tgt_sym(tgt, src)
        call create_sym_pool(nI, tgt, msCur, pool, i, tgt_sym, cc_unocc)
        call get_orb_from_pool(tgt, pool, i, pgen_pick)
        pool_sizes(i + 1) = size(pool)

        ! There might be no symmetry-allowed picks for the last orbital - in that case
        ! abort the excitation
        if (any(tgt == 0)) then
            tgt = 0
            return
        end if

        call add_permutations_to_pgen(pgen, pgen_pick, pool_sizes, ms, tgt, cc_unocc)

        ! sort the target orbitals for further usage
        tgt = sort_unique(tgt)
    end subroutine pick_three_orbs_sym