SymGenExcitIt_GetNextPair Subroutine

Subroutine SymGenExcitIt_GetNextPair(K, iTo, iLooped, iTo1, iTo2, tDebugPrint)
    use SymData, only: tStoreStateList
    use SymData, only: SymLabelCounts, SymStatePairs
    use SymData, only: SymLabelList, SymPairProds
    use SymData, only: SymClasses
    use constants, only: stdout
    IMPLICIT NONE
    INTEGER K, iTo
    INTEGER iLooped
    INTEGER iTo1, iTo2

    INTEGER iLabelPairIndexOff
    INTEGER iLabelPairBase
    INTEGER iLabel1Index, iLabel2Index
    INTEGER iLabel1, iLabel2
    LOGICAL tDebugPrint
    iLooped = 0
    if (.not. tStoreStateList) then
! We must generate the pair as we don't have a stored list
! SymStatePairs holds a list of SymLabels.
! SymPairProds(iTo) tells is which set of pairs we're using
! iTo is the symmetry of the 'to' spin-orb pair
! we split K into 3 integers as we need to store 3 variables rather than 1
!   (LabelPairIndexOff, Label1Index, Label2Index)
!bits      30-24         23-12           11-0
!The index of the current state of the 1st symmetry of the pair (Label1Index)
!The index of the current state of the 2nd symmetry of the pair (Label2Index)

!we loop over these with Label2Index changing fastest, then Label1Index, then iLabelPairIndexOff
! LabelPairIndexOff:  Which pair of symmetries (in SymStatePairs)
!                       This is the index for a given iTo.  The base is iLabelPairBase
! Label1Index:        The index of the current first state within its symmetry class
! Label2Index:        The index of the current second state within it symmetry class

!  the symmetries of the states are determined by the LAbelPairIndexOff.
!   We iterate through all possible [LabelPair,(state1,state2)], in order state2, state1, LabelPair.

        if (K < 0) then
            iLabelPairIndexOff = 0
            iLabel1Index = 0
            iLabel2Index = -1
            iLooped = -1 !indicate we've reset from nothing
        else
            iLabelPairIndexOff = ishft(iand(K, 2130706432), -24)
            iLabel1Index = ishft(iand(K, 16773120), -12)
            iLabel2Index = iand(K, 4095)
        end if
        iLabelPairBase = SymPairProds(iTo)%nIndex
        iLabel1 = SymStatePairs(1, iLabelPairBase + iLabelPairIndexOff)
        iLabel2 = SymStatePairs(2, iLabelPairBase + iLabelPairIndexOff)
! We cycle through these with Label2Index fastest
        iLabel2Index = iLabel2Index + 1
! This is a do rather than an if because if iLabel1==iLabel2, it's possible
!  that there are no valid combinations after an increment, so we have to go round again.
        if (iLabel2Index >= SymLabelCounts(2, iLabel2)) then
! Increment label 1
            iLabel1Index = iLabel1Index + 1
            if (iLabel1Index >= SymLabelCounts(2, iLabel1)) then
! Increment the Pair Index
                iLabelPairIndexOff = iLabelPairIndexOff + 1
                if (iLabelPairIndexOff >= SymPairProds(iTo)%nPairs) then
                    iLooped = 1 ! indicate we're looping
                    iLabelPairIndexOff = 0
                end if
                iLabel1 = SymStatePairs(1, iLabelPairBase + iLabelPairIndexOff)
                iLabel2 = SymStatePairs(2, iLabelPairBase + iLabelPairIndexOff)
                iLabel1Index = 0
            end if
!  If the two are from the same symlabel, we enforce 2's index>1
            if (iLabel1 == iLabel2) then
                iLabel2Index = iLabel1Index
            else
                iLabel2Index = 0
            end if
        end if
! Now extract the states and turn them into orbitals.
        iTo1 = (SymLabelList(SymLabelCounts(1, iLabel1) + iLabel1Index)) * 2 - 1
        iTo2 = (SymLabelList(SymLabelCounts(1, iLabel2) + iLabel2Index)) * 2 - 1
! Put the K index back together
        K = IShft(iLabelPairIndexOff, 24) + IShft(iLabel1Index, 12) + iLabel2Index
        if (tDebugPrint) THEN
            Write(stdout, *) "To: ", iTo, "LabelPair:", iLabelPairBase, "+", iLabelPairIndexOff, "/", SymPairProds(iTo)%nPairs
            Write(stdout, "(A10,I4,A1,I4,A1,I4,I4)") "Label1: (", iLabel1, ")", iLabel1Index, "/", SymLabelCounts(2, iLabel1), iTo1
            Write(stdout, "(A10,I4,A1,I4,A1,I4,I4)") "Label2: (", iLabel2, ")", iLabel2Index, "/", SymLabelCounts(2, iLabel2), iTo2
            WRITE(stdout, "(A,Z10,A,I3)") "K:", K, " Looped:", iLooped
        end if
    else
! We have a stored list, so just get the next from there.
        K = K + 1
        IF (K < 0 .or. K >= SymPairProds(ITO)%nIndex + SymPairProds(ITO)%nPairs) THEN
            iLooped = 1 !we've looped
            if (K < 0) iLooped = -1 ! we've reset
            K = SymPairProds(ITO)%nIndex
        END IF
        iTo1 = SymStatePairs(1, K) * 2 - 1
        iTo2 = SymStatePairs(2, K) * 2 - 1
        if (tDebugPrint) THEN
            Write(stdout, *) "To: ", iTo
            WRITE(stdout, "(A10,I4,A1,I4)") "Label1: (", SymClasses((iTo1 + 1) / 2), ")", iTo1
            WRITE(stdout, "(A10,I4,A1,I4)") "Label2: (", SymClasses((iTo2 + 1) / 2), ")", iTo2
            WRITE(stdout, "(A,I10,A,I3)") "K:", K, " Looped:", iLooped
        end if
    END IF
    RETURN
End Subroutine SymGenExcitIt_GetNextPair