Short answer: By the uncertainty principle, the harmonic oscillator can't be localized at the minimum value of potential energy, i.e., $x=0$, because, by the uncertainty principle, it's momentum would become large (strictly speaking, the expectation value of $p^2$, and thereby it's kinetic energy, becomes large). The lowest energy state of the harmonic oscillator is a compromise between minimizing potential energy (i.e., $x^2$) and kinetic energy (i.e., $p^2$), which cannot be done simultaneously, because $\langle x^2\rangle \langle p^2 \rangle \ge \frac{\hbar^2}{4}$.
I have to say that I don't fully understand your question, and there still seem to be inconsistent factors of $4\pi$. But I want to get at the heart of the Green function for your boundary conditions, and then you can always ask another question about how to handle it for other boundary conditions.
$$G(r,r')=\frac{1}{\sqrt{(x-x')^2+(y-y')^2+(z-a)^2}}$$ is not your Green function, and neither is
$$G(r,r')=\frac{1}{\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}}$$, but I think
$$G(r,r')=\frac{1}{\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}}-\frac{1}{\sqrt{(x-x')^2+(y-y')^2+(z+z')^2}}$$
is the function you are looking for. It tells you the potential at $r$ due to a charge at $r'$, and when $z=0$ it gives $G=0$. Now for a charge distribution $\rho$ you can integrate to get $\int G(r,r')\rho(r')dx'dy'dz'$
The point of a Green function is that if you can find the solution at $r$ due to a single unit charge at an arbitrary point $r'$ that meets your boundary conditions, and call that function $G(r,r')$ then the work you did to get $G$ now allows you to solve for any charge distribution $\rho$ by doing an integral to get $V(r)=\int G(r,r')\rho(r')dx'dy'dz'$.
One Green function $G(r,r')$ that allows you to get $V(r)$ for many different $\rho$ is the whole point.
Best Answer
Let me continue on from BeastRaban's exposition: In 1866, Mehler figured out how to carry out the sum for the $u_n(x)$ eigenfunctions of the quantum harmonic oscillator, Hermite polynomials, adding them all to an elegant and compact eponymous Mehler kernel, or equivalently.
The Green's function (propagator) for the quantum harmonic oscillator is then: $$ K(x,x';t)=\left(\frac{m\omega}{2\pi i\hbar \sin \omega t}\right)^{\frac{1}{2}}\exp\left(-\frac{m\omega((x^2+x'^2)\cos\omega t-2xx')}{2i\hbar \sin\omega t}\right) ~. $$ Consequently, given a configuration at t=0, in general a superposition of lots of components with characteristic frequency each and every one of them, you now know how the quantum oscillator hamiltonian will propagate it to arbitrary time t: $$ \psi(x,t) = \int_{-\infty}^\infty K(x,x';t) \psi(x',0) ~ dx'.$$
So the takeaway lesson of all such elaborate eigenfunction expansions is the Mehler kernel and one integral to perform, which accounts for everything.