$$g^{\mu\nu}\delta R_{\mu\nu}=\nabla_{\sigma}(g^{\mu\nu}\delta\Gamma^{\sigma}_{\nu\mu}-g^{\mu\sigma}\Gamma^{\rho}_{\rho\mu}),$$

and this term can be converted to a total derivative when it's multiplied by ##\sqrt{-g}##:

$$\sqrt{-g}\nabla_{\sigma}(g^{\mu\nu}\delta\Gamma^{\sigma}_{\nu\mu}-g^{\mu\sigma}\Gamma^{\rho}_{\rho\mu})=\partial_{\sigma}[\sqrt{-g}(g^{\mu\nu}\delta\Gamma^{\sigma}_{\nu\mu}-g^{\mu\sigma}\Gamma^{\rho}_{\rho\mu})].$$

When we integrate it over the whole spacetime and

**, it will contribute nothing to the variation of the action ##\int dx^{4}\sqrt{-g}f(R)##.**

__assume that the variation of the metric ##\delta g^{\mu\nu}## vanishes at infinity__However i'm considering things in the spatially flat Friedmann universe ##ds^{2}=dt^{2}-a^2(t)\delta_{ik}dx^{i}dx^{k}.## The assumption above can't be true any more. To make my question specific, could you compute ##\delta R/\delta g^{\mu\nu}## in the spatially flat Friedmann universe? Thanks a lot.