gen_excit_back_spawn Subroutine

    subroutine gen_excit_back_spawn(nI, ilutI, nJ, ilutJ, exFlag, ic, &
                                    ExcitMat, tParity, pgen, HelGen, store, part_type)
        integer, intent(in) :: nI(nel), exFlag
        integer(n_int), intent(in) :: ilutI(0:niftot)
        integer, intent(out) :: nJ(nel), ic, ExcitMat(2, maxExcit)
        integer(n_int), intent(out) :: ilutJ(0:niftot)
        logical, intent(out) :: tParity
        real(dp), intent(out) :: pgen
        HElement_t(dp), intent(out) :: HElGen
        type(excit_gen_store_type), intent(inout), target :: store
        integer, intent(in), optional :: part_type
        character(*), parameter :: this_routine = "gen_excit_back_spawn"

#ifdef DEBUG_
        HElement_t(dp) :: temp_hel
        real(dp) :: pgen2
#endif
        unused_var(exFlag); unused_var(store)

        HElGen = 0.0_dp
        ! check the non-initiator criteria beforehand
        ! i also have to consider that back-spawn gets turned on later on
        ! so i have to check if back-spawn is active already or not..

        if ((t_back_spawn_flex .or. t_back_spawn) .and. &
            .not. test_flag(ilutI, get_initiator_flag(part_type))) then

            ! otherwise use custom made ones
            if (genrand_real2_dSFMT() < pSingles) then

                ic = 1
                call gen_single_back_spawn(nI, ilutI, part_type, nJ, ilutJ, ExcitMat, &
                                           tParity, pgen)
                pgen = pgen * pSingles

            else

                ic = 2
                call gen_double_back_spawn(nI, ilutI, part_type, nJ, ilutJ, ExcitMat, &
                                           tParity, pgen)
                pgen = pgen * pDoubles

            end if

#ifdef DEBUG_
            if (.not. IsNullDet(nJ)) then
                pgen2 = calc_pgen_back_spawn(nI, ilutI, ExcitMat, ic, part_type)
                if (abs(pgen - pgen2) > 1.0e-6_dp) then
                    if (tHPHF) then
                        print *, "due to circular dependence, no matrix element calc possible!"
!                         temp_hel = hphf_off_diag_helement(nI,nJ,ilutI,ilutJ)
                        temp_hel = 0.0_dp
                    else
                        temp_hel = get_helement(nI, nJ, ilutI, ilutJ)
                    end if

                    write(stdout, *) 'Calculated and actual pgens differ. for non-initiator'
                    write(stdout, *) 'This will break HPHF calculations'
                    write(stdout, *) "reference determinant: "
                    call write_det(stdout, projedet(:, part_type_to_run(part_type)), .true.)
                    call write_det(stdout, nI, .false.)
                    write(stdout, '(" --> ")', advance='no')
                    call write_det(stdout, nJ, .true.)
                    write(stdout, *) 'Excitation matrix: ', ExcitMat(1, 1:ic), '-->', &
                        ExcitMat(2, 1:ic)
                    write(stdout, *) 'Generated pGen:  ', pgen
                    write(stdout, *) 'Calculated pGen: ', pgen2
                    write(stdout, *) 'matrix element: ', temp_hel
                    call stop_all(this_routine, "Invalid pGen")
                end if
            end if
#endif
        else

            ! do the "normal" excitation type if it is an initiator
            if (genrand_real2_dSFMT() < pSingles) then

                ic = 1
                call gen_single_4ind_ex(nI, ilutI, nJ, ilutJ, ExcitMat, &
                                        tParity, pGen)
                pgen = pgen * pSingles

            else

                ic = 2
                call gen_double_4ind_ex2(nI, ilutI, nJ, ilutJ, ExcitMat, &
                                         tParity, pGen)
                pgen = pgen * pDoubles

            end if
#ifdef DEBUG_
            if (.not. IsNullDet(nJ)) then
                pgen2 = calc_pgen_4ind_weighted2(nI, ilutI, ExcitMat, ic)
                if (abs(pgen - pgen2) > 1.0e-6_dp) then
                    if (tHPHF) then
                        print *, "due to circular dependence, no matrix element calc possible!"
!                         temp_hel = hphf_off_diag_helement(nI,nJ,ilutI,ilutJ)
                        temp_hel = 0.0_dp
                    else
                        temp_hel = get_helement(nI, nJ, ilutI, ilutJ)
                    end if

                    write(stdout, *) 'Calculated and actual pgens differ. initiator'
                    write(stdout, *) 'This will break HPHF calculations'
                    call write_det(stdout, nI, .false.)
                    write(stdout, '(" --> ")', advance='no')
                    call write_det(stdout, nJ, .true.)
                    write(stdout, *) 'Excitation matrix: ', ExcitMat(1, 1:ic), '-->', &
                        ExcitMat(2, 1:ic)
                    write(stdout, *) 'Generated pGen:  ', pgen
                    write(stdout, *) 'Calculated pGen: ', pgen2
                    write(stdout, *) 'matrix element: ', temp_hel
                    call stop_all(this_routine, "Invalid pGen")
                end if
            end if
#endif
        end if
    end subroutine gen_excit_back_spawn