Linear algebra libraries(線形幾何学・行列計算ライブラリ)

Fortran, C, C++, etc
LAPACK (Linear Algebra PACKage)
ScaLAPACK (Scalable LAPACK)Intel Math Kernel Library (MKL)
Python: numpy.linalg, scipy.linalg
matrix.py
Product of matrixes
Inner productOuter productInverse matrix
Determinant
Eigen values/vectors
Solve simul. linear eqs.
LU decomposition
Cholesky decomposition
QR decomposition

AB
V1･V2

:C
=A@ B
: inner = numpy.dot(V1, V2)
inner = numpy.inner(V1, V2)
V1 × V2 : V3
= numpy.cross(V1, V2): Ai= numpy.linalg.inv(A): det= numpy.linalg.det(A): lA, vA = numpy.linalg.eig(A)AX = B : X= numpy.linalg.solve(A, B): P, L, U = numpy.linalg.lu(A)= numpy.linalg.cholesky(A)A=LLT : L
A=QR : Q, R
= scypy.linalg.qr(A)一般座標系 (general coordinate system)
一般座標/非直交系 (Non-Cartesian)

直交座標系 (Orthogonal)
デカルト座標系 (Cartesian)P r = xc,1e1+ xc,2e2

e2

xn,2

P r = xg,1a1+ xg,2a2

a2
e1

a1 xn,1

正規直交系 (orthonormal system)

𝒆𝒆𝒊𝒊 � 𝒆𝒆𝒋𝒋 = 𝛿𝛿𝑖𝑖𝑖𝑖
𝒆𝒆𝒊𝒊 = 1

一般座標系 (general coordinate system)

𝒂𝒂𝒊𝒊 � 𝒂𝒂𝒋𝒋 ≠ 𝛿𝛿𝑖𝑖𝑖𝑖

𝒆𝒆𝒊𝒊 , 𝒂𝒂𝒊𝒊 : 基底ベクトル (base vecor)

Cartesian – general coord. Conversion(直交系 － 一般座標系変換)

r = xc,1e1+ xc,2e2= xg,1a1+ xg,2a2
xc,1= xg,1 𝒂𝒂𝟏𝟏 � 𝒆𝒆𝟏𝟏 + xg,2 𝒂𝒂𝟐𝟐 � 𝒆𝒆𝟏𝟏
xc,2= xg,1 𝒂𝒂𝟏𝟏 � 𝒆𝒆𝟐𝟐 + xg,2 𝒂𝒂𝟐𝟐 � 𝒆𝒆𝟐𝟐

If 𝒂𝒂𝟏𝟏 = a11𝒆𝒆𝟏𝟏 + a12𝒆𝒆𝟐𝟐
𝒂𝒂𝟐𝟐 = a21𝒆𝒆𝟏𝟏 + a22𝒆𝒆𝟐𝟐
are given,
xc,1= 𝑥𝑥𝑔𝑔,1 𝑎𝑎11 + 𝑥𝑥𝑔𝑔,2 𝑎𝑎21
xc,2= 𝑥𝑥𝑔𝑔,1 𝑎𝑎12 + 𝑥𝑥𝑔𝑔,2 𝑎𝑎22

𝑎𝑎11
𝒂𝒂𝟏𝟏
𝒂𝒂𝟐𝟐 = 𝑎𝑎21

𝑥𝑥𝑐𝑐,1𝑎𝑎11
𝑥𝑥𝑐𝑐,2 = 𝑎𝑎12

𝑎𝑎12
𝑎𝑎22

𝑎𝑎21
𝑎𝑎22

𝒆𝒆𝟏𝟏
𝒆𝒆𝟐𝟐

𝑥𝑥𝑔𝑔,1
𝑥𝑥𝑔𝑔,2

Fractional coordinates in crystal(結晶の内部座標)

Lattice parameters: 𝑎𝑎, 𝑏𝑏, 𝑐𝑐 (= 𝑎𝑎1 , 𝑎𝑎2 , 𝑎𝑎3 ), 𝛼𝛼, 𝛽𝛽, 𝛾𝛾(= 𝛼𝛼23 , 𝛼𝛼13 , 𝛼𝛼12 )Lattice vectors: 𝒂𝒂𝟏𝟏 , 𝒂𝒂𝟐𝟐 , 𝒂𝒂𝟑𝟑 = 𝒂𝒂, 𝒃𝒃, 𝒄𝒄

