NEWTMatInd Function

    INTEGER FUNCTION NEWTMatInd(I, J)
        ! In:
        !    i,j: spin orbitals.
        ! Return the index of the <i|h|j> element in TMatSym2.
        ! See notes for TMatInd. Used post-freezing.
        IMPLICIT NONE
        INTEGER I, J, 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
            NEWTMatInd = -1
            RETURN
        end if
        !To convert to spatial orbitals
        K = (I + 1) / 2
        L = (J + 1) / 2

        symI = SYMCLASSES2(K)
        symJ = SYMCLASSES2(L)

        IF (symI /= symJ) THEN
            NEWTMatInd = -1
            RETURN
        ELSE
            IF (symI == 1) THEN
                Block = 0
            ELSE
                Block = SYMLABELINTSCUM2(symI - 1)
            end if
            ! Convert K and L into their index within their symmetry class.
            K = StateSymMap(K)
            L = StateSymMap(L)
            IF (K >= L) THEN
                ind = (K * (K - 1)) / 2 + L
            ELSE
                ind = (L * (L - 1)) / 2 + K
            end if
            NEWTMatInd = Block + ind
            RETURN
        end if
    END FUNCTION NewTMatInd