gen_uniform_double_para Subroutine

    subroutine gen_uniform_double_para(nI, ilutI, nJ, ilutJ, ex, tParity, pgen)
        ! routine to do a uniform paralles spin double 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_double_para"
#endif
        integer :: elecs(2), ispn, spin, i, a, b
        integer, parameter :: max_trials = 1000
        real(dp) :: p_elec, p_orb

        nJ(1) = 0
        ! pick two spin-parallel electrons
        call pick_spin_par_elecs(nI, elecs, p_elec, ispn)

        ! i have to figure out the probabilty of two spin-parallel
        if (ispn == 1) then
            spin = 1
            p_orb = 1.0_dp / real(nbasis / 2 - nOccBeta, dp)
        else if (ispn == 3) then
            spin = 2
            p_orb = 1.0_dp / real(nbasis / 2 - nOccAlpha, dp)
#ifdef DEBUG_
        else
            call stop_all(this_routine, "no parallel spin!")
#endif
        end if

        ! pick 2 holes now
        do i = 1, max_trials

            a = 2 * int(nbasis / 2 * genrand_real2_dsfmt()) + spin

            if (IsOcc(ilutI, a)) cycle

            b = get_orb_from_kpoints(nI(elecs(1)), nI(elecs(2)), a)

            ! do we have to reject or can we cycle if not fitting?
            ! a == b test has to be here for the spin-parallel
            ! excitations!
            ! try not returning but cycling
            if (IsOcc(ilutI, b) .or. a == b) then
                nJ(1) = 0
                return
            end if

            call make_double(nI, nJ, elecs(1), elecs(2), a, b, ex, tParity)

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

        ! times 2, since both ab, ba orders are possible
        pgen = p_elec * p_orb * 2.0_dp

    end subroutine gen_uniform_double_para