Assume that $V$ is a vector field on a 
Riemannian  manifold $(M,g)$ with natural  volum form $\Omega$ arising from $g$.
Assume that the solution curves of $V$ are parametrized geodesics of the Riemannian meteic $g$.

>Is it true to say that the space of harmonic functions are invariant under the derivational operator $D(f)=V.f=df(V)$?