Note that due to the time dependency of $\psi(\vec r, t)$ the retarted potential has to be considered. So we're looking for
$$\begin{align*}
\langle\phi(\vec r, t)\rangle &= \int_{\mathbb R^3}d^3x \frac q{4\pi\epsilon_0}\frac{\Big|\psi(\vec x, t - \frac1c|\vec x - \vec r|))\Big|^2}{|\vec x - \vec r|} \quad\Bigg|\quad \vec x \to \vec x + \vec r
\\ &= \int_{\mathbb R^3}d^3x \frac q{4\pi\epsilon_0}\frac{\Big|\psi(\vec x + \vec r, \overbrace{t - \frac xc}^{=:t_c(x)}))\Big|^2}{x}
\\ &\stackrel{(B)}= \int_{\mathbb R^3} d^3x \underbrace{
\frac q{4\pi\epsilon_0}\overbrace{\Bigg|{a \over a + i\hbar t_c/m}}^{=\frac{a}{\sqrt{a^2+\hbar^2t_c^2/m^2}}=\frac{\sqrt{2a}}{\sigma}}\Bigg|^3}_{=:N}
\frac{\exp\overbrace{\left(- \Re \frac{(\vec x + \vec r)^2}{(a + i\hbar t_c/m)} \right)}^{=-\frac{a(\vec x + \vec r)^2}{a^2+\hbar^2t_c^2/m^2} =: -\frac{(\vec x + \vec r)^2}{2\sigma^2(t_c)}}}{x}
\end{align*}$$
The $x$-dependency of $t_c$ makes this a very nasty integral, so let's take the non-relativistic limit $c\to\infty$ such that $t_c(x) \to t$ (or alternatively, assume $m\to\infty$). Then we're looking for
$$\begin{align*}
\langle\phi(\vec r,t)\rangle &= N\int_{\mathbb R^3}d^3x \frac{\exp\left(-\frac{(\vec x + \vec r)^2}{2\sigma^2}\right)}{x}. \tag{tk.C}\label{tk.C}
\end{align*}$$
Note how for $\vec r=0$ we can use spherical coordinates to obtain
$$\begin{align*}
\langle\phi(0,t)\rangle &= 4\pi N\int_0^\infty \rho^2 d\rho \frac{\exp\left(-\frac{\rho^2}{2\sigma^2}\right)}{\rho}
\\ &= -4\pi\sigma^2 N\int_0^\infty d\rho\, \partial_\rho \exp\left(-\frac{\rho^2}{2\sigma^2}\right)
\\ &= 4\pi\sigma^2 N = \frac{q\sqrt{2a}^3}{\epsilon_0\sigma}
\end{align*}$$
which for $\sigma\neq0$ (that would yield a truly localized particle which would however diffuse infinitely the very next instant) is finite, thus the answer to the question is, considering uncertainty, there is no longer a singularity.
On paper, I actually calculated $\eqref{tk.C}$ to end up with a correction factor $\mathrm{erf}\left(\frac r{\sqrt{2\sigma^2}}\right)$ to the classical Coulomb potential (by using the Fourier transform w.r.t. $\vec r$, swapping integrals and integrating over $\sigma^2$), but that's too tedious to typeset for now, and there's a sign error in it somewhere...
