Simple calculation on weighted Riemannian manifolds gone wrong

Question

Let us consider a weighted Riemannian Manifold $M$ with measure $\nu(dx) = e^{-V(x)} dx$ such that $V\in C^2$. Let also $L = \Delta - \nabla V\cdot\nabla$, where $\Delta$ is the Laplace-Beltrami operator, and take $f\in C_c^\infty(M)$. The Reilly formula states that \begin{align*} (Lf)^2 = \lVert \text{Hess}(f)\rVert^2_{HS} + \nabla f\cdot\text{Ric}_{\nu}\nabla f ,\end{align*} so the following computation should be wrong, since it leads to a Hessian of $V$ instead of the Ricci curvature. I think I did not understand some basic thing on manifolds and did something I was not allowed to. \begin{align*} \int(Lf)^2 d\nu \quad = \quad & \int (\Delta f)^2 d\nu -2\int \Delta f(\nabla V\cdot\nabla f)d\nu + \int (\nabla V\cdot\nabla f)^2 d\nu \\ \stackrel{\text{IbP.}}{=} \quad & \int \nabla f\cdot(\Delta f\nabla V -\nabla\Delta f) d\nu - 2\int \Delta f(\nabla V\cdot\nabla f) d\nu \\ & + \int \Delta f\nabla V\cdot\nabla f+\nabla f\cdot\nabla(\nabla V\cdot\nabla f) d\nu \\ = \quad & - \int\nabla f\cdot\Delta\nabla f d\nu + \int\nabla f\cdot \text{Hess}(f)\nabla V d\nu + \int \nabla f\cdot\text{Hess}(V)\nabla f d\nu \end{align*} and then we have \begin{align*} \int \lVert \text{Hess}(f)\rVert^2_{HS}d\nu \quad = \quad & \sum_{i,j=1}^n\int(\partial_{ij}f)^2e^{-V}d x \\ \stackrel{\text{IbP.}}{=} \quad & \sum_{i,j=1}^n\int\partial_{j}f(\partial_{iij}f e^{-V}-\partial_{ij}f\partial_iVe^{-V})d x \\ = \quad & \int\nabla f\cdot\text{Hess}(f)\nabla V d\nu-\int \nabla f\cdot\Delta\nabla f d\nu .\end{align*} so together it gives \begin{align*} \int(Lf)^2 d\nu = \int \lVert \text{Hess}(f)\rVert^2_{HS} + \nabla f\cdot \text{Hess}(V)\nabla f d\nu ,\end{align*} which differs slightly from what should be true by the Reilly formula. What is my mistake? Any comment greatly appreciated! And sorry if it is a trivial mistake, I am not very experienced with manifolds.