Background: If we have given some open, bounded subset $U$ of $\mathbb{R}^n $, smooth functions $a^{ij}(x)$ on $U$ and smooth functions $u, v: U \to \mathbb{R}$, where $v$ and its derivatives vanish on the boundary of $U$, then $$\int_U (a^{ij} u_{x_i})_{x_j} v dx = - \int _U a^{ij} u_{x_i} v_{x_j} dx, $$ where the index $_{x_i}$ means differentiation with respect to $x_i$.
My goal is to convert this setting to compact manifolds with no boundary.
Given a smooth, compact $n$-dimensional Riemannian manifold $(M,g)$ with no boundary. We also have given (in some coordinate chart) $Q^{ij}:M \to \mathbb{R}$ smooth functions. Does it follow that for $f_1,f_2 \in C^1(M)$ we have $$\int_M \nabla_j (\nabla _i ( Q^{ij} f_1)) f_2 dv(g) = - \int _M \nabla_i(Q^{ij} f_1) (\nabla_j f_2) dv(g)?$$ If so, how could I prove it? My attempt at proving was to bring both integrals over to one side and use the product rule for the covariant derivative, i.e. use that $$\nabla_j (\nabla_i(Q^{ij} f_1) f_2) = \nabla_j (\nabla_i(Q^{ij} f_1)) f_2 + \nabla_i(Q^{ij} f_1) \nabla_j(f_2).$$ However I'm not sure whether that's actually valid in this case and neither do I know how to proceed from here.