get_lmat_el Function

    function get_lmat_el(a, b, c, i, j, k) result(matel)
        ! Input: a,b,c - indices of orbitals to excite to
        !        i,j,k - indices of orbitals to excite from
        ! Output: matel - matrix element of this excitation, including all exchange terms
        use UMatCache, only: gtID
        use SystemData, only: t_exclude_pure_parallel
        ! Gets an entry of the 3-body tensor L:
        ! L_{abc}^{ijk} - triple excitation from abc to ijk
        integer, value :: a, b, c
        integer, value :: i, j, k
        HElement_t(dp) :: matel
        integer(int64) :: ida, idb, idc, idi, idj, idk

        ! initialize spin-correlator check: if all spins are the same, use LMat
        ! without spin-dependent correlator, always use LMat

        matel = h_cast(0.0_dp)

        if (t_exclude_pure_parallel) then
            if (     (G1(a)%ms == G1(b)%ms .and. G1(a)%ms == G1(c)%ms) &
                .or. (G1(i)%ms == G1(j)%ms .and. G1(i)%ms == G1(k)%ms)) then
                return
            end if
        end if

        ! convert to spatial orbs if required
        ida = gtID(a)
        idb = gtID(b)
        idc = gtID(c)
        idi = gtID(i)
        idj = gtID(j)
        idk = gtID(k)

        ! only add the contribution if the spins match

        ! here, we add up all the exchange terms
        call addMatelContribution(i, j, k, idi, idj, idk, 1)
        call addMatelContribution(j, k, i, idj, idk, idi, 1)
        call addMatelContribution(k, i, j, idk, idi, idj, 1)
        call addMatelContribution(j, i, k, idj, idi, idk, -1)
        call addMatelContribution(i, k, j, idi, idk, idj, -1)
        call addMatelContribution(k, j, i, idk, idj, idi, -1)

    contains

        subroutine addMatelContribution(p, q, r, idp, idq, idr, sgn)
            ! get a single entry of the LMat array and add it to the matrix element
            integer(int64), value :: idp, idq, idr
            integer, value :: p, q, r
            integer, intent(in) :: sgn
            integer(int64) :: index
            integer :: spinPos
            class(lMat_t), pointer :: lMatPtr
            real(dp) :: lMatVal
            !     integer(int64) :: ai,bj,ck

            if (G1(p)%ms == G1(a)%ms .and. G1(q)%ms == G1(b)%ms .and. G1(r)%ms == G1(c)%ms) then

                if (tContact) then
                    lMatVal = get_lmat_ueg(ida, idb, idc, idp, idq, idr)
                else if (tLMatCalc .or. t_rs_factors) then
                    lMatVal = lMatCalc(ida, idb, idc, idp, idq, idr)

                else
                    ! the indexing function is contained in the lMat object
                    index = lMat%indexFunc(ida, idb, idc, idp, idq, idr)
                    lMatVal = real(lMat%get_elem(index), dp)
                end if
                matel = matel + sgn * lMatVal
            end if

        end subroutine addMatelContribution

    end function get_lmat_el