An efficient way to evaluate this expression with deeply nested Do loops

Question

Is there an efficient way to evaluate something like this?

Np=10;
xx = 0.;
Do[
 Do[
  Do[
   Do[
    Do[
     Do[
      If[(1. ix - 1. jx)^2 + (1. iy - 1. jy)^2 + (1. iz - 1. jz)^2 == 
        0., xx += 0., 
       xx += (1. iz - 1. jz)^2/(Sqrt[(1. ix - 1. jx)^2 + (1. iy - 
             1. jy)^2 + (1. iz - 1. jz)^2])^3
      , {jz, Np}],
     {iz, Np}],
    {jy, Np}],
   {iy, Np}],
  {jx, Np}],
 {ix, Np}]

The problem is not the complexity, but the fact that there are six independent variables and the test in that If. You can imagine what happens if Np is increased. The time it take to evaluate this increases quite a bit.

I tried to evaluate the expression that can become zero separately and put the results in a six way nested Table. The plan was to then Flatten the Table and take out the values equal to 0 and then perform the summation with something like Total on the inverse values. But this way I bump into a memory problem, because if I increase Np then all that huge Table has to be written to memory.

So the question is: Is there a way to efficiently evaluate the sum in my example, for Np=100?

Edit @ciao

It appears that I oversimplified the question too much. i thought that only the denominator had issues.

So the terms look something like this in reality:

$\frac{(iz-jz)^2}{((ix-jx)^2+(iy-jy)^2+(iz-jz)^2)^{3/2}}$, where ix, jx, iy, jy, iz, jz go from 1 to N. I didn't think that the nominator should be included in the original question. How does this modify your answer?

Answer 1

UPDATE: After reposting the below function after modifications from the deleted post to make it work after OP question was changed, I determined a way to utilize ideas from the original deleted post's second, faster function. This is much faster than the second function below for larger arguments (~100 X faster for Np argument of 500) and has less memory pressure. It will comfortably handle Np arguments into four digit range.

sumx[np_] := Module[{t, c2, u, r, v, p, x, c},
 t = Tally[(Subtract @@@ Tuples[Range@np, {2}])^2];
 p = Tr[x^t[[All, 1]]*t[[All, 2]]];
 c = CoefficientList[p, x];
 c2 = ListConvolve[c, c, {1, -1}, 0];
 r = Pick[Range[0, Length[c2] - 1], u = Unitize[c2], 1];
 v = Pick[c2, u, 1];
 Tr[Rest[Join @@ Outer[Times, v, t[[All, 2]]]] Rest[Join @@ ConstantArray[t[[All, 1]], Length@r]]/
    Rest[Join @@ Outer[Plus, r, t[[All, 1]]]]^N@(3/2)]]

==================================================

Original reposted answer:

Per OP comment request, modified original function to calculate this sum:

sum[np_] := 
  Module[{r2 = Tuples[Tally[(Subtract @@@ Tuples[N@Range@np,{2}])^2], {3}]},
          Tr[Times @@@ Transpose[{r2[[2 ;;, 3, 1]]/Total[r2[[2 ;;, All, 1]], {2}]^(3/2), 
             Times @@@ r2[[2 ;;, All, 2]]}]]];

This should be usefully quicker, and reasonable time to execute up to low 100s arguments:

sum[100]//AbsoluteTiming

{2.97555,4.67618*10^7}