setup_lattice_symmetry Subroutine

    subroutine setup_lattice_symmetry
        ! since i need it also in the real-space lattice for the
        ! hopping transcorrelation move the symmetry setup for the
        ! k-spae hubbard model into the lattice_mod
#ifdef DEBUG_
        character(*), parameter :: this_routine = "setup_lattice_symmetry"
#endif
        integer :: i, kmin(3), kmax(3), j, k_i(3), k, l, ind
        ASSERT(associated(lat))

        if (allocated(lat%k_to_sym)) deallocate(lat%k_to_sym)
        if (allocated(lat%sym_to_k)) deallocate(lat%sym_to_k)
        if (allocated(lat%mult_table)) deallocate(lat%mult_table)
        if (allocated(lat%inv_table)) deallocate(lat%inv_table)

        allocate(lat%sym_to_k(lat%get_nsites(), 3))
        allocate(lat%mult_table(lat%get_nsites(), lat%get_nsites()))
        allocate(lat%inv_table(lat%get_nsites()))

        ! i have to setup the symlabels first ofc..
        do i = 1, lat%get_nsites()
            ! and also just encode the symmetry labels as integers, instead of
            ! 2^(k-1), to be able to treat more than 64 orbitals (in the old
            ! implementation, an integer overflow happened in this case!)
            ind = get_spatial(brr(2 * i))

            call lat%set_sym(ind, i)

        end do

        kmin = 0
        kmax = 0
        do i = 1, lat%get_nsites()
            k_i = lat%get_k_vec(i)
            do j = 1, lat%get_ndim()
                if (k_i(j) < kmin(j)) kmin(j) = k_i(j)
                if (k_i(j) > kmax(j)) kmax(j) = k_i(j)
            end do
        end do

        allocate(lat%k_to_sym(kmin(1):kmax(1), kmin(2):kmax(2), kmin(3):kmax(3)))

        lat%k_to_sym = 0

        ! now find the inverses:
        do i = 1, lat%get_nsites()

            ! find the orbital of -k
            j = lat%get_orb_from_k_vec(-lat%get_k_vec(i))

            lat%inv_table(lat%get_sym(i)) = lat%get_sym(j)

            lat%sym_to_k(lat%get_sym(i), :) = lat%get_k_vec(i)

            k_i = lat%get_k_vec(i)

            lat%k_to_sym(k_i(1), k_i(2), k_i(3)) = lat%get_sym(i)

            ! and create the symmetry product of (i) with every other symmetry
            do k = 1, lat%get_nsites()
                ! i just have to add the momenta and map it to the first BZ
                l = lat%get_orb_from_k_vec(lat%get_k_vec(i) + lat%get_k_vec(k))

                lat%mult_table(lat%get_sym(i), lat%get_sym(k)) = lat%get_sym(l)

            end do
        end do

    end subroutine setup_lattice_symmetry