1. 概要
さらに原子の種類を増やして、シミュレーションする原子数も増やしてみます。
可視化も行います。
2. パラメータ
こんな感じで増やしました。
質量の設定がなんかダサいですね。改良を考えましょう。
LJ_params = {
# epsilon, sigma, m
'Ne': [0.50e-21, 0.274e-9, 20.18 / NA * 1e-3],
'Ar': [1.67e-21, 0.340e-9, 39.948 / NA * 1e-3],
'Kr': [2.25e-21, 0.365e-9, 83.80 / NA * 1e-3],
'Xe': [3.20e-21, 0.398e-9, 131.30 / NA * 1e-3],
}
3. Atomクラス
運動エネルギーの計算をできるよう改良しました。
class Atom:
def __init__(self, kind, r, v):
self.kind = kind
self.r = r
self.r_before = r
self.v_half_after = np.array([0.0, 0.0, 0.0])
self.v_half_before = v
self.f = np.array([0.0, 0.0, 0.0])
self.m = LJ_params[kind][2]
def __str__(self):
return f'kind: {self.kind}\n\tr = {self.r}\n\tf = {self.f}'
def clear_force(self):
self.f = np.array([0.0, 0.0, 0.0])
def update_position_and_calc_uk(self, initial=False):
a = self.f / self.m
if initial:
self.v_half_after = self.v_half_before + 0.5 * a * dt
else:
self.v_half_after = self.v_half_before + a * dt
self.r = self.r_before + self.v_half_after * dt
# 更新前に運動エネルギーの計算
uk = self.calc_uk()
# 値の更新
self.v_half_before = self.v_half_after
self.r_before = self.r
return uk
def calc_uk(self):
v = np.linalg.norm((self.v_half_after + self.v_half_before) / 2)
return 0.5 * self.m * v ** 2
4. calc_potential
ポテンシャルエネルギーの計算をする関数も作りました。
def calc_potential(a1: Atom, a2: Atom) -> float:
disp = a2.r - a1.r
abs_r = np.linalg.norm(disp)
A = LJ_pair_params[a1.kind][a2.kind]['A']
B = LJ_pair_params[a1.kind][a2.kind]['B']
return A / abs_r ** 12 - B / abs_r ** 6
5. do_step
エネルギーの計算を入れました。他は前回と同じです。
def do_step(atom_list, initial=False):
up, uk = 0, 0
for a in atom_list:
a.clear_force()
for i in range(len(atom_list)):
a1 = atom_list[i]
for j in range(i+1, len(atom_list)):
a2 = atom_list[j]
f = calc_force(a1, a2)
a1.f -= f
a2.f += f
for a in atom_list:
uk += a.update_position_and_calc_uk(initial=initial)
for i in range(len(atom_list)):
a1 = atom_list[i]
for j in range(i+1, len(atom_list)):
a2 = atom_list[j]
up += calc_potential(a1, a2)
return up, uk
6. mainと結果
原子の種類、数を増やしました。
def main():
set_params()
a1 = Atom('Ar', np.array([0, 0, 0]), np.array([0, 0, 0]))
a2 = Atom('Ar', np.array([5e-10, 0, 0]), np.array([0, 0, 0]))
a3 = Atom('Ne', np.array([0, 5e-10, 0]), np.array([0, 0, 0]))
a4 = Atom('Xe', np.array([0, 0, 5e-10]), np.array([0, 0, 0]))
atom_list = [a1, a2, a3, a4]
do_step(atom_list, initial=True)
up_list = []
uk_list = []
total_list = []
for _ in range(10000):
up, uk = do_step(atom_list)
up_list.append(up)
uk_list.append(uk)
total_list.append(up + uk)
plt.plot(up_list, color='r', label='Up')
plt.plot(uk_list, color='b', label='Uk')
plt.plot(total_list, color='k', label='Total')
plt.legend()
plt.show()
出力
ちゃんと保存していますねえ。
7. GIFにしてみる
10000フレーム全て可視化していたら日が暮れるので、100ステップごとに位置を出力するようにします。
traj = []
for step in range(10000):
do_step(atom_list)
if step % 100 == 0:
data = [a.r for a in atom_list]
traj.append(data)
保存したデータをGIFにします。
こちらの二つの記事を参考にさせていただきました。
import io
from PIL import Image
images = []
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for data in traj:
ax.scatter(*np.array(data).T)
ax.set_xlim(-5e-10, 5e-10)
ax.set_ylim(-5e-10, 5e-10)
ax.set_zlim(-5e-10, 5e-10)
buf = io.BytesIO()
fig.savefig(buf, format="png", dpi=240) # 解像度は調整
plt.close(fig)
buf.seek(0) # バッファの先頭に移動
img = Image.open(buf).convert('RGB') # RGBAになるので変換
images.append(img)
ax.cla()
images[0].save('./test.gif', save_all=True, append_images=images[1:], optimize=False, duration=40, loop=0)
出力
お〜〜いい感じですね。
-1. 次回予告
次回はdriverクラスを実装して拡張性を高めます。
リンクはこちら
前回はこちらhttps://qiita.com/PlusF/items/cef88e56bda7191af444