"The divergence of any electric field created by any surface charge distribution would be zero."
No it's not. Consider a charged conducting sphere with uniform surface charge density and a Gaussian sphere of radius greater than the original one. The electric field is diverging through the surface of the Gaussian sphere. So the divergence cannot be zero.
Being ρ(r) the volume charge density at that point.
No. ρ(r) is the volume charge density situated some where and we are trying to find the divergence of electric field through a surface enclosing that charge, no matter what's the size of the enclosing surface, the divergence you get will be the same, provided you have at least a volume to accommodate all the charges. This is what Gauss's law tells us.
What if there is zero volume charge density but non-zero surface charge density?
If you have only a surface charge density or linear charge density, use Gauss's law in integral form. Don't stick with the differential form as the differential form of Gauss's law is applicable to volume charges only. As Gauss's divergence theorem states,
$∫E.dS$ gives you the exact measure of divergence by measuring the flux across the surface.
A point charge is confined to a single point in space. Let's call it as $q\left(\vec{r}\right)$. This means that the charge has a magnitude only at $\vec{r}$, and at all other points it is zero.
The charge density for a point charge is given by charge per unit volume. Since the there is no charge except at the position $\vec{r}$, the charge density vanishes at all points except $\vec{r}$. Now, at $\vec{r}$, the volume gives off to the limit $V\rightarrow0$. So the charge density blows up to infinity at that point. It's the case for all point sources, not just restricted to point charges: for example, the mass density of a point mass blows up at the origin.
Does this work for all point charges? Yes. But, there is an exception when the magnitude of the point charge is zero. However if that's the case, why should one be concerned about all this? So, a point charge has a non-zero magnitude at the point it occupies. So, it works all time.
This particular property of the charge density of a point charge is exactly identical to the definition of the Dirac-delta function, which, for the point $\vec{r}$ can be defined as
$$
\delta^3\left(\vec{r}\right) =
\begin{cases}
\infty, & \text{at the point $\vec{r}$} \\[2ex]
0, & \text{at all other points}
\end{cases}
$$
So, it seems quite comfortable that we could use this function to represent the volume charge density of a point charge. It is to be noted that at all other points except at $\vec{r}$, the charge is zero (which contributes to the second case of the Dirac-delta function defined above, and at the point $\vec{r}$, the first case of Dirac- delta function is satisfied by the inverse volume.
However, your question is concerned by the definition of the charge density of a point charge at the "origin" (here $\vec{r}$). So, at the origin, the charge is non-zero, but the inverse volume blows up. So, we substitute the Dirac-delta function in the place of inverse volume as
$$\rho=q\delta^3\left(\vec{r}\right)$$
This definition is however valid at other points than origin also, since at all other points the delta function vanishes and so does the charge density. But for a point charge, the result is trivial that $q=0$ at all other points.
Best Answer
First equation is wrong, it should say $\rho(\vec{r}) = \delta(\vec{r} - \vec{r}')q$. (Note that you had two errors).
You treat it like a normal charge density $\rho(\vec{r})$, if you integrate the density over any volume you get the total charge within that volume.