1
$\begingroup$

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. enter image description here.

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)
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – abarnert May 27 '15 at 1:06
  • $\begingroup$ 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? $\endgroup$ – abarnert May 27 '15 at 1:07
  • $\begingroup$ @ 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 $\endgroup$ – FireFistAce May 27 '15 at 1:24
  • 1
    $\begingroup$ Well, I explained the difference between the two, but I don't know which one matches what you want. $\endgroup$ – abarnert May 27 '15 at 1:30
3
$\begingroup$

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

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ don't know how to make the code look pretty $\endgroup$ – Ryan Thomas Dec 28 '15 at 2:14
0
$\begingroup$

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.