init_heisenberg_model Subroutine

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

        root_print "initialising Heisenberg model with "
        root_print "J: ", exchange_j

        thub = .false.
        treal = .false.

        call check_real_space_hubbard_input()

        if (t_trans_corr_2body) then
            if (t_trans_corr) then
                call stop_all(this_routine, &
                              "1-body transcorrelation not allowed in the heisenberg model!")
            end if
        end if

        get_umat_el => get_umat_el_heisenberg

        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)

        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

        ! and i have to calculate the optimal time-step for the hubbard models.
        ! where i need the connectivity of the lattice i guess?
        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

        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

    end subroutine init_heisenberg_model