# Need help with the python code: Calculating Madelung constant CsCl crystal structure

Need help with the code to estimate the Madelung constant for CsCl lattice:

Cs at (0,0,0)

Cl at (0.5, 0.5, 0.5)

Answer: Converged value I am getting is 0.465. However, the correct answer is approximately 1.763.

What am I doing wrong? Code as below.

M = 0
L = 10

for i in range(-L,L+1):
for j in range(-L,L+1):
for k in range(-L,L+1):
if (i==0 and j==0 and k==0):
#negative ion
a = i+0.5
b = j+0.5
c = k+0.5
M += 1/np.sqrt(a*a + b*b + c*c)
else:
#postive ion
M -= 1/np.sqrt(i*i + j*j + k*k)
#negative ion
a = i+0.5
b = j+0.5
c = k+0.5
M += 1/np.sqrt(a*a + b*b + c*c)



Seems like you double count when accounting for offset.

import numpy as np

M = 0
L = 10
for i in range(-L, L + 1):
for j in range(-L, L + 1):
for k in range(-L, L + 1):
if i == 0 and j == 0 and k == 0:
continue
if (i+j+k) % 2 == 0:
sign = -1
else:
sign = 1
r = np.sqrt((i**2) + (j**2) + (k**2))
M += sign / r


which prints for $L=10$

Madelung constant: 1.6925789282594415


and for $L=100$

Madelung constant: 1.7418198158396654


You can also plot the convergence using something like this

import numpy as np
import matplotlib.pyplot as plt

M = 0
for i in range(-L, L + 1):
for j in range(-L, L + 1):
for k in range(-L, L + 1):
if i == 0 and j == 0 and k == 0:
continue
if (i+j+k) % 2 == 0:
sign = -1
else:
sign = 1
r = np.sqrt((i**2) + (j**2) + (k**2))
M += sign / r
return M

L_values = range(1, 50)
M_values = []

for L in L_values: