Title says it all. By now I have tried by hand and I think that it is indeed solvable, but I can't handle the very long terms.
I tried to run the thing through SAGEs integrator:
%var x, y, z, i, n
ai = function('ai')
xi = function('xi')
yi = function('yi')
zi = function('zi')
gauss(a,x,y,z) = a * e^(-(x^2+y^2+z^2)/2)
f(x,y,z) = (sum( gauss(ai(i),x-xi(i), y-yi(i), z-zi(i)) ,i,1,n))^2
integral(integral(integral(f,x),y),z)
Which results in:
(x, y, z) |--> integrate(integrate(integrate(sum(ai(i)*e^(-1/2*x^2 - 1/2*y^2 - 1/2*z^2 + x*xi(i) - 1/2*xi(i)^2 + y*yi(i) - 1/2*yi(i)^2 + z*zi(i) - 1/2*zi(i)^2), i, 1, n)^2, x), y), z)
So I think it is telling me that it is not possible?
Funny thing is that if I remove the outer ^2 it has no problem with coming up with a result. Problem is that I need that to calculate the absolute area under the integral. Is there a better way to achieve this?
edit just to avoid confusion: I asked the same question on StackOverflow but it seemed to be a bad fit there.

Compared to the numeric solution this is the correct solution:
$\iiint\limits_{\rm I\!R^3} \left( \sum_{i=0}^n a_i e^{-\frac{(x-b_i)^2+(y-c_i)^2+(z-d_i)^2}{2}} \right)^2 = \frac{1}{8}\pi^{\frac{3}{2}} \sum_{i=0}^n \left( a_i^2 erf(x-b_i) erf(y-c_i) erf(z-d_i) \right) + \\ \frac{1}{4}\pi^\frac{3}{2} \sum_{i=0}^n \sum_{j=i}^n \left( a_i a_j erf(x - \frac{b_i+b_j}{2}) erf(y-\frac{c_i+c_j}{2}) erf(z-\frac{d_i+d_j}{2}) e^{ \left( - \frac{1}{4}(b_i-b_j)^2 - \frac{1}{4}(c_i-c_j)^2 - \frac{1}{4}(d_i-d_j)^2 \right) } \right) $
Jacks answer is not giving me the results I expect, but put me on the right track to solve this myself.