I am new to Cadabra and I'm trying to it to play with infinitesimal variation of some gravitatoinal terms. 
I am defining a "derivative" \delta, that I would like to commute with partial derivatives. I am doing the following , based on the example about Einstein equations
{\mu,\nu, \rho,\alpha,\beta#}::Indices(spacetime);
g_{\mu\nu}::Metric;
g^{\mu\nu}::InverseMetric;
{g_{\mu}^{\nu},g^{\mu}_{\nu}}::KroneckerDelta;
{\partial{#}}::PartialDerivative;
{\nabla{#},\delta{#}}::Derivative;
In the long run I'm interested in infinitesimal Weyl variations, so I define the following list of variations:
 variations := {
    \delta{g^{\mu \nu}} -> +2 \sigma g^{\mu \nu},
    \delta{g_{\mu \nu}} -> -2 \sigma g_{\mu \nu}
};
As a warm up I tried to apply that to the Christoffel Symbol:
expr := \delta{
    1/2 g^{\rho\alpha} * (
        \partial_{\mu}{g_{\nu\alpha}}
        + \partial_{\nu}{g_{\mu\alpha}}
        - \partial_{\alpha}{g_{\mu\nu}}
    )
};
product_rule(_);
distribute(_);
substitute(_, $\delta{\partial_{\mu}{A}} -> \partial_{\mu}{\delta{A}}$ );
substitute(_, variations);
The first substitute gives me the same as distribute(), namely it does not commute the variation and the partial derivtive. I tried instead replacing A by g_{\nu\rho} explicitely, but that does not do anything either. Is there something I am missing?