GENSymStatePairs Subroutine

    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