setup_symmetry_table Subroutine

    subroutine setup_symmetry_table()
        ! implement a new symmetry setup to decouple it from the
        ! old hubbard.F code..

        character(*), parameter :: this_routine = "setup_symmetry_table"

        integer :: i, j, k, l, k_i(3), ind, kmin(3), kmax(3)

        ! the only problem could be that we reorderd the orbitals already or?
        ! so G1 has a different ordering than just 1, nBasis/2...
        ! yes it really is reordered already! hm..
        ASSERT(associated(lat))
        ASSERT(associated(G1))

        nsym = nBasis / 2
        nSymLabels = nsym

        ! copy the output from the old hubbard code:
        write(stdout, "(A,I3,A)") "Generating abelian symmetry table with", &
            nsym, " generators for Hubbard momentum"
        if (allocated(SymLabels)) then
            write(stdout, '(a/a)') &
                'Warning: symmetry info already allocated.', &
                'Deallocating and reallocating.'
            deallocate(SymLabels)
            call LogMemDealloc(this_routine, tagSymLabels)
        end if
        allocate(SymLabels(nSym))
        call LogMemAlloc('SymLabels', nSym, SymmetrySize, this_routine, tagSymLabels)
        if (associated(SymClasses)) then
            deallocate(SymClasses)
            call LogMemDealloc(this_routine, tagSymClasses)
        end if
        ! [W.D]
        ! why is symclasses allocated to nBasis? it is actually only used
        ! and filled up to nBasis/2, also in the rest of the code..
        allocate(SymClasses(nBasis / 2))
        call LogMemAlloc('SymClasses', nBasis, 4, this_routine, tagSymClasses)

        ! for some strange reason the (0,0,0) k-vector is treated
        ! special in the old implementation.. since it is its own
        ! symmetry inverse.. but this might be not the case in the
        ! general lattices or? there can be other states which are
        ! also its own inverse or?

        ! so try to change that here..
        ! ok i should still treat the gamma point special and set its
        ! sym label to 1 and maybe also order the symmetries by energy?

        if (allocated(SymTable)) then
            deallocate(SymTable)
            call LogMemDealloc(this_routine, tagSymTable)
        end if
        allocate(SymTable(nSym, nSym))
        call LogMemAlloc('SymTable', nSym**2, SymmetrySize, this_routine, tagSymTable)
        if (allocated(SymConjTab)) then
            deallocate(SymConjTab)
            call LogMemDealloc(this_routine, tagSymConjTab)
        end if
        allocate(SymConjTab(nSym))
        call LogMemAlloc('SymConjTable', nSym, 4, this_routine, tagSymConjTab)

        SYMTABLE = Symmetry(0)

        tAbelian = .false.

        ! also setup the stuff contained in the lattice class.
        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()
            ind = get_spatial(brr(2 * i))
            SymClasses(ind) = i
            ! 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!)
            SymLabels(ind)%s = i
            ! i also need it in G1:
            ! is G1 already ordered by energy??

            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()
            G1(2 * i - 1)%Sym = SymLabels(i)
            G1(2 * i)%Sym = SymLabels(i)

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

            ! since i have a linear encoding of the symmetries i do not need
            ! to use SymClasses here..
            SymConjTab(SymClasses(i)) = SymClasses(j)

            lat%inv_table(SymClasses(i)) = SymClasses(j)
            ! maybe also store the k_inverse of a orbital..

            lat%sym_to_k(SymClasses(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)) = SymClasses(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))

                SymTable(SymClasses(i), SymClasses(k)) = SymLabels(l)

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

            end do
        end do
#ifdef DEBUG_
        write(stdout, *) "Symmetry, Symmetry Conjugate"
        do i = 1, lat%get_nsites()
            print *, i, SymConjTab(i)
        end do

        print *, "lat%sym_to_k: "
        do i = 1, lat%get_nsites()
            print *, i, "|", SymClasses(i), "|", lat%sym_to_k(i, :)
        end do

        print *, "lat%inv_table: "
        do i = 1, lat%get_nsites()
            print *, i, "|", SymClasses(i), "|", lat%inv_table(i)
        end do

        print *, "lat%mult_table: "
        do i = 1, lat%get_nsites()
            print *, lat%mult_table(i, :)
        end do

        print *, "lat%k_to_sym: "
        do i = 1, lat%get_nsites()
            k_i = lat%get_k_vec(i)
            print *, k_i, "|", lat%k_to_sym(k_i(1), k_i(2), k_i(3))
        end do

#endif

    end subroutine setup_symmetry_table