Let us consider one-loop counterterms of the the vector field operator, which appears in the theory of the massive vector field: \[ D_\alpha{}^\beta = \delta_\alpha{}^\beta \Box - \lambda \, \nabla_\alpha \nabla^\beta + P_\alpha{}^\beta, \] where \(\Box =g^{\mu\nu}\nabla_\mu\nabla_\nu\) and \(\lambda = 1 + 1/\xi\). This is a second order operator, and in order to rewrite it in the symmetric form, it is necessary to symmetrize the second term by commutation of the covariant derivatives:
\[ D_\alpha{}^\beta = \left( g^{\mu\nu} \delta_\alpha^\beta - \frac{\lambda}{2} \left( g^{\mu\beta}\delta_\alpha^\nu + g^{\nu\beta}\delta_\alpha^\mu \right) \right) \nabla_\mu\nabla_\nu + P_\alpha{}^\beta + \frac{\lambda}{2}R_\alpha{}^\beta, \]
where $R_{\alpha\beta}$ is the Ricci tensor.
It can be easily found that $Kn$ and $(Kn)^{-1}$ are:
\begin{equation*} (Kn)_{\alpha}{}^\beta = \delta_{\alpha}^\beta - \lambda \, n_\alpha n^\beta, \qquad (Kn)^{-1}{}_{\alpha}{}^\beta = \delta_{\alpha}^\beta + \frac{\lambda}{1- \lambda} n_\alpha n^\beta. \end{equation*}
Hereby, at this point we have whole set of input tensors required by the algorithm: \begin{eqnarray*} &&K^{\mu\nu}{}_\alpha{}^\beta \,=\, g^{\mu\nu} \delta_\alpha^\beta - \frac{\lambda}{2}\left( g^{\mu\beta}\delta_\alpha^\nu + g^{\nu\beta}\delta_\alpha^\mu \right)\,, \quad S^\mu{}_\alpha{}^\beta \,=\, 0\,, \quad W_\alpha{}^\beta \,=\, P_\alpha{}^\beta + \frac{\lambda}{2}R_\alpha{}^\beta\,, \\&& (Kn)^{-1}{}_{\alpha}{}^\beta \,=\, \delta_{\alpha}^\beta + \frac{\lambda}{1- \lambda} n_\alpha n^\beta\,, \quad F_{\mu\nu}{}_\alpha{}^\beta \,=\, R_{\mu\nu}{}_\alpha{}^\beta\,. \end{eqnarray*}
In the further calculations we shall use the definition $\lambda =\gamma/(1+\gamma)$ for convenience (so $\gamma = \lambda/(1-\lambda)$).
The following code calculates one-loop counterterms of the vector field theory in curved space-time (here g used for $\gamma$):
//setup symmetries of Riemann tensor
addSymmetries 'R_abcd', -[[0, 1]].p, [[0, 2], [1, 3]].p
setSymmetric 'R_ab', 'P_ab'
//tensor (Kn)^{-1}
def iK = 'iK_a^b = d_a^b + g*n_a*n^b'.t
//tensor K
def K = '''K^{mn}_i^j =
g^{mn}*d_{i}^{j}
- g/(2*(1+g))*(g^{mj}*d_i^n + g^{nj}*d_i^m)'''.t
//tensor S
def S = 'S^p^l_m=0'.t
//tensor W
def W = 'W^{a}_{b}=P^{a}_{b} + g/(2*(1+g))*R^a_b'.t
//tensor F
def F = 'F_abcd = R_abcd'.t
//divergent part of one-loop effective action
def div = oneloopdiv2(iK, K, S, W, F)
//counterterms
def counterterms = div.counterterms[1];
//simplifying counterterms
counterterms = ('P^a_a = P'.t
& Collect['R', 'P', Factor[[FactorScalars: false]]]
) >> counterterms
println counterterms
> (1/120)*(10*g-32+5*g**2)*R^{lm}*R_{lm}+(1/24)*P*R*(2*g+g**2+4)
+(1/12)*g*(g+4)*P^{ml}*R_{lm}+(1/48)*g**2*P**2
+(1/240)*(20*g+28+5*g**2)*R**2
+(1/24)*(6*g+g**2+12)*P^{a_{5}}_{a}*P^{a}_{a_{5}}
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:
\begin{multline*} \Gamma^{(1)}_{\infty} = \frac{1}{16\pi(d-4)} \int d^4 x \sqrt{-g} \left( \frac{1}{120}(-32+5 \gamma^2+10 \gamma) R_{\epsilon\mu} R^{\epsilon\mu} +\frac{1}{48}\gamma^2 P^2 +\right.\\ +\frac{1}{240} R^2 (28+5 \gamma^2+20 \gamma)+\frac{1}{24} (\gamma^2+12+6 \gamma) P_{\beta\alpha} P^{\alpha\beta} +\\ \left.+\frac{1}{12}\gamma (4+\gamma) R_{\nu\epsilon} P^{\nu\epsilon} +\frac{1}{24} R (\gamma^2+4+2 \gamma) P\right) \end{multline*}