2
$\begingroup$

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)

print('Madelung constant:', M)
$\endgroup$

1 Answer 1

2
$\begingroup$

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
print('Madelung constant:', M)

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

def calculate_madelung(L):
    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
    print('Madelung constant:', M)
    return M

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

for L in L_values:
    M = calculate_madelung(L)
    M_values.append(M)
    print(f'L = {L}, Madelung constant = {M}')

plt.plot(L_values, M_values, marker='o')
plt.xlabel('L')
plt.ylabel('Madelung Constant')
plt.title('Convergence of Madelung Constant for CsCl')
plt.grid(True)
plt.show()

Convergence

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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