bHubIndexFunction Function

    elemental function bHubIndexFunction(i, j, k, l) result(bInd)
        ! this function gets the index of the matrix element of the breathing term
        ! corresponding to the four states i,j,k,l

        ! not sure if it is too expensive to do the index generation on the fly
        ! maybe storing all matrix elements ignoring redundancies might be worthwile
        ! even though memory cost will grow significantly (but still no noteworthy memory usage)
        implicit none
        integer, intent(in) :: i, j, k, l
        integer :: bInd, buf(3), li(3), tgt
        ! to get the breathing effect, we need to know the exchanged momentum
        ! therefore, compare the momentum of i to that of k/l with the matching
        ! spin
        unused_var(j)

        unused_var(j)

        if (G1(i)%MS == G1(k)%MS) then
            tgt = k
        else
            tgt = l
        end if
        buf = G1(i)%k - G1(tgt)%k
        call MomPbcSym(buf, nBasisMax)
        bInd = momIndexTable(buf(1), buf(2), buf(3), G1(i)%MS)
    end function bHubIndexFunction