subroutine setup_k_space_hub_sym(in_lat)
class(lattice), intent(in), optional :: in_lat
character(*), parameter :: this_routine = "setup_k_space_hub_sym"
INTEGER :: i, j, SymInd, Lab, spin, sym0
INTEGER, ALLOCATABLE :: Temp(:)
! These are for the hubbard and UEG model look-up table
type(Symmetry) :: SymProduct, SymI, SymJ
if (present(in_lat)) then
if (.not. associated(G1)) then
call setup_g1(in_lat)
end if
if (all(nBasisMax == 0)) then
call setup_nbasismax(in_lat)
end if
! do only the necessary setup here!
! this is essentially from symrandexcit2.F90 SpinOrbSymSetup()
ElecPairs = (NEl * (NEl - 1)) / 2
MaxABPairs = (nBasis * (nBasis - 1) / 2)
ScratchSize = 2 * nSymLabels
if (allocated(SpinOrbSymLabel)) deallocate(SpinOrbSymLabel)
allocate(SpinOrbSymLabel(nBasis))
do i = 1, nBasis
!This ensures that the symmetry labels go from 0 -> nSymLabels-1
SpinOrbSymLabel(i) = SymClasses(((i + 1) / 2)) - 1
end do
#ifdef DEBUG_
write(stdout, *) "SpinOrbSymLabel: "
do i = 1, nBasis
write(stdout, *) i, SpinOrbSymLabel(i)
end do
#endif
if (allocated(SymTableLabels)) deallocate(SymTableLabels)
allocate(SymTableLabels(0:nSymLabels - 1, 0:nSymLabels - 1))
SymTableLabels(:, :) = -9000 !To make it easier to track bugs
do i = 0, nSymLabels - 1
do j = 0, i
SymI = SymLabels(i + 1) !Convert to the other symlabel convention to use SymLabels -
!TODO: I will fix this to make them consistent when working (ghb24)!
SymJ = SymLabels(j + 1)
SymProduct = SymProd(SymI, SymJ)
!Now, we need to find the label according to this symmetry!
!Run through all symmetries to make working (this could be far more efficient, but its only once, so sod it...
do Lab = 1, nSymLabels
if (SymLabels(Lab)%S == SymProduct%S) then
EXIT
end if
end do
if (Lab == nSymLabels + 1) then
call stop_all("SpinOrbSymSetup", "Cannot find symmetry label")
end if
SymTableLabels(i, j) = Lab - 1
SymTableLabels(j, i) = Lab - 1
end do
end do
#ifdef DEBUG_
write(stdout, *) "SymTable:"
do i = 0, nSymLabels - 1
do j = 0, nSymLabels - 1
write(stdout, "(I6)", advance='no') SymTableLabels(i, j)
end do
write(stdout, *) ""
end do
#endif
!SymInvLabel takes the label (0 -> nSymLabels-1) of a spin orbital, and returns the inverse symmetry label, suitable for
!use in ClassCountInd.
if (allocated(SymInvLabel)) deallocate(SymInvLabel)
allocate(SymInvLabel(0:nSymLabels - 1))
SymInvLabel = -999
! Dongxia changes the gamma point away from center.
! SDS: Provide a default sym0 for cases where this doesn't apply
sym0 = 0
do i = 1, nsymlabels
if (symlabels(i)%s == 0) sym0 = i - 1
end do
do i = 0, nSymLabels - 1
! Change the sym label back to the representation used by the rest
! of the code, use SymConjTab, then change back to other rep of
! labels SymConjTab only works when all irreps are self-inverse.
! Therefore, instead, we will calculate the inverses by just
! finding the symmetry which will give A1.
do j = 0, nSymLabels - 1
! Run through all labels to find what gives totally symmetric
! rep
if (SymTableLabels(i, j) == sym0) then
if (SymInvLabel(i) /= -999) then
write(stdout, *) "SymLabel: ", i
call stop_all(this_routine, &
"Multiple inverse irreps found - error")
end if
! This is the inverse
SymInvLabel(i) = j
end if
end do
if (SymInvLabel(i) == -999) then
write(stdout, *) "SymLabel: ", i
call stop_all(this_routine, "No inverse symmetry found - error")
end if
end do
#ifdef DEBUG_
write(stdout, *) "SymInvLabel: "
do i = 0, nSymLabels - 1
write(stdout, *) i, SymInvLabel(i)
end do
#endif
!SymLabelList2 and SymLabelCounts2 are now organised differently, so that it is more efficient, and easier to add new symmetries.
!SymLabelCounts is of size (2,ScratchSize), where 1,x gives the index in SymlabelList2 where the orbitals of symmetry x start.
!SymLabelCounts(2,x) tells us the number of orbitals of spin & symmetry x there are.
!Therefore, if you want to run over all orbitals of a specific symmetry, you want to run over
!SymLabelList from SymLabelCounts(1,sym) to SymLabelCounts(1,sym)+SymLabelCounts(2,sym)-1
if (allocated(SymLabelList2)) deallocate(SymLabelList2)
if (allocated(SymLabelCounts2)) deallocate(SymLabelCounts2)
allocate(SymLabelList2(nBasis))
allocate(SymLabelCounts2(2, ScratchSize))
SymLabelList2(:) = 0 !Indices: spin-orbital number
SymLabelCounts2(:, :) = 0 !Indices: index/Number , symmetry(inc. spin)
allocate(Temp(ScratchSize))
do j = 1, nBasis
IF (G1(j)%Ms == 1) THEN
Spin = 1
ELSE
Spin = 2
end if
SymInd = ClassCountInd(Spin, SpinOrbSymLabel(j), G1(j)%Ml)
SymLabelCounts2(2, SymInd) = SymLabelCounts2(2, SymInd) + 1
end do
SymLabelCounts2(1, 1) = 1
do j = 2, ScratchSize
SymLabelCounts2(1, j) = SymLabelCounts2(1, j - 1) + SymLabelCounts2(2, j - 1)
end do
Temp(:) = SymLabelCounts2(1, :)
do j = 1, nBasis
IF (G1(j)%Ms == 1) THEN
Spin = 1
ELSE
Spin = 2
end if
SymInd = ClassCountInd(Spin, SpinOrbSymLabel(j), G1(j)%Ml)
SymLabelList2(Temp(SymInd)) = j
Temp(SymInd) = Temp(SymInd) + 1
end do
Deallocate(Temp)
if (allocated(OrbClassCount)) deallocate(OrbClassCount)
allocate(OrbClassCount(ScratchSize))
OrbClassCount(:) = 0
do i = 1, nBasis
IF (G1(i)%Ms == 1) THEN
OrbClassCount(ClassCountInd(1, SpinOrbSymLabel(i), G1(i)%Ml)) = &
& OrbClassCount(ClassCountInd(1, SpinOrbSymLabel(i), G1(i)%Ml)) + 1
ELSE
OrbClassCount(ClassCountInd(2, SpinOrbSymLabel(i), G1(i)%Ml)) = &
& OrbClassCount(ClassCountInd(2, SpinOrbSymLabel(i), G1(i)%Ml)) + 1
end if
end do
call setup_kPointToBasisFn(in_lat)
call gensymstatepairs(nbasis / 2, .false.)
else
call Stop_All(this_routine, "not yet implemented")
end if
end subroutine setup_k_space_hub_sym