Let us consider gravitational field in $\lambda$-family gauge conditions (for details see Sec. 5.5 in Nucl.Phys. B485 (1997) 517-544). The action with gauge-fixing term can be written as follows: \[ S = \int d^4 x \sqrt{-g} \left( R \,-\, \frac{1}{2}g_{\mu\nu} \chi^\mu \chi^\nu \right), \] where \[ \chi^\mu = \frac{1}{\sqrt{1 + \lambda}}\left( g^{\mu\alpha} \nabla^\beta h_{\alpha\beta} - \frac{1}{2}g^{\alpha\beta}\nabla^\mu h_{\alpha\beta} \right) \]
Calculating the second variation and symmetrising over indices by commuting covariant derivatives, it is easy to obtain all required input for the algorithm:
\begin{eqnarray*} &&K^{\mu\nu}{}_{\alpha\beta}{}^{\gamma\delta} = g^{\mu\nu} \delta_{\alpha\beta}{}^{\gamma\delta} -\ \frac{\lambda}{4(1+\lambda)} \left( \delta_\alpha{}^\gamma \delta_\beta{}^\mu g^{\delta\nu} + \delta_\alpha{}^\gamma \delta_\beta{}^\nu g^{\delta\mu} + \delta_\alpha{}^\delta \delta_\beta{}^\mu g^{\gamma\nu} \right. \nonumber\\ && \vphantom{\frac{1}{2}} + \delta_\alpha{}^\delta \delta_\beta{}^\nu g^{\gamma\mu} + \delta_\beta{}^\gamma \delta_\alpha{}^\mu g^{\delta\nu} + \delta_\beta{}^\gamma \delta_\alpha{}^\nu g^{\delta\mu} + \delta_\beta{}^\delta \delta_\alpha{}^\mu g^{\gamma\nu} \left. + \delta_\beta{}^\delta \delta_\alpha{}^\nu g^{\gamma\mu} \right) \nonumber\\ &&+ \frac{\lambda}{2(1+\lambda)} g^{\gamma\delta} \left(\delta_\alpha{}^\mu \delta_\beta{}^\nu + \delta_\alpha{}^\nu \delta_\beta{}^\mu \right),\\ \nonumber\\ &&W_{\alpha\beta}{}^{\gamma\delta} = P_{\alpha\beta}{}^{\gamma\delta} -\ \frac{\lambda}{2(1+\lambda)} \left(R_\alpha{}^\gamma{}_\beta{}^\delta + R_\alpha{}^\delta{}_\beta{}^\gamma\right) +\ \frac{\lambda}{4(1+\lambda)} \left( \delta_\alpha{}^\gamma R_\beta{}^\delta \right. \nonumber\\ && \vphantom{\frac{1}{2}} \left. + \delta_\alpha{}^\delta R_\beta{}^\gamma + \delta_\beta{}^\gamma R_\alpha{}^\delta + \delta_\beta{}^\delta R_\alpha{}^\gamma \right),\\ \nonumber\\ &&(Kn)^{-1}{}_{\alpha\beta}{}^{\gamma\delta} = \delta_{\alpha\beta}^{\gamma\delta} + \frac{\lambda}{2} \left( \delta_\alpha{}^\gamma n_\beta n^\delta + \delta_\alpha{}^\delta n_\beta n^\gamma + \delta_\beta{}^\gamma n_\alpha n^\delta + \delta_\beta{}^\delta n_\alpha n^\gamma \right) \nonumber\\ &&- \lambda g^{\gamma\delta} n_\alpha n_\beta, \end{eqnarray*} where \begin{eqnarray*} &&\delta^{\mu\nu}_{\alpha\beta} = \frac{1}{2} \left( \delta_\alpha^\mu \delta_\beta^\nu + \delta_\alpha^\nu \delta_\beta^\mu \right), \nonumber\\ &&P_{\gamma\delta}{}^{\mu\nu} = R_\gamma{}^\mu{}_\delta{}^\nu + R_\gamma{}^\nu{}_\delta{}^\mu + \frac{1}{2} \left(\delta_\gamma^\mu R_\delta{}^\nu + \delta_\gamma^\nu R_\delta{}^\mu + \delta_\delta^\mu R_\gamma{}^\nu + \delta_\delta^\nu R_\gamma{}^\mu\right) - g^{\mu\nu} R_{\gamma\delta} \nonumber\\ && - g_{\gamma\delta} R^{\mu\nu} - \delta^{\mu\nu}_{\gamma\delta} R + \frac{1}{2} g_{\gamma\delta} g^{\mu\nu} R. \end{eqnarray*}
Having all this, we can do in Redberry:
//setup symmetries of Riemann tensor
addSymmetries 'R_abcd', -[[0, 1]].p, [[0, 2], [1, 3]].p
setSymmetric 'R_ab'
//tensor (Kn)^{-1}
def iK = '''iK_ab^cd =
(d_a^c*d_b^d+d_b^c*d_a^d)/2
+ la/2*( d_a^c*n_b*n^d + d_a^d*n_b*n^c
+ d_b^c*n_a*n^d + d_b^d*n_a*n^c)
- la*g^cd*n_a*n_b'''.t
//tensor K
def K = '''K^lm_ab^cd =
g^lm*(d_a^c*d_b^d + d_b^c*d_a^d)/2
-la/4/(1 + la)*(
d_a^c*d_b^l*g^dm + d_a^c*d_b^m*g^dl
+ d_a^d*d_b^l*g^cm + d_a^d*d_b^m*g^cl
+ d_b^c*d_a^l*g^dm + d_b^c*d_a^m*g^dl
+ d_b^d*d_a^l*g^cm + d_b^d*d_a^m*g^cl)
+la/2/(1 + la)*g^cd*(d_a^l*d_b^m + d_a^m*d_b^l)'''.t
//tensor S
def S = '''S^p_{ab}^{cd} = 0'''.t
//tensor W
def W = '''W_{ab}^{cd} = P_ab^cd
-la/2/(1 + la)*(R_a^c_b^d + R_a^d_b^c)
+la/4/(1 + la)*( d_a^c*R_b^d + d_a^d*R_b^c
+d_b^c*R_a^d + d_b^d*R_a^c)'''.t
//auxiliary tensor
def P = '''P_cd^lm =
R_c^l_d^m+R_c^m_d^l
+ 1/2*( d_c^l*R_d^m + d_c^m*R_d^l
+d_d^l*R_c^m + d_d^m*R_c^l)
- g^lm*R_cd - R^lm*g_cd
+ (-d_c^l*d_d^m - d_c^m*d_d^l + g^lm*g_cd)*R/2'''.t
//substitute P in W
W <<= P
//curvature on principle bundle
def F = '''F_lm^kd_pr = R^k_plm*d^d_r + R^d_rlm*d^k_p'''.t
//main calculation
def counterterms = oneloopdiv2(iK, K, S, W, F).counterterms
//simplifying the result
counterterms <<= ExpandAll & Factor[[FactorScalars: false]]
println counterterms
> counterterms = (1/12)*R**2*(7+4*la**2)+(1/6)*(4*la+7+4*la**2)*R_{ab}*R^{ab}
In order to obtain one-loop counterterms in the dimensional regularization, one should multiply the result produced by Redberry by $1/16\pi(d-4)$ and integrate it over the space-time volume:
\[ \frac{1}{16\pi(d-4)} \int d^4 x \sqrt{-g} \left( \frac{1}{6} (4 \lambda^2 + 4 \lambda + 7) R^{\mu\nu} R_{\mu\nu} + \frac{1}{12}(4\lambda^2 + 7) R^2 \right) \]
To obtain full effective action counterterms, one should add contribution from Faddeev-Popov ghosts with the following Lagrangian: \[ \mathcal L_{ghosts} = \bar c^\alpha \left(\delta^\alpha_\beta \nabla^\mu \nabla_\mu + R^\alpha{}_\beta \right) c^\beta. \] The corresponding operator is just minimal second order operator, so substituting $W_{\mu\nu} = R_{\mu\nu}$ and $F_{\mu\nu\alpha\beta} = R_{\mu\nu\alpha\beta}$ to expressions obtained for minimal operator, we obtain ghosts counterterms: \[ -2 \times \frac{1}{16\pi(d-4)} \int d^4 x \sqrt{-g} \left( \frac{17}{60}R^2 + \frac{7}{30}R^{\mu\nu} R_{\mu\nu} \right). \]
Finally, summing the above results, we obtain complete counterterms: \[ \Gamma^{(1)}_{\infty} = \frac{1}{16\pi(d-4)} \int d^4 x \sqrt{-g} \left( \frac{1}{30} (20 \lambda^2 + 20 \lambda + 21) R^{\mu\nu} R_{\mu\nu} + \frac{1}{60}(20\lambda^2 + 1) R^2 \right) \]