# Calculating Madelung constant via Ewald summation

I am attempting to calculate the Madelung constant for NaCl using an Ewald sum derived by Nijboer. There are other methods and published codes to do this, but I am specifically interested in working within Nijboer's framework for a more complicated problem unrelated to a Madelung constant.

Ref: Nijboer, B. R. A., De Wette, F. W. (1957). On the calculation of lattice sums. Physica, 23(1–5), 309–321.

My work so far, Nijboer starts with the sum

$$\alpha_m=\sum_{\lambda=- \infty}^{\infty '}\dfrac{(-1)^{\lambda_1+\lambda_2+\lambda_3}}{r_\lambda}$$

Where $$\lambda=( \lambda_1 ,\lambda_2 ,\lambda_3)$$ and $$\lambda_i$$ take all integer values. The prime on the sum indicates that the index corresponding to $$\lambda_1=\lambda_2=\lambda_3=0$$ is omitted.

Then Nijboer rewrites the series as an Ewald sum

$$\alpha_m=\sum_{\lambda=- \infty}^{\infty '}(-1)^{\lambda_1+\lambda_2+\lambda_3}{r_\lambda} \dfrac{\Phi(\sqrt{\pi}r_\lambda)}{r_\lambda}+ \dfrac{1}{\pi v_a} \sum_{\lambda=- \infty}^{\infty}\dfrac{e^{\pi |\boldsymbol{h}_\lambda - \boldsymbol{k}_{1/2}|^2}}{|\boldsymbol{h}_\lambda - \boldsymbol{k}_{1/2}|^2}$$

The value $$r_\lambda$$ is the Euclidean distance of $$\boldsymbol{r_\lambda}=\lambda_1 \boldsymbol{a}_1 +\lambda_2 \boldsymbol{a}_2+\lambda_3 \boldsymbol{a}_3$$ where $$\boldsymbol{a}_i$$ corresponds to the primitive lattice vectors. For NaCl, these vectors are: $$\boldsymbol{a}_1=[0,1,1]$$, $$\boldsymbol{a}_2=[1,0,1]$$, $$\boldsymbol{a}_3=[1,1,0]$$

The reciprocal lattice vectors are given by $$\boldsymbol{h}_\lambda=\lambda_1 \boldsymbol{b}_1 + \lambda_2 \boldsymbol{b}_2 + \lambda_3 \boldsymbol{b}_3$$

Where $$\boldsymbol{b}_1=2\pi \dfrac{\boldsymbol{a}_2 \times \boldsymbol{a}_3}{\boldsymbol{a}_1 \cdot [\boldsymbol{a}_2 \times \boldsymbol{a}_3]}$$, $$\boldsymbol{b}_2=2\pi \dfrac{\boldsymbol{a}_3 \times \boldsymbol{a}_1}{\boldsymbol{a}_1 \cdot [\boldsymbol{a}_2 \times \boldsymbol{a}_3]}$$, $$\boldsymbol{b}_3=2\pi \dfrac{\boldsymbol{a}_1 \times \boldsymbol{a}_3}{\boldsymbol{a}_1 \cdot [\boldsymbol{a}_2 \times \boldsymbol{a}_3]}$$ Ref: Kittel

and $$\boldsymbol{k}_{1/2}=\dfrac{1}{2} (\boldsymbol{b}_1 + \boldsymbol{b}_2 + \boldsymbol{b}_3)$$

The expression $$|\boldsymbol{h}_\lambda - \boldsymbol{k}_{1/2}|^2$$ is the Euclidean distance between the reciprocal vector and the reciprocal lattice vector squared. The cell volume $$v_a$$ is defined as: $$v_a=\boldsymbol{a}_1 \cdot [\boldsymbol{a}_2 \times \boldsymbol{a}_3]$$

Finally, the function $$\Phi$$ is the complementary error function.

I believe(d) this correctly defined all parameters needed to compute the sum, unfortunately I do not get a meaningful answer at all (both sums converge rapidly but to no effect). I essentially get a Madelung constant $$\alpha = -2$$ versus the correct answer $$\alpha = 1.75$$.

import numpy as np
from scipy.special import gammainc, erfc

#Here primitive vectors are defined
a1 = np.array([0,1,1])
a2 = np.array([1,0,1])
a3 = np.array([1,1,0])

#Now the cell volume is calculated
A = np.array([a1,a2,a3])
va = np.abs(np.linalg.det(A))

#Here the reciprocal vectors are
b1 = 2*np.pi*np.divide(np.cross(a2,a3),np.dot(a1,np.cross(a2,a3)))
b2 = 2*np.pi*np.divide(np.cross(a3,a1),np.dot(a1,np.cross(a2,a3)))
b3 = 2*np.pi*np.divide(np.cross(a1,a2),np.dot(a1,np.cross(a2,a3)))

#Initiate the real and reciprocal space variables
a_real = 0
a_recip = 0

#Iterate across lambdas, hopping across lattice points
for lambda1 in range(-4,4):
for lambda2 in range(-4,4):
for lambda3 in range(-4,4):

r_lambda = (lambda1*a1+lambda2*a2+lambda3*a3) #Real lattice vector
h_lambda = (lambda1*b1+lambda2*b2+lambda3*b3) #Reciprocal lattice vector

k_half=1/2*(b1+b2+b3)

recipNorm = np.linalg.norm(h_lambda-k_half)**2
realNorm = np.linalg.norm(r_lambda)**2 #Assume that the basis atom is at (0,0,0)

if lambda1 == 0 and lambda2 == 0 and lambda3 == 0: #do not count the real sum at the singularity
a_recip += 1/(np.pi*va)*np.exp(-np.pi*recipNorm)/recipNorm
else:
a_real += (-1)**(lambda1+lambda2+lambda3) * erfc(np.sqrt(np.pi)*realNorm)/realNorm
a_recip += 1/(np.pi*va)*np.exp(-np.pi*recipNorm)/recipNorm



I think something is wrong with how I am handling lattice vectors, but I can't figure out what. I am defining everything as per definitions found in Kittel's textbook and Nijboer's paper. Any help at all would be so appreciated.

• Unfortunately there are two common definitions of the reciprocal lattice vectors which differ by a factor of 2*pi and it looks as though Kittel and the paper use different forms - see equation (4) in the paper and compare that with what you have above (I prefer the version in the paper). Basically drop the 2*pi in front of the recip lat vecs and see what happens. I haven't checked if this fixes it or looked for other errors. Sep 22, 2023 at 5:47
• Also note your equation for alpha as typeset in the image is incorrect - in the numerator of the recip space term the term with the lattice vectors should be inside the exponential - as it is an image I can't correct it. Sep 22, 2023 at 5:49
• OK, just about to go out but other errors are 1) You are using the square of the norm in real space when it should just be the norm 2) The condition should be if all of the lambdas are zero, not the sum of them, i.e. just skip one term at r=0 3) You should be using the simple cubic unit cell. Correcting those I get -1.7475645946331826 (with a range set to 8). I don't use the Madelung constant in my own work, so I assume there's a sign convention issue as well. Sep 22, 2023 at 6:36
• Thank you so much! I am able to replicate the result. This mostly makes sense, but I am confused by the use of the SC unit cell. NaCl has an FCC structure, no? I replicate the same result either by using the SC unit cell or by dividing the basis vectors of the FCC cell by two. Sep 22, 2023 at 13:43
• I need to think on the cell (when I will answer properly) but it could be because you aren't using unit vectors for the lattice. Sep 22, 2023 at 14:16

