Unit harmonic vector field

The answer is 'no'. Here is a counterexample:

Let $S$ be a compact oriented surface endowed with a nonvanishing $2$-form $\Omega$. Let $f:S\to S$ be a diffeomorphism preserving $\Omega$ with a hyperbolic fixed point $s\in S$. Let $\mathbb{Z}$ act freely on $\mathbb{R}\times S$ by the rule $$ n\cdot(t,p) = \bigl(t+n, f^n(p)\bigr) $$ for all $n\in\mathbb{Z}$. Since $\pi_1^*\mathrm{d}t$ and $\pi_2^*\Omega$ are invariant under this action (where $\pi_1:\mathbb{R}\times S\to\mathbb{R}$ and $\pi_2:\mathbb{R}\times S\to S$ are the projections), they are well-defined on the quotient $M = (\mathbb{R}\times S)/\mathbb{Z}$. (I will drop the projection pullback notation henceforth.)

Using a partition of unity on the circle $S^1 = \mathbb{R}/\mathbb{Z}$, we can construct a smooth Riemannian metric $g$ on $M$ such that $g = \mathrm{d}t^2 + h(t)$, where $h(t)=h(t+1)$ is, for each $t\in\mathbb{R}$, a positive definite Riemannian metric on $S$ whose oriented area form is $\Omega$. Let $X$ be the vector field on $M$ that is null for $h(t)$ and satisfies $\mathrm{d}t(X)=1$. Then, with respect to $g$, $\mathrm{d}t$ is the $1$-form dual of $X$ and $\ast\mathrm{d} t = \Omega$. Since both $\mathrm{d}t$ and $\Omega$ are closed on $M$, it follows that $X$ is a unit harmonic vector field on $M$.

However, because the closed $\mathbb{Z}$-periodic $X$-orbit $\gamma(t) = [t,s]$ is hyperbolic, there cannot be any Riemannian metric on $M$ that is invariant under the flow of $X$. I.e., $X$ is not a Killing field for any Riemannian metric on $M$.