PythonによるNaCl型結晶のマーデルング定数の実装

授業の課題で計算面倒くさかったので、NaCl型結晶のマーデルング定数(Madelung constant)をPythonで実装してみました。

マーデルング定数とは

マーデルング定数
イオン結晶におけるあるイオン(原子)の周りにあるイオン達の静電エネルギーの和

$$M = \frac{d}{Z_i Z_j}\sum_{j \neq i} \frac{Z_i Z_j}{r_{ij}}$$

参考：　初歩から学ぶ固体物理学(KS物理専門書) p.89

何の変数かはまた今度更新します。

アルゴリズム

1. 一対のイオン-イオン間のマーデルング定数を計算する
2. 上の計算を全てのイオン-イオン間について、順番にやる
3. 全ての和を計算する

1次元マーデルング定数の実装

m: list[float] = []
outermost: int = int(input())
i: int = 0
x: int = 0
r: float = 0

for x in range(-outermost, outermost + 1):
   if x == 0:
      continue
   r = abs(x)
   if x % 2 == 0:
      r *= -1
   m.append(1 / r)


print(sum(m))

2次元マーデルング定数の実装

m: list[float] = []
outermost: int = int(input())
i: int = 0
x: float = 0
y: float = 0


for x in range(-outermost, outermost + 1):
   for y in range(-outermost, outermost + 1):
      if x == y == 0:
         continue
      r: float = (x**2 + y**2) ** 0.5
      if (x + y ) % 2 == 0:
         r *= -1
      if abs(x) == abs(y) == outermost:
         #最外殻の四隅
         contribution: float = 0.25
      if abs(x) != abs(y) and (abs(x) == outermost or abs(y) == outermost):
         #最外殻の辺
         contribution = 0.5
      if abs(x) != outermost and abs(y) != outermost:
         #内殻
         contribution = 1.0
      m.append(contribution / r)


print(sum(m))

result

 outermost = 50
 1.6155415661410888

3次元マーデルング定数の実装

m: list[float] = []     #1原子ごとのマーデリング定数
outermost: int = int(input())   #最外殻　
x: int = 0      #x座標
y: int = 0      #y座標
z: int = 0      #z座標
r: float = 0    #中心原子との距離
test = outermost ** 2

#頂点の条件
is_vertex = (
    abs(x) == outermost and
    abs(y) == outermost and
    abs(z) == outermost
)
#辺の条件
is_edge = (
    (abs(x) == outermost and abs(y) == outermost) or
    (abs(x) == outermost and abs(z) == outermost) or
    (abs(y) == outermost and abs(z) == outermost)
)
#面の条件
is_face = (
    (abs(x) == outermost) or
    (abs(y) == outermost) or
    (abs(z) == outermost)
)
#内部の条件
is_inside = (
    abs(x) < outermost and
    abs(y) < outermost and
    abs(z) < outermost
)

for x in range(-outermost, outermost + 1):
    for y in range(-outermost, outermost + 1):
        for z in range(-outermost, outermost + 1):
            if x == y == z == 0: 
                #中心は計算しない
                continue
            r = (x**2 + y**2 + z**2) ** 0.5
            if (x + y + z) % 2 == 0:
                r *= -1
            if is_face:
                #立方体の面
                contribution = 0.5
            elif is_edge:
                #立方体の辺
                contribution = 0.25
            elif is_vertex:
                #立方体の頂点
                contribution: float = 0.125
            elif is_inside:
                #内殻
                contribution = 1.0
            m.append(contribution / r)
print(sum(m))

result

outermost = 100
1.7418198158396654

一言

３次元マーデルング定数は結構計算時間長い。NaCl型結晶は対称性があるので、8割くらい計算量を削減できそう。
本当はもっと丁寧な説明とか図とか作って入れたいけど、続かない気がするのでとりあえず投稿します。