Feb 1, 2000

Notes for Lecture Notes for Lecture #4

Electrostatic energy

Section 1.11 in Jackson's text

The total electrostatic potential energy of interaction between point charges {qi} at the positions {ri} is given by

W = 1
4pe0

å
i < j 
qiqj
|ri-rj|
= 1
8pe0

å
i ¹ j 
qi qj
|ri-rj|
.
(1)

In this expression, the first form explicitly counts all pairs, while the second form counts all interactions and divides by 2 to compensate for double counting.

For a finite system of charges, this expression can be evaluated directly, however, for a large or infinite system, the expression 1 does not converge and numerical tricks must be used to evaluate the energy. In fact, for the infinite system, one has an infinite amount of charge and the energy of interaction is undefined. If the system is neutral, it is possible to define a meaningful interaction energy by use of an Ewald transformation.

The basic idea of the Ewald approach is as follows. The error function erf(x) and its complement erfc(x) are defined as:

erf(x) = 2
Ö
p
 
ó
õ
x

0 
e-t2 dt     \and       erfc(x) = 1 - erf(x) = 2
Ö
p
 
ó
õ
¥

x 
e-t2 dt.
(2)

Ewald noted that

1
r
=
erf( 1
2
Ö
h
 
r)

r
+
erfc( 1
2
Ö
h
 
r)

r
.
(3)
In this expression, the first term goes to a constant (Ö{h/p}) as r ® 0, but has a long range tail as r ® ¥. The second term has a singular behavior as r ® 0, but vanishes exponentially as r ® ¥. Thus, Ewald's idea is to replace a single divergent summation with two convergent summations. The first summation has a convergent summation in the form of its Fourier transform and the second has a convergent direct summation. Thus the calculation of the electrostatic energy would be evaluated using:

W = 1
8pe0

å
i ¹ j 
qiqj
|ri-rj|
= 1
8pe0
æ
ç
ç
ç
ç
è

å
i ¹ j 
qi qjerf( 1
2
Ö
h
 
|ri-rj|)

|ri-rj|
+
å
i ¹ j 
qi qjerfc( 1
2
Ö
h
 
|ri-rj|)

|ri-rj|
ö
÷
÷
÷
÷
ø
.
(4)

For an appropriate choice of the parameter h, the second summation in Eq. 4 converges quickly and can be evaluated directly. The first term in the summation of Eq. 4 must be transformed into Fourier space.

In order to described these summations explicitly, we assume that we have a periodic lattice so that every ion can be located by ri = ta + T a location ta within a unit cell and a periodic translation vector T, In this way, the summation becomes:1


å
ij 
= N
å
abT 
,
(5)
where N denotes the number of unit cells in the system. Since we have a periodic system, N is infinite, but the energy per unit cell W/N is well defined. The other identity that we must use is that a sum over lattice translations T may be transformed into an equivalent sum over ``reciprocal lattice translations'' G according to the identity:

å
T 
d3(r - T) = 1
W

å
G 
ei G ·r,
(6)
where W denotes the unit cell volume. The proof of this relation is given in the appendix.

The first term of Eq. 4 thus becomes


å
i ¹ j 
qi qjerf( 1
2
Ö
h
 
|ri-rj|))

|ri-rj|
= N æ
ç
ç
ç
ç
è

å
ab 
qa qb
å
T 
erf( 1
2
Ö
h
 
|ta-tb+ T |)

|ta-tb + T|
-
å
a 
q2a   æ
 ú
Ö

h
p
 
ö
÷
÷
÷
÷
ø
,
(7)
where the last term in Eq. 7 comes from subtracting out the self-interaction (i = j) term from the complete lattice sum. Using the short hand notation, tab º ta-tb, the lattice sum can be evaluated:

å
T 
erf( 1
2
Ö
h
 
|tab + T|)

|tab + T|
= ó
õ
  d3r
å
T 
d3(r-T)
erf( 1
2
Ö
h
 
|tab + r|)

|tab + r|
.
(8)
This becomes,
1
W

å
G 
ó
õ
  d3r ei G ·r
erf( 1
2
Ö
h
 
|tab + r|)

|tab + r|
= 4 p
W
æ
ç
è

å
G ¹ 0 
e- i G ·tab   e- G2/h2
G2
+ 1
2
ó
õ
1/2Ö{h}

0 
du
u3
ö
÷
ø
,
(9)
where the last term, which is infinite, comes from the G = 0 contribution.

If the last term of Eq. 9 cannot be eliminated, it is clear that the electrostatic energy is infinite. The term can be eliminated if and only if the system is neutral. Thus, it is only meaningful to calculate the electrostatic energy of a neutral periodic system. If the actual system has a net charge of Q º åa qa, we could calculate a meanful energy if we subtract the interaction energy due to a uniform density of Q/W. Taking all of these terms into account, we find the final Ewald expression to be:

W
N
=
å
ab 
qa qb
8 pe0
æ
ç
ç
ç
ç
è
4 p
W

å
G ¹ 0 
e- i G ·tab   e-G2/h
G2
-   æ
 ú
Ö

h
p
 
dab + ¢
å
T 
erfc( 1
2
Ö
h
 
|tab+T|)

