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$.