# Calculating a 2D Ewald sum for a multipolar expansion

I am attempting to calculate the potential of a particle at the center of an infinite two-dimensional lattice as per the following reference:

The authors represent the multipolar expansion as a series of spherical harmonics which, for certain degrees/order of harmonics, do not converge (or converge only slowly). They provide an Ewald sum to ease the calculation. I am trying to replicate Table I, which gives some well-known results regarding the self-potentials of a particle situated at the (0,0,0) position in a 2D hexagonal lattice for multipole moments ranging from two to six. I implement both the Ewald sum given in (5) and the direct sum given in (11). The solution I calculate to the series converges rapidly, but to the wrong answer. I have verified:

1. The mpmath convention for the spherical harmonics matches that referenced by the author
2. The polar and azimuthal angles are not being mixed up
3. My convergence functions match the simplifications given for n=0 in the description of (7) (e.g. $$R_0(X)=\sqrt(\pi)erfc(X)/X$$)
4. My definition of the primitive lattice vectors are correct, or at least they successfully lead to the correct number in the direct sum employed at the end of the script
5. The mpmath implementation of the upper incomplete gamma function evaluates conventionally for negative "n" values (verified against the Mathematica upper incomplete gamma function).

I am at my wits end with this. It seems simple but I can't get it right. I do not think the author's did anything wrong because these results are used for future work and match prior literature well. If anyone has experience with this, I would love some feedback.

# -*- coding: utf-8 -*-
"""
This code computes two-dimensional Ewald sums for periodic lattices per the
method prescribed by Senet and Lambin in "Ewald Summation of Multipolar
Interactions at an Arbitrary Order on a Two-Dimensional Lattice" and later
referenced in Stone's "Lattice sums and their derivatives for surface adlayers"
"""

import numpy as np
from mpmath import gammainc, inf, gamma, spherharm

#Define the real-space convergence function
def R(X,n):
return X**(2*n-1)*gammainc(-n+0.5,X**2,inf,regularized=False)

#Define the reciprocal-sapce convergence function
def S(X,n):
return gammainc(n+0.5,X**2,inf,regularized=False)/gamma(n+0.5)

#Convert from Cartesian to spherical coordinates
def cart2sph(vector):
x, y, z = vector
hxy = np.hypot(x, y)
r = np.hypot(hxy, z)
theta = np.arctan2(z, hxy)
phi = np.arctan2(y, x)
return r, theta,phi

#Kronecker delta functions (unused at the moment)
def kdelta(x,y):
if x == y:
out = 1
else:
out = 0
return out

#Define unit cell vectors (simple cubic)
a = 1 #Lattice constant
c = a
a1 = np.array([a,0,0])
a2 = np.array([a/2,np.sqrt(3)/2*a,0])
a3 = np.array([0,0,c])

A = np.linalg.norm(np.cross(a1,a2)) #Unit cell area

#Define the reciprocal space vectors
Vc = abs(np.dot(a1,np.cross(a2,a3))) #Unit cell volume
b1 = 2*np.pi/Vc*np.cross(a2,a3)
b2 = 2*np.pi/Vc*np.cross(a3,a1)
b3 = 2*np.pi/Vc*np.cross(a1,a2)

#Define the degree and order of multipole to be inspected
l = 6
m = 0

#Initialize the parts of the sum
recipSum = 0
realSum = 0
const = 0

#Displacement at which to inspect to the multipoles - just use the (0,0,0)
#particle as a referenece
rho = np.array([0,0,0])

#Wavevector at which to the inspect the multipoles - just using Q = 0 for now,
#but according to the author, any vector in the first Brillouin zone works?
Q = np.array([0,0,0])

#Ewald convergence parameter
xi = 2*np.sqrt(np.pi)*a/((np.linalg.norm(a1)+np.linalg.norm(a2)))

#Number of terms to consider - usually converges to a few sig figs after
#a couple iterations
N=20

for lambda1 in range(-N,N):
for lambda2 in range(-N,N):

G =  lambda1*b1 + lambda2*b2 #Define the translation vector
g = np.linalg.norm(G) #Calculate the length of the translation vector
phi = np.arctan2(G[1],G[0]) #Compute the angle in reciprocal space
#(irrelevant for m = 0)

#Compute the contribution to the reciprocal part of the sum
recipSum += np.pi*a**2*xi**(l-1) \
*np.exp(np.dot(-(1j)*(Q+G),rho)) \
*(a*g/(2*xi))**abs(m) \
*R(a*g/(2*xi),0.5*(l-abs(m))) \
*(1j)**abs(m)*np.conjugate(spherharm(l,m,np.pi/2,phi))

for lambda1 in range(-N,N):
for lambda2 in range(-N,N):

T =  lambda1*a1 + lambda2*a2 #Define the reciprocal translation vector
t = np.linalg.norm(T-rho) #Calculate the length of the translation vector
phi = np.arctan2(T[1],T[0]) #Compute the real space angle

if lambda1 == 0 and lambda2 == 0: #Valid if rho = (0,0,0), stop the sum and add the constant

const += (-xi)**(l+1) \
/gamma(0.5*(l+3)) \
*np.exp(1j*np.dot(Q,T)) \
*np.conjugate(spherharm(l,m,np.pi/2,0))

continue #Skip the next addition of the real part

#Compute the contribution to the real part of the sum
realSum += np.exp((1j)*np.dot(Q,T)) \
*(a/t)**(l+1) \
*S(xi*t/a,0.5*(abs(m)+l)) \
*np.conjugate(spherharm(l,m,np.pi/2,phi))

#Convert the results to floats
recipSum = float(recipSum.real)
realSum = float(realSum.real)
const = float(const.real)

total = const + realSum + recipSum

#Print the outputs
print("The reciprocal sum evaluates to: " + str(recipSum))
print("The real sum evaluates to: " + str(realSum))
print("The constant is: " + str(const))
print("The total is: " + str((realSum+recipSum+const)))

#Check this against the direct summation
directSum = 0
N=40

for lambda1 in range(-N,N):
for lambda2 in range(-N,N):
T = lambda1*a1 + lambda2*a2

if lambda1 == 0 and lambda2 == 0:
continue
else:
realMag, theta,phi = cart2sph(T)
directSum += (a/realMag)**(l+1)

print('The direct sum is: '  + str(directSum))  #Known not to converge for l <= 2, but is usually close anyway
print('The literature value is: '  + str(6.19524))

print(float(directSum.real))

• What have you already tried to do to find the problem? Apr 3 at 22:11
• I saw four possibilities for hiccups. 1: the implementation needs to be symbolically faithful to the author's equation, there are no typos. I have quintuple checked that the keystrokes are referring to the correct quantities. 2: the definition of the translation vectors must be correct, you need the right lattice. I confirm this by finding the correct direct summation value. 3: the definitions of the spherical harmonic/gamma functions much match what the authors used. I confirmed this. 4: the vectors $\rho$ and $Q$ must be physically correct. I am directly using values they gave under Eq (11) Apr 3 at 23:54
• Additionally I have read and reread several papers within the citation web of the one I linked, spanning from Nijboer in 1957 to Stone in 2005, and after all of this I feel I understand the physics and mathematics well. In the past, I have used a similar Ewald sum to correctly calculate Madelung constants, so I think that my overall framework applying these methods numerically is correct. This problem seems so simple (I am just copying someone!) and I am at a loss for what isn't working and why... Apr 3 at 23:57