init_guga_tJ_model Subroutine

    subroutine init_guga_tJ_model
        ! write a specific subroutine to init the spin-free tJ model
        character(*), parameter :: this_routine = "init_guga_tJ_model"

        real(dp) :: tau_opt

        root_print "initializing the spin-free t-J model, with parameter: "
        root_print "t:", bhub
        root_print "J:", exchange_j

        tHub = .false.
        treal = .false.
        t_new_real_space_hubbard = .true.

        call check_real_space_hubbard_input()

        nSpatOrbs = nBasis .div. 2
        nSpatOrbs_i8 = int(nBasis .div. 2, int64)
        ij_pairs = fuse_symm_idx(nSpatOrbs_i8, nSpatOrbs_i8)

        pSingles = 0.1_dp
        pDoubles = 1.0_dp - pSingles

        get_umat_el => get_umat_heisenberg_spin_free

        if (trim(adjustl(lattice_type)) == 'read') then
            ! then i have to construct tmat first
            call stop_all(this_routine, "starting from fcidump not yet implemented!")
            ! 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
            ! otherwise i have to do it the other way around
            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, t_bipartite_order = t_bipartite_order)
            ! if nbasis was not yet provided:
            if (nbasis <= 0) then
                nbasis = 2 * lat%get_nsites()
            end if

            call init_tmat(lat)
            call init_spin_free_tmat(lat)
            call setup_exchange_matrix(lat)
            call setup_spin_free_exchange(lat)

        end if

        if (nel >= nbasis / 2) then
            call stop_all(this_routine, &
                          " too many electrons for the tJ model! nel >= nbasis/2")
        end if

        ! and also check the double occupancy in the starting det,
        ! no double occupancy allowed!
        if (count_double_orbs(ilutRef) > 0) then
            call stop_all(this_routine, &
                          "incorrect starting state for tJ model: there is an doubly occupied site!")
        end if

        ! i guess i have to setup G1 also.. argh.. i hate this!
        allocate(G1(nbasis))
        G1(1:nbasis - 1:2)%ms = -1
        G1(2:nbasis:2)%ms = 1

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

        tau_opt = determine_optimal_time_step()
        if (tau < EPS) then
            call assign_value_to_tau(&
                lat_tau_factor * tau_opt, &
                'Initialization with optimal tau value')
        else
            root_print "optimal time-step would be: ", tau_opt
            root_print "but tau specified in input!"
        end if

    end subroutine init_guga_tJ_model