init_real_space_hubbard Subroutine

    subroutine init_real_space_hubbard()
        use SystemData, only: tExch, thub, treal
        ! routine, which does all of the necessary initialisation
        character(*), parameter :: this_routine = "init_real_space_hubbard"
        ! especially do some stop_all here so no wrong input is used
        real(dp) :: tau_opt

        real(dp), allocatable :: Works(:)
        real(dp) :: T(nbasis, nbasis)

        integer :: i, j, lwork
        real(dp) :: Energie_ground(nbasis)
        integer :: Info

        root_print "using new real-space hubbard implementation: "

        ! i do not need exchange integrals in the real-space hubbard model
        if (.not. t_trans_corr_hop) then
            tExch = .false.
        end if
        ! after the whole setup i can set thub to false or?
        thub = .false.
        ! and treal i can also set to false or?
        treal = .false.

        ! just to be save swithc of Brillouins
        tNoBrillouin = .true.
        tUseBrillouin = .false.

        lnosymmetry = .true.

        ! first assert all the right input!
        call check_real_space_hubbard_input()

        ! which stuff do i need to initialize here?
        ! for now also use the umat in the spin-dependent transcorr
        ! although it is not
        if (t_trans_corr_hop .or. t_spin_dependent_transcorr) then
            get_umat_el => get_umat_rs_hub_trans
        else
            get_umat_el => get_umat_el_hub
        end if

        ! also use the new lattice matrix elements

        ! i have to check if the lattice should be constructed from an fcidump
        ! or created internally..
        if (trim(adjustl(lattice_type)) == 'read') then
            ! then i have to construct tmat first
            call init_tmat()
            ! and then construct the lattice
            lat => lattice(lattice_type, length_x, length_y, length_z,.not. t_open_bc_x, &
                           .not. t_open_bc_y,.not. t_open_bc_z)
        else
            lat => lattice(lattice_type, length_x, length_y, length_z,.not. t_open_bc_x, &
                       .not. t_open_bc_y,.not. t_open_bc_z, 'real-space', t_bipartite_order = t_bipartite_order)

            ! if nbaiss was not yet provided:
            if (nbasis <= 0) then
                nbasis = 2 * lat%get_nsites()
            end if

            call init_tmat(lat)

        end if

        allocate(G1(nbasis))
        G1(1::2) = BasisFn(Symmetry(0), [0, 0, 0], -1, 0, 0)
        G1(2::2) = BasisFn(Symmetry(0), [0, 0, 0], 1, 0, 0)

        ! Ecore should default to 0, but be sure anyway!
        ecore = 0.0_dp

        if (t_trans_corr_hop) then
            ! we have double excitations with the hopping correlation
            ! but only anti-parallel excitations!
            if (allocated(pSinglesIn) .and. allocated(pDoublesIn)) then
                if (.not. (pSinglesIn + pDoublesIn .isclose. 1.0_dp)) then
                    call stop_all(this_routine, "pSinglesIn + pDoublesIn /= 1.0!")
                else
                    pSingles = pSinglesIn
                    pDoubles = pDoublesIn
                end if
            else if (allocated(pSinglesIn) .and. (.not. allocated(pDoublesIn))) then
                pSingles = pSinglesIn
                pDoubles = 1.0_dp - pSingles
            else if (allocated(pDoublesIn) .and. (.not. allocated(pSinglesIn))) then
                pDoubles = pDoublesIn
                pSingles = 1.0_dp - pDoubles

            ! For consistency pParallelIn should be taken as well or error out
            else
                pSingles = 0.8_dp
                pDoubles = 1.0_dp - pSingles
            end if
        else
            ! and i have to point to the new hubbard excitation generator
            pSingles = 1.0_dp
            pDoubles = 0.0_dp
        end if

        ! and i have to calculate the optimal time-step for the hubbard models.
        ! where i need the connectivity of the lattice i guess?
        if (t_trans_corr_hop .and. .not. tHPHF) then
            if (t_twisted_bc) then
                call stop_all(this_routine, &
                              "twisted BC + Transcorr not yet implemented!")
            end if
        end if
        ! i have to calculate the optimal time-step
        ! and maybe i have to be a bit more safe here and not be too near to
        ! the optimal time-step
        tau_opt = determine_optimal_time_step()
        if (tau < EPS) then
            root_print "setting time-step to optimally determined time-step: ", tau_opt
            root_print "times: ", lat_tau_factor
            call assign_value_to_tau(&
                clamp(lat_tau_factor * tau_opt, min_tau, max_tau), &
                'Initialization with optimal real-space Hubbard value')

        else
            root_print "optimal time-step would be: ", tau_opt
            root_print "but tau specified in input!"
        end if

        ! re-enable tau-search if we have transcorrelation
        if (.not. (t_trans_corr_2body .or. t_trans_corr .or. t_trans_corr_hop &
                   .or. t_spin_dependent_transcorr)) then

            if (tau_search_method /= possible_tau_search_methods%OFF) then
                call stop_all(this_routine, "tau-search should be switched off")
            end if
            t_scale_tau_to_death = .true.
        end if

        if (t_start_neel_state) then

            root_print "starting from the Neel state: "
            if (nel > nbasis / 2) then
                call stop_all(this_routine, &
                              "more than half-filling! does neel state make sense?")
            end if

        end if

        if (tHartreeGutzwiller) then
            do i=1, nbasis
                do j=1, nbasis
                    T(i, j) = real(getTMatEl(i, j), dp)
                end do
            end do

            lwork = max(1,3*nbasis-1)
            allocate(Works(lwork))

            if (.not. allocated(U_matrix)) allocate(U_matrix(nbasis, nbasis))
            ! Calculate eigenvectors of the T matrix
            ! This then is the U_matrix
            call dsyev('V', 'U', nbasis, T, nbasis, Energie_ground, Works, lwork, Info)
            do i=1, nbasis
                do j=1, nbasis
                    U_matrix(i,j) = T(i,j)
                end do
            end do
        end if

        ! i need to setup the necessary stuff for the new hopping
        ! transcorrelated real-space hubbard!
        call init_get_helement_hubbard()

    end subroutine init_real_space_hubbard