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, 2015 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, 2015 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 May 27, 2015 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, 2015 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)


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.")


Wikipedia cites ±1.747565 for NaCl. The following code gives For L = 200 M = -1.7446850421707383

# madelung.py: compute the Madelung constant for sodium chloride
from math import sqrt
L = int(input("Enter L: "))
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 i == j == k == 0:
continue
M += (-1)**(i+j+k) / sqrt( i*i + j*j + k*k )
print("For L =", L, "M =", M)
`