With the Liouville's theorem $$\frac{{d\rho }}{{dt}} = \frac{{\partial \rho }}{{\partial t}} + \sum\limits_{a = 1}^{3N} {(\frac{{\partial \rho }}{{\partial {p_a}}}\frac{{d{p_a}}}{{dt}} + \frac{{\partial \rho }}{{\partial {q_a}}}\frac{{d{q_a}}}{{dt}})} = 0$$ when we calculate the time evolution of the ensemble average of a quantity $O(p,q)$ we have $$\frac{{d\left\langle O \right\rangle }}{{dt}} = \int {d\Gamma \frac{{\partial \rho (p,q,t)}}{{\partial t}}O(p,q)} = \sum\limits_{a = 1}^{3N} {\int {d\Gamma } } O(p,q)(\frac{{\partial \rho }}{{\partial {p_a}}}\frac{{\partial H}}{{\partial {q_a}}} - \frac{{\partial \rho }}{{\partial {q_a}}}\frac{{\partial H}}{{\partial {p_a}}})$$ here $p,q$ represents a bunch of generalized coordinates and momentum ${p_a},{q_a},a = 1,...,3N$. Then by using the method of integration by parts, the above integration becomes $$\frac{{d\left\langle O \right\rangle }}{{dt}} = - \sum\limits_{a = 1}^{3N} {\int {d\Gamma } } \rho [(\frac{{\partial O}}{{\partial {p_a}}}\frac{{\partial H}}{{\partial {q_a}}} - \frac{{\partial O}}{{\partial {q_a}}}\frac{{\partial H}}{{\partial {p_a}}}) + O(\frac{{{\partial ^2}H}}{{\partial {p_a}\partial {q_a}}} - \frac{{{\partial ^2}H}}{{\partial {q_a}\partial {p_a}}})$$
Here comes my questions, I think the integration by parts uses $\int {d{p_a}} \frac{{\partial \rho }}{{\partial {p_a}}} = \int {d\rho }$. However as $\rho$ ~$\rho (p,q,t)$, should we have $\frac{{d\rho }}{{d{p_a}}} = \frac{{\partial \rho }}{{\partial {p_a}}} + \frac{{\partial \rho }}{{\partial t}}\frac{{dt}}{{d{p_a}}}$.
I take that last relation for granted because when we calculate $\frac{{dy}}{{dx}}$, if $y=y(x)$ determines implicitly by some relation $F(x,y)=0$, we use $$\frac{{\partial F(x,y)}}{{\partial x}} + \frac{{\partial F(x,y)}}{{\partial y}}\frac{{dy}}{{dx}} = 0$$ and get $$\frac{{dy}}{{dx}} = - \frac{{\frac{{\partial F(x,y)}}{{\partial x}}}}{{\frac{{\partial F(x,y)}}{{\partial y}}}}$$ If we differentiate $F(x,y)=0$ with y, in order to get the same value of $\frac{{dy}}{{dx}}$, we need to have $$\frac{{\partial F(x,y)}}{{\partial x}}\frac{{dx}}{{dy}} + \frac{{\partial F(x,y)}}{{\partial y}} = 0$$ where we consider $x$~$x(y)$. So back to the $\rho$ case, in calculating $\frac{{d\rho }}{{d{p_a}}}$, I think $\frac{{\partial \rho }}{{\partial t}}\frac{{dt}}{{d{p_a}}}$ term should be taken into account since $p_a$~$p_a(t)$. Then the integration by parts seems to be wrong, so I think I have made a mistake.
My second question concerns Liouville's theorem itself, if we have $\frac{{d\rho }}{{dt}}=0$, then can we have $\frac{{d\rho }}{{d{q_a}}} = \frac{{d\rho }}{{dt}}\frac{{d{q_a}}}{{dt}} = 0$? It sounds ridiculous as it indicates that the probability density is everywhere the same in the phase space, regardless of what the system is. Then where did I make a mistake? Thanks for your patience reading my long questions!