hubbard_find_tau_from_refdet_conn Subroutine

    subroutine hubbard_find_tau_from_refdet_conn()

        ! Routine to find an upper bound to tau, by consideration of the
        ! singles and doubles connected to the reference determinant
        !
        ! Obviously, this make assumptions about the possible range of pgen,
        ! so may actually give a tau that is too SMALL for the latest
        ! excitation generators, which is exciting!

        character(len=*), parameter :: this_routine = "find_tau_from_refdet_conn"
        integer :: ex(2, maxExcit), ic, nJ(nel), n_excits, i, ex_3(2, 3)
        real(dp) :: nAddFac, MagHel, pGen, pGenFac
        logical :: tParity
        integer(n_int), allocatable :: det_list(:, :)
        real(dp) :: new_tau

        ASSERT(t_k_space_hubbard)
        ASSERT(.not. tGUGA)

        ! NOTE: test if the new real-space implementation works with this
        ! function! maybe i also have to use a specific routine for this !
        ! since it might be necessary in the transcorrelated approach to
        ! the real-space hubbard

        new_tau = huge(new_tau)

        nAddFac = MaxWalkerBloom

        if (tHPHF) then
            call Stop_All(this_routine, &
                          "not yet implemented with HPHF, since gen_all_excits not atapted to it!")
        end if

        call gen_all_excits_k_space_hubbard(ProjEDet(:, 1), n_excits, det_list)

        ! now loop over all of them and determine the worst case H_ij/pgen ratio
        do i = 1, n_excits
            call decode_bit_det(nJ, det_list(:, i))
            ! i have to take the right direction in the case of the
            ! transcorrelated, due to non-hermiticity..
            ic = FindBitExcitlevel(det_list(:, i), ilutRef(:, 1))
            ASSERT(ic == 2 .or. ic == 3)
            if (ic == 2) then
                call GetBitExcitation(ilutRef(:, 1), det_list(:, i), ex, tParity)
            else if (ic == 3) then
                call GetBitExcitation(ilutRef(:, 1), det_list(:, i), ex_3, tParity)
            end if

            MagHel = abs(get_helement_lattice(nJ, ProjEDet(:, 1)))
            ! and also get the generation probability
            if (t_trans_corr_2body) then
                if (t_uniform_excits) then
                    ! i have to setup pDoubles and the other quantities
                    ! before i call this functionality!
                    pgen = calc_pgen_k_space_hubbard_uniform_transcorr(ex_3, ic)
                else
                    pgen = calc_pgen_k_space_hubbard_transcorr( &
                           ProjEDet(:, 1), ilutRef(:, 1), ex_3, ic)
                end if
            else
                pgen = calc_pgen_k_space_hubbard( &
                       ProjEDet(:, 1), ilutRef(:, 1), ex, ic)
            end if

            if (MagHel > EPS) then
                pGenFac = pgen * nAddFac / MagHel
                new_tau = clamp(&
                    merge(new_tau, min(pGenFac, new_tau), near_zero(pGenFac)), &
                    min_tau, max_tau)
            end if
        end do

        if (new_tau > 0.075_dp) then
            new_tau = clamp(0.075_dp, min_tau, max_tau)
            write(stdout, "(A,F8.5,A)") "Small system. Setting initial timestep to be ", Tau, " although this &
                                            &may be inappropriate. Care needed"
        end if
        call assign_value_to_tau(new_tau, this_routine)
    end subroutine hubbard_find_tau_from_refdet_conn