#!/usr/bin/env cadabra2
{i,j,m,n,k,p,q,l,r,r#}::Indices(vector).
C_{m n p q}::WeylTensor.
\nabla{#}::Derivative.
\nabla_{r}{ C_{m n p q} }::SatisfiesBianchi.
\delta_{i j}::KroneckerDelta.

Eij:=- C_{i m k l} C_{j p k q} C_{l p m q} + 1/4 C_{i m k l} C_{j m p q} C_{k l p q}
     - 1/2 C_{i k j l} C_{k m p q} C_{l m p q}:

E:=  C_{j m n k} C_{m p q n} C_{p j k q} + 1/2 C_{j k m n} C_{p q m n} C_{j k p q}:

exp:= \nabla_{i}{\nabla_{j}{ @(Eij) }} - 1/6 \nabla_{i}{\nabla_{i}{ @(E) }};
distribute(_); product_rule(_); distribute(_); product_rule(_)
sort_product(_); canonicalise(_)
rename_dummies(_);
substitute(_, $\nabla_{i}{ C_{k i l m} } -> 0, \nabla_{i}{ C_{k m l i} } -> 0$);
young_project_product(_);