Divergence of Yamabe soliton In the article 
Ma, Li; Cheng, Liang, Properties of complete non-compact Yamabe solitons, Ann. Global Anal. Geom. 40, No. 3, 379-387 (2011). ZBL1225.53038. at page 382, there is a calculation. It says that if take divergence of both sides of the equation
$$\nabla^2f=Rg,$$ where $\nabla^2$ is the Hessian operator, $R$ is the scalar curvature and $f\in C^\infty(M)$, we get
$$\nabla_jR+\frac{1}{n-1}R_{jk}\nabla^kf=0.$$ Please show me the intermediate calculation. 

First, $\nabla^k(Rg_{jk})=\nabla_jR$ because $\nabla g=0$.
Second, commuting derivatives using the definition of the Ricci tensor implies that
$$ \nabla^k(f_{jk}) = \nabla^k\nabla_j\nabla_k f= \nabla_j\Delta f + R_{jk}\nabla^k f . $$
Third, the trace of the equation gives $\nabla_j\Delta f = n\nabla_j R$.
Combining these three equations gives the result.