There are a number of issues here (as I noted in the comments)

Firstly you have the maths wrong in a number of places, both in the question and the implementation. The correct final formula should be

$$\alpha_m=\sum_{\lambda=- \infty}^{\infty '}(-1)^{\lambda_1+\lambda_2+\lambda_3} \dfrac{\Phi(\sqrt{\pi}r_\lambda)}{r_\lambda}+ \dfrac{1}{\pi v_a} \sum_{\lambda=- \infty}^{\infty}\dfrac{e^{\pi |\boldsymbol{h}_\lambda - \boldsymbol{k}_{1/2}|^2}}{|\boldsymbol{h}_\lambda - \boldsymbol{k}_{1/2}|^2}$$

Note I have dropped an $$r_\lambda$$ in the first, real space term - and also note your original form can not be correct as the first and second terms have different dimensions - dimensionless for the real space term and $$\bf{L}^{-1}$$ for the reciprocal space one.

In the implementation you have two main "mathematical" mistakes

• You use the square of the length of the real space vector, so you are missing a square root
• The sum should be over all terms except when all the $$\lambda$$s are zero, not when the sum of the $$\lambda$$s is zero (so you could say the missing term is when the product of the $$\lambda$$s is zero, but I think that is less clear). What you are doing physically is omitting the r=0 term, for which the Coulomb potential diverges
• (I don't know Python, but I would also avoid writing 1/2 - in other languages such as Fortran or C this would be an integer division, which would evaluate to zero, a common bug, so it's a good habit to avoid in general IMO)

There's also a couple of problems with the lattice vectors

• In real space the primitive vectors are too large by a factor of 2 - see here Effectively you are using different units for the summation.
• In reciprocal space you have run into the problem of there being two common conventions for the reciprocal lattice vectors, which differ by a factor of $$2\pi$$. The paper uses the convention that $$\boldsymbol{b}_1= \dfrac{\boldsymbol{a}_2 \times \boldsymbol{a}_3}{\boldsymbol{a}_1 \cdot [\boldsymbol{a}_2 \times \boldsymbol{a}_3]}$$ et cetera. Unfortunately your reference book, Kittel, uses the other convention. Ultimately this is related to whether you define a Fourier transform using $$e^{2\pi i {\bf f}\cdot.{\bf r}}$$ or $$e^{i {\bf k}\cdot.{\bf r}}$$

Correcting all this gives (Note I am not a python programmer, this si all by analogy to other languages I know, so please check!!):

ijb@ijb-Latitude-5410:~/work/stack\$ cat mad.py
import numpy as np
from scipy.special import gammainc, erfc

#Here primitive vectors are defined
a1 = np.array([0,1,1]) * 0.5 # IJB, correct primitve vectors
a2 = np.array([1,0,1]) * 0.5
a3 = np.array([1,1,0]) * 0.5

#Now the cell volume is calculated
A = np.array([a1,a2,a3])
va = np.abs(np.linalg.det(A))

#Here the reciprocal vectors are
#b1 = 2*np.pi*np.divide(np.cross(a2,a3),np.dot(a1,np.cross(a2,a3)))
#b2 = 2*np.pi*np.divide(np.cross(a3,a1),np.dot(a1,np.cross(a2,a3)))
#b3 = 2*np.pi*np.divide(np.cross(a1,a2),np.dot(a1,np.cross(a2,a3)))
# IJB, correct recip vecs to the correct convention
b1 = np.divide(np.cross(a2,a3),np.dot(a1,np.cross(a2,a3)))
b2 = np.divide(np.cross(a3,a1),np.dot(a1,np.cross(a2,a3)))
b3 = np.divide(np.cross(a1,a2),np.dot(a1,np.cross(a2,a3)))

#Initiate the real and reciprocal space variables
a_real = 0
a_recip = 0

#Iterate across lambdas, hopping across lattice points
for lambda1 in range(-4,4):
for lambda2 in range(-4,4):
for lambda3 in range(-4,4):

r_lambda = (lambda1*a1+lambda2*a2+lambda3*a3) #Real lattice vector
h_lambda = (lambda1*b1+lambda2*b2+lambda3*b3) #Reciprocal lattice vector

#            k_half=1/2*(b1+b2+b3)
k_half=0.5*(b1+b2+b3)

# IJB, give it a better name to reflect what it is
#            recipNorm = np.linalg.norm(h_lambda-k_half)**2
recipNorm_sq = np.linalg.norm(h_lambda-k_half)**2
# IJB, correct to not use square
#            realNorm     = np.linalg.norm(r_lambda)**2 #Assume that the basis atom is at (0,0,0)
realNorm     = np.linalg.norm(r_lambda) #Assume that the basis atom is at (0,0,0)

if lambda1 == 0 and lambda2 == 0 and lambda3 == 0: #do not count the real sum at the singularity
a_recip += 1/(np.pi*va)*np.exp(-np.pi*recipNorm_sq)/recipNorm_sq
else:
a_real += (-1)**(lambda1+lambda2+lambda3) * erfc(np.sqrt(np.pi)*realNorm)/realNorm
a_recip += 1/(np.pi*va)*np.exp(-np.pi*recipNorm_sq)/recipNorm_sq

ijb@ijb-Latitude-5410:~/work/stack$$python3 mad.py Calculated: -1.7475645919728273 Reference : 1.747565 AbsDiff : -4.080271727158191e-07 Frac error: -2.3348325968751897e-07 ijb@ijb-Latitude-5410:~/work/stack$$