PythonでMDことはじめ Vol. 3

1. 概要

数値積分の手法の一つである蛙飛び(Leap-frog)法をPythonで実装します。

2. Leap-frog法とは

例のごとく、Wikipediaを信じきっている私です。
質量$m$の質点に加わる力$F$が位置$r$によってのみ決まる変数だとします。
加速度$a$は

a=\frac{F(r)}{m}

で決まります。あるステップ$i$の位置、速度、加速度を$r_i$, $v_i$, $a_i$と表すと、

\begin{align}
r_i&=r_{i-1}+v_{i-1/2}\Delta t\\
a_i&=\frac{F(r_i)}{m}\\
v_{i+1/2}&=v_{i-1/2}+a_i\Delta t
\end{align}

という流れで値を更新していくのがleap-frog法です。この手法は時間可逆性、エネルギーの保存性が優れているそうです。
他にもたくさんの数値積分法がありますが、エネルギーの保存性に優れているというのが魅力的だったのと、ぱっと見シンプルだったので採用しました。

3. えっ、初期値ってどうするの？

前の値がわかっていれば次の値を決められるのはわかったけど・・・初期値はどう与えるのが正解なの？？
と思ったあなたは頭がいいです。そうなんです。ここで妥協します。
数値積分の基本的な手法の一つにEuler法というものがあり、これを使って最初のステップだけ進めてあげる、というのが解決策です。初期値$r_0$, $v_0$は既知として、次のとおりになります。

\begin{align}
a_0&=\frac{F(r_0)}{m}\\
v_{1/2}&=v_0+\frac{1}{2}a_0\Delta t\\
r_1&=r_0+v_{1/2}\Delta t\\
\end{align}

$v_{1/2}$だけはEuler法で求めてあげれば、あとはleap-frogで進められます。

4. Pythonで実装

まずは数式の通りベタ打ちで書きます。後で綺麗にするので安心してください。
原点にAr原子が固定されていて、もう一個のAr原子をある程度離れた場所に置いてあげるようなイメージです。

コード

NA = 6.022e23  # アボガドロ数
m = 39.948 / NA * 1e-3  # Arの原子量をアボガドロ数で割り、kgに単位を直す
dt = 1e-14  # 時間ステップ幅
r_0, v_0 = 5e-10, 0  # 初期値

# 最初の1ステップ
a = F(r_0) / m
v_half_after = v_0 + 0.5 * a * dt
r_after = r_0 + v_half_after * dt
# 値の更新
v_half_before = v_half_after
r_before = r_after

r_traj = []
for _ in range(1000):
    # ステップを進める
    r_after = r_before + v_half_before * dt
    a = F(r_after) / m
    v_half_after = v_half_before + a * dt
    # 値の更新
    r_before = r_after
    v_half_before = v_half_after

    r_traj.append(r_after)

plt.plot(r_traj)
plt.show()

出力

何やら振動していますね。前回のポテンシャルのグラフと照らし合わせてみましょう。

$r$の最大値と最小値に線を引いてみました。
初期値からスタートして、エネルギーが最小になるところに向かって進んで、ちょっと行き過ぎて、また帰ってくるという振動をしているのですね。

5. エネルギーの保存性

数値積分で気になってくるのはエネルギーの保存性です。
孤立系であれば、次の等式が成り立つはずです。

U_p + U_k = const.

$U_p$がポテンシャルエネルギー、$U_k$が運動エネルギーです。
実際に計算して保存されているか見てみます。運動エネルギーの計算に$v_i$が必要になるので、

v_i=\frac{v_{i-1/2}+v_{i+1/2}}{2}

とします。

コード

Up = []
Uk = []
Utotal = []
for _ in range(1000):
    # ステップを進める
    r_after = r_before + v_half_before * dt
    a = F(r_after) / m
    v_half_after = v_half_before + a * dt
    #エネルギーを計算
    up_tmp = U(r_after)
    v = (v_half_before + v_half_after) / 2
    uk_tmp = 0.5 * m * v ** 2
    Up.append(up_tmp)
    Uk.append(uk_tmp)
    Utotal.append(up_tmp + uk_tmp)
    # 値の更新
    r_before = r_after
    v_half_before = v_half_after

plt.plot(Up, color='r', label='Up')
plt.plot(Uk, color='b', label='Uk')
plt.plot(Utotal, color='k', label='Total')
plt.legend()
plt.show()

出力

面白い形をしていますが、合計すると真っ直ぐですね拡大すると合計の値も振動していますが、増えていったり減っていったりなどはしていません。エネルギー保存性は分子数が増えたときにまた確認します。

6. 次回予告

次回は分子数を増やしてみます。
リンクはこちらhttps://qiita.com/PlusF/items/cef88e56bda7191af444
前回はこちらhttps://qiita.com/PlusF/items/e356ce5ebcbad530c0d5