gen_excit_impurity_model Subroutine

    subroutine gen_excit_impurity_model(nI, ilut, nJ, ilutnJ, exFlag, IC, ExcitMat, &
                                        tParity, pGen, HElGen, store, part_type)
        use FciMCData, only: excit_gen_store_type
        ! generate a random excitation for an impurity-type system, that is,
        ! there is significant part of the system which is non-interacting (called bath)
        ! for the one-electron integrals, we have the connections cached
        integer, intent(in) :: nI(nel)
        integer, intent(in) :: exFlag
        integer(n_int), intent(in) :: iLut(0:NIfTot)
        integer, intent(out) :: nJ(nel), IC, ExcitMat(2, maxExcit)
        logical, intent(out) :: tParity
        type(excit_gen_store_type), intent(inout), target :: store
        real(dp), intent(out) :: pgen
        !type(excit_gen_store_type), intent(inout), target :: store

        integer(n_int), intent(out) :: ilutnJ(0:NIfTot)
        HElement_t(dp), intent(out) :: HElGen

        integer, intent(in), optional :: part_type

        real(dp) :: r
        integer :: nImpEls
        character(*), parameter :: this_routine = 'gen_rand_excit'

        unused_var(store)
        unused_var(part_type)
        unused_var(exFlag)

        HElGen = HEl_zero

        ! first, determine if a single or a double is created
        ! if there are enoguh holes and electrons in the impurity, pick randomly
        nImpEls = binary_search_first_ge(nI, nImp + 1) - 1
        if (nImpEls > 1 .and. (nImp - nImpEls) > 1) then
            r = genrand_real2_dSFMT()
            if (r < pSingles) then
                IC = 1
                pGen = pSingles
                call generate_imp_single_excitation(nI, ilut, nJ, ilutnJ, ExcitMat, tParity, pGen)
            else
                IC = 2
                pGen = 1.0_dp - pSingles
                call generate_imp_double_excitation(nI, ilut, nJ, ilutnJ, ExcitMat, tParity, pGen)
            end if
        else
            ! else, we target a single
            IC = 1
            pGen = 1.0_dp
            call generate_imp_single_excitation(nI, ilut, nJ, ilutnJ, ExcitMat, tParity, pGen)
        end if

    end subroutine gen_excit_impurity_model