Evaluating the Madelung Constant

Question

The Madelung constant is used in determining the electrostatic potential of a single ion in a crystal by approximating the ions by point charges:

$$ V_i = \frac{e}{4\pi\epsilon_0r_0}\sum_j\frac{z_jr_0}{r_{ij}} = \frac{e}{4\pi\epsilon_0r_0}M_i $$

where $r_0$ is the nearest neighbour distance, and the Madelung constant of the ith ion is given by the sum over all the other ions in the crystal:

$$ M_i = \sum_j\frac{z_j}{r_{ij}/r_0}. $$

For the cubic crystal NaCl, the summation may be performed over three orthogonal co-ordinates:

$$ M_\mathrm{Na} = \sum_{j,k,l=-\infty}^{\infty}' \frac{(-1)^{j+k+l}}{\sqrt{j^2+k^2+l^2}}, $$

where the prime indicates that the term $(0,0,0)$ is excluded. Benson's formula provides a practical and efficient way of evaluating this sum as:

$$ M_\mathrm{Na^+} = -12\pi\sum_{m,n=1,3,...}^{\infty}\mathrm{sech}^2\left(\frac{1}{2}\pi\sqrt{m^2+n^2}\right), $$

where the summation is performed over all positive odd integers $m$ and $n$.

Write a program to calculate the Madelung constant for $\mathrm{Na^+}$ ions in NaCl, -1.74756...

Solution

import math

nmax = 13
M = 0.
for n in range(1,nmax,2):
    for m in range(n,nmax,2):
        fac = 2
        if n == m:
            fac = 1
        M += fac * math.cosh(0.5*math.pi*math.hypot(n, m))**-2
print(-12*math.pi*M)

import math

M = 0.
n = 1
while True:
    oldM = M
    for m in range(1,n+1,2):
            fac = 2
            if n == m:
                fac = 1
            M += fac * math.cosh(0.5*math.pi*math.hypot(n, m))**-2
    if abs(M - oldM) < 0.0001:
        break
    n += 2
print(-12*math.pi*M)