The problem I'm having is how to find the gradient of delta.
Does the gradient do the following?
$$\nabla \delta(r)=\frac{\partial\delta(r)}{\partial x}\hat{x}+\frac{\partial\delta(r)}{\partial y}\hat{y}+\frac{\partial\delta(r)}{\partial z}\hat{z}$$
If that is correct, how do I derive delta of $r$ in respect to $x, y, z$?
I mean, that answer may be unsatisfying. The $\delta$-function is defined so that for a smooth functions $f$, we have that: $$\int f(x)\delta(x)dx=f(0)$$ So if you want to understand the derivative of the $\delta$-function, one can note that by integration by parts, and assuming that $f$ vanishes at infinity, you should have: $$\int f(x)\delta'(x)dx=-\int f'(x)\delta(x)dx=-f'(0)$$ In general, you just define $\delta'(x)$ to be the function with this property on all smooth functions.
Edit: I'll compute $\nabla \delta$. Note that we have that: $$\nabla \delta=\delta'(r)\hat{r}$$ This is a badly-behaved distribution in the following sense. Note that if $f(0,0,0)\neq 0$ then we formally have that: $$\int\delta'(r) f(r,\theta,\phi)drd\theta d\phi= \delta(0)f(0,0,0)-\int \int\delta(r) \partial_r f(r,\theta,\phi)drd\theta d\phi$$ The latter is well-defined by the prior is not. On the other hand, we may write: $$\int \delta'(r) f(r,\theta,\phi)drd\theta d\phi=-\int\delta(r) f(r,\theta,\phi)drd\theta d\phi-\int\delta(r) r\partial_r f(r,\theta,\phi)drd\theta d\phi=-f(0,0,0)$$
So you achieve the relation: $$r\delta'(r)=-\delta(r)$$
Edit2: To answer a point in the comments, there are two distributions with similar notation. One is defined by: $$\int\int \int \delta(r)f(r,\theta,\phi)drd\theta d\phi=f(0,0,0)$$ The other is defined by: $$\int\int \int \delta({\bf r})f(x,y,z)dxdy dz=f(0,0,0)$$
One has the formal relation: $$4\pi r^2\delta({\bf r})=\delta(r)$$ I believe the the asker was interested in the prior, rather than the former.