Does Hamilton Mechanics give a general phase-space conserving flux? Hamiltonian dynamics fulfil the Liouville's theorem, which means that one can imagine the flux of a phase space volume under a Hamiltonian theory like the flux of an ideal fluid, which doesn't change it's volume. But is it possible to reproduce every phase space conserving flux with an appropriate Hamiltonian?
So can I simply imagine the entity of all possible Hamiltonian dynamics as all possible phase space conserving fluxes? Or are Hamiltonian dynamics a special case for phase space conserving fluxes? If they are a special case, what would be an example for a phase space conserving flux for which there is no Hamiltonian that can produce it?
 A: First, let's take a look at one-dimensional systems with phase space dimension $2$.
The volume form is just the symplectic one, ie any volume-preserving flow is symplectic and thus at least locally Hamiltonian (but not necessarily globally so).
Now, consider an arbitrary phase space of dimension $2n\geq4$ with canonical coordinates $q^i,p^i$.
Up to a constant factor, the volume form is
$$
\Omega = dq^1\wedge\cdots\wedge dq^n\wedge dp^1\wedge\cdots\wedge dp^n
$$
and the symplectic form
$$
\omega = \sum_i dp^i\wedge dq^i
$$
Let's take a look at the vector field $X$ given by
$$
X=q^1\frac{\partial}{\partial q^2}
$$
As
$$
\mathcal L_X\Omega = 0
$$
phase space volume will be preserved, but as
$$
\mathcal L_X\omega \not= 0
$$
the vector field is not symplectic and thus also not Hamiltonian.
A: I) More generally, OP is essentially pondering: 

Let $X\in \Gamma(TM)$ be a given vector field on a $2n$-dimensional manifold $M$. Under what conditions is the evolution equation 
  $$\tag{1} \frac{df}{dt} ~=~ X[f]+\frac{\partial f}{\partial t}  $$
  a Hamiltonian system? In other words, under what conditions is $X$ a Hamiltonian vector field?

Globally, there can be topological obstructions, so let us from now on in this answer only consider local conditions. Locally, the Hamiltonian vector fields $X$ (wrt. to a symplectic two-form $\omega$) are precisely the vector fields $X$ that preserve
$$\tag{2} {\cal L}_X \omega ~=~0$$ 
the symplectic two-form $\omega$. 
II) Discarding global issues, OP's actual question might be interpreted as follows:

Is a divergence-free vector field $X\in \Gamma(TM)$ on a given $2n$-dimensional symplectic manifold $(M,\omega)$ a local Hamiltonian vector field? Here the volume form is 
  $$\tag{3} \Omega~=~\frac{1}{n!}\omega^{\wedge n},$$ 
  and a divergence-free vector field satisfies 
  $$\tag{4} {\cal L}_X \Omega ~=~0.$$

Answer: No. It is indeed always true in two dimensions. However, it is not necessarily true in dimensions $\geq 4$. It is easy to produce counterexamples, cf. e.g. the answer by Christoph. 
III) Alternatively, OP's actual question might be interpreted as follows:

Let there be given a $2n$-dimensional manifold $(M,\Omega)$ with a volume form $\Omega$ and a non-vanishing$^1$ vector field $X\in \Gamma(TM)$ that is divergence-free (4). Does there locally exist a symplectic two-form $\omega$ such that eqs. (2) and (3) are satisfied?

Answer: Yes! 
Sketched proof: Given a point $p\in M$ with $X_p\neq 0$. One may show that there exists a local chart $U\subseteq M$ with coordinates $z=(z^1, \ldots z^{2n})$ such that $X$ is locally of the form 
$$\tag{5} X~=~\frac{\partial }{\partial z^2}. $$
We may assume that the fixed point $p\in M$ corresponds to $z=0$. The volume form will locally be of the form  
$$\tag{6}\Omega~=~\rho(z)~ \mathrm{d}z^1\wedge \ldots \wedge \mathrm{d}z^{2n}, \qquad \rho(z)~\neq~ 0. $$
Eqs. (4) and (5) imply that $\rho(z)$ does not depend on $z^2$,
$$\tag{7} \frac{\partial \rho(z)}{\partial z^2}~=~0. $$
Locally there exists a function $f=f(z)$ also independent of $z^2$, such that  
$$\tag{8} \rho(z)~=~\frac{\partial f(z)}{\partial z^1}, \qquad 
\frac{\partial f(z)}{\partial z^2}~=~0. $$
Now change coordinates 
$$\tag{9} w^1 ~:=~ f(z), \quad w^2 ~:=~z^2,\quad w^3 ~:=~z^3, \quad \ldots,\quad w^{2n} ~:=~z^{2n}.$$
The Jacobian will be 
$$\tag{10} J~:=~\det\frac{\partial w}{\partial z}~=~\frac{\partial f}{\partial z^1}~=~\rho~\neq 0,$$
so in the new coordinates $w^I$, the volume form (6) is 
$$\tag{11}\Omega~=~ \mathrm{d}w^1\wedge \ldots \wedge \mathrm{d}w^{2n},$$
while the vector field (5) becomes
$$\tag{12} X~=~\frac{\partial }{\partial z^2}
~=~\sum_{I=1}^{2n}\frac{\partial  w^I}{\partial z^2}
\frac{\partial }{\partial w^I}
~=~\frac{\partial f(z)}{\partial z^2}\frac{\partial }{\partial w^1}
+\frac{\partial }{\partial w^2}
~\stackrel{(8)}{=}~\frac{\partial }{\partial w^2} . $$
Next we rename the new coordinates
$$\tag{13} p_1 ~:=~ w^1, \quad q^1 ~:=~w^2,\quad p_2 ~:=~w^3, \quad q^2 ~:=~w^4, \quad \ldots,\quad q^{n} ~:=~w^{2n},$$
and define a symplectic two-form
$$ \tag{14} \omega ~:=~\sum_{i=1}^n\mathrm{d}p_i\wedge  \mathrm{d}q^{i}~=~\sum_{i=1}^n\mathrm{d}w^{2i-1}\wedge  \mathrm{d}w^{2i}. $$
The vector field (12) is locally Hamiltonian
$$\tag{15} X~=~\frac{\partial }{\partial w^2}~=~\frac{\partial }{\partial q^1}~=~ -\{p_1, \cdot\}_{PB}. $$
It is straightforward to see that eqs. (2) and (3) are satisfied. End of proof.
IV) The proof generalizes to the following question:

Let there be given a $2n$-dimensional manifold $M$ with a non-vanishing$^1$ vector field $X\in \Gamma(TM)$. Does there locally exist a symplectic two-form $\omega$, such that $X$ is a Hamiltonian vector field?

Answer: Yes! 
V) The two-dimensional problem is also considered in this Phys.SE post.
--
$^1$ The vector field $X$ is for technical reasons assumed to be non-vanishing, which is the generic case. We leave it to the reader to ponder about the special case where the vector field $X$ vanishes $X_p=0$ in a point $p\in M$.  
