module sym_mod use constants, only: dp, int64, sizeof_int, stdout use SymExcitDataMod, only: SymTableLabels, SymLabelList2, SymLabelCounts2, & OrbClassCount use SystemData, only: tKpntSym, tNoSymGenRandExcits, tHub, t_new_hubbard, & t_k_space_hubbard, Symmetry, BasisFN, SymmetrySize, & tUEG, treal, arr, brr, tSymIgnoreEnergies, & nBasis, G1, tKPntSym, tFixLz, lNoSymmetry, nBasisMax, nel, & NullBasisFn, Tperiodicinmom use lattice_mod, only: lat ! use SymData, only: SymTable,nProp,tAbelian,TwoCycleSymGens, SymConjTab, nSym, & ! SymClasses, nSymLabels, SymLabels, tagSymLabels, tagSymClasses, & ! tagSymConjTab, SymReps, SymClasses2, tagSymClasses2, & ! tagSymReps use SymData use util_mod, only: int_fmt, stop_all, neci_flush use global_utilities use sort_mod implicit none contains ! Symmetries can be unset(=0), or have bits set for the different irreps included ! bit 0 corresponds to totally symmetric. ! JSS: if Abelian symmetry, then ! symmetry=0 corresponds to the symmetric representation. It is not ! possible to unset symmetries for k-point jobs. ! To multiply symmetries, for each bit set in each of the two symmetries, we look up the ! product in the symmetry table, and OR that with the product. ! JSS: In Abelian (k-point) symmetry, a representation can be described by "quantum" numbers. ! The multiplication of two irreps is equal to the sum of such vectors ! describing the irreps, subject to a modulo of the periodic conditions ! due to the size of the symmetry cell (which corresponds to the k-point ! mesh). ! Note that many symmetry parameters for CPMD-NECI jobs are set in ! kpntrep.F in CPMD source. elemental FUNCTION SYMPROD(ISYM1, ISYM2) TYPE(Symmetry), intent(in) :: ISYM1, ISYM2 TYPE(Symmetry) SYMPROD TYPE(Symmetry) IS1, IS2 INTEGER I, J, Abel1(3), Abel2(3) character(*), parameter :: this_routine = 'SYMPROD' if (TAbelian) then IF (TwoCycleSymGens) THEN !For molecular systems, we can only have a maximum of 8 irreps, and so can do a simple !xor to get the symmetry. SymProd%s = IEOR(ISym1%s, ISym2%s) ELSE call DecomposeAbelianSym(ISym1%s, Abel1) call DecomposeAbelianSym(ISym2%s, Abel2) !Slightly faster when calling a lot to do it in an array operation Abel1 = modulo(Abel1 + Abel2, NProp) SymProd%s = ComposeAbelianSym(Abel1) end if else IF (ISYM1%s == 0 .OR. ISYM2%s == 0) THEN SYMPROD%s = 0 RETURN end if ! fix if all symmetries are set to 1 ! if (isym1%s == isym2%s .and. isym1%s == 1) then ! symprod%s = 1 ! return ! end if ! change for the real-space hubbard IF (.not. allocated(SYMTABLE)) call stop_all(this_routine, 'SYMMETRY TABLE NOT ALLOCATED') IS1 = ISYM1 I = 1 SYMPROD%s = 0 if (t_k_space_hubbard) then symprod = SYMTABLE(isym1%s, isym2%s) else DO WHILE (IS1%s /= 0) IF (BTEST(IS1%s, 0)) THEN IS2 = ISYM2 J = 1 DO WHILE (IS2%s /= 0) IF (BTEST(IS2%s, 0)) THEN SYMPROD%s = IOR(SYMPROD%s, SYMTABLE(I, J)%s) end if ! RSHIFT(,1) IS2%s = ISHFT(IS2%s, -1) J = J + 1 end do end if !RSHIFT(,1) IS1%s = ISHFT(IS1%s, -1) I = I + 1 end do end if end if RETURN END FUNCTION SYMPROD elemental FUNCTION SymConj(s2) type(symmetry), intent(in) :: s2 TYPE(Symmetry) s, SymConj INTEGER i, AbelConj(3) if (TAbelian) then IF (TwoCycleSymGens) THEN !For molecular systems, we only have symmetry generators which are two cycles, and so the inverse of !a symmetry is simply itself. SymConj%s = s2%s RETURN ELSE ! K-point symmetry has k_-i=k_i. We store k-vectors from 0 ! to N rather than -N/2 to N. Hence k_-i=mod(-k_i+N,N) for ! each component of the vector. ! This also works for abelian groups made out of symmetry ! generators which are 2-cycles call DecomposeAbelianSym(s2%s, AbelConj) do i = 1, 3 AbelConj(i) = modulo(-AbelConj(i) + NProp(i), NProp(i)) end do SymConj%s = ComposeAbelianSym(AbelConj) end if else SymConj%s = 0 s = s2 i = 1 DO WHILE (s%s /= 0) !LSHIFT(,) ! 1_8 is 1 in integer(int64): we need to have consistent ! kinds for bit-wise operations. IF (BTEST(s%s, 0)) SymConj%s = IOR(SymConj%s, ISHFT(1_int64, int(SymConjTab(I) - 1, int64))) !RSHIFT(,1) s%s = ISHFT(s%s, -1) i = i + 1 end do end if END FUNCTION SYMCONJ SUBROUTINE WRITESYMTABLE(IUNIT) IMPLICIT NONE INTEGER IUNIT, I, J DO I = 0, NSYM - 1 DO J = 0, NSYM - 1 CALL WRITESYM(IUNIT, SYMTABLE(I + 1, J + 1), .FALSE.) end do write(IUNIT, *) end do IF (NSYM == 0) THEN write(stdout, *) "No Symmetry table found." end if END SUBROUTINE WRITESYMTABLE LOGICAL elemental FUNCTION LSYMSYM(SYM) Type(Symmetry), intent(in) :: SYM if (TAbelian) then LSymSym = Sym%s == 0 else LSYMSYM = (SYM%s == 0 .OR. BTEST(SYM%s, 0)) end if RETURN END FUNCTION LSYMSYM ! Generate a symmetry table with molpro symmetries. ! irreps are simply made out of up to three generators. ! MOLPRO classifies irreps with bits 0-3 corresponding to those generators. ! symmetry products are merely exclusive ors of the molpro irrep numbers ! We set each of the MOLPRO irreps to a bit in our symmetry specifier. ! A1 corresponds to bit 0 (i.e. irrep 1) SUBROUTINE GENMOLPSYMTABLE(NSYMMAX, G1, NBASIS) IMPLICIT NONE INTEGER NSYMMAX, nSymGen INTEGER I, ILABEL TYPE(BasisFN) G1(*) INTEGER NBASIS character(*), parameter :: this_routine = 'GenMolPSymTable' TAbelian = .true. nSymGen = INT(log(NSYMMAX + 0.0_dp) / log(2.0_dp) + .4_dp) write(stdout, "(A,I3,A)") " Generating abelian symmetry table with", & nSymGen, " generators" write(stdout, '(A,'//int_fmt(nSymMax)//')') & " Number of symmetry classes: ", nSymMax ! We actually use momentum conservation directly for the UEG and ! Hubbard mode so just fake the symmetry information here. ! WARNING: do *not* use SymConj etc for these systems without fixing ! this---functions which rely upon the wavevectors being encoded into ! a symmetry integer will not work. if (TwoCycleSymGens .or. tUEG .or. tHUB) then ! Set propogation information. ! If not TwoCycleSymGens we assume the user has already ! done so... nprop = 1 nprop(1:min(nSymGen, 3)) = 2 end if ! Now generate a list of sym labels. NSYMLABELS = NSYMMAX if (.not. allocated(SymLabels)) allocate(SymLabels(nSymLabels)) call LogMemAlloc('SymLabels', nSymLabels, SymmetrySize, this_routine, tagSymLabels) if (.not. associated(SymClasses)) allocate(SymClasses(nBasis)) call LogMemAlloc('SymClasses', nBasis, 4, this_routine, tagSymClasses) if (.not. allocated(SymConjTab)) allocate(SymConjTab(nSymlabels)) call LogMemAlloc('SymConjTab', nSymlabels, 4, this_routine, tagSymConjTab) if (TwoCycleSymGens .or. tUEG) then DO I = 1, NBASIS, 2 ! place the sym label of each state in SymClasses(ISTATE). For molp sym, this is ! the log_2 of the symmetry bit string IF (G1(I)%Sym%s == 0) THEN ! we don't have symmetry, so fake it. SymClasses((I + 1) / 2) = 1 ELSE SymClasses((I + 1) / 2) = int(G1(I)%Sym%s) + 1 end if end do ! list the symmetry string of each sym label DO I = 1, NSYMLABELS SYMLABELS(I)%s = I - 1 ! Abelian representations are self-inverses if the group is ! real. SymConjTab(I) = I end do #ifdef DEBUG_ write(stdout, *) "Label, Sym, SymConjLabel, SymConj, SymProd" do i = 1, nsymlabels write(stdout, "(5I12)") i, symlabels(i), SymConjTab(i), symlabels(SymConjTab(i)), & SYMPROD(symlabels(i), symlabels(SymConjTab(i))) end do #endif else if (.not. tHUB .or. treal) then ! Hubbard symmetry info set up in GenHubMomIrrepsSymTable. ! except for the real-space lattice! symlabels(:)%s = -1 do i = 1, nbasis, 2 do ilabel = 1, nsymlabels if (symlabels(ilabel)%s == g1(i)%sym%s) then ! Already found this symmetry label. symclasses((i + 1) / 2) = ilabel exit else if (symlabels(ilabel)%s == -1) then ! Have not found this label... symclasses((i + 1) / 2) = ilabel symlabels(ilabel)%s = g1(i)%sym%s exit end if end do end do ! Find inverses. do ilabel = 1, nsymlabels do i = 1, nsymlabels if (SymEq(symlabels(i), SymConj(symlabels(ilabel)))) then SymConjTab(ilabel) = i exit end if end do end do #ifdef DEBUG_ write(stdout, *) "Label, Sym, SymConjLabel, SymConj, SymProd" do i = 1, nsymlabels write(stdout, "(5I12)") i, symlabels(i), SymConjTab(i), symlabels(SymConjTab(i)), & SYMPROD(symlabels(i), symlabels(SymConjTab(i))) end do #endif end if END SUBROUTINE GENMOLPSYMTABLE ! Freeze the SYM LABELS and reps ! NHG is the old number of orbs ! NBASIS is the new number of orbs ! GG(I) is the new index of the (old) orb I SUBROUTINE FREEZESYMLABELS(NHG, NBASIS, GG, FRZ) IMPLICIT NONE INTEGER NHG, NBASIS, GG(NHG) INTEGER I INTEGER NSL(NBASIS) LOGICAL FRZ character(*), parameter :: this_routine = 'FreezeSymLabels' ! SYMLABELS is used to classify all states which transform with the same symmetry ! for the excitation generation routines ! Each state's symmetry falls into a class SymClasses(ISTATE). ! The symmetry bit string, decomposing the sym label into its component irreps is in ! SymLabels(ISYMLABEL) ! The characters of this class are stored in SYMLABELCHARS(1:NROT, SymClasses(ISTATE)) ! The total number of symmetry labels is NSYMLABELS !.. SYMREPS(1,IBASISFN) contains the numnber of the representation !.. of which IBASISFN is a part. IF (.NOT. FRZ) THEN DO I = 1, NHG, 2 IF (GG(I) /= 0) THEN NSL((GG(I) + 1) / 2) = SymClasses((I + 1) / 2) end if end do DO I = 1, NBASIS / 2 SymClasses(I) = NSL(I) ! write(stdout,*) "SL",I,SymClasses(I) end do DO i = 1, nhg IF (GG(i) /= 0) NSL(GG(i)) = Symreps(1, i) end do DO i = 1, nbasis Symreps(1, i) = NSL(i) end do ELSE IF (associated(SYMCLASSES2)) call stop_all(this_routine, 'Problem in freezing') allocate(SymClasses2(nBasis / 2)) call LogMemAlloc('SymClasses2', nBasis / 2, 4, this_routine, tagSymClasses2) DO I = 1, NHG, 2 IF (GG(I) /= 0) THEN NSL((GG(I) + 1) / 2) = SymClasses((I + 1) / 2) end if end do DO I = 1, NBASIS / 2 SymClasses2(I) = NSL(I) ! write(stdout,*) "SL",I,SymClasses(I) end do DO i = 1, nhg IF (GG(i) /= 0) NSL(GG(i)) = Symreps(1, i) end do DO i = 1, nbasis Symreps(1, i) = NSL(i) end do end if ! write(stdout,*) "Sym Reps after Freezing" ! DO i=1,nbasis ! write(stdout,*) i,Symreps(1,i),Symreps(2,i) ! end do END SUBROUTINE FREEZESYMLABELS SUBROUTINE GENMOLPSYMREPS() IMPLICIT NONE INTEGER I, J ! TYPE(BasisFN) G1(NBASIS) ! INTEGER NBASIS,BRR(NBASIS) ! real(dp) ARR(NBASIS) character(*), parameter :: this_routine = 'GENMOLPSYMREPS' LOGICAL tNew if (tKPntSym) then !These symmetry routines only work for cases where all irreps are their !own inverse. In systems with multiple kpoints, this will not be the !case. Setup the symreps for non-abelian symmetries. CALL GENSYMREPS(G1, NBASIS, ARR, 1.e-6_dp) return end if ! now work out which reps are degenerate and label them allocate(SymReps(2, nBasis)) call LogMemAlloc('SymReps', 2 * nBasis, 4, this_routine, tagSymReps) SymReps(:, :) = 0 J = 0 DO I = 1, NBASIS ! write(stdout,*) I,nbasis tNew = .true. IF (tSymIgnoreEnergies .AND. MOD(I, 2) == 0) THEN !Pair even orbs up with the odd ones. SYMREPS(2, J) = SYMREPS(2, J) + 1 tNew = .false. else if (I > 1) THEN IF ((ABS(ARR(I, 1) - ARR(I - 1, 1)) < 1.0e-5_dp) & .AND. (G1(BRR(I))%Sym%s == G1(BRR(I - 1))%Sym%s)) THEN ! we have the same degenerate rep as the previous entry SYMREPS(2, J) = SYMREPS(2, J) + 1 tNew = .false. end if end if IF (tNew) THEN ! we have a new rep J = J + 1 SYMREPS(2, J) = 1 end if SYMREPS(1, BRR(I)) = J end do ! write(stdout,*) "Sym Reps MOLPRO" ! DO i=1,nbasis ! write(stdout,*) i,Symreps(1,i),Symreps(2,i) ! end do END SUBROUTINE GENMOLPSYMREPS ! delete a symmetry table if one existed. SUBROUTINE ENDSYM() IMPLICIT NONE character(*), parameter :: this_routine = 'EndSym' if (allocated(SymTable)) then deallocate(SymTable) call LogMemDealloc(this_routine, tagSymTable) end if if (allocated(SymConjTab)) then deallocate(SymConjTab) call LogMemDealloc(this_routine, tagSymConjTab) end if if (allocated(SymReps)) then deallocate(SymReps) call LogMemDealloc(this_routine, tagSymReps) end if if (associated(SymClasses)) then deallocate(SymClasses) call LogMemDealloc(this_routine, tagSymClasses) end if nullify (SymClasses) if (associated(SymClasses2)) then deallocate(SymClasses2) call LogMemDealloc(this_routine, tagSymClasses2) end if nullify (SymClasses2) if (allocated(SymLabels)) then deallocate(SymLabels) call LogMemDealloc(this_routine, tagSymLabels) end if if (allocated(SymLabelChars)) then deallocate(SymLabelChars) call LogMemDealloc(this_routine, tagSymLabelChars) end if if (allocated(IRREPCHARS)) then deallocate(IRREPCHARS) call LogMemDealloc(this_routine, tagIRREPCHARS) end if if (allocated(SymStatePairs)) then deallocate(SymStatePairs) call LogMemDealloc(this_routine, tagSymStatePairs) end if if (allocated(SymLabelList)) then deallocate(SymLabelList) call LogMemDealloc(this_routine, tagSymLabelList) end if if (allocated(SymLabelCounts)) then deallocate(SymLabelCounts) call LogMemDealloc(this_routine, tagSymLabelCounts) end if if (associated(SymIndex)) then deallocate(SymIndex) call LogMemDealloc(this_routine, tagSymIndex) end if if (associated(SymIndex2)) then deallocate(SymIndex2) call LogMemDealloc(this_routine, tagSymIndex2) end if if (allocated(SymLabelList2)) then deallocate(SymLabelList2) end if if (allocated(SymLabelCounts2)) then deallocate(SymLabelCounts2) end if if (allocated(OrbClassCount)) then deallocate(OrbClassCount) end if if (allocated(SymPairProds)) then deallocate(SymPairProds) call LogMemDealloc(this_routine, tagSymPairProds) end if END SUBROUTINE ENDSYM ! Precompute a list of the symmetry product of all pairs of symmetry labels SUBROUTINE GENSymStatePairs(NSTATES, FRZ) IMPLICIT NONE INTEGER I, TOT, NSTATES INTEGER TEMPLIST(NSTATES), lo, hi LOGICAL FRZ character(*), parameter :: this_routine = 'GenSymStatePairs' if (tAbelianFastExcitGen .AND. .NOT. tAbelian) THEN write(stdout, *) "Fast Abelian excitation generators specified,", & "but abelian symmetry not in use. Using slow generators." tAbelianFastExcitGen = .false. end if if (.not. tAbelianFastExcitGen) THEN !We are not in abelian fast excitgen - we are always storing state pairs, whether specified or not. tStoreStateList = .true. end if !May need to deallocate, since this info is allocated in storage of UMAT before freezing if (allocated(SymPairProds)) then deallocate(SymPairProds) call LogMemDealloc(this_routine, tagSymPairProds) end if if (allocated(SymLabelList)) then deallocate(SymLabelList) call LogMemDealloc(this_routine, tagSymLabelList) end if if (allocated(SymLabelCounts)) then deallocate(SymLabelCounts) call LogMemDealloc(this_routine, tagSymLabelCounts) end if if (allocated(SymStatePairs)) then deallocate(SymStatePairs) call LogMemDealloc(this_routine, tagSymStatePairs) end if allocate(SymLabelList(nStates)) call LogMemAlloc('SymLabelList', nStates, 4, this_routine, tagSymLabelList) allocate(SymLabelCounts(2, nSymLabels)) call LogMemAlloc('SymLabelCounts', 2 * nSymLabels, 4, this_routine, tagSymLabelCounts) ! First deal with listing single states DO I = 1, NSTATES SYMLABELLIST(I) = I IF (FRZ) THEN TEMPLIST(I) = SymClasses2(I) ELSE TEMPLIST(I) = SymClasses(I) end if end do ! order according to sym label, so SYMLABELLIST gets a list of states grouped under SYMLABEL call sort(tempList, symLabelList) ! 1:nStates SYMLABELCOUNTS(:, :) = 0 SYMLABELCOUNTS(1, TEMPLIST(1)) = 1 SYMLABELCOUNTS(2, TEMPLIST(1)) = 1 DO I = 2, NSTATES IF (TEMPLIST(I) /= TEMPLIST(I - 1)) THEN ! add a new sym label SYMLABELCOUNTS(2, TEMPLIST(I)) = 1 SYMLABELCOUNTS(1, TEMPLIST(I)) = I ! sort the symlabellist lo = symLabelcounts(1, tempList(i - 1)) hi = lo + symLabelCounts(2, tempList(i - 1)) - 1 call sort(symLabelList(lo:hi)) ELSE SYMLABELCOUNTS(2, TEMPLIST(I)) = SYMLABELCOUNTS(2, TEMPLIST(I)) + 1 end if end do lo = symLabelcounts(1, tempList(i - 1)) hi = lo + symLabelCounts(2, tempList(i - 1)) - 1 call sort(symLabelList(lo:hi)) ! DO I=1,NSYMLABELS ! write(stdout,*) "NSL",I,SYMLABELCOUNTS(1,I),SYMLABELCOUNTS(2,I) ! & SymLabels(I) ! end do ! Now deal with pairs of states if (.not. tStoreStateList) then !.. We don't bother listing all pairs of orbs, because we can calculate the number !.. and they're easy to generate. !.. Instead of listing all pairs of states, we can list all pairs of sym classes (labels), !.. ordered according to their sym prod. allocate(SymPairProds(nSymLabels**2)) call LogMemAlloc('SymPairProds', nSymLabels**2, & SymPairProdSize * 8, this_routine, tagSymPairProds) SymPairProds(:) = SymPairProd(Symmetry(0), 0, 0, 0, 0) ! Now enumerate all pairs, and classify their product, but don't store them. nSymPairProds = 0 CALL GenSymPairs(nSymLabels, 0) ! Now sort the SymPairProds into order call sort(symPairProds(1:nSymPairProds)) TOT = 0 DO I = 1, nSymPairProds ! write(stdout,"(I4,Z8,4I4)") ! & I,SymPairProds(I)%Sym,SymPairProds(I)%nPairs, ! & SymPairProds(I)%nIndex, SymPairProds(I)%nPairsStateSS, ! & SymPairProds(I)%nPairsStateOS SymPairProds(I)%nIndex = TOT TOT = TOT + SymPairProds(I)%nPairs SymPairProds(I)%nPairs = 0 SymPairProds(I)%nPairsStateSS = 0 SymPairProds(I)%nPairsStateOS = 0 end do write(stdout, *) TOT, " Symmetry PAIRS" write(stdout, *) NSYMPAIRPRODS, " DISTINCT ORBITAL PAIR PRODUCT SYMS" allocate(SymStatePairs(2, 0:TOT - 1)) call LogMemAlloc('SymStatePairs', 2 * TOT, 4, this_routine, tagSymStatePairs) SymStatePairs(:, :) = 0 CALL GenSymPairs(nSymLabels, 1) ! write(stdout,*) "Sym State Pairs" ! DO I=0,TOT-1 ! write(stdout,*) I,SymStatePairs(1:2,I) ! end do ! write(stdout,*) "SymLabelList: ",SymLabelList(:) ! write(stdout,*) "***","SymLabelCounts..." ! write(stdout,*) SymLabelCounts(1,:) ! write(stdout,*) "***" ! write(stdout,*) SymLabelCounts(2,:) else !.. Non-abelian symmetry requires us to go through and work out all the possible pairs of orbs. allocate(SymPairProds(nSymLabels**2)) call LogMemAlloc('SymPairProds', nSymLabels**2, SymPairProdSize * 8, this_routine, tagSymPairProds) SymPairProds = SymPairProd(Symmetry(0), 0, 0, 0, 0) ! Now enumerate all pairs, and classify their product, but don't store them. nSymPairProds = 0 CALL GENALLSymStatePairs(NSTATES, .FALSE., FRZ) ! Now sort the SymPairProds into order call sort(symPairProds(1:nSymPairProds)) TOT = 0 ! write(stdout,*) "SymPairs",nSymPairProds DO I = 1, nSymPairProds SymPairProds(I)%nIndex = TOT TOT = TOT + SymPairProds(I)%nPairs SymPairProds(I)%nPairs = 0 end do write(stdout, *) TOT, " STATE PAIRS" write(stdout, *) NSYMPAIRPRODS, " DISTINCT ORBITAL PAIR PRODUCT SYMS" allocate(SymStatePairs(2, 0:TOT - 1)) call LogMemAlloc('SymStatePairs', 2 * TOT, 4, this_routine, tagSymStatePairs) SymStatePairs(:, :) = 0 CALL GENALLSymStatePairs(NSTATES, .TRUE., FRZ) end if END SUBROUTINE GENSymStatePairs != Generates Symmetry pairs in three passes: != iPass != 0 Count number of pairs of symmetries for each possible symmetry product != 1 Store each pair of symmetries for each symmetry product != 2 Count the number of pairs of STATES for each pair of symmetries. != != NB This differs from GENAllSymStatePairs which goes through every possible != pair of states (and so is O(N^2)), and eventually stores them all. != Here we only store pairs of symmetries (but calculate the number of pairs of states) != This will only work for Abelian symmetries. SUBROUTINE GenSymPairs(nSymLabels, iPass) IMPLICIT NONE INTEGER iPass INTEGER I, J TYPE(Symmetry) PROD INTEGER nSymLabels, iProd INTEGER iSS, iOS DO I = 1, nSymLabels DO J = I, nSymLabels ! write(stdout,*) I,J PROD = SYMPROD(SymLabels(I), SymLabels(J)) CALL FindSymProd(Prod, SymPairProds, nSymPairProds, iProd) IF (iProd == nSymPairProds + 1) THEN nSymPairProds = nSymPairProds + 1 SymPairProds(iProd)%Sym = Prod SymPairProds(iProd)%nIndex = 0 SymPairProds(iProd)%nPairs = 0 SymPairProds(iProd)%nPairsStateSS = 0 SymPairProds(iProd)%nPairsStateOS = 0 end if !.. iOS counts the number of pairs of spin-orbitals with the opposite spin, which give rise to the !.. given symmetry product. iSS is for same spin orbital pairs. iOS = SymLabelCounts(2, I) * SymLabelCounts(2, J) if (i /= j) iOS = iOS * 2 !.. Same spin has n(n-1)/2 if same state if (i /= J) then iSS = SymLabelCounts(2, I) * SymLabelCounts(2, J) else iSS = (SymLabelCounts(2, I) * (SymLabelCounts(2, J) - 1)) / 2 end if if (iOS > 0 .or. iSS > 0) THEN IF (iPass == 1) THEN ! put the pair into the list of pairs. SymStatePairs(1, SymPairProds(iProd)%nIndex + SymPairProds(iProd)%nPairs) = I SymStatePairs(2, SymPairProds(iProd)%nIndex + SymPairProds(iProd)%nPairs) = J ! write(stdout,"(3I5,Z10,3I5)") ! & iProd,I,J,PROD,SymPairProds(iProd)%nIndex ! & +SymPairProds(iProd)%nPairs, ! & SymPairProds(iProd)%nIndex, ! & SymPairProds(iProd)%nPairs end if ! increment the counter in the pairlist SymPairProds(iProd)%nPairs = SymPairProds(iProd)%nPairs + 1 SymPairProds(iProd)%nPairsStateOS = SymPairProds(iProd)%nPairsStateOS + iOS SymPairProds(iProd)%nPairsStateSS = SymPairProds(iProd)%nPairsStateSS + iSS ! write(stdout,*) "NN",SymLabelCounts(2,I),SymLabelCounts(2,J), ! & SymPairProds(iProd)%nPairsStateSS, ! & SymPairProds(iProd)%nPairsStateOS end if end do end do END SUBROUTINE GenSymPairs SUBROUTINE GENALLSymStatePairs(NSTATES, TSTORE, FRZ) IMPLICIT NONE LOGICAL TSTORE, FRZ INTEGER I, J TYPE(Symmetry) PROD INTEGER NSTATES, iProd DO I = 1, NSTATES DO J = I, NSTATES ! write(stdout,*) I,J,SymClasses(I),SymClasses(J) IF (FRZ) THEN PROD = SYMPROD(SymLabels(SymClasses2(I)), SymLabels(SymClasses2(J))) ELSE PROD = SYMPROD(SymLabels(SymClasses(I)), SymLabels(SymClasses(J))) end if CALL FindSymProd(Prod, SymPairProds, nSymPairProds, iProd) IF (TSTORE) THEN ! put the pair into the list of pairs. SymStatePairs(1, SymPairProds(iProd)%nIndex + SymPairProds(iProd)%nPairs) = I SymStatePairs(2, SymPairProds(iProd)%nIndex + SymPairProds(iProd)%nPairs) = J end if IF (iProd == nSymPairProds + 1) THEN nSymPairProds = nSymPairProds + 1 SymPairProds(iProd)%Sym = Prod SymPairProds(iProd)%nIndex = 0 SymPairProds(iProd)%nPairs = 0 end if ! incrememnt the counter in the pairlist SymPairProds(iProd)%nPairs = SymPairProds(iProd)%nPairs + 1 end do end do END SUBROUTINE GENALLSymStatePairs SUBROUTINE FindSymProd(Prod, SymPairProds, nSymPairProds, iProd) implicit none INTEGER nSymPairProds, iProd TYPE(SymPairProd) SymPairProds(nSymPairProds) TYPE(Symmetry) Prod DO iProd = 1, nSymPairProds IF (SYMEQ(SymPairProds(iProd)%Sym, Prod)) EXIT end do END SUBROUTINE FindSymProd !.. SYMREPS is used to group together degenerate sets of orbitals of the same sym !.. (e.g. the six orbitals which might make up a T2g set), and is used for working !.. out the symmetry of a determinant in GETSYM !.. It uses that fact that even for non-abelian groups a completely filled degenerate symmetry set is totally symmetric. !.. Thus each member of a set of states which when completely filled gives a totally !symmetric det should be labelled with the same symrep ! SYMREPS(2,*) has two sets of data: ! SYMREPS(1,IBASISFN) contains the numnber of the representation ! of which IBASISFN is a part. ! SYMPREPS(2,IREP) contains the degeneracy of the rep IREP ! The new method does the following: ! Identify all the completely filled reps. ! use ADDELECSYM to add together the momenta of these. ! These together are totally symmetric ! Identify all part-filled degenerate non-reduced representations ! Use SYMPROD and ADDELECSYM to generate the resultant symmetry of these subroutine getsym_wrapper(det, sym) integer, intent(in) :: det(nel) type(basisfn), intent(out) :: sym call getsym(det, nel, G1, nBasisMax, sym) end subroutine SUBROUTINE GETSYM(NI2, NEL, G1, NBASISMAX, ISYM) IMPLICIT NONE INTEGER NEL, NI(NEL), nBasisMax(5, *) TYPE(BasisFn) G1(*), ISym INTEGER I, J, NI2(NEL) INTEGER NREPS(NEL), NELECS(NEL), SSYM I = 1 NREPS(1:NEL) = 0 CALL SETUPSYM(ISYM) IF (tFixLz) THEN CALL GetLz(NI2, NEL, ISYM%Ml) ELSE ISYM%Ml = 0 end if NI(1:NEL) = NI2(1:NEL) IF (tAbelian) THEN !For Abelian symmetry we don't need the symreps malarky. DO I = 1, NEL ISYM%Sym = SYMPROD(ISYM%Sym, G1(NI(I))%Sym) ! add the momentum CALL ADDELECSYM(NI(I), G1, NBASISMAX, ISYM) end do ELSE DO I = 1, NEL ! Count all electrons in each rep ! NREPS(J) is the rep, and NELECS(J) is the number of electrons in that rep J = 1 DO WHILE (J < NEL) IF (NREPS(J) == 0) exit IF (NREPS(J) == SYMREPS(1, NI(I))) THEN ! We've found the slot for the rep. increment it and leave. NELECS(J) = NELECS(J) + 1 J = NEL end if J = J + 1 end do IF (J <= NEL) THEN ! need to put the new rep in a new space NREPS(J) = SYMREPS(1, NI(I)) NELECS(J) = 1 end if END DO ! now go through and see which are closed and which open DO I = 1, NEL J = 1 DO WHILE (NREPS(J) /= SYMREPS(1, NI(I))) J = J + 1 end do ! electron NI(I) is in rep NREPS(J) IF (NELECS(J) /= SYMREPS(2, NREPS(J))) THEN ! we don't have a closed shell ! add the sym product ! print *, "i, ni(i): ", i, ni(i) ! print *, "s1,s2: ", isym%sym%s, G1(ni(i))%sym%s ISYM%Sym = SYMPROD(ISYM%Sym, G1(NI(I))%Sym) ! print *, "isym after: ", isym%sym%s end if ! add the momentum CALL ADDELECSYM(NI(I), G1, NBASISMAX, ISYM) end do end if ! round the momentum CALL ROUNDSYM(ISYM, NBASISMAX) RETURN END SUBROUTINE GETSYM SUBROUTINE GetLz(nI, NElec, Lz) INTEGER :: NElec INTEGER :: nI(NElec), Lz, i Lz = 0 do i = 1, NElec Lz = Lz + G1(nI(i))%Ml end do END SUBROUTINE GetLz ! Given a set of characters of states, generate all relevant irreps which span the set of characters. SUBROUTINE GENIRREPS(TKP, IMPROPER_OP, NROTOP) IMPLICIT NONE INTEGER I, J, K LOGICAL LDO, LDO2 TYPE(Symmetry) iDecomp INTEGER NEXTSYMLAB complex(dp) REPCHARS(NROT, NSYMLABELS * 10) INTEGER NREPS, NROTOP real(dp) NORM LOGICAL TKP, INV, IMPROPER_OP(NROTOP) character(*), parameter :: this_routine = 'GENIRREPS' NREPS = 0 ! Initialize the table with the totally symmetric rep. INV = .FALSE. DO I = 1, NROT IRREPCHARS(I, 1) = 1 IF (IMPROPER_OP(MOD(I - 1, NROTOP) + 1) .and. .not. TKP) INV = .TRUE. end do NSYM = 1 IF (INV) THEN write(stdout, *) "Inversion centre detected" NSYM = NSYM + 1 ! There's an inversion centre, so we can immediately create an A1u irrep DO I = 1, NROT IF (IMPROPER_OP(MOD(I - 1, NROTOP) + 1)) THEN IRREPCHARS(I, NSYM) = -1 ELSE IRREPCHARS(I, NSYM) = 1 end if end do ! CALL writeirreptab(stdout,IRREPCHARS,NROT,NSYM) end if LDO = .TRUE. NEXTSYMLAB = 1 LDO2 = .TRUE. DO WHILE (LDO .OR. LDO2) ! CALL writeirreptab(stdout,IRREPCHARS,NROT,NSYM) ! write(stdout,*) NREPS," non-reducible" ! CALL writeirreptab(stdout,REPCHARS,NROT,NREPS) ! First see if all the products of chars are decomposable LDO = .FALSE. NREPS = 0 lp1: DO I = 1, NSYM DO J = I, NSYM NREPS = NREPS + 1 IF (NREPS > NSYMLABELS * 10) call stop_all(this_routine, 'TOO MANY REPS') DO K = 1, NROT REPCHARS(K, NREPS) = CONJG(IRREPCHARS(K, I)) * IRREPCHARS(K, J) end do ! write(stdout,*) NREPS,"PROD",I,J ! CALL N_MEMORY_CHECK ! CALL writechars(stdout,REPCHARS(1,NREPS),NROT,"ADDPRD") IF (GETIRREPDECOMP(REPCHARS(1, NREPS), IRREPCHARS, NSYM, NROT, IDECOMP, NORM, TAbelian)) THEN ! CHARWORK now contains the remainder, which will be a new irrep (or combination or irreps), which we need to add IF (ABS(NORM - NROT) <= 1.0e-2_dp) THEN ! if it's an irrep NSYM = NSYM + 1 IF (NSYM > 64) call stop_all(this_routine, "MORE than 64 irreps") DO K = 1, NROT IRREPCHARS(K, NSYM) = REPCHARS(K, NREPS) end do ! CALL writeirreptab(stdout,IRREPCHARS,NROT,NSYM) NREPS = NREPS - 1 LDO = .TRUE. EXIT lp1 ELSE ! write(stdout,*) "IDECOMP:", IDECOMP,NORM,"SYMS:",NSYM ! CALL writechars(stdout,REPCHARS(1,NREPS),NROT,"REMAIN") ! It's not an irrep, but we cannot reduce it. Store only if we think we've got all the irreps. ! write(stdout,*) "NR",NREPS,LDO2 IF (LDO2) NREPS = NREPS - 1 ! NREPS=NREPS-1 end if ELSE ! write(stdout,*) "IDECOMP:", IDECOMP NREPS = NREPS - 1 end if END DO END DO lp1 ! write(stdout,*) LDO,NEXTSYMLAB,NSYMLABELS IF (LDO) CYCLE ! Check to see if the next symlabel's char is decomposable lp2: DO WHILE (NEXTSYMLAB <= NSYMLABELS) NREPS = NREPS + 1 IF (NREPS > NSYMLABELS * 10) call stop_all(this_routine, 'TOO MANY REPS') DO I = 1, NROT REPCHARS(I, NREPS) = SYMLABELCHARS(I, NEXTSYMLAB) end do ! CALL writechars(stdout,REPCHARS(1,NREPS),NROT,"ADDST ") IF (GETIRREPDECOMP(REPCHARS(1, NREPS), IRREPCHARS, NSYM, NROT, IDECOMP, NORM, TAbelian)) THEN ! CHARWORK now contains the remainder, which will be a new irrep (or combination or irreps), which we need to add IF (ABS(NORM - NROT) <= 1.0e-2_dp) THEN ! if it's an irrep NSYM = NSYM + 1 IF (NSYM > 64) call stop_all(this_routine, "MORE than 64 irreps") DO I = 1, NROT IRREPCHARS(I, NSYM) = REPCHARS(I, NREPS) end do ! CALL writeirreptab(stdout,IRREPCHARS,NROT,NSYM) NREPS = NREPS - 1 LDO = .TRUE. EXIT lp2 ELSE ! write(stdout,*) "IDECOMP:", IDECOMP,NORM,"SYMS:",NSYM ! CALL writechars(stdout,REPCHARS(1,NREPS),NROT,"REMAIN") ! It's not an irrep, but we cannot reduce it. Store only if we think we've got all the irreps. IF (LDO2) NREPS = NREPS - 1 end if ELSE ! write(stdout,*) "IDECOMP:", IDECOMP NREPS = NREPS - 1 end if NEXTSYMLAB = NEXTSYMLAB + 1 IF (.NOT. LDO) THEN ! We've not manage to add any more irreps, so we have achieved self-consistency. !Do one more pass to check, saving all C.. non-reducible reps LDO = .TRUE. LDO2 = .FALSE. NREPS = 0 end if END DO lp2 end do ! write(stdout, *) "IRREP TABLE" CALL writeirreptab(stdout, IRREPCHARS, NROT, NSYM) IF (NREPS > 0) THEN write(stdout, *) NREPS, " non-reducible" CALL writeirreptab(stdout, REPCHARS, NROT, NREPS) ! IF(NREPS.GT.1) THEN call stop_all(this_routine, "More than 1 non-reducible reps found.") ! end if ! we can cope with a single reducible rep. ! NSYM=NSYM+1 ! DO I=1,NROT ! IRREPCHARS(I,NSYM)=REPCHARS(I,NREPS) ! end do end if ! Classify each of the symlabels with its decomposition into irreps DO I = 1, NSYMLABELS CALL DECOMPOSEREP(SYMLABELCHARS(1, I), IDECOMP) SymLabels(I) = IDECOMP end do END SUBROUTINE GENIRREPS ! Display irrep table SUBROUTINE WRITEIRREPTAB(IUNIT, CHARS, NROT, NSYM) IMPLICIT NONE INTEGER IUNIT, NROT, NSYM complex(dp) CHARS(NROT, NSYM) CHARACTER(6) STR INTEGER I, J LOGICAL LCOMP, LREAL LCOMP = .FALSE. LREAL = .FALSE. DO I = 1, NSYM DO J = 1, NROT IF (ABS(REAL(CHARS(J, I))) > 1.0e-2_dp .AND. ABS(AIMAG(CHARS(J, I))) > 1.0e-2_dp) LCOMP = .TRUE. IF (ABS(REAL(CHARS(J, I)) - NINT(REAL(CHARS(J, I)))) > 1.0e-2_dp) LREAL = .TRUE. IF (ABS(AIMAG(CHARS(J, I)) - NINT(AIMAG(CHARS(J, I)))) > 1.0e-2_dp) LREAL = .TRUE. end do end do DO I = 1, NSYM write(STR, "(A,I3)") "SYM", I CALL WRITECHARSF(IUNIT, CHARS(1, I), NROT, STR, LCOMP, LREAL) end do write(IUNIT, *) END SUBROUTINE WRITEIRREPTAB ! Display a line of characters SUBROUTINE WRITECHARS(IUNIT, CHARS, NROT, STR) IMPLICIT NONE INTEGER IUNIT, NROT complex(dp) CHARS(NROT) INTEGER J CHARACTER(6) STR LOGICAL LCOMP, LREAL ! First do a check for the format LCOMP = .FALSE. LREAL = .FALSE. DO J = 1, NROT IF (ABS(REAL(CHARS(J))) > 1.0e-2_dp .AND. ABS(AIMAG(CHARS(J))) > 1.0e-2_dp) LCOMP = .TRUE. IF (ABS(REAL(CHARS(J)) - NINT(REAL(CHARS(J)))) > 1.0e-2_dp) LREAL = .TRUE. IF (ABS(AIMAG(CHARS(J)) - NINT(AIMAG(CHARS(J)))) > 1.0e-2_dp) LREAL = .TRUE. end do CALL WRITECHARSF(IUNIT, CHARS, NROT, STR, LCOMP, LREAL) END SUBROUTINE WRITECHARS SUBROUTINE WRITECHARSF(IUNIT, CHARS, NROT, STR, LCOMP, LREAL) IMPLICIT NONE INTEGER IUNIT, NROT complex(dp) CHARS(NROT) INTEGER J CHARACTER(6) STR LOGICAL LCOMP, LREAL write(IUNIT, "(A6,A)", advance='no') STR, ": " DO J = 1, NROT IF (LCOMP) THEN IF (LREAL) THEN write(IUNIT, "(A,2G16.9,A)", advance='no') "(", & NINT(REAL(CHARS(J)) * 1000) / 1000.0_dp, & NINT(AIMAG(CHARS(J)) * 1000) / 1000.0_dp & , ")" ELSE write(IUNIT, "(A,2F6.3,A)", advance='no') "(", CHARS(J), ")" end if ELSE IF (ABS(AIMAG(CHARS(J))) > 1.0e-2_dp) THEN ! write in terms of I. IF (LREAL) THEN write(IUNIT, "(G14.9,A)", advance='no') CHARS(J), " " ELSE IF (ABS(AIMAG(CHARS(J)) + 1.0_dp) < 1.0e-2_dp) THEN write(IUNIT, "(A)", advance='no') " -I " else if (ABS(AIMAG(CHARS(J)) - 1.0_dp) < 1.0e-2_dp) THEN write(IUNIT, "(A)", advance='no') " I " ELSE write(IUNIT, "(I2,A)", advance='no') NINT(AIMAG(CHARS(J))), "I " end if end if ELSE IF (LREAL) THEN write(IUNIT, "(G21.9)", advance='no') REAL(CHARS(J)), " " ELSE write(IUNIT, "(I3)", advance='no') NINT(REAL(CHARS(J))) end if end if end if end do write(IUNIT, *) END SUBROUTINE WRITECHARSF ! Decompose rep CHARS into irreps in IRREPCHARS. Bit 0 in IDECOMP corresponds to the first irrep etc. ! CHARS at the end contains the remainder after the decomposition. SUBROUTINE DECOMPOSEREP(CHARSIN, IDECOMP) IMPLICIT NONE TYPE(Symmetry) IDECOMP complex(dp) CHARS(NROT), CHARSIN(NROT), TOT real(dp) CNORM INTEGER I, J real(dp) NORM, DIFF character(*), parameter :: this_routine = 'DECOMPOSEREP' if (TAbelian) then ! We shouldn't be here! Using symmetry "quantum" numbers ! rather than irreps. call stop_all(this_routine, "Should not be decomposing irreps with Abelian sym") end if IDECOMP%s = 0 CALL DCOPY(NROT * 2, CHARSIN, 1, CHARS, 1) ! write(stdout,*) "Decompose Rep" ! CALL writechars(stdout,CHARS,NROT,"REP ") !,. First check norm of this state CNORM = 0 DO J = 1, NROT CNORM = CNORM + real(CHARS(J) * CHARS(J), dp) end do DO I = 1, NSYM TOT = 0 ! CALL writechars(stdout,IRREPCHARS(1,I),NROT,"IR") ! CALL writechars(stdout,CHARS(1),NROT,"CH") DO J = 1, NROT TOT = TOT + CONJG(IRREPCHARS(J, I)) * CHARS(J) end do ! write(stdout,*) I,TOT IF (abs(TOT) > 1.0e-12_dp) THEN ! Calculate the normalization of the state I which matches (if it's an irrep, this will be 1) NORM = 0 DO J = 1, NROT NORM = NORM + real(CONJG(IRREPCHARS(J, I)) * IRREPCHARS(J, I), dp) end do ! write(stdout,*) "IRREP ",I,(TOT+0.0_dp)/NORM DIFF = ABS(TOT - NINT(ABS(TOT / NORM)) * NORM) IF (DIFF > 1.0e-2_dp) THEN write(stdout, *) 'Symmetry decomposition not complete' CALL writechars(stdout, IRREPCHARS(1, I), NROT, "IRREP ") CALL writechars(stdout, CHARS, NROT, "CHARS ") write(stdout, *) "Dot product: ", (TOT + 0.0_dp) / NORM, TOT, NORM call stop_all(this_routine, 'Incomplete symmetry decomposition') ! The given representation CHARS has fewer irreps in it than the one in IRREPCHARS, and is an irrep ! Hurrah! Remove it from the one in IRREPCHARS, and keep on going) else if (ABS(TOT) > 1.0e-2_dp) THEN ! We've found an (ir)rep which is wholly in CHARS IDECOMP%s = IBSET(IDECOMP%s, I - 1) CNORM = 0 ! write(stdout,*) I,DIFF,TOT,TOT/NORM DO J = 1, NROT CHARS(J) = CHARS(J) - (IRREPCHARS(J, I) * TOT) / NORM CNORM = CNORM + real(CONJG(CHARS(J)) * CHARS(J), dp) end do ! CALL writechars(stdout,IRREPCHARS(1,I),NROT,"DIRREP") ! CALL writechars(stdout,CHARS,NROT,"DCHARS") end if end if end do END SUBROUTINE DECOMPOSEREP ! Decompose rep CHARS into irreps in IRREPCHARS. Bit 0 in IDECOMP corresponds to the first irrep etc. ! CHARS at the end contains the remainder after the decomposition. ! Return .FALSE. if the decomposition is complete and CHARS contains only 0. ! This is used internally in the symmetry routine and destroys CHARS. For general decomposition, !,, use DECOMPOSEREP LOGICAL FUNCTION GETIRREPDECOMP(CHARS, IRREPCHARS, NIRREPS, NROT, IDECOMP, CNORM, TAbelian) IMPLICIT NONE INTEGER NIRREPS, NROT TYPE(Symmetry) IDECOMP complex(dp) IRREPCHARS(NROT, NIRREPS), CHARS(NROT) real(dp) CNORM, NORM, DIFF complex(dp) TOT INTEGER I, J logical TAbelian character(*), parameter :: this_routine = 'GETIRREPDECOMP' if (TAbelian) then ! We shouldn't be here! Using symmetry "quantum" numbers ! rather than irreps. call stop_all(this_routine, "Should not be decomposing irreps with Abelian sym") end if IDECOMP%s = 0 !,. First check norm of this state CNORM = 0 DO J = 1, NROT CNORM = CNORM + real(CONJG(CHARS(J)) * CHARS(J), dp) end do DO I = 1, NIRREPS TOT = 0 DO J = 1, NROT TOT = TOT + CONJG(IRREPCHARS(J, I)) * CHARS(J) end do IF (ABS(TOT) >= 1.0e-2_dp) THEN ! Calculate the normalization of the state I which matches (if it's an irrep, this will be 1) NORM = 0 DO J = 1, NROT NORM = NORM + real(CONJG(IRREPCHARS(J, I)) * IRREPCHARS(J, I), dp) end do ! write(stdout,*) "IRREP ",I,(TOT+0.0_dp)/NORM ! CALL writechars(stdout,CHARS,NROT,"REP ") ! CALL writechars(stdout,IRREPCHARS(1,I),NROT,"IRREP ") DIFF = ABS(TOT - NINT(ABS(TOT / NORM)) * NORM) IF (DIFF >= 1.0e-2_dp .AND. abs(CNORM - NROT) < 1.0e-12_dp) THEN ! The given representation CHARS has fewer irreps in it than the one in IRREPCHARS, and is an irrep ! Hurrah! Remove it from the one in IRREPCHARS, and keep on going) ! DO J=1,NROT ! IRREPCHARS(J,I)=IRREPCHARS(J,I)-CHARS(J)*TOT/CNORM ! end do ! CALL writechars(stdout,IRREPCHARS(1,I),NROT,"NOW ") else if (DIFF < 1.0e-2_dp) THEN ! We've found an (ir)rep which is wholly in CHARS IDECOMP%s = IBSET(IDECOMP%s, I - 1) CNORM = 0 DO J = 1, NROT CHARS(J) = CHARS(J) - (IRREPCHARS(J, I) * TOT) / NORM CNORM = CNORM + real(CONJG(CHARS(J)) * CHARS(J), dp) end do end if end if end do GETIRREPDECOMP = .FALSE. DO J = 1, NROT IF (ABS(CHARS(J)) > 1.0e-2_dp) GETIRREPDECOMP = .TRUE. end do END FUNCTION GETIRREPDECOMP SUBROUTINE GENSYMTABLE IMPLICIT NONE INTEGER I, J, K complex(dp) CHARS(NROT) TYPE(Symmetry) IDECOMP real(dp) CNORM character(*), parameter :: this_routine = 'GENSYMTABLE' allocate(SymTable(nSym, nSym)) call LogMemAlloc('SymTable', nSym**2, SymmetrySize, this_routine, tagSymTable) allocate(SymConjTab(nSym)) call LogMemAlloc('SymConjTab', nSym, 4, this_routine, tagSymConjTab) DO I = 1, NSYM DO K = 1, NROT CHARS(K) = CONJG(IRREPCHARS(K, I)) end do IF (GETIRREPDECOMP(CHARS, IRREPCHARS, NSYM, NROT, IDECOMP, CNORM, TAbelian)) THEN write(stdout, *) "Conjugate of SYM ", I, " not reducible," CALL writechars(stdout, CHARS, NROT, "REMAIN") call stop_all(this_routine, "Symmetry table element not conjugable") end if K = 0 DO WHILE (.NOT. BTEST(IDECOMP%s, 0)) K = K + 1 !RSHIFT(,1) IDECOMP%s = ISHFT(IDECOMP%s, -1) end do IF (IDECOMP%s /= 1) THEN write(stdout, *) "Conjugate of SYM ", I, " not a single SYM," call stop_all(this_routine, 'Incorrect sym conjugate') end if SymConjTab(I) = K + 1 DO J = I, NSYM DO K = 1, NROT CHARS(K) = IRREPCHARS(K, I) * IRREPCHARS(K, J) end do IF (GETIRREPDECOMP(CHARS, IRREPCHARS, NSYM, NROT, IDECOMP, CNORM, TAbelian)) THEN write(stdout, *) "Multiplication of SYMS ", I, J, " not reducible," CALL writechars(stdout, CHARS, NROT, "REMAIN") call stop_all(this_routine, "Symmetry table element not reducible") end if SYMTABLE(I, J) = IDECOMP SYMTABLE(J, I) = IDECOMP ! write(stdout,"(2I3,B12)") I,J,IDECOMP end do end do write(stdout, *) "Symmetry, Symmetry Conjugate" DO I = 1, NSYM write(stdout, *) I, SymConjTab(I) end do END SUBROUTINE GENSYMTABLE SUBROUTINE GENSYMREPS(G1, NBASIS, ARR, DEGENTOL) IMPLICIT NONE INTEGER I, J INTEGER NBASIS TYPE(BasisFN) G1(nBasis) real(dp) ARR(NBASIS, 2) real(dp) DEGENTOL logical lTmp character(*), parameter :: this_routine = 'GenSymReps' ! now work out which reps are degenerate and label them allocate(SymReps(2, nBasis)) call LogMemAlloc('SymReps', 2 * nBasis, 4, this_routine, tagSymReps) symreps(2,1) = 1 symreps(1,1) = 1 J = 1 DO I = 2, NBASIS ltmp = .false. if (abs(arr(i, 2) - arr(i - 1, 2)) < degentol .and. & (tAbelian .or. symeq(G1(i)%sym, G1(i - 1)%sym))) then ! We have the same degenerate rep as the previous entry symreps(2, J) = symreps(2, J) + 1 lTmp = .true. end if if (.not. lTmp) then ! We have a new rep J = J + 1 symreps(2, J) = 1 end if SYMREPS(1, I) = J end do END SUBROUTINE GENSYMREPS !. Irrep symmetries are specified in SYM(5). ! if SYM(5)=0, we assume it's totally symmetric ! Other irreps contributing to the symmetry have bits set in ! SYM. ! e.g. if irreps are a1,a2,b1,b2 LOGICAL elemental FUNCTION LCHKSYM(ISYM, JSYM) type(basisfn), intent(in) :: isym, jsym INTEGER I LCHKSYM = .TRUE. DO I = 1, 3 IF (ISYM%K(I) /= JSYM%K(I)) LCHKSYM = .FALSE. end do IF (ISYM%Ms /= JSYM%Ms) LCHKSYM = .FALSE. IF (ISYM%Ml /= JSYM%Ml) LCHKSYM = .FALSE. ! if the symmetry product of I and J doesn't contain the totally ! symmetric irrep, we set sym to .FALSE. ! [W.D]: does this still work with my new hubbard implementation? LCHKSYM = LCHKSYM .AND. LSYMSYM(SYMPROD(SymConj(ISYM%SYM), JSYM%SYM)) RETURN END FUNCTION LCHKSYM LOGICAL FUNCTION LCHKSYMD(NI, NJ, NEL, G1, NBASISMAX) IMPLICIT NONE TYPE(BASISFN) ISYM, JSYM, G1(*) INTEGER NEL, NI(NEL), NJ(NEL), nBasisMax(5, *) CALL GETSYM(NI, NEL, G1, NBASISMAX, ISYM) CALL GETSYM(NJ, NEL, G1, NBASISMAX, JSYM) LCHKSYMD = LCHKSYM(ISYM, JSYM) RETURN END FUNCTION LCHKSYMD ! NBASISMAX descriptor (1,3) ! ! HUBBARD: ! BITS !. 0 Tilted !. 1 non-pbc !. 2 real-space ! which effects to values ! MOM SPACE ! 0 Non-Tilted Lattice - pbc ! 1 Tilted Lattice - pbc ! 2 Non-Tilted lattice - no pbc ! 3 Tilted Lattice - no pbc ! four following are REAL ! 4 Non-Tilted Lattice - pbc ! 5 Tilted Lattice - pbc ! 6 Non-Tilted lattice - no pbc ! 7 Tilted Lattice - no pbc ! ! (3,3) ! -2 Particle in a box ! -1 UEG ! 0 Hubbard ! 1 Generic spatial ! This only works for momentum variables - 1-4 pure SUBROUTINE ADDELECSYM(IEL, G1, NBASISMAX, ISYM) integer, intent(in) :: iEl, nBasisMax(5, *) TYPE(BASISFN), intent(in) :: G1(*) TYPE(BASISFN), intent(inout) :: ISYM integer :: ielec, i, ssym ielec = iel ssym = 0 IF (NBASISMAX(1, 3) < 4) THEN ! Momentum space if (t_k_space_hubbard) then isym%k = lat%add_k_vec(isym%k, G1(ielec)%k) else ISYM%k = ISYM%k + G1(IELEC)%k end if ! Symmetry space else if (NBASISMAX(3, 3) == 0 .AND. NBASISMAX(1, 3) >= 4) THEN ! We have no symmetries, so do nothing. (we're in real space) ! except Ms else if (NBASISMAX(3, 3) == 1) THEN ! deal with momentum if (t_k_space_hubbard) then isym%k = lat%add_k_vec(isym%k, G1(ielec)%k) else ISYM%k = ISYM%k + G1(IELEC)%k end if end if ISYM%MS = ISYM%MS + G1(IELEC)%MS ISYM%Ml = ISYM%Ml + G1(IELEC)%Ml ! SSYM keeps track of the total S change on adding this electron ! (it is +/-CSF_NSBASIS) I = ISYM%MS + 0 ISYM%Ms = I + SSYM END SUBROUTINE ADDELECSYM pure subroutine roundsym(isym, nbasismax) type(basisfn), intent(inout) :: isym integer, intent(in) :: nbasismax(5, *) integer :: i if (t_new_hubbard .and. t_k_space_hubbard) then ! deal differently with the new k-space hubbard ! use the lattice intrinsic function ! also do something in the real-space case!! ! maybe there i have to set k to 0.. isym = lat%round_sym(isym) return end if IF (NBASISMAX(3, 3) == -2) THEN ! particle in a box ! parity symmetries DO I = 1, 3 ISYM%k(I) = MOD(ISYM%k(I), 2) end do else if (NBASISMAX(3, 3) == -1) THEN ! UEG (can't remember the symmetries of that ! probably momentum conservation) if (tperiodicinmom) CALL MOMPBCSYM(ISYM%k, NBASISMAX) else if (NBASISMAX(3, 3) == 0) THEN ! Hubbard model IF (NBASISMAX(1, 3) < 2) THEN ! momentum conservation - various PBC issues !ALEX PLEASE CHECK. CALL MOMPBCSYM(ISYM%k, NBASISMAX) else if (NBASISMAX(1, 3) >= 2) THEN ! we're in real space so no sym DO I = 1, 3 ISYM%k(I) = 0 end do end if else if (NBASISMAX(3, 3) == 1) THEN ! Generic spatial symmetries ! We need do nothing. ! However, there is still momentum conservation - various PBC issues !ALEX PLEASE CHECK. CALL MOMPBCSYM(ISYM%k, NBASISMAX) end if RETURN END SUBROUTINE ROUNDSYM ! NBASISMAX descriptor (1,3) ! ! HUBBARD: ! BITS !. 0 Tilted !. 1 non-pbc !. 2 real-space ! which effects to values ! MOM SPACE ! 0 Non-Tilted Lattice - pbc ! 1 Tilted Lattice - pbc ! 2 Non-Tilted lattice - no pbc ! 3 Tilted Lattice - no pbc ! four following are REAL ! 4 Non-Tilted Lattice - pbc ! 5 Tilted Lattice - pbc ! 6 Non-Tilted lattice - no pbc ! 7 Tilted Lattice - no pbc ! ! (3,3) ! -2 Particle in a box ! -1 UEG ! 0 Hubbard ! 1 Generic spatial ! This function imposes periodic boundary conditions for the Hubbard Model ! K1(3) is a K vector which is then mapped back into the Brillouin Zone of the Full cell. ! The PBC are such that, a point, (KX,KY) (in terms of the cell lattice vectors), has bounds ! -LENX/2 < KX <= LENX/2 and -LENY/2 <= KY < LENY/2 ! My ASCII art isn't up to drawing a picture unfortunately. pure SUBROUTINE MOMPBCSYM(K1, NBASISMAX) ! NB the third column of NBASISMAX tells us whether it is tilted INTEGER, intent(inout) :: K1(3) INTEGER, intent(in) :: nBasisMax(5, *) INTEGER J, LDIM, AX, AY, LENX, LENY, KK2, T1, T2 real(dp) R1, R2, NORM ! [W.D] ! in case of the new k-space hubbard implementation use the ! built-in k-vec mapping if (t_k_space_hubbard) then ! this should actually never be called anymore! k1 = lat%map_k_vec(k1) return end if ! (AX,AY) is the tilt of the lattice, and corresponds to the lattice vector of the cell. The other lattice vector is (-AY,AX). !These are expressed in terms of the primitive Hubbard lattice vectors AX = NBASISMAX(1, 4) AY = NBASISMAX(2, 4) ! LENX is the length (i.e. number of lattice vectors) in direction (AX,AY). LENY is in the other lattice vector direction LENX = NBASISMAX(1, 5) LENY = NBASISMAX(2, 5) IF (NBASISMAX(1, 3) == 0 .OR. NBASISMAX(1, 3) == 0) THEN ! A non-tilted lattice with PBC DO J = 1, 3 ! non-tilted KK2 = K1(J) LDIM = NBASISMAX(J, 2) - NBASISMAX(J, 1) + 1 KK2 = MOD(KK2, LDIM) IF (KK2 < NBASISMAX(J, 1)) KK2 = KK2 + LDIM IF (KK2 > NBASISMAX(J, 2)) KK2 = KK2 - LDIM K1(J) = KK2 end do else if (NBASISMAX(1, 3) == 1) THEN ! we have a tilted lattice with PBC ! we want the a1,a2 components of k NORM = AX * AX + AY * AY R1 = (AX * K1(1) + AY * K1(2)) / NORM R2 = (AX * K1(2) - AY * K1(1)) / NORM R1 = R1 / LENX + 0.5_dp R2 = R2 / LENY + 0.5_dp ! R1 is now in terms of the shifted extended unit cell. We want 0<R1<=1 ! R2 is now in terms of the shifted extended unit cell. We want 0<=R2<1 ! T1 and T2 will be the vectors to translate the point K1 back into our cell T1 = INT(ABS(R1)) IF (R1 < 0.0_dp) THEN T1 = -T1 IF (abs(t1 - r1) > 1.0e-12_dp) T1 = T1 - 1 end if ! Now T1= highest integer less than R1 (i.e. FLOOR) ! We want to include R1=1, so we have one fewer translations in that case. IF (abs(r1 - t1) < 1.0e-12_dp) T1 = T1 - 1 T2 = INT(ABS(R2)) IF (R2 < 0.0_dp) THEN T2 = -T2 IF (abs(t2 - r2) > 1.0e-12_dp) T2 = T2 - 1 end if ! The following line conserves the top-right rule for edge effects IF (abs(r1 - t1) < 1.0e-12_dp) T1 = T1 - 1 ! Now T2= highest integer less than R2 (i.e. FLOOR) ! Do the translation IF (R1 > 1.0_dp .OR. R1 <= 0.0_dp) R1 = R1 - T1 IF (R2 >= 1.0_dp .OR. R2 < 0.0_dp) R2 = R2 - T2 ! Now convert back to the the Coords in the extended brillouin zone (ish) R1 = (R1 - 0.5_dp) * LENX R2 = (R2 - 0.5_dp) * LENY ! Now K1 is defined as the re-scaled (R1,R2) vector K1(1) = NINT(R1 * AX - R2 * AY) K1(2) = NINT(R1 * AY + R2 * AX) else if (NBASISMAX(1, 3) == -1) THEN if ((abs(k1(1)) > NBASISMAX(1, 2)) .and. tperiodicinmom) then if (k1(1) > 0) then k1(1) = k1(1) - (2 * NBASISMAX(1, 2) + 1) else k1(1) = k1(1) + (2 * NBASISMAX(1, 2) + 1) end if end if end if END SUBROUTINE MOMPBCSYM LOGICAL FUNCTION SYMLT(A, B) IMPLICIT NONE TYPE(Symmetry) A, B IF (A%s >= 0) THEN IF (B%s >= 0) THEN SYMLT = A%s < B%s ELSE SYMLT = .TRUE. end if ELSE IF (B%s >= 0) THEN SYMLT = .FALSE. ELSE SYMLT = A%s < B%s end if end if RETURN END FUNCTION SYMLT LOGICAL FUNCTION SYMNE(A, B) IMPLICIT NONE TYPE(Symmetry) A, B SYMNE = A%s /= B%s RETURN END FUNCTION SYMNE LOGICAL FUNCTION SYMEQ(A, B) use SystemData, only: Symmetry IMPLICIT NONE TYPE(Symmetry) A, B SYMEQ = A%s == B%s RETURN !Need to cope with 'unsigned integers' END FUNCTION SYMEQ LOGICAL FUNCTION SYMGT(A, B) IMPLICIT NONE TYPE(Symmetry) A, B IF (A%s >= 0) THEN IF (B%s >= 0) THEN SYMGT = A%s > B%s ELSE SYMGT = .FALSE. end if ELSE IF (B%s >= 0) THEN SYMGT = .TRUE. ELSE SYMGT = A%s > B%s end if end if RETURN END FUNCTION SYMGT integer Function FindSymLabel(s) IMPLICIT NONE Type(Symmetry) s integer i do i = 1, nSymLabels if (symeq(SymLabels(i), s)) exit end do if (i > nSymLabels) i = 0 FindSymLabel = i return end Function FindSymLabel ! A binary search to find VAL in TAB. TAB is sorted SUBROUTINE BINARYSEARCHSYM(VAL, TAB, LEN, LOC) IMPLICIT NONE TYPE(Symmetry) VAL INTEGER LOC, LEN type(Symmetry) TAB(LEN) INTEGER I, J, IFIRST, N, ILAST I = 1 J = LEN IFIRST = I ILAST = J DO WHILE (J - I >= 1) N = (I + J) / 2 ! write(stdout,"(3I4)",advance='no') I,J,N ! CALL writesym(stdout,TAB(1,I),.FALSE.) ! CALL writesym(stdout,TAB(1,J),.FALSE.) ! CALL writesym(stdout,TAB(1,N),.TRUE.) IF (SYMLT(TAB(N), VAL) .AND. I /= N) THEN IF (SYMNE(TAB(N), TAB(IFIRST))) IFIRST = N ! reset the lower limit I = N else if (SYMGT(TAB(N), VAL)) THEN IF (SYMNE(TAB(N), TAB(ILAST))) ILAST = N ! reset the upper limit J = N else if (SYMEQ(TAB(N), VAL)) THEN ! bingo, we've got it! LOC = N RETURN ELSE ! we've reached a situation where I and J's entries have the same value, and it's ! not the one we want. Leave the loop. I = J end if end do IF (SYMEQ(TAB(I), VAL)) THEN LOC = I else if (SYMEQ(TAB(J), VAL)) THEN LOC = J ELSE ! Failure LOC = 0 end if END SUBROUTINE BINARYSEARCHSYM !This function when called repeatedly, generates all allowed symmmetries. !It appears to not have been modified for Lz symmetry, so will only generate Lz=0 ! AJWT 20110121 SUBROUTINE GENNEXTSYM(NEL, NBASISMAX, TSPN, LMS, TPARITY, IPARITY, TSETUP, TDONE, IMAX, ISYM) IMPLICIT NONE INTEGER NEL, nBasisMax(5, *) INTEGER LMS TYPE(BasisFN) IPARITY, ISYM, IMax(2) LOGICAL TSPN, TPARITY, TSETUP, TMORE, TDONE, TMORE2 INTEGER ILEV IF (TSETUP) THEN IMAX(1) = NullBasisFn IMAX(2) = NullBasisFn DO ILEV = 1, 3 IF (TPARITY) THEN IMAX(1)%k(iLev) = IPARITY%k(ILEV) IMAX(2)%k(ILEV) = IPARITY%k(ILEV) ELSE IMAX(1)%k(iLev) = NBASISMAX(ILEV, 1) IMAX(2)%k(iLev) = NBASISMAX(ILEV, 2) ! IF(NBASISMAX(1,3).EQ.2) THEN ! hubbard non-pbc mom space ! IMAX(ILEV,1)=IMAX(ILEV,1)*NEL ! IMAX(ILEV,2)=IMAX(ILEV,2)*NEL ! end if end if end do IF (TSPN) THEN IMAX(1)%Ms = LMS IMAX(2)%Ms = LMS ELSE IMAX(1)%Ms = NBASISMAX(4, 1) * NEL IMAX(2)%Ms = NBASISMAX(4, 2) * NEL end if ! If we're specifying a sym (TPARITY) in IPARITY(5), and ! we have a system with all 1D reducible orbs, then we put ! that into IMAX IF (NBASISMAX(5, 2) /= 0 .OR. TAbelian) THEN IF (TPARITY) THEN IMAX(1)%Sym%s = IPARITY%Sym%s IMAX(2)%Sym%s = IMAX(1)%Sym%s ELSE IMAX(1)%Sym%s = MinSymRep() IMAX(2)%Sym%s = MaxSymRep() end if ELSE ! we've got a sym system with polydimensional irreps, which leads to ! dets with combinations of irreps, so we cannot put sym into blocks ! JSS: if only 1D symmetries, then a determinant can only interact with ! other determinants of the same symmetry. This applies to Abelian ! groups. If there are multi-dimensional irreps, then this is no ! longer the case, so we set the symmetries to be 0 (i.e. ignore ! symmetry when generating determinants which interact). This is not ! equivalent to setting %s=0 if the Abelian case (which corresponds to ! the totally symmetric irrep). IMAX(1)%Sym%s = 0 IMAX(2)%Sym%s = 0 end if TDONE = .FALSE. CALL DOSYMLIMDEGEN(IMAX, NBASISMAX) ISym = IMax(1) end if IF (TSETUP .AND. KALLOWED(ISYM, NBASISMAX)) RETURN ! Go to the next sym. TMORE2 = .TRUE. TMORE = .TRUE. ILEV = 5 DO WHILE (TMORE2) DO WHILE (ILEV > 0) IF (ILEV == 5) THEN IF (IMAX(1)%Sym%s /= 0) THEN ! symmetry specifiers are incremented by multiplying*2 (unless there are no syms counted) ISYM%Sym%s = ISYM%Sym%s * 2 ELSE Call IncrSym(ISym%Sym) end if IF (ISYM%Sym%s == IMAX(1)%Sym%s) THEN ILEV = ILEV - 1 IF (ILEV == 0) THEN TMORE2 = .FALSE. ! If we've run out of syms, we give up TMORE = .FALSE. end if else if (KALLOWED(ISYM, NBASISMAX)) THEN TMORE2 = .FALSE. ILEV = 0 end if ELSE ISYM%k(ILEV) = ISYM%k(ILEV) + 1 IF (ISYM%k(ILEV) > IMAX(2)%k(ILEV)) THEN ISYM%k(ILEV) = IMAX(1)%k(ILEV) ILEV = ILEV - 1 IF (ILEV == 0) THEN TMORE2 = .FALSE. ! If we've run out of syms, we give up TMORE = .FALSE. end if else if (ILEV < 4) THEN ! We've just incremented one of the higher columns, now go down to the ! lower ones. ILEV = ILEV + 1 ISYM%k(ILEV) = IMAX(1)%k(ILEV) - 1 else if (KALLOWED(ISYM, NBASISMAX)) THEN TMORE2 = .FALSE. ILEV = 0 end if end if end do end do TDONE = .NOT. TMORE END SUBROUTINE GENNEXTSYM SUBROUTINE DOSYMLIMDEGEN(IMAX, NBASISMAX) IMPLICIT NONE TYPE(BasisFN) IMax(2) INTEGER nBasisMax(5, *), I IF (NBASISMAX(3, 3) == 0) THEN DO I = 1, 3 IF (IMax(2)%k(I) /= IMAX(1)%k(I)) IMAX(1)%k(I) = 0 end do end if ! always a spin symmetry IF (IMAX(1)%Ms /= IMAX(2)%Ms) IMAX(1)%Ms = 0 END SUBROUTINE DOSYMLIMDEGEN SUBROUTINE GETSYMDEGEN(ISYM, NBASISMAX, IDEGEN) IMPLICIT NONE TYPE(BasisFN) ISym, ISym2 INTEGER nBasisMax(5, *), IDEGEN, I, J LOGICAL TDO IDEGEN = 0 IF (NBASISMAX(3, 3) == 0) THEN ! Hubbard DO I = 0, 7 TDO = .TRUE. DO J = 1, 3 IF (.NOT. BTEST(I, J - 1)) THEN ISYM2%k(J) = ISYM%k(J) ELSE ISYM2%k(J) = -ISYM%k(J) IF (ISYM%k(J) == 0) TDO = .FALSE. end if end do IF (TDO .AND. KALLOWED(ISYM2, NBASISMAX)) IDEGEN = IDEGEN + 1 end do ELSE IDEGEN = 1 end if ! Spin IF (ISYM%Ms /= 0) IDEGEN = IDEGEN * 2 END SUBROUTINE GETSYMDEGEN ! Initialize symmetry to take into account the core electrons elemental subroutine setupsym(isym) type(basisfn), intent(inout) :: isym isym = frozensym end subroutine setupsym SUBROUTINE WRITEALLSYM(IUNIT, SYM, LTERM) IMPLICIT NONE INTEGER IUNIT TYPE(BASISFN) SYM LOGICAL LTERM write(IUNIT, "(4I5)", advance='no') SYM%K(1), SYM%K(2), SYM%K(3), SYM%MS CALL WRITESYM(IUNIT, SYM%SYM, LTERM) END SUBROUTINE WRITEALLSYM SUBROUTINE WRITESYM(IUNIT, SYM, LTERM) IMPLICIT NONE INTEGER IUNIT TYPE(SYMMETRY) SYM LOGICAL LTERM INTEGER Abel(3) if (t_k_space_hubbard) then write(iunit, "(I4)", advance='no') Sym return end if IF (TAbelian) THEN CALL DecomposeAbelianSym(SYM%s, Abel) if (TwoCycleSymGens) then write(IUNIT, '(" (",I2,",",I2,",",I2,")",I2)', advance='no') Abel(1:3), SYM%s else write(IUNIT, '(" (",I2,",",I2,",",I2,")",I2)', advance='no') Abel(1:3) end if else if (NSYM <= 16) THEN write(IUNIT, "(Z5)", advance='no') SYM else if (NSYM <= 24) THEN write(IUNIT, "(Z7)", advance='no') SYM else if (NSYM <= 32) THEN write(IUNIT, "(Z9)", advance='no') SYM else if (NSYM <= 40) THEN write(IUNIT, "(Z11)", advance='no') SYM else if (NSYM <= 48) THEN write(IUNIT, "(Z13)", advance='no') SYM else if (NSYM <= 56) THEN write(IUNIT, "(Z15)", advance='no') SYM ELSE write(IUNIT, "(Z17)", advance='no') SYM end if IF (LTERM) write(IUNIT, *) END SUBROUTINE WRITESYM SUBROUTINE SetupFreezeAllSym(Sym) IMPLICIT NONE TYPE(BasisFN) Sym ! Set to be totally symmetric FrozenSym = Sym END SUBROUTINE SetupFreezeAllSym SUBROUTINE SetupFreezeSym(Sym) IMPLICIT NONE TYPE(BasisFN) Sym ! FrozenSym=NullBasisFn ! Set to be totally symmetric FrozenSym = Sym END SUBROUTINE SetupFreezeSym !Deal with K-point symmetries, using translational symmetry operations. ! JSS: use Abelian symmetry formulation (allows us to go beyond 64 ! symmetry operations, and hence deal with larger k-point meshes). SUBROUTINE GenKPtIrreps(nTranslat, nKps, KpntInd, nStates) IMPLICIT NONE INTEGER I, nStates INTEGER nTranslat, nKps, KpntInd(nStates) character(*), parameter :: this_routine = 'GenKPtIrreps' nSymLabels = nKps nRot = nTranslat allocate(SymLabelChars(nRot, nSymLabels)) call LogMemAlloc('SymLabelChars', nSymLabels * nRot, 16, this_routine, tagSymLabelChars) allocate(SymLabels(nSymLabels)) call LogMemAlloc('SymLabels', nSymLabels, 4, this_routine, tagSymLabels) allocate(SymClasses(nStates)) call LogMemAlloc('SymClasses', nStates, 4, this_routine, tagSymClasses) SYMLABELCHARS = 0.0_dp DO I = 1, nStates SymClasses(I) = KpntInd(I) SymLabels(KPntInd(I))%s = ComposeAbelianSym(KpntSym(:, KPntInd(I))) END DO write(stdout, *) write(stdout, '(a11," |",a13,"|",a10)') ' K-vector', ' Label ', 'Conjugate' write(stdout, '(39("-"))') do i = 1, nSymLabels write(stdout, '("(",3i3,")"," | ")', advance='no') KpntSym(:, I) call writesym(stdout, SymLabels(I), .false.) write(stdout, '(A)', advance='no') " | " call writesym(stdout, SymConj(SymLabels(I)), .true.) end do ! write(stdout,'(/,a)') 'Symmetry Multiplication Table' ! do i=1,nSymLabels ! do j=1,nSymLabels ! write(stdout,'(z12)',advance='no') SymProd(SymLabels(I),SymLabels(J)) ! end do ! write(stdout,*) ! end do ! write(stdout,'(/)') ! write(stdout,*) "SYMMETRY CLASSES" ! CALL writeirreptab(stdout, SYMLABELCHARS,NROT,NSYMLABELS) !.. Allocate memory gor irreps. !.. Assume there will be no more than 64 irreps ! CALL N_MEMORY(IP_IRREPCHARS,NROT*64*2,"IRREPCH") END SUBROUTINE GenKPtIrreps pure subroutine DecomposeAbelianSym(ISym, AbelSym) ! Store the symmetry index as integer(int64). For Abelian symmetry ! we need to have 3 numbers stored in this index. We store ! according to isym=\sum_i AbelSym(i)*32768**(i-1). ! This allows point groups with more than 64 irreps to be used in ! the point group is Abelian (as all translational/k-point ! symmetries are). ! Decompose the symmetry label back into the appropriate "quantum" ! numbers. ! Store the symmetry index as integer(int64). For Abelian symmetry ! we need to have 3 numbers stored in this index. We store ! according to isym=1+\sum_i AbelSym(i)*32768**(i-1). ! Note that many symmetry parameters for CPMD-NECI jobs are set in ! kpntrep.F in CPMD source. ! Decompose the symmetry label back into the appropriate ! numbers... integer(int64), intent(in) :: Isym integer, intent(out) :: AbelSym(3) !RShift AbelSym(3) = int(IShft(Isym, -(PropBitLen * 2))) !RShift AbelSym(2) = int(Iand(IShft(ISym, -PropBitLen), 2_int64**PropBitLen - 1)) AbelSym(1) = int(Iand(Isym, 2_int64**PropBitLen - 1)) end subroutine DecomposeAbelianSym integer(int64) pure function ComposeAbelianSym(AbelSym) integer, intent(in) :: AbelSym(3) integer(int64) :: TempVar TempVar = AbelSym(3) !LShift ComposeAbelianSym = IShft(Tempvar, PropBitLen * 2) & + IShft(AbelSym(2), PropBitLen) & + AbelSym(1) end function ComposeAbelianSym pure function TotSymRep() ! Our definition of the totally symmetric representation ! changes according to whether we're using Abelian/k-point ! symmetry or the standard symmetry. It's just a matter of ! convenience, rather than some deep theoretical insight! Type(Symmetry) TotSymRep if (TAbelian .or. tUEG) then TotSymRep%s = 0 else TotSymRep%s = 1 end if end function TotSymRep ! nBasisMax might well be needed in the future in these functions. integer(int64) function MinSymRep() implicit none if (TAbelian) then MinSymRep = 0 else MinSymRep = 0 end if end function MinSymRep integer(int64) function MaxSymRep() implicit none integer abel(3) abel(:) = nprop(:) - 1 if (TAbelian) then MaxSymRep = ComposeAbelianSym(abel) else MaxSymRep = 0 end if end function MaxSymRep subroutine IncrSym(Sym) implicit none type(Symmetry) Sym integer abel(3), i logical lcont if (TAbelian) then call DecomposeAbelianSym(Sym%s, abel) i = 1 lcont = .true. do while (i < 4 .and. lcont) abel(i) = mod(abel(i) + 1, nprop(i)) lcont = abel(i) == 0 i = i + 1 end do Sym%s = ComposeAbelianSym(abel) else Sym%s = 0 end if end subroutine IncrSym SUBROUTINE GETSYMTMATSIZE(Nirrep, nBasis, iSS, iSize) implicit none integer Nirrep, nBasis, iSS, nBi, i, basirrep, t integer(int64) iSize character(*), parameter :: this_routine = 'GetSymTMATSize' nBi = nBasis / iSS iSize = 0 allocate(SymLabelIntsCum(nIrrep)) call LogMemAlloc('SymLabelIntsCum', nIrrep, 4, this_routine, tagSymLabelIntsCum) allocate(SymLabelCountsCum(nIrrep)) call LogMemAlloc('SymLabelCountsCum', nIrrep, 4, this_routine, tagSymLabelCountsCum) SYMLABELINTSCUM(1:Nirrep) = 0 SYMLABELCOUNTSCUM(1:Nirrep) = 0 do i = 1, Nirrep basirrep = SYMLABELCOUNTS(2, i) iSize = iSize + (basirrep * (basirrep + 1)) / 2 SYMLABELINTSCUM(i) = int(iSize) IF (i == 1) THEN SYMLABELCOUNTSCUM(i) = 0 ELSE DO t = 1, (i - 1) SYMLABELCOUNTSCUM(i) = SYMLABELCOUNTSCUM(i) + SYMLABELCOUNTS(2, t) end do end if write(stdout, *) basirrep, SYMLABELINTSCUM(i), SYMLABELCOUNTSCUM(i) call neci_flush(stdout) end do iSize = iSize + 2 !This is to allow the index of '-1' in the array to give a zero value END SUBROUTINE GETSYMTMATSIZE ! This function returns the label (0 -> nSymlabels-1) of the symmetry ! product of two symmetry labels. PURE INTEGER FUNCTION RandExcitSymLabelProd(SymLabel1, SymLabel2) IMPLICIT NONE INTEGER, INTENT(IN) :: SymLabel1, SymLabel2 IF (tNoSymGenRandExcits) THEN RandExcitSymLabelProd = 0 else if (tKPntSym) THEN !Look up the symmetry in the product table for labels (returning labels, not syms) RandExcitSymLabelProd = SymTableLabels(SymLabel1, SymLabel2) ELSE RandExcitSymLabelProd = IEOR(SymLabel1, SymLabel2) ! write(stdout,*) "***",SymLabel1,SymLabel2,RandExcitSymLabelProd end if END FUNCTION RandExcitSymLabelProd ! this function checks whether a determinant nI can appear if the ! total momentum is to be targetK and determinants with momentum ! cK have already been picked. rEls is the number of electrons that can ! still be distributed recursive function checkMomentumInvalidity(nI, cK, targetK, rElsUp, & rElsDown) result(momcheck) implicit none integer, intent(in) :: nI, rElsUp, rElsDown integer, dimension(3), intent(in) :: cK integer, dimension(3), intent(in) :: targetK integer, dimension(3) :: bufK logical :: momcheck integer :: lI, rElsUpNew, rElsDownNew ! returns true if nI can not appear momcheck = .true. if (G1(brr(nI))%Ms == -1) then rElsUpNew = rElsUp - 1 rElsDownNew = rElsDown else rElsDownNew = rElsDown - 1 rElsUpNew = rElsUp end if ! in case there are no electrons left matching the spin of nI if (rElsUpNew < 0 .or. rElsDownNew < 0) return ! this is the momentum from which we want to reach targetK bufK = cK + G1(brr(nI))%k if (tHub .or. Tperiodicinmom) then if (t_k_space_hubbard) then ! add in a way to never leave the first BZ instead of mapping! bufK = lat%add_k_vec(cK, G1(brr(nI))%k) else call MomPbcSym(bufK, nBasisMax) end if end if ! for the last electron, the total momentum has to be hit if ((rElsUpNew + rElsDownNew) == 0) then if (all(abs(bufK - targetK) == 0)) momcheck = .false. return end if ! try if the target momentum can still be reached ! this requires some computational effort, but the number of orbitals is not that large do lI = nI + 1, nbasis momcheck = checkMomentumInvalidity(lI, bufK, targetK, rElsUpNew, rElsDownNew) if (.not. momcheck) return end do end function checkMomentumInvalidity LOGICAL FUNCTION KALLOWED(G, NBASISMAX) ! See if a given G vector is within a (possibly tilted) unit cell. ! Used to generate the basis functions for the hubbard model (or perhaps electrons in boxes) Type(BasisFN), intent(in) :: G INTEGER nBasisMax(5, *), NMAXX, I, J, AX, AY INTEGER KX, KY real(dp) MX, MY, XX, YY LOGICAL TALLOW TALLOW = .TRUE. IF (NBASISMAX(3, 3) == 1) THEN !.. spatial symmetries IF (G%k(1) /= 0) TALLOW = .FALSE. ELSEIF (NBASISMAX(3, 3) == 0) THEN !.. Hubbard IF (NBASISMAX(1, 3) == 1) THEN !.. Tilted hubbard NMAXX = NBASISMAX(1, 5) ! NMAXY= AX = NBASISMAX(1, 4) AY = NBASISMAX(2, 4) !.. (XX,YY) is the position of the bottom right corner of the unit cell XX = ((AX + AY) / 2.0_dp) * NMAXX YY = ((AY - AX) / 2.0_dp) * NMAXX MX = XX * AX + YY * AY MY = XX * AY - YY * AX I = G%k(1) J = G%k(2) KX = I * AX + J * AY KY = I * AY - J * AX IF (KX > MX) TALLOW = .FALSE. IF (KY > MY) TALLOW = .FALSE. IF (KX <= -MX) TALLOW = .FALSE. IF (KY <= -MY) TALLOW = .FALSE. ELSEIF (NBASISMAX(1, 3) >= 4 .OR. NBASISMAX(1, 3) == 2) THEN !.. Real space Hubbard IF (G%k(1) == 0 .AND. G%k(2) == 0 .AND. G%k(3) == 0) THEN TALLOW = .TRUE. ELSE TALLOW = .FALSE. END IF ! ELSEIF(NBASISMAX(1,3).EQ.2) THEN !.. mom space non-pbc non-tilt hub - parity sym ! IF( (G(1).EQ.0.OR.G(1).EQ.1) ! & .AND.(G(2).EQ.0.OR.G(2).EQ.1) ! & .AND.(G(3).EQ.0.OR.G(3).EQ.1)) THEN ! TALLOW=.TRUE. ! ELSE ! TALLOW=.FALSE. ! ENDIF ! ELSEIF(NBASISMAX(1,3).EQ.2) THEN !.. non-pbc hubbard ! TALLOW=.TRUE. ! IF(G(1).GT.NBASISMAX(1,2).OR.G(1).LT.NBASISMAX(1,1)) ! & TALLOW=.FALSE. ! IF(G(2).GT.NBASISMAX(2,2).OR.G(2).LT.NBASISMAX(2,1)) ! & TALLOW=.FALSE. ! IF(G(3).GT.NBASISMAX(3,2).OR.G(3).LT.NBASISMAX(3,1)) ! & TALLOW=.FALSE. ELSE !.. Normal Hubbard TALLOW = .TRUE. IF (G%k(1) > NBASISMAX(1, 2) .OR. G%k(1) < NBASISMAX(1, 1)) TALLOW = .FALSE. IF (G%k(2) > NBASISMAX(2, 2) .OR. G%k(2) < NBASISMAX(2, 1)) TALLOW = .FALSE. IF (G%k(3) > NBASISMAX(3, 2) .OR. G%k(3) < NBASISMAX(3, 1)) TALLOW = .FALSE. END IF END IF KALLOWED = TALLOW RETURN END FUNCTION KALLOWED end module sym_mod