generate_imp_double_excitation Subroutine

    subroutine generate_imp_double_excitation(nI, ilut, nJ, ilutnJ, ex, tParity, pGen)
        integer, intent(in) :: nI(nel)
        integer, intent(out) :: nJ(nel), ex(2, 2)
        logical, intent(out) :: tParity
        real(dp), intent(inout) :: pGen
        integer(n_int), intent(in) :: ilut(0:NIfTot)
        integer(n_int), intent(out) :: ilutnJ(0:NIfTot)

        integer :: i, j, k, l, nOccImp
        procedure(sumFunc), pointer :: p_getImpImpMatElDoubs

        ! pick the first electron from the impurity at random
        i = pick_random_occ_impurity(nI, nOccImp)
        ! If an excitation is not possible, abort
        if (nOccImp == 1 .or. (nImp - nOccImp) < 2) then
            call invalidate()
            return
        end if
        ! then pick the second electron randomly from the remaining ones
        j = pick_random_occ_impurity(nI, nOccImp)
        if (i == j) then
            call invalidate()
            return
        end if
        ! pgen is 1.0/nOccImp per orbital, but they can be picked in either order
        pgen = pgen * 2.0 / (nOccImp**2)

        ! we now have two electrons from the impurity
        ! so get two holes weighted with the matrix elements
        k = pick_random_unocc_impurity(nI, G1(i)%Ms, pgen)
        l = pick_random_unocc_impurity(nI, G1(j)%Ms, pgen)
        ! Since the spin is fixed, the holes cannot be picked in any order if they have
        ! different spin, only if it is the same
        if (G1(i)%Ms == G1(j)%Ms) then
            if (k == l) then
                call invalidate()
                return
            end if
            ! Now they could have been chosen in any order
            pgen = pgen * 2.0_dp
        end if
        ! and write the resulting det to output
        call make_double(nI, nJ, i, j, k, l, ex, tparity)
        call assign_output_ilut(ilut, ilutnJ, i, k, j, l)

    contains

        subroutine invalidate()
            nJ = 0
            ilutnJ = 0
            pgen = 1.0
        end subroutine invalidate
    end subroutine generate_imp_double_excitation