r = xf,1a1+ xf,2a2 + xf,3a3 = xc,1e1+ xc,2e2 + xc,3e3(xf,1, xf,1, xf,3): Fractional coordinate (部分座標)
Internal coordinate (内部座標)
𝒂𝒂𝒊𝒊 = 𝑎𝑎𝑖𝑖
𝒂𝒂𝒊𝒊 � 𝒂𝒂𝒊𝒊 = 𝑎𝑎𝑖𝑖 𝑎𝑎𝑖𝑖 cos 𝛼𝛼𝑖𝑖𝑖𝑖 (𝑖𝑖 ≠ 𝑗𝑗)

𝑎𝑎11
𝒂𝒂𝟏𝟏
𝒂𝒂𝟐𝟐 = 𝑎𝑎21
𝑎𝑎31
𝒂𝒂𝟑𝟑

𝑎𝑎12
𝑎𝑎22
𝑎𝑎32

𝑎𝑎13
𝑎𝑎23
𝑎𝑎33

𝒆𝒆𝟏𝟏
𝒆𝒆𝟐𝟐
𝒆𝒆𝟑𝟑

Fractional coordinate to Cartesian coordinate

𝑥𝑥𝑐𝑐,1𝑎𝑎11
𝑥𝑥𝑐𝑐,2 = 𝑎𝑎12𝑥𝑥𝑐𝑐,3
𝑎𝑎13

𝑎𝑎21
𝑎𝑎22
𝑎𝑎23

𝑎𝑎31
𝑎𝑎32
𝑎𝑎33

𝑥𝑥𝑓𝑓,1𝑥𝑥𝑓𝑓,2𝑥𝑥𝑓𝑓,3𝑎𝑎11
𝒂𝒂𝟏𝟏
𝒂𝒂𝟐𝟐 = 𝑎𝑎21
𝑎𝑎31
𝒂𝒂𝟑𝟑

Conversion matrix

𝑎𝑎12
𝑎𝑎22
𝑎𝑎32

𝑎𝑎13
𝑎𝑎23
𝑎𝑎33

𝒂𝒂𝒊𝒊 = 𝑎𝑎𝑖𝑖
𝒂𝒂𝒊𝒊 � 𝒂𝒂𝒋𝒋 = cos 𝛼𝛼𝑖𝑖𝑖𝑖 (𝑖𝑖 ≠ 𝑗𝑗)

𝒆𝒆𝟏𝟏
𝒆𝒆𝟐𝟐
𝒆𝒆𝟑𝟑

tkcrystalbase.py
cal_lattice_vectors()𝒂𝒂𝟏𝟏
𝑎𝑎𝒂𝒂𝟐𝟐 = 𝑏𝑏 cos 𝛾𝛾
𝒂𝒂𝟑𝟑
𝑐𝑐 cos 𝛽𝛽𝑎𝑎, 𝑏𝑏, 𝑐𝑐 (= 𝑎𝑎1 , 𝑎𝑎2 , 𝑎𝑎3 )𝛼𝛼, 𝛽𝛽, 𝛾𝛾 (= 𝛼𝛼23 , 𝛼𝛼13 , 𝛼𝛼12 )0𝑏𝑏 sin 𝛾𝛾
𝑐𝑐 cos 𝛽𝛽 − 𝑐𝑐 cos 𝛽𝛽 cos 𝛾𝛾
𝑎𝑎33 =

00
𝑎𝑎33

𝑐𝑐 2 − 𝑎𝑎31 2 − 𝑎𝑎32 2

𝒆𝒆𝟏𝟏
𝒆𝒆𝟐𝟐
𝒆𝒆𝟑𝟑

Lattice properties

Unit cell volume
tkcrystalbase.cal_volume ()𝑉𝑉 = 𝒂𝒂𝟏𝟏 � 𝒂𝒂𝟐𝟐 × 𝒂𝒂𝟑𝟑
Distance rkl = rk – rl
tkcrystalbase.distance2() / .distance()
𝑟𝑟𝑘𝑘𝑘𝑘 2 = 𝒓𝒓𝒌𝒌𝒌𝒌 2 = ∑𝟐𝟐𝒊𝒊=𝟎𝟎 ∑𝟐𝟐𝒋𝒋=𝟎𝟎 𝒂𝒂𝒊𝒊 � 𝒂𝒂𝒋𝒋 𝑥𝑥𝑘𝑘𝑙𝑙,𝑖𝑖 𝑥𝑥𝑘𝑘𝑘𝑘,𝑗𝑗 = ∑𝑖𝑖,𝑗𝑗 𝑔𝑔𝑖𝑖𝑖𝑖 𝑥𝑥𝑘𝑘𝑘𝑘,𝑖𝑖 𝑥𝑥𝑘𝑘𝑘𝑘,𝑗𝑗
𝑔𝑔𝑖𝑖𝑖𝑖 = 𝒂𝒂𝒊𝒊 � 𝒂𝒂𝒋𝒋 : Metric tensor (計量テンソル)tkcrystalbase.cal_metrics()

