More than 1 year has passed since last update.

PythonでMDことはじめ vol. 5

Last updated at 2022-10-16Posted at 2022-10-06

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)