get_lmat_el_ua Function

    pure function get_lmat_el_ua(a, b, c, i, j, k) result(matel)
        use SystemData, only: G1
        use UMatCache, only: gtID
        ! Gets an entry of the 3-body tensor L:
        ! L_{abc}^{ijk} - triple excitation from abc to ijk
        integer, value :: a, b, c
        integer :: a2, b2, c2
        integer, intent(in) :: i, j, k
        HElement_t(dp) :: matel

        ! convert to spatial orbs if required

        matel = 0

        if (G1(a)%ms == G1(b)%ms .and. G1(a)%ms /= G1(c)%ms) then
            a2 = a
            b2 = b
            c2 = c
        else if (G1(a)%ms == G1(c)%ms .and. G1(a)%ms /= G1(b)%ms) then
            a2 = c
            b2 = a
            c2 = b
        else if (G1(b)%ms == G1(c)%ms .and. G1(a)%ms /= G1(b)%ms) then
            a2 = b
            b2 = c
            c2 = a
        else
            return
        end if

        ! only add the contribution if the spins match
        call addMatelContribution_ua(i, j, k, 1, matel)
        call addMatelContribution_ua(j, k, i, 1, matel)
        call addMatelContribution_ua(k, i, j, 1, matel)
        call addMatelContribution_ua(j, i, k, -1, matel)
        call addMatelContribution_ua(i, k, j, -1, matel)
        call addMatelContribution_ua(k, j, i, -1, matel)
    contains

        pure subroutine addMatelContribution_ua(p, q, r, sgn, matel)
            integer, value :: p, q, r
            integer, intent(in) :: sgn
            HElement_t(dp), intent(inout) :: matel

            if (G1(p)%ms == G1(a2)%ms .and. G1(q)%ms == G1(b2)%ms .and. G1(r)%ms == G1(c2)%ms) then
                matel = matel + 2.d0 * sgn * get_lmat_ua(a2, b2, c2, p, q, r)
            end if

        end subroutine addMatelContribution_ua

    end function get_lmat_el_ua