determine_optimal_time_step Function

    function determine_optimal_time_step(time_step_death) result(time_step)
        use SystemData, only: nel, bhub, uhub, t_new_real_space_hubbard, &
                              t_tJ_model, t_heisenberg_model, exchange_j, &
                              nOccAlpha, nOccBeta, t_k_space_hubbard, omega, &
                              nbasis

        real(dp), optional, intent(out) :: time_step_death
        real(dp) :: time_step
        ! move this time-step determination to this routine for the real
        ! space hubbard to have it fully conained

        real(dp) :: p_elec, p_hole, mat_ele, max_diag

        ! determine the optimal hubbard time-step for an optimized
        ! hubbard excitation generation
        ! the first electron is chosen at random
        if (t_k_space_hubbard) then
            p_elec = 1.0_dp / real(nOccAlpha * nOccBeta, dp)
        else
            p_elec = 1.0_dp / real(nel, dp)
        end if

        ! and for a picked electron the possible neighbors are looked for
        ! open holes so the lowest probability is determined by the
        ! maximum numbers of connections
        ! do this in general in the same way for all types of
        ! lattice models. thats not exact, but a good enough estimate
        if (t_k_space_hubbard) then
            p_hole = 2.0_dp / real(nbasis - nel, dp)
        else
            p_hole = 1.0_dp / real(lat%get_nconnect_max(), dp)
        end if

        if (t_new_real_space_hubbard) then

            mat_ele = real(abs(bhub), dp)

        else if (t_tJ_model) then

            ! for the t-J take the maximum of hopping or exchange
            mat_ele = real(max(abs(bhub), abs(exchange_j)), dp)

        else if (t_heisenberg_model) then

            mat_ele = real(abs(exchange_j), dp)

        else if (t_k_space_hubbard) then
            mat_ele = abs(real(uhub, dp) / real(omega, dp))

        end if

        ! so the time-step is
        time_step = p_elec * p_hole / mat_ele

        if (present(time_step_death)) then
            ! the maximum death contribution is max(H_ii) - shift
            ! and the maximum diagonal matrix element we know it is
            ! the maximum possible number of doubly occupied sites times U
            ! but this might be a too hard limit, since especially for high U
            ! such an amount of double occupations is probably never reached
            ! and the shift must also be included in the calculation..
            if (t_new_real_space_hubbard) then
                max_diag = real(abs(uhub) * min(nOccAlpha, nOccBeta), dp)
            else if (t_tJ_model) then
                print *, "todo!"
            else if (t_heisenberg_model) then
                ! we have to find the maximum number of non-frustrated
                ! nearest neighbor spins!
                print *, "todo: "
            else if (t_k_space_hubbard) then
                print *, "todo"

            end if

            time_step_death = 1.0_dp / max_diag
        end if

    end function determine_optimal_time_step