Lennard-Jonesポテンシャルを用いたMDシミュレーション

力の計算

原子間相互作用ポテンシャル\(V\)の下で\(i\)番目の原子に加わる力\(\boldsymbol{F}_i\)は

$$ \boldsymbol{F}_i = -\nabla_i V = – \frac{\partial}{\partial \boldsymbol{r}_i }V $$

で与えられる。ペアポテンシャルの場合は

$$ \boldsymbol{F}_i = – \frac{\partial}{\partial \boldsymbol{r}_i }\sum_a \sum_{b>a} v(r_{ab}) = -\sum_a \sum_{b>a} \frac{\partial v(r_{ab})}{\partial \boldsymbol{r}_i} $$

ここで、\(a=i\)もしくは\(b=i\)の時だけ和の中身がゼロにならずに残るから

$$ \boldsymbol{F}_i = -\sum_{b>i} \frac{\partial v(r_{ib})}{\partial \boldsymbol{r}_i}-\sum_{i>a}\frac{\partial v(r_{ai})}{\partial \boldsymbol{r}_i} = -\sum_{j\neq i} \frac{\partial v(r_{ij})}{\partial \boldsymbol{r}_i} $$

ただしここで\(r_{ij}=r_{ji}\)の関係を用いた。さらに\(\sum\)の中身を変形すると

$$ \boldsymbol{F}_i = -\sum_{j\neq i} \frac{d v(r_{ij})}{dr_{ij}}\frac{\partial r_{ij}}{\partial \boldsymbol{r}_i} = -\sum_{j\neq i} \frac{d v(r_{ij})}{dr_{ij}}\frac{\boldsymbol{r}_i-\boldsymbol{r}_j}{r_{ij}}
$$

と書ける。「力はポテンシャルの微分」ということで、初心者はとかく\(dv(r_{ij})/dr_{ij}\)の計算に目が向きがちなようだが、そんなことはMathematicaにやってもらえばよろしい。クラスター展開した原子間ポテンシャルの微分では、\(\sum\)の範囲に注意を払う必要がある。

ところで

$$ \frac{\partial r_{ij}}{\partial \boldsymbol{r}_i} = \frac{\boldsymbol{r}_i-\boldsymbol{r}_j}{r_{ij}} $$

となるのはお分かりだろうか? 一応ここで導出しておこう。

$$ r_{ij} = \sqrt{ (r_{ix} – r_{jx})^2 + (r_{iy} – r_{jy})^2 + (r_{iz} – r_{jz})^2} $$

であるから、

$$ \begin{eqnarray} \frac{\partial r_{ij}}{\partial r_{ix}} & = & \frac{1}{2} {\left \{ (r_{ix} – r_{jx})^2 + (r_{iy} – r_{jy})^2 + (r_{iz} – r_{jz})^2 \right \} }^{-1/2} 2 (r_{ix} – r_{jx}) \\
& = & \frac{r_{ix} – r_{jx}}{r_{ij}} \end{eqnarray} $$

となる。同様に

$$ \frac{\partial r_{ij}}{\partial r_{iy}} = \frac{r_{iy} – r_{jy}}{r_{ij}} $$

$$ \frac{\partial r_{ij}}{\partial r_{iz}} = \frac{r_{iz} – r_{jz}}{r_{ij}} $$

なので、\(x,y,z\)成分をまとめて

$$ \frac{\partial r_{ij}}{\partial \boldsymbol{r}_i} = \frac{\boldsymbol{r}_i-\boldsymbol{r}_j}{r_{ij}} $$

と書けることがわかる。

Lennard-Jonesポテンシャルの場合の力は

$$ \boldsymbol{F}_i = -\sum_{j\neq i} 4 \varepsilon \left ( -12 \frac{\sigma^{12}}{r_{ij}
^{13}} + 6 \frac{\sigma^6}{r_{ij}^7} \right)\frac{\boldsymbol{r}_i-\boldsymbol{r}_j}{r_{ij}}
$$

となる。