Assume each summation is over all possible values of an index.
Assume the two indices of each Kronecker delta have the same domain.
Simplify Kronecker delta sum:
test = Sum[A[j, k, l]*KroneckerDelta[j, a]*KroneckerDelta[k, b]*
B[j, k, l, m] +
A[j, k, m]*KroneckerDelta[l, c]*KroneckerDelta[j, a]*
KroneckerDelta[k, b]*
B[j, k, l, m], {j, 1, J}, {k, 1, K}, {l, 1, L}, {m, 1, L}]
f[exp_] :=
exp /. {(x1_) + (x2__) :> f[x1] + f[Total[{x2}]],
Sum[(x1_) + (x2__), z__] :> f[Sum[x1, z]] + f[Sum[x2, z]],
Sum[c1___*(x1_ + x2__)*c2___, z__] :>
f[Sum[c1*c2*x1, z]] + f[Sum[c1*c2*x2, z]], (c__)*
Sum[x1__, z__]*(c2___) :> c*c2*f[Sum[x1, z]],
Sum[Sum[x1__, z1__], z2__] :> f[Sum[x1, z1, z2]],
Sum[(c1___)*KroneckerDelta[r_, s_]*(c2___), z1__, {s_, 1, p_},
z2___] :> f[Sum[c1*c2 /. s -> r, z1, z2]],
Sum[(c1___)*KroneckerDelta[r_, s_]*(c2___), z1___, {s_, 1, p_},
z2__] :> f[Sum[c1*c2 /. s -> r, z1, z2]],
Sum[(c1___)*KroneckerDelta[r_, s_]*(c2___), {s_, 1, p_}] :>
f[c1*c2 /. s -> r]}
f[test]