はじめに
分子動力学について勉強したことをまとめます！
前提知識はテイラー展開くらいで、誰でも読めるように書いていると思います。
コードはすべてFortranで書いていますが,簡単なアルゴリズムなので書き換えは容易です。
Fortranのメモリ順序が列順であることに注意してください。
また、私自身、分子動力学初学者なため誤っていることを書いてるかもしれません。
マサカリ大歓迎です。
背景
分子集団系における連立微分方程式を考えます。我々は多体問題を解析的に解くとこができませんから、時間を細かく刻んで数値的に解くほかありません。計算する際、如何に誤差を小さくするか、また時間可逆性を満足するかを考えた計算法が幾つも提案されています。ここでは、ベルレ法と速度ベルレ法について紹介したいと思います。
前提条件として、
\ddot{r}_i(t) =\frac{F_i(t)}{m_i}
と簡単な微分方程式で書けるモデルについてのみ考えます。（単原子分子やイオンなど）
ベルレ法
まず、$r_i(t±Δt)$を$t=t$周りでテイラー展開し、足し合わせます。
r_i(t+\Delta t)=r_i(t) +\Delta t \dot{r}_i(t) +\frac{\Delta t^2}{2!}\ddot{r}_i(t) + \frac{\Delta t^3}{3!}\dddot{r}_i(t) +\mathcal{O}(\Delta t^4) \\
r_i(t-\Delta t)=r_i(t) -\Delta t \dot{r}_i(t) +\frac{\Delta t^2}{2!}\ddot{r}_i(t) - \frac{\Delta t^3}{3!}\dddot{r}_i(t) +\mathcal{O}(\Delta t^4)
より、
r_i(t+\Delta t) +r_i(t-\Delta t) = 2r_i(t) +\Delta t^2 \ddot{r}_i(t)+\mathcal{O}(\Delta t^4)
となります。一点のみをテイラー展開すると誤差が$\mathcal{O}(\Delta t^3)$になりますが、上記のように変形することで誤差のオーダーを小さくすることができます。変形すると次のようになります。
r_i(t+\Delta t) = 2r_i(t) - r_i(t-\Delta t) + \Delta t^2\frac{F_i(t)}{m_i} + \mathcal{O}(\Delta t^4)
先ほどテイラー展開した２式を差し引いて$\dot{r}_i(t)$について解くと
\dot{r}_i(t) =\frac{r_i(t+\Delta t) -r_i(t -\Delta t)}{2\Delta t} + \mathcal{O}(\Delta t^2)
となり、運動エネルギーや力場を与えてやればポテンシャルエネルギーが計算可能になります。
ベルレ法のアルゴリズムを図式すると以下のようになります。初期値は$t=t$と1ステップ前の位置を与えます。
最初のサイクルでは位置が決まるとポテンシャルが定まり、そこから力を求まることを表しています。なお、力の計算については今回の記事で紹介しませんので悪しからず。
一見、問題無さそうなアルゴリズムに思えますが、ベルレ法には桁落ちする罠が潜んでいます。
各ステップごとに
r_i(t)-r_i(t-\Delta t)+\Delta t^2\frac{F_i(t)}{m_i}
を足していることになりますが、分子の動きを追うためには$\Delta t$のオーダーがfsとなってくるため、三項目が非常に小さくなります。そのため、大 ー 小となり桁落ちが発生します。以下サンプルコードです。
subroutine verlet
use MDPARAM,only : dt,MASS
use VARIABLES,only : R,ROLD,RNEW,V,F,T
implicit none
integer :: i,j
T=0.0d0
do i=1,n !粒子数
do j=1,3 !xyz
RNEW(j,i)=2.0d0*R(j,i)-ROLD(j,i)+dT**2*F(j,i)/MASS
V(j,i)=(RNEW(j,i)-ROLD(j,i))/(2.0d0*dt)
T=T+V(j,i)**2
enddo
enddo
T=T*MASS*0.5d0
return
end
subroutine NEXT
use VARIABLES,only : R,ROLD,RNEW
implicit none
ROLD=R
R=RNEW
return
end
速度ベルレ法 (VV法)
ベルレ法は桁落ちが起こったりと数値計算上あまり良いアルゴリズムということができません。
ここで紹介する速度ベルレ法が与える式は、ベルレ法の式と等価な式ですが、数値的には桁落ちを回避することができたり、メモリを節約できたりとメリットが沢山あります。
それでは早速式を求めていきます。ベルレ法の式から、
\begin{align}
r_i(t+\Delta t)&=r_i(t)+\frac{1}{2}r_i(t)+\frac{1}{2}r_i(t)-r_i(t-\Delta t)
+\frac{1}{2}r_i(t-2\Delta t)-\frac{1}{2}r_i(t-2\Delta t)+\Delta t^2\frac{F_i(t)}{m_i}+\mathcal{O}(\Delta t^4)\\
&=r_i(t)+\Delta t\frac{r_i(t)-r_i(t-2\Delta t)}{2\Delta t}+\frac{1}{2}\{2r_i(t-\Delta t)-r_i(t-2\Delta t)+\frac{\Delta t^2}{m_i}F_i(t-\Delta t)\}
-r_i(t-\Delta t)+\frac{1}{2}r_i(t-2\Delta t)+\frac{\Delta t^2}{2}\frac{F_i(t)}{m_i}+\frac{\Delta t^2}{2}\frac{F_i(t)}{m_i}+\mathcal{O}(\Delta t^4) \\
&= r_i(t)+\Delta t \Big\{ \dot{r}_i(t-\Delta t)+\frac{\Delta t}{2}\frac{F_i(t)+F_i(t-\Delta t)}{m_i} \Big\}+\frac{\Delta t^2}{2}\frac{F_i(t)}{m_i}+\mathcal{O}(\Delta t^4)
\end{align}
二行目の{}内はベルレ法の$r_i(t)$の式を代入しました。
ここで速度を
v_i(t)=v_i(t-\Delta t)+\frac{\Delta t}{2}\frac{F_i(t)+F_i(t-\Delta t)}{m_i}
で求めることとすると、
r_i(t+\Delta t) = r_i(t)+\Delta tv_i(t) +\frac{\Delta t^2}{2}\frac{F_i(t)}{m_i}+\mathcal{O}(\Delta t^4)
少しテクニカルな式変形でしたが、結果を見てみればベルレ法の式よりも自然な形に思えます。
ここで桁落ちが解消されたことを確認してみましょう。
nステップ後において位置の式は以下のようになります。
r_i(t+n\Delta t) =r_i(t) +\Delta t\sum_{k=1}^n v_i[t+(k-1)\Delta t]+\frac{\Delta t^2}{2m_i}\sum_{k=1}^n F_i[t+(k-1)\Delta t]
上記の式をそのまま実装しても問題ないのですが、まだ改良の余地があります。
速度の式をよく見ると、力を半分ずつ計算してから速度を求めています。この操作を俗称でそれぞれMOVEA、MOVEBと言い、メモリを節約することができます。
以下、速度ベルレ法のアルゴリズムを図示したものです。初期値は$t=t$での位置と速度を与えます。
MOVEAとMOVEBは次のように書くことができます。
subroutine MOVEA
use MDPARAM,only : WMASS,dt
use VARIABLES,only : R,V,F
implicit none
R=R+dt*V+(dt**2)*0.5d0*F/WMASS
V=V+0.5d0*dt*F/WMASS
return
end
subroutine MOVEB
use MDPARAM,only : MASS,dt
use VARIABLES,only : V,F,T
implicit none
integer :: i,j
T=0.0d0
do i=1,n !粒子数
do j=1,3 !xyz
V(j,i)=V(j,i)+0.5d0*dt*F(j,i)/MASS
T=T+V(j,i)**2
enddo
enddo
T=T*MASS*0.5d0
return
end
速度ベルレ法のメリット
- 桁落ちを防ぐことが出来る
- １ステップ前の情報のみ必要なためメモリを節約できる
- 速度を半分づつ動かすことでメモリを節約できる
おわりに
勘の良い人はもう気づいているかもしれませんが、図で示した動かし方とは別の方法も勿論存在します。気になる人は位置ベルレ法やleapfrogで調べてみてください。
次回はエワルド法について書く予定です。
誤りやコードの改善点等ありましたらコメントください。
参考文献
ALLEN and TILDESLEY『Computer Simulation of Liquids』
岡崎進・吉井範行『コンピューター・シミュレーションの基礎』