TMatInd Function

    pure INTEGER FUNCTION TMatInd(I, J)
        ! In:
        !    i,j: spin orbitals.
        ! Return the index of the <i|h|j> element in TMatSym2.
        ! We assume a restricted calculation.  We note that TMat is a Hermitian matrix.
        ! For the TMat(i,j) to be non-zero, i and j have to belong to the same symmetry, as we're working in Abelian groups.
        ! We store the non-zero elements in TMatSym(:).
        !
        ! The indexing scheme is:
        !   Arrange basis into blocks of the same symmetry, ie so TMat is of a block
        !   (but not necessarily block diagonal) nature.  Order the blocks by their
        !   symmetry index.
        !   Within each block, convert i and j to a spatial index and label by the
        !   (energy) order in which they come within that symmetry: i->k,j->l.
        !   We only need to store the upper diagonal of each block:
        !        blockind=k*(k-1)/2+l, l<k.
        !   The overall index depends on the number of non-zero integrals in the blocks
        !   preceeding the block i and j are in.  This is given by SYMLABELINTSCUM(symI-1).
        !        TMatInd=k*(k-1)/2+l + SymLabelIntsCum(symI-1).
        !   If the element is zero by symmetry, return -1 (TMatSym(-1) is set to 0).
        use SymData, only: SymClasses, StateSymMap, SymLabelIntsCum
        IMPLICIT NONE
        integer, intent(in) :: i, j
        integer A, B, symI, symJ, Block, ind, K, L
        A = mod(I, 2)
        B = mod(J, 2)
        !If TMatInd = -1, then the spin-orbitals have different spins, or are symmetry disallowed
!therefore have a zero integral (apart from in UHF - might cause problem if we want this)
        IF (A /= B) THEN
            TMatInd = -1
            RETURN
        end if
        !To convert to spatial orbitals
        K = (I + 1) / 2
        L = (J + 1) / 2

        symI = SYMCLASSES(K)
        symJ = SYMCLASSES(L)

        IF (symI /= symJ) THEN
            ! <i|t|j> is symmetry forbidden.
            TMatInd = -1
            RETURN
        ELSE
            ! Convert K and L into their index within their symmetry class.
            K = StateSymMap(K)
            L = StateSymMap(L)
            ! Block index: how many symmetry allowed integrals are there in the
            ! preceeding blocks?
            IF (symI == 1) THEN
                Block = 0
            ELSE
                Block = SYMLABELINTSCUM(symI - 1)
            end if
            ! Index within symmetry block.
            IF (K >= L) THEN
                ind = (K * (K - 1)) / 2 + L
            ELSE
                ind = (L * (L - 1)) / 2 + K
            end if
            TMatInd = Block + ind
            RETURN
        end if
    END FUNCTION TMatInd