I took the time to give a clean version of the answer and to eliminate the dependence on really hard theorems.

**Construction of geodesible, volume-preserving flows on $S^3$ that are not Reeb flows for any contact form.**

Ingredients:

A contact form $\alpha$ whose Reeb vector field $X$ has no nonconstant integrals of motion and which carries a closed characteristic. Such a form can be obtained from a Finsler or Riemannian metric on the two-sphere with ergodic geodesic flow (examples by Katok and Donnay).

A smooth function $f$ that is nowhere zero and such that for some closed characteristic $\gamma$ in $(S^3, \alpha)$ we have that
$$
{ \int_{S^3} f^{-1} \ \alpha \wedge d\alpha \over \int_{S^3} \alpha \wedge d\alpha}
\neq
{\int_\gamma f^{-1} \ \alpha \over \int_\gamma \alpha}.
$$

**Theorem.**
*The vector field $fX$ preserves the volume form $f^{-1} \alpha \wedge d\alpha$, but it is not the Reeb vector field of any contact form on $S^3$.*

The proof is by contradiction and we will need the following

**Claim.**
*If $fX$ is the Reeb vector field of a contact form $\beta$, then for some nonzero constant $\lambda$ and some smooth function $g$ we have that $\beta = \lambda \alpha + dg$.*

*Proof of the claim.* Since the kernels of the $2$-forms
$d\alpha$ and $d\beta$ coincide and the space is three-dimensional, the forms are multiples:
$d\beta = \lambda d\alpha$, with $\lambda$ some smooth function. Moreover
$\lambda$ must be an integral of motion of $X$,
$$
0 = {\cal L}_X d\beta = {\cal L}_X \lambda d\alpha = X(\lambda) \alpha .
$$
By hypothesis $\lambda$ must be a constant and, therefore,
$d(\beta - \lambda \alpha) = 0$ or $\beta = \lambda \alpha + dg$.

Now we proceed with the proof of the theorem. First we show that $\lambda$ is the
average of the function $f^{-1}$ over $S^3$. Indeed,
$$
1 = \beta(fX) = \lambda \alpha (fX) + dg(fX) = \lambda f + fX(g),
$$
and so
$$
\int_{S^3} f^{-1} \ \alpha \wedge d\alpha =
\int_{S^3} (\lambda f + fX(g))f^{-1} \ \alpha \wedge d\alpha.
$$
Since $X$ preserves the form $\alpha \wedge d\alpha$, the integral of
$X(g)\alpha \wedge d\alpha$ is zero and we obtain
$$
\int_{S^3} f^{-1} \ \alpha \wedge d\alpha = \lambda \int_{S^3} \alpha \wedge d\alpha .
$$

Now we notice that $\lambda$ is also the average of $f^{-1}$ over *every* closed leaf of the characteritic foliation of $(S^3, \alpha)$. If $\gamma$ is a closed leaf, then
$$
\int_\gamma f^{-1} \ \alpha = \int_\gamma (\lambda f + fX(g)) f^{-1} \ \alpha.
$$
As before, using that $X$ preserves $\alpha$, we obtain that the integral of
$X(g)\alpha$ along $\gamma$ is zero, and hence
$$
\int_\gamma f^{-1} \ \alpha = \lambda \int_\gamma \alpha.
$$

This contradicts our assumption of the existence of a closed characteristic for which
$$
{ \int_{S^3} f^{-1} \ \alpha \wedge d\alpha \over \int_{S^3} \alpha \wedge d\alpha}
\neq
{\int_\gamma f^{-1} \ \alpha \over \int_\gamma \alpha}.
$$

%%%%%%%%%%%%%%%%%% **Old, messy answer** %%%%%%%%%%%%%

Edit. Long answer and a bit "thinking aloud" or "typing while thinking" in its organization (I'll reorganize it when I'll have the time), but at the end you have the following result:

**Theorem.** *Consider a contact form $\alpha$ on the three-sphere and so that its Reeb vector field $X$ does not admit any smooth integral of motion except constants (e.g., the lift of some ergodic Finslerian or Riemannian geodesic flow on the two-sphere). If $f$ is a nowhere zero smooth function on $S^3$ such that it is not constant, and such there exist two closed Reeb orbits of $\alpha$ over which the averages of $1/f$ are distinct, then the vector field $fX$ preserves the volume form $(1/f) \alpha \wedge d\alpha$, but it is not the Reeb vector field of any contact form on $S^3$.*

