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>Since you've noted in a comment that you may use <code>numpy</code>, I suggest doing that. You can construct a 3d grid for your integers, and compute each term simultaneously, thereby vectorizing your calculation. You only need to watch out for the singular case where every integer is 0, for instance using <a href="https://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html" rel="nofollow"><code>numpy.where</code></a>:</p>

<pre><code>import numpy as np

ran = np.arange(-L,L+1)
i,j,k = np.meshgrid(ran,ran,ran)
M = np.where((i!=0) | (j!=0) | (k!=0),
             (-1)**(i+j+k+1)/np.sqrt(i**2+j**2+k**2),
             0).sum()
</code></pre>

<p><code>ran</code> is a numpy array with the same elements as in <code>range()</code> (if cast to a list in python 3). <code>meshgrid</code> then constructs three 3d arrays, which together span the 3d space where you need to perform your sum.</p>

<p>Note that for large domains this approach has much higher memory need. This is ususal with vectorization: you can spare CPU time at the cost of increased memory need.</p>

python performance numpy vectorization