PythonでMDことはじめ Vol. 3

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

# 2. Leap-frog法とは
例のごとく、[Wikipedia](https://en.wikipedia.org/wiki/Leapfrog_integration)を信じきっている私です。:pray:
質量$m$の質点に加わる力$F$が位置$r$によってのみ決まる変数だとします。
加速度$a$は
```math
a=\frac{F(r)}{m}
```
で決まります。あるステップ$i$の位置、速度、加速度を$r_i$, $v_i$, $a_i$と表すと、
```math
\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法です。この手法は時間可逆性、エネルギーの保存性が優れているそうです。:thinking:
他にもたくさんの数値積分法がありますが、エネルギーの保存性に優れているというのが魅力的だったのと、ぱっと見シンプルだったので採用しました。

# 3. えっ、初期値ってどうするの？
前の値がわかっていれば次の値を決められるのはわかったけど・・・初期値はどう与えるのが正解なの？？
と思ったあなたは頭がいいです。そうなんです。ここで妥協します。
数値積分の基本的な手法の一つに**Euler法**というものがあり、これを使って最初のステップだけ進めてあげる、というのが解決策です。初期値$r_0$, $v_0$は既知として、次のとおりになります。
```math
\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で実装
まずは数式の通りベタ打ちで書きます。後で綺麗にするので安心してください。:hugging:
原点にAr原子が固定されていて、もう一個のAr原子をある程度離れた場所に置いてあげるようなイメージです。

**コード**
```python
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()
```
**出力**
![Figure_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2702367/92de5c31-0321-5ceb-9752-3e27268c3e0c.png)
何やら振動していますね。前回のポテンシャルのグラフと照らし合わせてみましょう。
![Figure_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2702367/ec11f7a2-1597-a192-56be-3cef4d1421ed.png)
$r$の最大値と最小値に線を引いてみました。
初期値からスタートして、エネルギーが最小になるところに向かって進んで、ちょっと行き過ぎて、また帰ってくるという振動をしているのですね。

# 5. エネルギーの保存性
数値積分で気になってくるのはエネルギーの保存性です。
孤立系であれば、次の等式が成り立つはずです。
```math
U_p + U_k = const.
```
$U_p$がポテンシャルエネルギー、$U_k$が運動エネルギーです。
実際に計算して保存されているか見てみます。運動エネルギーの計算に$v_i$が必要になるので、
```math
v_i=\frac{v_{i-1/2}+v_{i+1/2}}{2}
```
とします。

**コード**
```python
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()
```
**出力**
![Figure_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2702367/5b81b1de-8057-1ff4-a682-eb7e2ee4202d.png)
面白い形をしていますが、合計すると真っ直ぐですね:hotdog:拡大すると合計の値も振動していますが、増えていったり減っていったりなどはしていません。エネルギー保存性は分子数が増えたときにまた確認します。

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