How to speed up three dimensional sum for Madelung Constant?

Question

For my Computational Physics class we have to compute the Madelung Constant 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

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)

Yes, but I don't understand how to use it for this problem.The professor has barely discussed how to use arrays other than to say that they can be used here. — Tom, Sep 29 '16 at 20:12

score 2 · Accepted Answer · answered Sep 29 '16 at 20:22

Since you've noted in a comment that you may use numpy, 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 numpy.where:

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()

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

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.

@Tom I suggest that you take a look at Divakar's new answer. Using `ogrid` should be much more efficient, as those huge 3d arrays are not stored in memory. His solution to fix the single singular point is also much simpler, and surely more efficient. — Andras Deak -- Слава Україні, Sep 30 '16 at 17:42

score 2 · Answer 2 · edited May 23 '17 at 10:27

2

Here's an approach using open meshes with np.ogrid instead of creating actual meshes, which based on this other post must be pretty efficient -

# 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()

edited May 23 '17 at 10:27

Community

1
1

answered Sep 30 '16 at 08:25

Divakar

218,885
19
262
358

I keep forgetting that `ogrid` exists...very nice answer:) – Andras Deak -- Слава Україні Sep 30 '16 at 17:42
1

@AndrasDeak Yeah even I haven't used it much. But those open meshes are just screaming for leveraging `broadcasting` and as such I should be looking more into it ;) – Divakar Sep 30 '16 at 18:20

How to speed up three dimensional sum for Madelung Constant?

2 Answers2