# Computing the Madelung constant

I am self teaching myself python and computational physics via Mark Newmans book Computational Physics the exercise is 2.9 of Computational Physics I have to compute the Madelung constant. .

I have two different solutions worked out I wanted to know which code is correct and is there a way to exclude 0 in my range so I don't have to write two separate for loops. I think my second answer is correct can somebody verify.

M = 0

for i in range(-200,0):
for j in range(-200,0):
for k in range(-200,0):

M += ((-1)**(i+j+k))/((i**2 + j**2 + k**2)**(1/2))

print(M)
for i in range(1,200):
for j in range(1,200):
for k in range(1,200):

M += ((-1)**(i+j+k))/((i**2 + j**2 + k**2)**(1/2))

print(M)


M = 0
i = 0
j = 0
k = 0
for l in range(-200,0):
i = l
j = l
k = l
M += ((-1)**(i+j+k))/((i**2 + j**2 + k**2)**(1/2))

print(M)

for l in range(1,200):
i = l
j = l
k = l
M += ((-1)**(i+j+k))/((i**2 + j**2 + k**2)**(1/2))

print(M)

• This seems like a question for either CodeReview or… is there a computational physics StackExchange? Or does it just belong on Physics? I'm not sure. Read their help, especially on-topic; if it fits better there, you can flag your own question to ask a moderator to migrate it. – abarnert May 27 '15 at 1:06
• Anyway, if this is a StackOverflow question, then it has to explain what you mean by "correct" in some way that makes sense to a programming expert, not a physics expert. Where do the two differ? – abarnert May 27 '15 at 1:07
• @ abarnert I was not sure since it involved programing if it was a purely physics based question and would be relevant there, for your second part comment I would like someone to explain which of the two is correct in the sense why one answer is correct not merely if it is right or wrong – FireFistAce May 27 '15 at 1:24
• Well, I explained the difference between the two, but I don't know which one matches what you want. – abarnert May 27 '15 at 1:30

Below is the code (in C, though) for the kind of loop you can write:

  for (int n = 1; n <= L; n++) {
for (int i = -n; i <= n; i++) {
for (int j = -n; j <= n; j++) {
for (int k = -n; k <= n; k++) {
if (abs(i) != n && abs(j) != n && abs(k) != n)
continue;

// Your update to the potential energy goes here.
}
}
}
}


Observe that this loop iterates over cubic shells of radius $n$. This is the way to properly compute the Madelung constant (see Section III of [1] for details).

[1] Borwein, D., Borwein, J. M., & Taylor, K. F. (1985). Convergence of lattice sums and Madelung’s constant. Journal of Mathematical Physics, 26(11), 2999. doi:10.1063/1.526675

There is a huge, obvious difference between the two.

The first one loops over every combination of negative i, j, k values, and every combination of positive i, j, k values. The second only loops over i == j == k values. So, for example, the first one has -200, -200, -200, then -200, -200, -199, and so on to -200, -200, -1, then -200, -199, -200, and so on; the second just goes right from -200, -200, -200 to -199, -199, -199.

But I have no idea which one is correct or how to verify it, because that's a physics question, not a Python question.

Meanwhile, there is a Python question hidden in there:

is there a way to exclude 0 in my range so I don't have to write two separate for loops.

Sure. You want an iterable over all the values in range(-200, 0) and all the values in range(1, 200). You can do that by just chaining the two iterables together:

import itertools

for i in itertools.chain(range(-200, 0), range(1, 200)):
print(i)


Or, if you prefer, just merge them manually and explicitly skip the 0:

for i in range(-200, 200):
if i == 0:
continue
print(i)


By the way, are you sure you wanted range(1, 200), not range(1, 201)? It seems strange that -200 is inside the set of values, but 200 isn't. But maybe that's correct.

I am working through the same book to teach myself Python. My working code is:

from math import *
from numpy import *

M = 0.0 #set a float variable for Madelung Constant
L = int(input("Enter the size of the crystal lattice L/2: "))

for x in range(-L, L+1):
for y in range(-L, L+1):
for z in range(-L, L+1):
if x == y == z == 0:
continue
r = (x**2 + y**2 + z**2)**-0.5
if x + y + z%2 == 0:
M = M+r
else:
M = M-r
print(x, y, z)
print(M)

• don't know how to make the code look pretty – Ryan Thomas Dec 28 '15 at 2:14

Try this :)

from math import sqrt,pi

e = 1.602e-19 #Coloumbs
epsNaught = 8.85e-12 #F/m
M = 0 #Madelung constant, total potential felt by origin sodium atom
n = 0 #Number of atoms

#lattice size
L = int(input("Enter lattice size: "))

#spacing of atoms
a = int(input("Distance between atoms: "))

for i in range(-L,L+1):
for j in range(-L,L+1):
for k in range(-L,L+1):
n += 1
distance = a * sqrt(i**2 + j**2 + k**2)

if (i == j == k == 0): #Case for soidium atom at origin
continue

potential = e / (4 * pi * epsNaught * distance)

if (i+j+k)%2 == 1: #Odd, chlorine atoms, attraction = neg potential
potential *= -1

M += potential

print("Madelung constant:",M,"due to",n-1,"atoms.")
`