# Improve Mandelung constant code

I'm learning and improving my Python skills.

I did a program in Python about Mandelung constant. But, I'm having kind of a problem. The Mangelung constant is calculated using this sum:

$$V_{total} = \sum_{\substack{i,j,k = -L\\ i, j, k \neq 0}}^L V(i,j,k)= \frac e{4\pi\epsilon_0} M$$

Or

$$M =\sum_{\substack{i,j,k = -L\\ i, j, k \neq 0}}^L \frac{(-1)^{i + j + k}}{\sqrt {i^2 + j^2 + k^2}}$$

When I run it, it takes a long time, when I put a huge number. So, I need to improve my code, to run it faster. Can someone explain me a way to do it? (importing others libs, using other stuff)

The code that I did:

import time

start_time = time.time()
L = int(input("Put the number of L:")) # size of the lattice
L = L+1 # this is for the vector (0,0,0)
# n = 0 # number of atoms
M = 0 # Madelung constant
for i in range(-L,L+1):
for j in range(-L,L+1):
for k in range(-L,L+1):
# n += 1 #counter for number of atoms
if i == j == k == 0: # doesn't count the origin (0,0,0)
continue
r = (i**2 + j**2 + k**2)**(-0.5)
if (i + j + k) % 2 == 1: # odd number
r *= -1
M += r
print ("Mandelung Constant is::", M)
print("It takes %s seconds" % (time.time() - start_time))


When I put a $$L = 300$$, it takes more than 7 minutes. This is why I'm trying to improve it.

• Python for loops are slow. Trade memory for time for building a cube in numpy and correcting for 0,0,0 entry. – Richard Aug 29 at 17:19

As mentioned by @Richard, loops in Python are slow. Two solutions come to my mind:

• Use NumPy and vectorize the operations. This will speed up your calculations at the cost of storing your intermediate arrays in memory.

• Wrap your calculation in a function and use Numba's jit decorator to (magically) accelerate your calculation. This implies refactoring the code sometimes.

Using NumPy capabilities you can rewrite to have something like the following:

import numpy as np
import time

start_time = time.time()
L = 300
i = np.array(range(-L, L + 1), dtype=np.float)
I, J, K = np.meshgrid(i, i, i)

M = (-1)**(I + J + K)/np.sqrt(I**2 + J**2 + K**2)
M[(I == 0)*(J == 0)*(K == 0)] = 0
M = np.sum(M)
print ("Madelung Constant is::", M)
print("It takes %s seconds" % (time.time() - start_time))


And this gives a result

Madelung Constant is:: -1.7456432959005515
It takes 17.750624418258667 seconds


Compared with

Madelung Constant is:: -1.745643295911936
It takes 310.843811750412 seconds


You could also use Numba.

import time
from numba import jit

@jit
M = 0 # Madelung constant
for i in range(-L , L + 1):
for j in range(-L, L + 1):
for k in range(-L, L + 1):
if i == j == k == 0: # doesn't count the origin (0,0,0)
continue
r = (i**2 + j**2 + k**2)**(-0.5)
if (i + j + k) % 2 == 1: # odd number
r *= -1
M += r
return M

start_time = time.time()
print ("Madelung Constant is::", M)
print("It takes %s seconds" % (time.time() - start_time))


with the result

Madelung Constant is:: -1.7456432959126562
It takes 2.8461902141571045 seconds

• I barely do python, but ... The last calculated value is markedly different from the other 2, are you sure you are doing the same calculation? Further the Madelung constant for NaCl structures is 1.748 (e.g. en.wikipedia.org/wiki/Madelung_constant) which disagrees with both sets of numbers, but a little Fortran implementation I did of the above gets 1.7475645946509435 (in about 0.7s). I just want to be sure that the various codes are doing what they are supposed to be. Is python single precision unless otherwise stated? – Ian Bush Aug 30 at 9:04
• @IanBush, I just reran all the snippets and updated the values. Regarding the value of 1.748, I just assumed that the summation from the OP was right. – nicoguaro Aug 30 at 16:02
• OK. Just really curious I can't reproduce the values you get with the language I know best in either single or double precision, and as far as I understand what is going on the Fortran version gets the "right" answer. But we are summing conditionally convergent series here so their could be issues over and above simple floating point maths. – Ian Bush Aug 30 at 16:08
• @IanBush, did you use the same summation? How many terms did you use? – nicoguaro Aug 30 at 16:13
• Yes, as I understand things, and exactly the same number of terms, even down to the lack of symmetry in the range (which is probably technically wrong). And the numpy version gets the same value on my machine as you show above, and changing it to double doesn't seem to affect things much. But I don't really speak python. Let me play over the weekend and see if I can work out what is going on. – Ian Bush Aug 30 at 16:26