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