init_guga_heisenberg_model Subroutine

    subroutine init_guga_heisenberg_model()
        character(*), parameter :: this_routine = "init_guga_heisenberg_model"
        real(dp) :: tau_opt

        root_print "initializing spin-free Heisenberg model, with "
        root_print "J: ", exchange_j

!         call init_guga()

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

        call check_real_space_hubbard_input()

        pSingles = 0.0_dp
        pDoubles = 1.0_dp

        get_umat_el => get_umat_heisenberg_spin_free

        if (associated(tmat2d)) deallocate(tmat2d)
        allocate(tmat2d(nbasis, nbasis), source=h_cast(0.0_dp))

        if (trim(adjustl(lattice_type)) == 'read') then
            ! then i have to construct tmat first
            ! no need for tmat in the heisenberg model
            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
            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 setup_exchange_matrix(lat)
            call setup_spin_free_exchange(lat)

        end if

        if (nel /= nbasis / 2) then
            call stop_all(this_routine, &
                          "heisenberg model need half filling nel == nbasis/2")
        end if

        if (count_double_orbs(ilutref) > 0) then
            call stop_all(this_routine, &
                          " no double occupancies allowed in the heisenberg model")
        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

        ! TODO: maybe I need the tau-search for the GUGA..
        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

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

    end subroutine init_guga_heisenberg_model