Computing the Madelung constant

Question

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.

My first answer:

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)

My second answer:

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)

Answer 1

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

Answer 2

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.

Answer 3

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)

Answer 4

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