pick_three_orbs_ueg Subroutine

    subroutine pick_three_orbs_ueg(nI, src, tgt, pgen, ms)
        use SymExcitDataMod, only: kPointToBasisFn
        use SystemData, only: nBasisMax, tOrbECutoff, OrbECutoff, nmaxx, nmaxy, nmaxz
        use sym_mod, only: mompbcsym

        ! 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, k1(3), k2(3), k3(3), p1(3), p2(3), p3(3), j
        integer :: k, l
        real(dp) :: testE
        logical :: i_occup, is_allowed

        msCur = ms
        do i = 0, 1
            ! get the ms of this orb
            ! we take ms = 1 until the total leftover ms is negative
            if (msCur > 0) then
                msOrb = 1
            else
                ! then we start taking ms = -1
                msOrb = -1
            end if
            call get_rand_orb(nI, tgt, msOrb, i, pgen)
            ! the remaining ms
            msCur = msCur - msOrb
        end do

        if (msCur > 0) then
            msOrb = 1
        else
            ! then we start taking ms = -1
            msOrb = -1
        end if

        do i = 1, 3
            k1(i) = G1(src(1))%k(i)
            k2(i) = G1(src(2))%k(i)
            k3(i) = G1(src(3))%k(i)
        end do

        do i = 1, 3
            p1(i) = G1(tgt(1))%k(i)
            p2(i) = G1(tgt(2))%k(i)
        end do

        p3 = k1 + k2 + k3 - p1 - p2

        ! Is p3 allowed by the size of the space?
        testE = sum(p3**2)
        if (abs(p3(1)) <= nmaxx .and. abs(p3(2)) <= nmaxy .and. &
            abs(p3(3)) <= nmaxz .and. &
            (.not. (tOrbECutoff .and. (testE > (OrbECutoff + 1.d-12))))) then
            is_allowed = .true.
        else
            is_allowed = .false.
        end if

        if (.not. is_allowed) then
            tgt(3) = 0
            pgen = 0.0_dp
            return
        end if

        i = kPointToBasisFn(p3(1), p3(2), p3(3), (msOrb + 1) / 2 + 1)

        if (i > nBasis .or. i < 1) then
            print *, 'bug kPointToBasisFn', p3, msOrb, i
            do i = -1, 1
            do j = -1, 1
            do k = -1, 1
            do l = -1, 1
                print *, i, j, k, l, kPointToBasisFn(i, j, k, (l + 1) / 2 + 1)
            end do
            end do
            end do
            end do
            call stop_all('pick_three_orbs_ueg', "kPointToBasisFn")
        end if

        ! i occupied? check:
        i_occup = .false.
        if (i == tgt(1) .or. i == tgt(2)) then
            i_occup = .true.
        else
            do j = 1, nel
                if (i == nI(j)) i_occup = .true.
            end do
        end if

        if (i_occup) then
            tgt(3) = 0
            pgen = 0.0_dp
            return
        else
            tgt(3) = i
            tgt = sort_unique(tgt)
        end if

        ! adjust the probability by taking permutations into account
        pgen = pgen * 2 * abs(ms)

    end subroutine pick_three_orbs_ueg