Calculate surface integral of point charge located outside the surface Problem
Given:
$$\vec r = r(\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta)$$
$S: | \vec r  | = a$, with $\hat n$ outwards
$$\vec r_0 = \frac{3a}{5}(\hat x + \hat y - \hat z) = \frac{3a}{5}(1,1,-1)$$
$$\vec F = k \frac{\vec r - \vec r_0}{ {| \vec r - \vec r_0 |}^3 }$$
Calculate:
$$\int_S \vec F \cdot \mathrm d \vec S $$
Solution (added after accepted answer)
Thanks to the accepted answer which confirmed that the explicit calculation on/in $S$ is messy, and the suggestion of using translation invariance, I have noted down the solution below (let me know if you have further suggestions).
An argument is used which concludes that the integrand is zero inside the sphere $S$. 
Using Gauss' theorem: $\int_S \vec F \cdot \mathrm d \vec S = \int_V \nabla \cdot \vec F \mathrm dV $. 
Since the calculation is not easy to explicitly calculate in (or on) $S$, instead it is shown that $\nabla \cdot \vec F = 0$ exterior to $S_{\epsilon}$, which is a sphere containing the singularity, and since $S$ is exterior to $S_{\epsilon}$ the integrand must be $0$ and so $\int_V \nabla \cdot \vec F \mathrm dV = 0$.
$$S_{\epsilon}: |\vec r - \vec r_0 | = \epsilon, \epsilon > 0$$
translate the coordinate system so that $S_{\epsilon}$ is the origin in the translated system:
$$\vec R = \vec r - \vec r_0 = R (\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta)$$
$$F(\vec R) = k \frac{\vec R}{|\vec R|^3} = k \frac{1}{R^2} \hat R$$
$$S_{\epsilon}: |\vec R| = \epsilon$$
$$\nabla \cdot \vec F (\vec R) = \frac{1}{R^2} \frac{\partial}{\partial R} \left(R^2 \vec F_R \right) = \frac{1}{R^2} \frac{\partial}{\partial R} \left(R^2 k \frac{1}{R^2} \right) = 0$$, ($k$ is a constant)
This is true for any value of $R$ except at the singularity ($\epsilon$ can be made as small as required), since $S$ does not contain the singularity, the integrand ($\nabla \cdot \vec F(\vec r)$ is 0 and the integral is therefore zero.
My question (updated)
I am having difficulties in explicitly calculating the value of this integral. 
Specifically, the divergence of the field becomes messy, I am not able to see how I can use the symmetry of $S$ due to $\vec r_0$. (When applying the Gauss' theorem, I am stuck in evaluating the divergence of the field). 
I am able to argue for that this integral is indeed zero (the field only has one singularity and it is exterior to the sphere $S$). However, I am not able to explicitly show (by calculation) that this integral is zero. Thankful for any help.
Please note that this is not homework, I am studying for an exam.
My question (original)
What is an easy way to calculate this integral? Any suggestions on the approaches below?
My apologies if this is due to lack of some basic knowledge (I am back studying after 2.5 years)
Intuitively I understand the integral is zero (the point charge is located outside the sphere, anything flowing into the sphere will also flow out), however, I have issues with the calculation.
Calculate directly
$$\int_S \vec F \cdot \mathrm d \vec S = \int_0^{\pi} \mathrm d \theta \int_0^{2\pi} \mathrm d \phi r^2 \sin {\theta} \vec F_r \cdot \hat r = \int_0^{\pi} \mathrm d \theta \int_0^{2\pi} \mathrm d \phi r^2 \sin \theta k \frac{a - \frac{3a}{5} \left(  \sin \theta \cos \phi + \sin \theta \sin \phi - \cos \theta) \right)}{ {\left| \vec r - \vec r_0 \right|}^3 }$$
Using Gauss theorem
I am tempted to use Gauss theorem, $\int_S \vec F \cdot \mathrm d \vec S = \int_V \nabla \cdot \vec F \mathrm dV $, in spherical coordinates. However $\vec F$ has components also in $\hat \theta$ and $\hat \phi$ (due to $\vec r_0$) and ${\left| \vec r - \vec r_0 \right|}^3$ is not that nice to derivate.
$$\vec F_r = k \frac{\vec r - \vec r_0}{ {\left| \vec r - \vec r_0 \right|}^3 } \cdot \hat r = k \frac{r - \frac{3a}{5} \left(  \sin \theta \cos \phi + \sin \theta \sin \phi - \cos \theta) \right)}{ {\left| \vec r - \vec r_0 \right|}^3 }$$
$$\vec F_{\theta} = k \frac{\vec r - \vec r_0}{ {\left| \vec r - \vec r_0 \right|}^3 } \cdot \hat \theta = k \frac{ \frac{3a}{5} \left( \cos \theta \cos \phi + \cos \theta \sin \phi - \sin \theta) \right)}{ {\left| \vec r - \vec r_0 \right|}^3 }$$
$$\vec F_{\phi} = k \frac{\vec r - \vec r_0}{ {\left| \vec r - \vec r_0 \right|}^3 } \cdot \hat \phi = k \frac{ \frac{3a}{5} \left( - \sin \phi + \cos \phi ) \right)}{ {\left| \vec r - \vec r_0 \right|}^3 }$$
$$\nabla \cdot \vec F = \frac{1}{r^2} \frac{\partial}{\partial r} \left(r^2 \vec F_r \right) + \frac{1}{r \sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta F_{\theta} \right) + \frac{1}{r \sin \theta} \frac{\partial F_{\phi}}{\partial \phi}$$
 A: (I don't know what you mean by $\hat x+\hat y -\hat z$, so I disregard the information about $r_0$.) Your field $\vec F$ is the gravitational field produced by a point mass at $r_0$. The divergence of this field is $\equiv0$ away from $r_0$. If the point $r_0$ is in the exterior of the sphere $S$ then it follows by Gauss' theorem that the integral in question is $0$. If $r_0$ lies in the interior of $S$ then you should draw a tiny sphere $S_\epsilon$ with center $r_0$ and apply Gauss' theorem to the region between $S_\epsilon$ and $S$. You are left with an integral over $S_\epsilon$ which you can almost do in your head.
A: EDIT: This reply assumes that the surface contains the singularity ("source") at its center, so it does not actually pertain to this problem.  For the sake of clarity, I have now changed the surface from $S$ to $S'$.
This is a lot like an electrostatics problem, and it is probably trying to mimic one due the field's specific form.  I am going to treat this problem like it is so; someone please point out if what follows is incorrect since I am not absolutely sure.
The field has rotational symmetry, and the $\vec{r_{0}}$ usually means the location of a "source" of the field.  Thus, no matter the location of the source, this field is still rotationally invariant.  Using that surface for the integral is called a "Gaussian surface" by physicists; it is a tool used to help solve a problem by exploiting its symmetry.
Now, I know you probably just want to solve the integral, so I will get to the point.  Your surface will be a sphere of radius $a$, and so will have a unit normal area vector $\hat{r}'$.  Let us write $\vec{r}'=\vec{r}-\vec{r_{0}}$ and $r=|\vec{r}|$.  If we rewrite the integrand, we get
$$k\frac{\vec{r}-\vec{r_{0}}}{|\vec{r}-\vec{r_{0}}|^{3}}=k\frac{\vec{r}'}{|\vec{r}'|^{3}}=k\frac{\hat{r}'}{|\vec{r}'|^{2}}~\text{since, for any vector }\vec{v},~\hat{v}=\frac{\vec{v}}{|\vec{v}|}.$$
The infinitesimal vector area element for the sphere will be
$$d\vec{S'}=\hat{r}'r'^{2}\sin\phi d\theta d\phi.$$
Thus, the integral becomes
$$\begin{aligned}
\int\vec{F}\cdot d\vec{S'}&=\int k\left(\frac{\hat{r}'}{|\vec{r}'|^{2}}\right)\cdot\left(\hat{r}'r'^{2}\sin\phi d\theta d\phi\right) \\
&=\int k\frac{r'^{2}}{r'^{2}}\hat{r}'^{2}\sin\phi d\theta d\phi \\
&=k\int_{0}^{2\pi}d\theta\int_{0}^{\pi}\sin\phi d\phi \\
&=4\pi k.
\end{aligned}$$
Notice that the radius of the surface is irrelevant.  As long as it contains the "source" point in the center, the value of the integral will always be the same for this field.
