generate_nGAS_double Subroutine

    subroutine generate_nGAS_double(this, nI, ilutI, nJ, ilutJ, ex, par, pgen)
        class(GAS_disc_ExcGenerator_t), intent(in) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(out) :: nJ(nel), ex(2, 2)
        integer(n_int), intent(out) :: ilutJ(0:NIfTot)
        logical, intent(out) :: par
        real(dp), intent(out) :: pgen

        integer :: elecs(2), src(2), sym_product, ispn, sum_ml, tgt(2)
        integer :: srcGAS(2)
        integer :: ms
        real(dp) :: r, pgen_pick1, pgen_pick2
        logical :: tExchange
        ! assuming that there are possible excitations within each active space,
        ! pick two random electrons (we would not include a full / empty space, so
        ! the assumption is very mild)
        call pick_biased_elecs(nI, elecs, src, sym_product, ispn, sum_ml, pgen)

        ! active spaces we consider
        srcGAS = this%GAS_table(src)

        tExchange = (ispn == 2) &
                    .and. (this%GAS_spec%recoupling() .or. srcGAS(1) == srcGAS(2))

        ! pick an empty orb from each active space chosen
        ! number of empty orbs in this active space
        r = genrand_real2_dSFMT()
        ! get the spin of the target orbital
        if (tExchange) then
            if (r > 0.5_dp) then
                ms = 1
                r = 2 * (r - 0.5_dp)
            else
                ms = 2
                r = 2 * (0.5_dp - r)
            end if
            pgen = pgen * 0.5_dp
        else
            ms = get_spin(src(1))
        end if
        tgt = 0
        ! the first hole is chosen randomly from the first active space
        tgt(1) = this%pick_hole_from_active_space(ilutI, nI, srcGAS(1), ms, genrand_real2_dSFMT(), pgen_pick1)
        if (tgt(1) == 0) then
            call zeroResult()
            return
        end if

        if (tExchange) then
            ! we picked the first spin randomly, now the second elec has the opposite spin
            ms = 3 - ms
        else
            ms = get_spin(src(2))
        end if
        ! the second hole is chosen in a weighted fashion
        tgt(2) = this%pick_weighted_hole(nI, Excite_2_t(src(1), tgt(1), src(2)), ms, srcGAS(2), pgen_pick1)
        if (any(tgt == 0) .or. tgt(1) == tgt(2)) then
            call zeroResult()
            return
        end if

        ! if both excits are in the same GAS, we could also
        ! have picked them the other way around
        if (srcGAS(1) == srcGAS(2)) then
            pgen_pick2 = &
                this%get_pgen_pick_weighted_hole(nI, Excite_2_t(src(1), tgt(2), src(2), tgt(1))) &
                * this%get_pgen_pick_hole_from_active_space(ilutI, srcGAS(2), get_spin(tgt(2)))
            pgen = pgen * (pgen_pick1 + pgen_pick2)
        else
            pgen = pgen * pgen_pick1
        end if

        call make_double(nI, nJ, elecs(1), elecs(2), tgt(1), tgt(2), ex, par)
        ilutJ = ilutI
        clr_orb(ilutJ, src(1))
        clr_orb(ilutJ, src(2))
        set_orb(ilutJ, tgt(1))
        set_orb(ilutJ, tgt(2))
    contains

        subroutine zeroResult()
            integer :: src_copy(2)

            pgen = pgen * pgen_pick1
            src_copy(:) = src(:)
            call sort(src_copy)
            ex(1, :) = src_copy
            ex(2, :) = tgt
            nJ(1) = 0
            ilutJ = 0_n_int
        end subroutine zeroResult

    end subroutine generate_nGAS_double