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

前回の記事で実装したマーデルング定数のプログラムに以下課題があったので、作り直しました。

1. 実行速度
   原因：for文の使用

2. なぜかis_insideが全部Trueになる。→ outermostが小さいほど、誤差が大きくなる
   原因：分からず!!悔しい

   ※ outermost:　中心イオンと任意の近接イオンとの直線上に何個イオンがあるか

   ...説明力皆無すぎる

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

import numpy as np

outermost = int(input())

# 座標を作成
x, y, z = np.meshgrid(
    np.arange(-outermost, outermost + 1),
    np.arange(-outermost, outermost + 1),
    np.arange(-outermost, outermost + 1),
    indexing='ij'
)

# 中心(0, 0, 0)を除外
mask_center = (x == 0) & (y == 0) & (z == 0)

# 距離の計算
r = np.sqrt(x**2 + y**2 + z**2)
r[mask_center] = np.inf  # 中心は計算しない

# 符号
sign = np.where((x + y + z) % 2 == 0, -1, 1)

# 頂点の条件
is_vertex = (np.abs(x) == outermost) & (np.abs(y) == outermost) & (np.abs(z) == outermost)

# 辺の条件
is_edge = ((np.abs(x) == outermost) & (np.abs(y) == outermost)) | \
          ((np.abs(x) == outermost) & (np.abs(z) == outermost)) | \
          ((np.abs(y) == outermost) & (np.abs(z) == outermost))

# 面の条件
is_face = (np.abs(x) == outermost) | (np.abs(y) == outermost) | (np.abs(z) == outermost)

# 内部の条件
is_inside = (np.abs(x) < outermost) & (np.abs(y) < outermost) & (np.abs(z) < outermost)

# 原子の位置による寄与(正方形の四隅は原子の1/4だけが含まれている)
contribution = np.ones_like(r)
contribution[is_face] = 0.5 # 立方体の面
contribution[is_edge] = 0.25 # 立方体の辺
contribution[is_vertex] = 0.125 # 立方体の頂点

m = sign * contribution / r
m[mask_center] = 0

print(f"outermost = {outermost}")
print(f"output = {np.sum(m)}")

実行速度比較

前回の3Dマーデルング定数との実行速度の比較
比較のためのプログラムは以下

import time

outermost: int = int(input())
start_time = time.time()

            #... マーデルング定数のプログラム ...

end_time = time.time()
print(f"実行時間: {end_time - start_time:.6f} 秒")

outermost		50	100	200
for文(前回)[s]		1.034154	7.722279	63.227078
Numpy(今回)[s]		0.073490	0.366250	7.722279

近接イオンが少ないとき(outermostが小さいとき)の誤差の解消

outermost = 1
output = 1.4560299256299833

outermost = 2
output = 1.7517691333369407

outermost = 3
output = 1.7470415645202513

outermost = 10
output = 1.7475686037717981

outermost = 100
output = 1.7475645950341636

前回の課題は解決できた