gen_uniform_triple Subroutine

    subroutine gen_uniform_triple(nI, ilutI, nJ, ilutJ, ex, tParity, pgen)
        ! routine to do a uniform triple excitation
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:niftot)
        integer, intent(out) :: nJ(nel), ex(2, 3)
        integer(n_int), intent(out) :: ilutJ(0:niftot)
        logical, intent(out) :: tParity
        real(dp), intent(out) :: pgen
#ifdef DEBUG_
        character(*), parameter :: this_routine = "gen_uniform_triple"
#endif
        integer :: a, b, c, elecs(3), i, ispn, src(3), sum_ms
        integer, parameter :: max_trials = 1000
        real(dp) :: p_elec, p_orb, p_orb_a

        nJ(1) = 0

        call pick_three_opp_elecs(nI, elecs, p_elec, sum_ms)
        src = nI(elecs)

        ASSERT(sum_ms == -1 .or. sum_ms == 1)

        call pick_a_orbital_hubbard(ilutI, a, p_orb_a, sum_ms)

        ! and now i have to pick orbital b and fitting c in a uniform
        ! way.. i hope this still works with the probabilities
        ! if A is beta, we need to pick a alpha B uniformly and vv.
        if (is_beta(a)) then
            p_orb = 1.0_dp / real(nbasis / 2 - nOccAlpha, dp)
            ! also use a spin to specify the spin-orbital
            ! is a is beta we want an alpha -> so add +1
            ispn = 0
        else
            p_orb = 1.0_dp / real(nBasis / 2 - nOccBeta, dp)
            ispn = 1
        end if

        do i = 1, max_trials

            b = 2 * (1 + int(genrand_real2_dsfmt() * nbasis / 2)) - ispn

            if (IsOcc(ilutI, b)) cycle

            c = get_orb_from_kpoints_three(src, a, b)

            if (IsOcc(ilutI, c) .or. b == c) then
                nJ(1) = 0
                return
            end if

            call make_triple(nI, nJ, elecs, [a, b, c], ex, tParity)

            ilutJ = make_ilutJ(ilutI, ex, 3)
            exit
        end do

        ! times 2 since BC <> CB is both possible
        pgen = p_elec * p_orb * p_orb_a * 2.0_dp

    end subroutine gen_uniform_triple