Bill Foster said:
\nabla=\hat{\textbf{r}}\frac{\partial}{\partial{r}}+\hat{\theta}\frac{1}{r}\frac{\partial}{\partial{\theta}}+\hat{\phi}\frac{1}{r\sin{\theta}}\frac{\partial}{\partial{\phi}}
\textbf{r}=r\hat{\textbf{r}}+\theta\hat{\theta}+\phi\hat{\phi}
\rho\left(\textbf{x}\right)=\rho\left(r,\theta,\phi\right)
Now taking into consideration just the middle term (the one that should cancel out), and ignoring the constants:
\int_0^{2\pi}\int_0^\pi \left(\textbf{r}\cdot\nabla\right)\rho\left(\textbf{x}\right)r^2\sin\theta d\theta d\phi
There are several problems here.
First, it is very important to distinguish \left(\textbf{r}\cdot\mathbf{\nabla}_{\overline{\textbf{x}}}\right)\rho(\overline{\textbf{x}})|_{\overline{\textbf{x}}=\textbf{x}} and \left(\textbf{r}\cdot\mathbf{\nabla}_{\overline{\textbf{x}}}\right)^2\rho(\overline{\textbf{x}})|_{\overline{\textbf{x}}=\textbf{x}}from \left(\textbf{r}\cdot\nabla\right)\rho\left(\textbf{x}\right) and \left(\textbf{r}\cdot\nabla\right)^2\rho\left(\textbf{x}\right). In the first two, differentiation is done with respect to the dummy variable \overline{\textbf{x}}, and so \textbf{r} is treated as a constant vector. In the latter two, \textbf{r} is not treated as a constant vector, and so your results will be quite different.
Second, by defintion, \textbf{r}=r\mathbf{\hat{r}}\neq r\mathbf{\hat{r}}+\theta\mathbf{\hat{\theta}}+\phi\mathbf{\hat{\phi}}...the units on your expression don't even make sense... r has unit of distance, while \theta and \phi are angles, with units of radians.
Third, unit vectors in curvilinear coordinates are position dependent. So, for example \textbf{r}\cdot\textbf{r}'=\left(r\mathbf{\hat{r}}\right)\cdot \left(r'\mathbf{\hat{r}'}\right)=r r' (\mathbf{\hat{r}}\cdot\mathbf{\hat{r}}') \neq rr' in general.
Finally, you only require that
\int_{0}^{\pi} \int_{0}^{2\pi}\left[-\left(\textbf{r}\cdot\mathbf{\nabla}_{\overline{\textbf{x}}}\right)\rho(\overline{\textbf{x}})|_{\overline{\textbf{x}}=\textbf{x}}+\frac{1}{2}\left(\textbf{r}\cdot\mathbf{\nabla}_{\overline{\textbf{x}}}\right)^2 \rho(\overline{\textbf{x}})|_{\overline{\textbf{x}}=\textbf{x}}\right]\sin\theta d\theta d\phi=\frac{2\pi r^2}{3}\nabla_{\overline{\textbf{x}}}^2\rho(\overline{\textbf{x}})|_{\overline{\textbf{x}}=\textbf{x}}
This does not necessarily mean that
\int_{0}^{\pi} \int_{0}^{2\pi}\left(\textbf{r}\cdot\mathbf{\nabla}_{\overline{\textbf{x}}}\right)\rho(\overline{\textbf{x}})|_{\overline{\textbf{x}}=\textbf{x}}\sin\theta d\theta d\phi=0
and
\int_{0}^{\pi} \int_{0}^{2\pi}\frac{1}{2}\left(\textbf{r}\cdot\mathbf{\nabla}_{\overline{\textbf{x}}}\right)^2 \rho(\overline{\textbf{x}})|_{\overline{\textbf{x}}=\textbf{x}}\right]\sin\theta d\theta d\phi=\frac{2\pi r^2}{3}\nabla_{\overline{\textbf{x}}}^2\rho(\overline{\textbf{x}})|_{\overline{\textbf{x}}=\textbf{x}}