Consider a contact form $\alpha$ on the three-sphere and so that its Reeb vector field $X$ does not admit any smooth integral of motion except constants (e.g., the lift of some ergodic Finslerian or Riemannian geodesic flow on the two-sphere). Let $f$ be a nowhere zero smooth function on $S^3$ that is **not** constant.

**Claim.** *The vector field $fX$ preserves the volume form $(1/f) \ \alpha \wedge d\alpha$, but it is not the Reeb vector field of any contact form DEFINING THE SAME CONTACT STRUCTURE AS $\alpha$*

The proof is by contradiction:

Assume $fX$ is the Reeb vector field of the contact form $\beta$. This means that the characteristic distributions of $d\alpha$ and $d\beta$ coincide. Since these are $2$-forms in a three-dimensional space, we conclude that these forms are multiples: $d\alpha = \lambda d\beta$. Moreover $\lambda$ must be an integral of motion of $X$,
$$
0= \mathcal{L}_X d\alpha = \mathcal{L}_X \lambda d\beta = X(\lambda) d\beta .
$$
By hypothesis, this means that $\lambda$ is a constant and, *since both forms define the same contact structure,* that $\beta = \lambda \alpha$. In this case, the Reeb vector field of $\beta$ is $fX = \lambda X$ and $f = \lambda$, contrary to the assumption that $f$ was not constant. $\square$

**Remark.** If $\alpha$ and $\beta$ do not define the same contact structure, the above proof still shows that $\alpha$ and $\beta$ differ only in a multiplicative constant and the addition of the differential of some function. I have to finish grading so I'll come back to this later if someone has not given some nice explicit example.

**Continuation.** If we drop the assumption that both forms define the same contact structure, the preceding proof give us that $\beta = \lambda \alpha + dg$, where $\lambda$ is a nonzero constant and $g$ is some smooth function on the three-sphere. The relation between $\lambda$ and $f$ is easy to establish. Indeed, if we note that
$$
1 = \beta(fX) = \lambda \alpha (fX) + dg(fX) = \lambda f + fX(g) ,
$$
multiply both sides by the volume form $(1/f) \ \alpha \wedge d\alpha$, and integrate over the three-sphere, we obtain that
$$
\int (1/f) \ \alpha \wedge d\alpha = \lambda \int \ \alpha \wedge d\alpha
$$.
Note that the "missing" integral is zero because $X$ preserves the form $ \alpha \wedge d\alpha$. *So $\lambda$ is the average of $1/f$ with respect to the volume form $ \alpha \wedge d\alpha$.*

Going back to the equation
$$
1 = \beta(fX) = \lambda f + fX(g) = f(\lambda + X(g)) \ \ {\rm or} \ 1/f = \lambda + X(g),
$$
observe that if we multiply both sides by $\alpha$ and integrate over any closed characteristic $\gamma$, we obtain
$$
\int_\gamma (1/f) \alpha = \lambda \int_\gamma \alpha
$$
Therefore $\lambda$ is also the average of $1/f$ over every single closed characteristic of $(S^3,\alpha)$. That is a bit too much to ask of $\lambda$!

Thanks to D. Cristofaro-Gardiner and M. Hutchings we know that every Reeb vector field in the sphere has alt least two closed characteristics. We can be more restrictive with our original choice of $f$ and require that the average of $1/f$ be different for two closed characteristics of $(S^3,\alpha)$. The form $\beta$ then does not exist.

We have then:

**Theorem.** *Consider a contact form $\alpha$ on the three-sphere and so that its Reeb vector field $X$ does not admit any smooth integral of motion except constants (e.g., the lift of some ergodic Finslerian or Riemannian geodesic flow on the two-sphere). If $f$ is a nowhere zero smooth function on $S^3$ such that it is not constant, and such there exist two closed Reeb orbits of $\alpha$ over which the averages of $1/f$ are distinct, then the vector field $fX$ preserves the volume form $(1/f) \alpha \wedge d\alpha$, but it is not the Reeb vector field of any contact form on $S^3$.*

1more comment