PythonでMDことはじめ vol. 4

# 1. 概要
前回まででMDの基礎はさらったと思うので、規模を大きくしていきましょう。
同時進行でコードも再利用性、保守性の高いものに変えていきます。

# 2. クラス化
分子(原子)に紐づいている情報は多いので、まずは分子（原子）のクラスを作りましょう。
とりあえず単原子分子だけ考えるのでAtomというクラス名にしましょう。
大きさのある分子を扱うのは剛体の力学を考えなければならず非常に難易度が高いです。挫折した経験があります。
```python
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 clear_force(self):
        self.f = np.array([0.0, 0.0, 0.0])

    def update_position(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
        # 値の更新
        self.v_half_before = self.v_half_after
        self.r_before = self.r
```

必要な関数は2つです。
- clear_force: 名前の通り、力をリセットします。複数の分子間で力を計算する際関さんが必要なため、事前にリセットしておきます。
- update_position: 名前の通り、位置を更新します。速度もちゃっかり更新してます。`initial`が`True`の時は前回解説したようにEuler法で更新を行います。

力の更新には他の分子の情報が必要です。外部の関数で計算してあげます。
```python
def calc_force(a1: Atom, a2:Atom) -> np.ndarray:
    disp = a2.r - a1.r
    abs_r = np.linalg.norm(disp)
    C = LJ_pair_params[a1.kind][a2.kind]['C']
    D = LJ_pair_params[a1.kind][a2.kind]['D']
    return (C / abs_r ** 14 - D / abs_r ** 8) * disp
```
`disp`はdisplacementの略です。`abs_r`は分子間距離です。
今回からちゃっかりベクトルにしているので、マイナーチェンジが行われています。
```math
\begin{align}
\boldsymbol{F}(\boldsymbol{r})&=\left(\frac{C}{|\boldsymbol{r}|^{13}}-\frac{D}{|\boldsymbol{r}|^{7}}\right)\frac{\boldsymbol{r}}{|\boldsymbol{r}|}\\
&=\left(\frac{C}{|\boldsymbol{r}|^{14}}-\frac{D}{|\boldsymbol{r}|^{8}}\right)\boldsymbol{r}
\end{align}
```
という簡略化です。冪乗で`np.power`の方が早いんでしょうか。検討しておきましょう。
検討した記事はこちら:point_right:https://qiita.com/PlusF/items/83bf1a1be6a82452d780

# 3. パラメータの設定
前節で`LJ_params`だとか`LJ_pair_params`だとかが出てきていました。こいつらの説明を加えておきます。
```python
dt = 1e-15

NA = 6.022e23  # アボガドロ数

LJ_params = {
    # epsilon, sigma, m
    'Ar': [1.67e-21, 3.40e-10, 39.948 / NA * 1e-3]
}

LJ_pair_params = {}

def set_params():
    for kind_1, (eps_1, sigma_1, _) in LJ_params.items():
        params_tmp = {}
        for kind_2, (eps_2, sigma_2, _) in LJ_params.items():
            eps = np.sqrt(eps_1 * eps_2)
            sigma = (sigma_1 + sigma_2) / 2
            params_tmp[kind_2] = {
                'A': 4 * eps * sigma ** 12,
                'B': 4 * eps * sigma ** 6,
                'C': 48 * eps * sigma ** 12,
                'D': 24 * eps * sigma ** 6
            }
        LJ_pair_params[kind_1] = params_tmp
```
`LJ_params`には辞書型で$\varepsilon$, $\sigma$, $m$が格納されています。これらは原子ごとの値です。`LJ_pair_params`には、ポテンシャルや力の計算で必要になる係数を計算して格納してあります。`LJ_params`に値を入れるだけで原子の種類を容易に増やすことができます。

# 4. `do_step`
ステップ一回の計算を行う`do_step`関数を紹介します。
```python
def do_step(atom_list, initial=False):
    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:
        a.update_position(initial=initial)
```
1. 力をリセットする
1. 各原子間で力を計算する
1. 力をもとに速度、位置の更新を行う

# 5. `main`
`set_params`の呼び出しを忘れないようにしましょう。
これで簡単に10000ステップのシミュレーションを行えました。
```python
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]))

    atom_list = [a1, a2]
    do_step(atom_list, initial=True)
    for _ in range(10000):
        do_step(atom_list)
```
エネルギーの保存を確認しておきましょう。
![Figure_6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2702367/29d4560c-b207-6a11-ef99-247a6f22dce8.png)
多少ブレがありますが、まあ良いでしょう。

# -1. 次回予告
次回はさらに原子を増やしていきます。
リンクはこちら:point_right:https://qiita.com/PlusF/items/64b7c8ecd81584a81f62
前回はこちら:point_right:https://qiita.com/PlusF/items/22ed49534b6fc6be174c