How to speed up three dimensional sum for Madelung Constant?
<p>For my Computational Physics class we have to compute the <a href="https://en.wikipedia.org/wiki/Madelung_constant" rel="nofollow">Madelung Constant</a> for NaCl. My code to do this uses three nested for loops and therefore runs very slowly. I was wondering if there was a way to use arrays or some other method to increase the speed of computation. Thanks </p>

<pre><code>from math import sqrt 
l= int(input("The lattice size is :"))
M = 0.0
for i in range(-L,L+1):
    for j in range(-L,L+1):
       for k in range(-L,L+1):
           if not (i==j==k==0):
              M += ((-1)**(i+j+k+1))/sqrt(i*i +j*j +k*k)

print("The Madelung constant for NaCl with lattice size",l,"is",M)
</code></pre>

<p>Here's an approach using open meshes with <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.ogrid.html" rel="nofollow"><code>np.ogrid</code></a> instead of creating actual meshes, which based on <a href="http://stackoverflow.com/a/39667342/3293881"><code>this other post</code></a> must be pretty efficient -</p>

<pre><code># Create open meshes for efficiency purposes
I,J,K = np.ogrid[-L:L+1,-L:L+1,-L:L+1]

# Perform all computations using those open meshes in a vectorized manner
all_vals = ((-1)**(I+J+K+1))/np.sqrt(I**2 +J**2+ K**2)

# Corresponding to "not (i==j==k==0)", which would be satisfied by 
# just one combination and that will occur at location (L,L,L), 
# let's set it as zero as a means to sum later on without adding for that elem
all_vals[L,L,L] = 0
M_out = all_vals.sum()
</code></pre>

python performance numpy vectorization