We have the Gauss-Bonnet term
$$L_{GB}=R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$$
where $R$, $R_{\mu\nu}$ and $R_{\mu\nu\rho\sigma}$ are the Ricci scalar, the Ricci tensor and the Riemann tensor respectively. The variation of $L_{GB}$ with respect to the metric $g_{\mu\nu}$ should give $$\frac{\delta L_{GB}}{\delta g^{\mu\nu}}=2(RR_{\mu\nu}-2R_{\mu\alpha}R^{\alpha}_{\nu}-2R^{\alpha\beta}R_{\mu\alpha\nu\beta}+R_{\mu}^{\alpha\beta\gamma}R_{\nu\alpha\beta\gamma})$$
My attempt:
$$\delta R^{2}=2R\delta R=2R\delta(g^{\mu\nu}R_{\mu\nu})=2RR_{\mu\nu}\delta g^{\mu\nu}+2Rg^{\mu\nu}\delta R_{\mu\nu} \\ \frac{\delta R^{2}}{\delta g^{\mu\nu}}=2RR_{\mu\nu}$$
$$ \delta(R_{\mu\nu}R^{\mu\nu})=R^{\mu\nu}\delta R_{\mu\nu}+R_{\mu\nu}\delta R^{\mu\nu} =R^{\mu\nu}\delta R_{\mu\nu}+R_{\mu\nu}\delta(g^{\mu\rho}g^{\nu\sigma}R_{\rho\sigma}) \\ =R^{\mu\nu}\delta R_{\mu\nu}+R_{\mu\nu}g^{\mu\rho}g^{\nu\sigma}\delta R_{\rho\sigma}+R_{\mu\nu}g^{\mu\rho}R_{\rho\sigma}\delta g^{\nu\sigma}+R_{\mu\nu}g^{\nu\sigma}R_{\rho\sigma}\delta g^{\mu\rho} \\ =2R^{\mu\nu}\delta R_{\mu\nu}+2R_{\mu\alpha}R^{\alpha}_{\nu}\delta g^{\mu\nu} =2R^{\mu\nu}\delta(R^{\alpha}_{\mu\beta\nu}g^{\beta}_{\alpha})+2R_{\mu\alpha}R^{\alpha}_{\nu}\delta g^{\mu\nu} \\ = 2R^{\alpha\beta}g^{\mu\nu}\delta R_{\mu\alpha\nu\beta}+2(R^{\alpha\beta}R_{\mu\alpha\nu\beta}+R_{\mu\alpha}R^{\alpha}_{\nu})\delta g^{\mu\nu} \\ \frac{\delta(R_{\mu\nu}R^{\mu\nu})}{\delta g^{\mu\nu}}=2(R^{\alpha\beta}R_{\mu\alpha\nu\beta}+R_{\mu\alpha}R^{\alpha}_{\nu})$$
$$\delta (R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma})=R_{\mu\nu\rho\sigma}\delta R^{\mu\nu\rho\sigma}+R^{\mu\nu\rho\sigma}\delta R_{\mu\nu\rho\sigma}=R^{\mu\nu\rho\sigma}\delta R_{\mu\nu\rho\sigma}+R_{\mu\nu\rho\sigma}\delta (g^{\mu\alpha}g^{\nu\beta}g^{\rho\gamma}g^{\sigma\delta}R^{\alpha\beta\gamma\delta})=...=2R^{\mu\nu\rho\sigma}\delta R_{\mu\nu\rho\sigma}+4R_{\mu}^{\alpha\beta\gamma}R_{\nu\alpha\beta\gamma}\delta g^{\mu\nu}=4R_{\mu}^{\alpha\beta\gamma}R_{\nu\alpha\beta\gamma}\delta g^{\mu\nu}+2R^{\mu\nu\rho\sigma}\delta (R^{\alpha}_{\nu\rho\sigma}g_{\alpha\mu})=2R_{\mu}^{\alpha\beta\gamma}R_{\nu\alpha\beta\gamma}\delta g_{\mu\nu}+6R^{\mu\alpha\rho\sigma}g_{\mu\nu}\delta R^{\nu}_{\alpha\rho\sigma} \\ \frac{\delta (R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma})}{\delta g^{\mu\nu}}=2R_{\mu}^{\alpha\beta\gamma}R_{\nu\alpha\beta\gamma}$$
where I took $g^{\mu\nu}R\frac{\delta R_{\mu\nu}}{\delta g^{\mu\nu}}=0$, $R^{\alpha\beta}g^{\mu\nu}\frac{\delta R_{\mu\alpha\nu\beta}}{\delta g^{\mu\nu}}=0$ and $R^{\mu\alpha\rho\sigma}g_{\mu\nu}\frac{\delta R^{\nu}_{\alpha\rho\sigma}}{\delta g^{\mu\nu}}=0$. Eventhough I obtain the terms I expect, the coefficient for $\delta(R_{\mu\nu}R^{\mu\nu})/\delta g^{\mu\nu}$ should be 1 instead of 2. What am I missing?
