I want to calculate the divergence of the Gravitational field: $$\nabla\cdot \vec{F}=\nabla\cdot\left( -\frac{GMm}{\lvert \vec{r} \rvert^2} \hat{e}_r\right )$$ in spherical coordinates.
I know that in spherical coordinates: $$\begin{aligned} & x=r \sin\theta \cos \phi \\&y=r\sin\theta \sin \phi \\& z=r\cos\theta \end{aligned}$$
and the unit vector are:
$$\begin{aligned} & e_r=\begin{pmatrix}\sin\theta\cos\phi\\\sin\theta \sin\phi\\\cos\theta \end{pmatrix} \\ & e_{\theta}=\begin{pmatrix}\cos\theta\cos\phi\\\cos\theta \sin\phi\\-\sin\theta \end{pmatrix}\\&e_{\phi}=\begin{pmatrix}-\sin\phi\\\cos\phi\\0\end{pmatrix}\end{aligned}$$
Now I need to convert my original vector field into spherical coordinates (this is the part I am not really sure about):
$$\vec{F}=-\frac{GMm}{x^2+y^2+z^2} \hat{e}_x-\frac{GMm}{x^2+y^2+z^2} \hat{e}_y-\frac{GMm}{x^2+y^2+z^2} \hat{e}_z $$
transforming the coordinates: $x^2+y^2+z^2=(r\sin\theta\cos\phi)^2+(r\sin\theta\sin\phi)^2+(r\cos\theta)^2=r^2$
$$\implies\vec{F}=\frac{-GMm}{r^2}\left(\hat{e}_x +\hat{e}_y +\hat{e}_z \right )$$
How can I transform the unit vectors now? Do I just replace them by the spherical unit verctors?
Is there a short really cool way to calculate the divergence of this vector field? I know that the answer should be zero except at $r=0$ the divergence should be undefined.