Reciprocal lattice vectors tkcrystalbase.cal_reciprocal_lattice_vectors()𝒂𝒂∗ 𝟏𝟏 = 𝒂𝒂𝟐𝟐 × 𝒂𝒂𝟑𝟑 /V𝒂𝒂∗ 𝟐𝟐 = 𝒂𝒂𝟑𝟑 × 𝒂𝒂𝟏𝟏 /V𝒂𝒂∗ 𝟑𝟑 = 𝒂𝒂𝟏𝟏 × 𝒂𝒂𝟐𝟐 /V
Reciprocal vector at (h k l)
𝑮𝑮𝒉𝒉𝒉𝒉𝒉𝒉 = ℎ𝒂𝒂∗ 𝟏𝟏 + 𝑘𝑘𝒂𝒂∗ 𝟐𝟐 + 𝑙𝑙𝒂𝒂∗ 𝟑𝟑
Lattice space
𝑑𝑑ℎ𝑘𝑘𝑘𝑘 −2 = 𝑮𝑮𝒉𝒉𝒉𝒉𝒉𝒉 2 = ∑𝟑𝟑𝒊𝒊=𝟎𝟎 ∑𝟑𝟑𝒋𝒋=𝟎𝟎 𝒂𝒂∗ 𝒊𝒊 � 𝒂𝒂∗ 𝒋𝒋 ℎ𝑖𝑖 ℎ𝑗𝑗 = ∑𝒊𝒊,𝒋𝒋 𝑅𝑅𝑅𝑅𝑖𝑖𝑖𝑖 ℎ𝑖𝑖 ℎ𝑗𝑗
Bragg angleℎ, 𝑘𝑘, 𝑙𝑙 (= ℎ1 , ℎ2 , ℎ3 )
2𝑑𝑑ℎ𝑘𝑘𝑘𝑘 sin 𝜃𝜃 = 𝜆𝜆𝑅𝑅𝑔𝑔𝑖𝑖𝑖𝑖 = 𝒂𝒂∗ 𝒊𝒊 � 𝒂𝒂∗ 𝒋𝒋

どうやって逆格子を描くか• 実格子のフーリエ変換
逆格子の長さ: 実格子の逆数
角度: α* = 180o – α• 逆格子の基本ベクトルは実格子の基本ベクトルに直交する• ブラベー格子の消滅側を考える(基本格子の基本ベクトルを考えるのと同じ結果が得られる)逆格子の基本ベクトル

実格子の基本ベクトル

a*2

a2
a1

a*1

a*1 = 2πa2 × a⊥a1 ⋅ (a 2 × a ⊥ )Fractional – Cartesian conversion

python crystal_draw_cell.py

Rhombohedral celland reciprocal unit cellInter-atomic distances

crystal_distance.py

NaCl

Source code:# Lattice parameters (angstrom and degree)
lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)sites = [
#
原子量
電荷
描画条件
内部座標
[‘Na’, ‘Na1’, 11, 22.98997, +1.0, 0.7, ‘red‘, np.array([0.0, 0.0, 0.0])],[‘Na’, ‘Na2’, 11, 22.98997, +1.0, 0.7, ‘red‘, np.array([0.0, 0.5, 0.5])],[‘Na’, ‘Na3’, 11, 22.98997, +1.0, 0.7, ‘red‘, np.array([0.5, 0.0, 0.5])],[‘Na’, ‘Na4’, 11, 22.98997, +1.0, 0.7, ‘red‘, np.array([0.5, 0.5, 0.0])],[‘Cl‘, ’Cl1’, 17, 35.4527, -1.0, 1.4, ‘blue‘, np.array([0.5, 0.0, 0.0])],[‘Cl‘, ’Cl2’, 17, 35.4527, -1.0, 1.4, ‘blue‘, np.array([0.5, 0.5, 0.5])],[‘Cl‘, ’Cl3’, 17, 35.4527, -1.0, 1.4, ‘blue‘, np.array([0.0, 0.0, 0.5])],[‘Cl‘, ’Cl4’, 17, 35.4527, -1.0, 1.4, ‘blue‘, np.array([0.0, 0.5, 0.0])]]# Distance range
rmin = 0.1 # angstrom. 原子間距離がrmin未満の場合、同一の原子とみなす
rmax = 4.5 # angstrom. rmaxまでの原子間距離を計算

