When does one consider the laplacian as a dirac delta function? If we consider that 
$$\nabla^2\left(\frac{1}{r}\right) = -4\pi\delta(\vec{r})$$
we can explain the dirac-delta function here via the results of Gauss' law
$$\int_V \nabla^2\left(\frac{1}{r}\right) = \int_S\nabla\left(\frac{1}{r}\right)\cdot d\vec{a} = \int_S \frac{-\hat{r}}{r^2}\cdot d\vec{a} = -4\pi\frac{R^2}{r^2} $$
where the surface integral has been taken over a sphere of radius $R$. Now, if we say $r>0$ and $R\to 0$, we get 
$$\int_V \nabla^2\left(\frac{1}{r}\right) = 0. $$
Furthermore, if we set $r=R$ and send $R\to0$ we get 
$$\int_V \nabla^2\left(\frac{1}{r}\right) = -4\pi $$
Thus we have demonstrated the validity of the first equation. However, if we were to directly calculate the Laplacian of $\frac 1r$ we obtain
$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial}{\partial r}\left(\frac{1}{r}\right)\right) = 0$$ 
So it seems like the -$4\pi\delta(\vec{r})$ has somehow been neglected by the direct calculation. That being said, what about the case where we have, for some $\epsilon\in\mathbb{R} > 0$
$$\nabla^2\left(\frac{1}{r^\epsilon}\right) = \frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial}{\partial r}\left(\frac{1}{r^\epsilon}\right)\right) = \frac{(\epsilon-1)\epsilon}{r^{\epsilon+2}}$$
Has there also been a dirac delta function neglected here? In general, when does one need to consider contributions due to the dirac-delta function? 
 A: What $\nabla^2\left(\frac{1}{r}\right) = -4\pi\delta(\vec{r})$ means, by the definition of $\delta$, is that $\int\frac{1}{r}\,\nabla^2\phi\,dV=-4\pi\phi(0)$ for any smooth function $\phi$ with compact support. (It is easily shown using Gauss theorem and the identity $f\nabla^2g-g\nabla^2 f=\nabla\cdot(f\nabla g-g\nabla f)$; the verification by Gauss's law that you mentioned is not quite sufficient - think about $\nabla^2(x/r^3)$).
As for $\nabla^2r^{-\epsilon}$: the function $r^{-\epsilon}$ is locally integrable for $\epsilon<3$, so only in this case it is a well-defined distribution (and thus $r^{-2-\epsilon}$ is a well-defined distribution for $\epsilon<1$). It can be extended to a distribution also for $\epsilon\geq3$, but the extension is not unique - it's up to $\delta$ and its derivatives. If we require the distribution to be even and homogeneous, it will be unique if $\epsilon$ is not an integer. If we require it to be invariant w.r.t. all rotations, it will be unique up to terms $(\nabla^2)^k\delta$ (which is homogeneous of degree $-3-2k$). So while $r^{-\epsilon}$ is uniquely specified by these two conditions for most $\epsilon$'s, it fails for $\epsilon$ an odd integer $\geq 3$.
To summarize - one needs to be a little bit careful when working with distributions.
A: While for $r \ne 0$, $\nabla^2\left(\frac1{4\pi r}\right)=0$, it is undefined at $r=0$.  In introducing the idea of the Dirac Delta, we need to provide a rigorous interpretation.
In THIS ANSWER, I showed how to interpret the expression 
$$\nabla^2\left(\frac{1}{ r}\right)=-4\pi \delta(0)$$ 
in terms of a regularization of the Dirac Delta.  In particular, I showed that for any test function $\phi(\vec r)$, the regularization $\vec \psi(\vec r;a)=\frac{\vec r}{(r^2+a^2)^{3/2}}$, of the Dirac Delta satisfies
$$\lim_{a\to 0}\int_V \nabla \cdot \vec \psi(\vec r;a) \phi(\vec r)dV=\lim_{a\to 0}\int_V \nabla^2 \left(\frac{-1}{\sqrt{r^2+a^2}}\right) \phi(\vec r)dV=
\begin{cases}4\pi \phi(0)&\vec r\in V\\\\
0&,\vec r \,\,\text{elsewhere}
\end{cases}$$
Note that 
$$\nabla \cdot \vec \psi(\vec r;a)=\nabla^2\left(\frac{-1}{\sqrt{r^2+a^2}}\right)=\frac{3a^2}{(r^2+a^2)^{5/2}}$$
is perfectly defined for all $r$ and $a\ne 0$.  And for $r\ne 0$, we have
$$\begin{align}
\lim_{a\to 0}\vec \psi(\vec r;a)&=\lim_{a\to 0}\nabla\left(\frac{-1}{\sqrt{r^2+a^2}}\right)\\\\
&=-\nabla\left(\frac1r\right)
\end{align}$$
Thus we can interpret the Dirac Delta as 
$$\delta(\vec r)=-\nabla^2\left(\frac1{4\pi r}\right)$$
in the sense that
$$\begin{align}
\lim_{a\to 0}\int_V \nabla \cdot \vec \psi(\vec r;a)\,\phi(\vec r)\,dV&=\lim_{a\to 0}\int_V \nabla^2\left(\frac{-1}{\sqrt{r^2+a^2}}\right)\,\phi(\vec r)
\,dV\\\\
&=4\pi \phi(0)
\end{align}$$
when the origin is contained in $V$.