|tab+T|
ö
÷
÷
÷
÷
ø
- 4 pQ2
8 pe0 Wh
,
(10)
where the ¢ in the summation over lattice translations T indicates that all self-interaction terms should be omitted.

Appendix I - comments on lattice vectors and reciprocal lattice vectors

In this discussion, will assume we have a 3-dimensional periodic system. It can be easily generalized to 1- or 2- dimensional systems. In general, a translation vector can be described a linear combination of the three primitive translation vectors T1, T2, and T3:

T = n1T1 + n2T2 + n3T3,
(11)
where {n1, n2,n3} are integers. Note that the unit cell volume W can be expressed in terms of the primitive translation vectors according to:
W = |T1 ·(T2 ×T3)|.
(12)
The reciprocal lattice vectors G can generally be written as a linear combination of the three primitive reciprocal lattice vectors G1, G2, and G3:
G = m1G1 + m2G2 + m3G3,
(13)
where {m1, m2,m3} are integers. The primitive reciprocal lattice vectors are determined from the primitive translation vectors according to the identities:
Gi ·Tj = 2 pdij.
(14)
Note that the ``volume'' of the primitive reciprocal lattice is given by
|G1 ·(G2 ×G3)| = (2 p)3
W
.
(15)

Some examples of this are given below.

``Proof'' of Eq. 6

Consider the geometric series

+M
å
k = -M 
ei k (G1 ·r) =
sin æ
ç
è
(M+ 1
2
)G1 ·r ö
÷
ø

sin( G1 ·r/2)
,
(16)
where the summation limit M will be taken in the limit M® ¥:

lim
M ® ¥ 
sin æ
ç
è
(M+ 1
2
)G1 ·r ö
÷
ø

sin( G1 ·r/2)
= 2 p
å
n1 
d(G1 ·(r-n1T1).
(17)
The summation over all lattice translations n1T1 is due to the fact that sin( G1 ·r/2) = 0 whenever r = n1T1. Carrying out the geometric summation in Eq. 6 for all three reciprocal lattice vectors and taking the limit as in Eq. 17,

å
G 
ei G ·r = (2 p)3
å
n1,n2,n3 
d(G1 · (r-n1T1)) d(G2 ·(r-n2T2)) d(G3 · (r-n3T3)).
(18)
Finally, the right hand side of Eq. 18 can be simplified by transforming the d functions into their spatial form:

å
G 
ei G ·r = (2 p)3
|G1 ·(G2 ×G3)|

å
T 
   d3(r-T),
(19)
which is consistent with Eq. 6.

Appendix II - examples

In these examples, denote the length of the unit cell by a.

CsCl structure

There are two kinds of sites - tCs = 0 and tCl = a/2 ([^(x)] + [^(y)] +[^(z)]).

T1 = a ^
x
 
      T2 = a ^
y
 
      T3 = a ^
z
 
.
(20)

G1 = 2 p
a
^
x
 
      G2 = 2 p
a
^
y
 
      G3 = 2 p
a
^
z
 
.
(21)

It is convenient to define h º m2/a2. Also we will denote by q the unit charge. In these terms, Eq. 10 becomes for this case:

W
N
= q2
8 pe0 a
(T1+T2 ),
(22)
where
T1 º 1
p

å
m1,m2,m3 ¹ 0 
2 (1-eip(m1+m2+m3)) e-[(4p2)/( m2)](m12+m22+m32)
m12+m22+m32
- 2 m
Ö
p
 
(23)
and
T2 º
å
n1,n2,n3 ¹ 0 
2 erfc( m
2
  ___________
Ön12+n22+n32
 
)

  ___________
Ön12+n22+n32
 
-
å
n1,n2,n3 
2erfc( m
2
  æ
 ú
Ö

( 1
2
+n1)2+( 1
2
+n2)2+( 1
2
+n3)2
 
)


  æ
 ú
Ö

( 1
2
+n1)2+( 1
2
+n2)2+( 1
2
+n3)2
 
.
(24)

In this expression, the sum over a and b has a total of 4 contributions - 2 pairs of identical contributions. For Na-Na or Cl-Cl interactions, tab = 0 and qa = qb resulting in repulsive contributions. For Na-Cl or Cl-Na interations, tab = ±a/2([^(x)]+[^(y)]+[^(z)]) and qa = qb resulting in attractive contributions. This expression can be evaluated using Maple, which gives the result

W
N
= q2
8 pe0 a
(-4.071).
(25)

NaCl structure

There are two kinds of sites - tNa = 0 and tCl = a/2 [^(x)].

T1 = a
2
( ^
x
 
+ ^
y
 
)    T2 = a
2
( ^
y
 
+ ^
z
 
)      T3 = a
2
( ^
x
 
+ ^
z
 
).
(26)

G1 = 2 p
a
( ^
x
 
+ ^
y
 
- ^
z
 
)     G2 = 2 p
a
(- ^
x
 
+ ^
y
 
+ ^
z
 
)     G3 = 2 p
a
( ^
x
 
- ^
y
 
+ ^
z
 
).
(27)


Footnotes:

1 Note that for any two lattice translations Ti and Tj, Ti-Tj = Tk, where Tk is also a lattice translation.


File translated from TEX by TTH, version 2.20.
On 1 Feb 2000, 12:15.