Inter-atomic distances

python crystal_distance.py
OUTPUT

NaCl

Lattice parameters: [5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
Lattice vectors:ax: ( 5.62,0,0) A
ay: ( 2.546e-10,
5.62,
0) A
az: ( 2.546e-10,0,5.62) A
Metric tensor:
gij: ( 31.58, 1.431e-09, 1.431e-09) A( 1.431e-09,
31.58, 6.48e-20) A( 1.431e-09, 6.48e-20,
31.58) A
Volume:177.5 A^3
Unit cell volume:177.5 A^3
Reciprocal lattice parameters: [0.17793594306049823, 0.17793594306049823, 0.17793594306049823, 90.00000000257246,
90.00000000516778, 90.00000000516778]
Reciprocal lattice vectors:
Rax: ( 0.1779, -8.06e-12, -8.06e-12) A^-1
Ray: (
0, 0.1779,0) A^-1
Raz: (
0,0, 0.1779) A^-1
Reciprocal lattice metric tensor:
Rgij: ( 0.03166, -1.422e-12, -1.422e-12) A^-1(-1.422e-12, 0.03166, 6.382e-23) A^-1(-1.422e-12, 6.382e-23, 0.03166) A^-1
Reciprocal unit cell volume: 0.005634 A^-3
nmax: 1 1 1Interatomic distances:
Cl1 ( 0.5,
0,(cut)
Na4 ( 0.5,
0.5,
Na4 ( 0.5,
0.5,
Na4 ( 0.5,
0.5,0) - Na4 (

0.5,

0.5,

0) + ( 0, -1, 0): dis =2.81 A0) - Na1 (0) - Na2 (0) - Na1 (

0,
0,
0,

0,
0.5,
0,

0) + ( 0, 1, 0): dis = 3.974 A
0.5) + ( 1, 0, -1): dis = 3.974 A0) + ( 1, 0, 0): dis = 3.974 A

三斜晶の面間隔 dhkl
2

−2

d hkl = G hkl = ha + kb + lc
*

*

* 2

1
1
2
2
2
=
+
+
+ 2 S12 hk + 2 S 23 kl + 2 S 31lh
S
h
S
k
S
l
11
22
33
2
2
d
V
S11 = b 2 c 2 sin 2 α(

S 22 = c 2 a 2 sin 2 β
S 33 = c 2 a 2 sin 2 γ

S12 = abc 2 (cos α cos β − cos γ )S 23 = a 2 bc(cos β cos γ − cos α )S 31 = ab 2c (cos γ cos α − cos β )V = abc 1 − cos 2 α − cos 2 β − cos 2 γ + 2 cos α cos β cos γ)python crystal_xrd.py

Bragg angles

NaCl

OUTPUT
Lattice parameters: [5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
Lattice vectors:ax: ( 5.62,0,0) A
ay: ( 2.546e-10,
5.62,
0) A
az: ( 2.546e-10,0,5.62) A
Metric tensor:
gij: ( 31.58, 1.431e-09, 1.431e-09) A( 1.431e-09,
31.58, 6.48e-20) A( 1.431e-09, 6.48e-20,
31.58) A
Volume:177.5 A^3
Unit cell volume:177.5 A^3
Reciprocal lattice parameters: [0.17793594306049823, 0.17793594306049823, 0.17793594306049823, 90.00000000257246,
90.00000000516778, 90.00000000516778]
Reciprocal lattice vectors:
Rax: ( 0.1779, -8.06e-12, -8.06e-12) A^-1
Ray: (
0, 0.1779,0) A^-1
Raz: (
0,0, 0.1779) A^-1
Reciprocal lattice metric tensor:
Rgij: ( 0.03166, -1.422e-12, -1.422e-12) A^-1(-1.422e-12, 0.03166, 6.382e-23) A^-1(-1.422e-12, 6.382e-23, 0.03166) A^-1
Reciprocal unit cell volume: 0.005634 A^-3
hkl range: 7 7 7
Diffraction angle, d, h, k, l:
2Q=
15.75 d=5.62 ( -1 0 0)
2Q=
15.75 d=5.62 ( 0 -1 0)
(cut)
2Q=
22.35 d= 3.97394 ( -1 -1 0)
2Q=
22.35 d= 3.97394 ( -1 0 -1)
2Q=
22.35 d= 3.97394 ( 1 0 1)

六方/三方格子 － 底心斜方格子変換
a'i = Σj (tij aj)(A' = TA)𝒂𝒂𝟏𝟏,𝐨𝐨
1 0 0 𝒂𝒂𝟏𝟏,𝐡𝐡
𝒂𝒂𝟐𝟐,𝐨𝐨 = 1 2 0 𝒂𝒂𝟐𝟐,𝐡𝐡
𝒂𝒂𝟑𝟑,𝐨𝐨0 0 1 𝒂𝒂𝟑𝟑,𝐡𝐡
R = Σi (xi ai ) = Σi (x'i a'i ) = Σi,j (x'i tij aj) から、
xi = Σ,j (x'j tji)
𝑥𝑥1,h
1𝑥𝑥2,h = 0𝑥𝑥3,h0x'i = Σ,j x'j (tji)-1𝑥𝑥1,o
1𝑥𝑥2,o = 0
𝑥𝑥3,o
0

a2,o
a2,h
yh

yo

P

O

xo xh

a1,h
a1,o

1
2
0

1
2
0

0
0
1

𝑥𝑥1,o
𝑥𝑥2,o
𝑥𝑥3,o

(X = tTX’)(X' = tT-1X)

0 −1 𝑥𝑥1,h
𝑥𝑥2,h
0
𝑥𝑥3,h
1

格子変換と行列実格子空間のベクトルの変換： (ai) => (a’i)
a'i = Σj (tij aj)( A' = TA )変換行列: (tij)R = Σi (xi ai ) = Σi (x'i a'i ) = Σi,j (x'i tij aj) から、
xi = Σ,j (x'j tji)
( X = tTX' )( X' = tT-1X )x'i = Σ,j x'j (tji)-1
R・Ghkl = hx + ky + lz はスカラー: R・Ghkl = R'・G'h'k'l'Σi (hixi) = Σi (h'ix'i) = Σi,j (hi x'j tji)h'j = Σi (hi tji)
( H' = TH )Ghkl = G'h'k'l' より、
Σi (h'ia*'i) = Σi (hia*i) = Σi,j (hi tji a*'i)
a*i = Σi,j (tji a*'j)( A* = tTA*' )a*'i = Σi,j ((tji)-1 a*j)
( A*' = tT-1A* )単位格子(ブラベー格子) と 基本格子
Siの構造 (室温)
空間群 Fd3m, No. 227 (立方晶系, ダイヤモンド構造)
ブラベー格子基本格子
aC = 0.5431 nm
aR = 0.3840 nm α=60oSia'i = Σj (tij aj)
𝒂𝒂𝟏𝟏,𝐩𝐩1/2
𝒂𝒂𝟐𝟐,𝐩𝐩 = 1/2
𝒂𝒂𝟑𝟑,𝐩𝐩
0(A' = TA)
1/2
0
1/2

0
1/2
1/2

𝒂𝒂𝟏𝟏,𝐅𝐅𝐅𝐅𝐅𝐅
𝒂𝒂𝟐𝟐,𝐅𝐅𝐅𝐅𝐅𝐅
𝒂𝒂𝟑𝟑,𝐅𝐅𝐅𝐅𝐅𝐅

R = Σi (xi ai ) = Σi (x'i a'i ) = Σi,j (x'i tij aj) から、
xi = Σ,j (x'j tji)

𝑥𝑥1,FCC
1/2
𝑥𝑥2,FCC = 1/2
𝑥𝑥3,FCC0x'i = Σ,j x'j (tji)-1

𝑥𝑥1,p
1/2
𝑥𝑥2,p = 1/2
𝑥𝑥3,p
0

1/2
0
1/2
1/2
0
1/2

(X = tTX’)𝑥𝑥1,p0
1/2 𝑥𝑥2,p1/2 𝑥𝑥3,p(X' = tT-1X)

0
1/2
1/2

−1

𝑥𝑥1,FCC
𝑥𝑥2,FCC
𝑥𝑥3,FCC

単位格子(ブラベー格子) と 基本格子
体心立方格子a'i = Σj (tij aj)
𝒂𝒂𝟏𝟏,𝐨𝐨−1/2
𝒂𝒂𝟐𝟐,𝐨𝐨 = 1/2
𝒂𝒂𝟑𝟑,𝐨𝐨−1(A' = TA)
1/2
−1/2
1/2

1/2
1/2
1/2

𝒂𝒂𝟏𝟏,𝐁𝐁𝐁𝐁𝐁𝐁
𝒂𝒂𝟐𝟐,𝐁𝐁𝐁𝐁𝐁𝐁
𝒂𝒂𝟑𝟑,𝐁𝐁𝐁𝐁𝐁𝐁

R = Σi (xi ai ) = Σi (x'i a'i ) = Σi,j (x'i tij aj) から、
xi = Σ,j (x'j tji)

𝑥𝑥1,BCC−1/2
𝑥𝑥2,BCC = 1/2
𝑥𝑥3,BCC−1x'i = Σ,j x'j (tji)-1

𝑥𝑥1,p
−1/2
𝑥𝑥2,p = 1/2
𝑥𝑥3,p
−1

(X = tTX’)
1/2 1/2 𝑥𝑥1,p−1/2 1/2 𝑥𝑥2,p1/2 1/2 𝑥𝑥3,p(X' = tT-1X)

1/2
−1/2
1/2

1/2
1/2
1/2

−1

𝑥𝑥1,BCC
𝑥𝑥2,BCC
𝑥𝑥3,BCC

六方格子－三方格子変換
菱面体晶: 六方格子軸と三方格子軸のどちらも取れる。
六方格子軸の基本ベクトル: a1(H), a2(H), a3(H)
三方格子軸の基本ベクトル: a1(R), a2(R), a3(R)a1(R) = (2a1(H) + a2(H) + a3(H)) / 3
a2(R) = (-a1(H) + a2(H) + a3(H)) / 3
a3(R) = (-a2(H) - 2a2(H) + a3(H)) /3
三方格子軸での逆格子座標 h k l
六方格子軸での逆格子座標 H K L
h = (2H + K + L) / 3
k = (-H + K + L) / 3
l = (-H - 2K + L) /3
H=h-k
K=k-l
L=h+k+l格子定数の関係
aR = √(3aH2 + c2)
sin(α/2) = 3/2 / √(3 + (c/aH)2)

a'i = Σj (tij aj)
𝒂𝒂𝟏𝟏,𝐡𝐡1
𝒂𝒂𝟐𝟐,𝐡𝐡 = 0
𝒂𝒂𝟑𝟑,𝐡𝐡1(A' = TA)

−1
1
1

0
−1
1

𝒂𝒂𝟏𝟏,𝐫𝐫
𝒂𝒂𝟐𝟐,𝐫𝐫
𝒂𝒂𝟑𝟑,𝐫𝐫Program: crystal_convert_cell.py
立方面心 => 基本格子変換
python crystal_convert_cell.py FCC FCCPrim

立方体心 => 基本格子変換
python crystal_convert_cell.py BCC BCCPrim

三方 => 六方格子変換
六方 => 直方格子変換
python crystal_convert_cell.py Rhomb RhombHex python crystal_convert_cell.py Hex HexOrthoMadelung potential
Sum of Coulomb potential in 3D is very slowly converging
Potential is proportional to r-1
Polarization potential due to +/- ions is to r -2
Number of ions on the sphere surface at radius r is to r2=> Contribution of ions from a surface region at r
to Coulomb sum is almost constant, independent of r

U ij (rij ) =Zi Z j e2 1
4πε 0 rij

+ U R ij (rij )Z 2e 2
1
U = ∑ U ij = − AM N A+UR
2 i≠ j4πε 0 R

Crystal structure
Rock salt type(NaCl)
CsCl type(CsCl)
Zinc blend (CuCl)
Wurzite(ZnO)
Cu2O type
Fluorite type (CaF2)

1
1
AM = ∑
Madelung constant2 i ≠ j rij / R

Ar
1.74761.76271.63801.6413
4.116
2.520

Madelung potential: Simple sum

python crystal_MP_simple.py

Coulomb sum in sphere with the radius r

Exact: -8.9 eVRock salt type

Efficient Coulomb sum: Evjen method
Sum up Coulomb potential inunits with zero net charge
Ion charges: Zi
On boundary plane : 1/2Zi
On boundary edge : 1/4Zi
On boundary corner : 1/8Zi

Madelung constant of Rock salt type structure∞1
nx + n y + nz()
1−
AM = −∑
2 nx ,n y ,nz = −∞ ,≠ ( 0, 0, 0 )

1
2

2

2

nx + n y + nz
1 1
1
1
1
1
AM = 6 × ×− 12 × ×+ 8× ×= 1.456
2
4
8 1+1+11
1+1

Madelung potential: Evjen method
Usage: python crystal_MP_Evjen.py ncell
ncell MP Madelung constant1-8.9766
1.7517691
2-8.95586
1.7477211
3-8.95521
1.7475955
4-8.95111.7475744
5-8.95508
1.7475686
6-8.95507
1.7475665
8-8.95506
1.7475652
10-8.95506
1.7475648
Exact (精確値)1.74756

Rock salt type

3D sum of Coulomb potential: Ewald method
Periodic calculation can be enhanced by FT?

Periodic positions of charge=> converted to the origin of FT data
But the charges are point charges=> converted to infinite in FT space=> Calculate for charges with finite width(拡がりのある電荷の周期配列として計算する)+

3D sum of Coulomb potential: Ewald method
The finite width charge distributions are converted by FT=> Take faster calculation parts in the real space and the reciprocal space
拡がった電荷のフーリエ変換を利用し、実空間和と逆空間和の計算の速い部分をとる
𝑍𝑍𝑗𝑗
𝑒𝑒 2
(𝐾𝐾𝐶𝐶 =)
Φ𝑖𝑖 = 𝐾𝐾𝐶𝐶 𝑍𝑍𝑖𝑖 �4𝜋𝜋𝜀𝜀0𝑟𝑟𝑖𝑖𝑖𝑖𝑗𝑗

Φ𝑖𝑖𝐼𝐼 = 𝐾𝐾𝐶𝐶 𝑍𝑍𝑖𝑖 � 𝑍𝑍𝑗𝑗
𝑗𝑗

erfc 𝛼𝛼 𝑟𝑟𝑖𝑖𝑖𝑖
𝑟𝑟𝑖𝑖𝑖𝑖

2
2 𝐆𝐆
𝑍𝑍
1
𝜋𝜋
𝑖𝑖
ℎ𝑘𝑘𝑘𝑘
�
Φ𝑖𝑖𝐼𝐼𝐼𝐼 = 𝐾𝐾𝐶𝐶
exp −
𝜋𝜋𝜋𝜋
𝐆𝐆ℎ𝑘𝑘𝑘𝑘 2
𝛼𝛼 2ℎ,𝑘𝑘,𝑙𝑙× cos 2𝜋𝜋𝐆𝐆ℎ𝑘𝑘𝑘𝑘 ⋅ 𝐫𝐫𝑖𝑖 ∑𝑗𝑗 𝑍𝑍𝑗𝑗 cos 2𝜋𝜋𝐆𝐆ℎ𝑘𝑘𝑘𝑘 ⋅ 𝐫𝐫𝑗𝑗 + sin 2𝜋𝜋𝐆𝐆ℎ𝑘𝑘𝑘𝑘 ⋅ 𝐫𝐫𝑖𝑖 ∑𝑗𝑗 𝑍𝑍𝑗𝑗 sin 2𝜋𝜋𝐆𝐆ℎ𝑘𝑘𝑘𝑘 ⋅ 𝐫𝐫𝑗𝑗
𝑮𝑮ℎ𝑘𝑘𝑘𝑘 ⋅ 𝒓𝒓𝒊𝒊 = ℎ𝑥𝑥𝑖𝑖 + 𝑘𝑘𝑦𝑦𝑖𝑖 + 𝑙𝑙𝑧𝑧𝑖𝑖

Φ𝑖𝑖𝐼𝐼𝐼𝐼𝐼𝐼 = 𝐾𝐾𝐶𝐶 𝑍𝑍𝑖𝑖

2𝛼𝛼𝑍𝑍𝑖𝑖
𝜋𝜋Φ𝑖𝑖 = Φ𝑖𝑖𝐼𝐼 + Φ𝑖𝑖𝐼𝐼𝐼𝐼 − Φ𝑖𝑖𝐼𝐼𝐼𝐼𝐼𝐼

Madelung potential: Ewald method
Usage: python crystal_MP_Ewald.py alpha prec
Alpha Precision
0.3
10-3
0.3
10-5
0.3
10-7
0.2
10-3
0.6
10-3
0.8
10-3
0.2
10-10
0.4
10-10
0.5
10-10
0.6
10-10
Exact (精確値)

MP
Madelung constant Range
Time (s)-8.95558
1.7476663 10.1/222 0.063 /222 0.016/0/0.016-8.95506
1.7475646 11.9/333 0.105 /222 0.031/0/0.031-8.95506
1.7475646 13.6/333 0.147 /333 0.047/0/0.047-8.95506
1.7475646 15.2/333 0.028 /111 0.042/0/0.042-8.95607
1.7477629 5.1/111 0.25 /333 0 /0.016 /0.016-8.95584
1.747718 3.8/111 0.45 /444 0 /0.016 /0.016-8.95506
1.7475646 24.3/555 0.093/222 0.16/0/0.16-8.95506
1.7475646 12.1/333 0.373/444 0.036/0.016/0.052-8.95506
1.7475646 9.7/222 0.58 /555 0.016/0.016/0.031-8.95506
1.7475646 8.1/222 0.84 /666 0.016/0.031/0.047
1.74756
Range: Rmax [Å]/nxmaxnymaxnzmax Gmax [Å-1]/hmaxkmaxlmax
Time: Real space sum / Reciprocal space sum / Total [s]Rock salt typeComparison: Evjen method
Rock salt type1
𝐴𝐴𝑀𝑀 = −2∞𝑛𝑛𝑥𝑥 ,𝑛𝑛𝑦𝑦 ,𝑛𝑛𝑧𝑧 =−∞,≠(0,0,0)

nx ny nzr
m Z S(mZ/r)
f
0
0
1
1
6 -1
-6 0.50
1
1 1.4142 12 1 8.48528 0.25
11
1 1.7321
8 -1 -4.6188 0.13

nx ny nz
0
0
1
0
1
1
1
1
1
0
0
2
0
1
2
0
2
2
1
1
2
1
2
2
2
2
2

r
1
1.4142
1.73212
2.23612.82842.44953
3.4641

m
6
12
8
6
24
12
24
24
8

�

−1 𝑛𝑛𝑥𝑥 +𝑛𝑛𝑦𝑦+𝑛𝑛𝑧𝑧

1𝑛𝑛𝑥𝑥 2 + 𝑛𝑛𝑦𝑦 2 + 𝑛𝑛𝑧𝑧 2

nx ny nzr
m Z S(mZ/ r)
S(mZf/r)
0
0
1
1
6 -1
-6
-3
0
1
1 1.4142 12
1 8.48528
2.12132034
1
1
11.73218-1-4.6188-0.5773503-2.13-1.456

Z S(mZ/r)
f-1-6
1
1 8.48528
1-1 -4.6188
1
1
3 0.5
-1 -10.733 0.5
1 4.24264 0.25
1 9.79796 0.5-1-8 0.2512.3094 0.13

S(mZf/r)-68.48528137-4.6188022
1.5-5.3665631
1.06066017
4.89897949-2
0.28867513-1.52-1.7518

0
0
0
1
1
2
0
0
0
0
1
1
1
2
2
3

0
1
2
1
2
2
0
1
2
3
1
2
3
2
3
3

2
2
2
2
2
2
3
3
3
3
3
3
3
3
3
3

2
2.23612.82842.44953
3.4641
3
3.16233.60564.24263.31663.74174.35894.12314.69045.1962

6
24
12
24
24
8
6
24
24
12
24
48
24
24
24
8

1
-1
1
1
-1
1
-1
1
-1
1
-1
1
-1
-1
1
-1

3
-10.733
4.242649.79796-82.3094-2
7.58947-6.65642.82843-7.2363
12.8285-5.506-5.82095.11682-1.5396-1.91

f

S(mZf/ r)

1
1
1
1
1
1
1
1
1
0.5
0.5
0.5
0.25
0.5
0.5
0.25
0.5
0.25
0.13

-6
8.485281374-4.61880215
3-10.7331263
4.242640687
9.797958971-8
2.309401077-1
3.794733192-3.32820118
0.707106781-3.61813613
6.414269806-1.3764944-2.9104275
1.279204298-0.19245009-1.7470

Exact value = 1.7476