初めに
今回、量子化学計算プログラムの作成に向けて、分子動力学での時間発展法のプログラム作成を行った。本稿では時間発展法として、陽的Euler法、4次の古典的Runge–Kutta法、LAI (Locally Analytic Integrator) 法を扱い、各々の手法のアルゴリズム、誤差の大きさ、計算コストの比較を扱う。
各手法の概要
数値解法の基礎
古典力学における Hamilton の運動方程式は、一般化座標 $\mathbf{q}={ q_{k} } $、一般化共役運動量 $\mathbf{p}={ p_{k}} $、Hamiltonian $H(\mathbf{q},\mathbf{p})$を用いて次式(\ref{new1})で表される。
$$
\dot{q_{k}} = \frac{\partial H}{\partial p_{k}}
$$
$$
\dot{p_{k}} = -\frac{\partial H}{\partial q_{k}}
$$
ここで、添字 $ k=1,2,\cdots ,n $ は運動の自由度を示す。保存系に対して質量荷重された直交座標系を用いた場合、$H(\mathbf{q},\mathbf{p})$ は次式の形になる。
$$
H(\mathbf{q},\mathbf{p}) = \frac{1}{2}|\mathbf{p}|^{2} + V(\mathbf{q})
$$
ここで、$V(\mathbf{q})$はポテンシャルであり、式は次式のようになる。
$$
\dot{ q_{k} } = p_{k} \nonumber
$$
$$
\dot{ p_{k} } = -\frac{ \partial V(\mathbf{q}) }{ \partial q_{k} }
$$
上式は $2n$ 個の連立常微分方程式である。一般には、式に解析的な解は存在しないため、数値積分法で近似解が求められる。
今回の取り組みでは、一次元ポテンシャル場での１個の粒子に対して、各々の時間発展法を適用して時間発展させ、運動量と一般化座標、エネルギーの変化を時間に対してプロットし、計算精度を手法ごとに比較する。計算精度の比較として、調和振動子ポテンシャルでは運動量・座標変位の値と解析解との誤差、Morse ポテンシャルではHamiltonianの収束値との誤差を基準とする。また、同じ時間発展法において、時間刻み幅を小さくすることで計算精度が向上するかの確認も行う。
一次元ポテンシャル $ V(\mathbf{q})$ の式は、調和振動子ポテンシャルでは式(\ref{osc})、Morseポテンシャルでは式(\ref{mor})の通りである。一般化座標を $q$、運動量を $p$ とすると、全ての時間発展法において、初めの時間 $t=0$ で調和振動子ポテンシャルでは質量座標 $q=0$, $p=1.0$、Morse ポテンシャルでは $q = 5.0$, $p = 1.0$ という同一の初期条件で数値積分計算を行う。
$$
V(q) =-\omega^{2} q
$$
$\omega$ : 角振動数。
$\omega=1.0$ とした。
$$
V(q) = D (1.0 – \mathrm{e}^{ -a(q – q_{e})} )^{2}
$$
$q_{e}$ : Morse ポテンシャルの極小点。$a$,$D$ : パラメータ。
$q_{e}=a=D=1.0$ とした。
Euler method
時間変数 $t$、初期時刻 $t_{0}$、未知関数 $y_{k}(t)$、初期値を $y_{0}$ を設定して、次式で与えられる常微分方程式を近似的に解くことを考える。
$$
\dot{y_{k}}(t) = f(t,y_{k}(t))
$$
$$
y_{k}(t_0) = y_0
$$
ある時刻 $t$ での勾配$f(t,y_{k})$は $t$ と $y(t)$ の関数として与えられる。初期時刻 $t_{0}$での初期値は $y_0$で与えられている。
時間刻み幅を $dt$ とすると、1ステップの時間区間[$t$,$t+dt$]であり、Euler法での解の近似計算式は次式のようになる。
$$
\dot{y_{k}}(t+dt) = y_{k}(t)+ f(t,y_{k}(t))dt
$$
"$y$"を 物理量"$p,q$"に置き換えると式は次式の形になる。
$$
q_{k}(t_0+dt) = q_{k}(t_0)+p_{k}(t_0)dt
$$
$$
p_{k}(t_0+dt) = p_{k}(t_0)-\frac{\partial V(q(t_0))}{\partial q_{k}(t_0)}dt
$$
式から分かるように、Euler法は各ステップあたりTaylor展開の打ち切り誤差として、座標に対して$O(dt^2)$のオーダーの誤差が蓄積するため、精度が低くなる。
Runge-Kutta method
Runge-Kutta法は、幾つかのEuler型のステップから得られる情報を使ってTaylor展開のある高次の次数まで一致させることで1つの区間での近似解を求める数値積分法である。今回は4次の古典的Runge-Kutta法を用いる。
2.1 節での常微分方程式(\ref{di})を用いると、1ステップの時間区間[$t$,$t+dt$]での4次Runge-Kutta法による解の近似計算式は次式(\ref{RK})のようになる。
$$
y_{k}(t_0+dt) = y_{k}(t_0)+\frac{k_1 + 2k_2+ 2k_3 +k_4}{6} dt
$$
$$
k_1 = f(t_0,y_{k}(t_0))
$$
$$
k_2 = f(t_{0}+\frac{1}{2}dt,y_{k}(t_0)+\frac{k_1}{2} dt )
$$
$$
k_3 = f(t_{0}+\frac{1}{2}dt,y_{k}(t_0)+\frac{k_2}{2} dt )
$$
$$
k_4 = f(t_{0}+dt,y_{k}(t_0)+k_3dt)
$$
4次のRunge-Kutta法の打ち切り誤差は$O(dt^{5})$であるが、時間区間$[a,b]$での累積誤差は(比較において区間の大きさ$b-a$を一定とすると)$\frac{b-a}{dt} \cdot O(dt^{5}) = O(dt^{4})$である。よって、数値積分においてEuler法より高い精度を持つ。
LAI method
$q(t)$ と $p(t)$ を $t_{0}$ において Taylor展開して次式を得る。
$$
q_{k}(t_0+dt) = q_{k}(t_0)+\dot{q_{k}}(t_0)dt + \frac{1}{2}\ddot{q_{k}}(t_0)dt^2 + \sum_{j=1}^{n} a^{(k)}_{j}(t_0) dt^{j+2}
$$
$$
p_{k}(t_0+dt) = p_{k}(t_0)+\dot{p_{k}}(t_0)dt + \sum_{j=1}^{n} b^{(k)}_{j}(t_0) dt^{j+1}
$$
ここで${a}^{(k)}{j}$と${b}^{(k)}{j}$は求める未知数として設定する。
さらに、式を式に代入して次式を得る。
$$
\dot{q_{k}}(t_0) = p_{k}(t_0)
$$
$$
\ddot{q_k}(t_0) = \dot{p_{k}}(t_0)
= -\frac{\partial V(q)}{\partial q_{k}}|_{t=t_0}
$$
$$
(j+2){a}^{(k)}{j}({t}{0})={b}^{(k)}{j}({t}{0})
$$
式を式に代入すると、${b}^{(k)}_{j}$ が消去されたより具体的な次式が得られる。
$$
q_{k}(t_0+dt) = q_{k}(t_0)+p_{k}(t_0)dt -\frac{1}{2}\frac{\partial V(q(t_0))}{\partial q_{k}(t_0)} dt^2+ \sum_{j=1}^{n} a^{(k)}{j}(t_0) dt^{j+2}
$$
$$
p{k}(t_0+dt) = p_{k}(t_0)-\frac{\partial V(q(t_0))}{\partial q_k(t_0)}dt+ \sum_{j=1}^{n} (j+2) a^{(k)}_{j}(t_0) dt^{j+1}
$$
さらに式の両辺を$dt$で微分すると、次式が得られる。
$$
\sum_{j=1}^{n} (j+2)(j+1)a^{(k)}{j}(t_0) dt^{j}
=-\frac{\partial V(q(t{0}+dt))}{\partial q_{k}(t_{0}+dt)}+\frac{\partial V(q(t_0))}{\partial q_{k}(t_0)}
$$
ここで、1ステップ区間[$t_0$, $t_0+dt$]を $n$ 分割し、異なる $n$ 個の時刻 $t_{0} < t_{1} < \cdots < t_{n} = t_{0}+dt$ を設定する。この各時刻を代入すると次式が得られる。
$$
\sum_{j=1}^{n} (j+2)(j+1)a^{(k)}{j}(t_0)(t{i}-t_{0})^{j}=-\frac{\partial V(q(t_{i}))}{\partial q_{k}(t_{i})}+\frac{\partial V(q(t_0))}{\partial q_{k}(t_0)}
$$
$$
{ i = 1, 2, \cdots, n } \nonumber
$$
ベクトル {$a^{(k)}(t_0)$} を解とする $n$ 個の連立方程式と見なすことができる。この連立方程式をプログラムで解き、解として得た新しい{$a^{(k)}(t_0)$}を再び式(\ref{LAI9})に代入して連立方程式を解く。この{$a^{(k)}(t_0)$}の解を求める過程を{$a^{(k)}(t_0)$}の値が収束するまで行う。そして、収束した {$a^{(k)}(t_0)$}を式に代入して、区間$[t_{0},t_{0}+dt]$での1ステップ分の数値積分を行う。今回の試行において各 $a^{(k)}(t_0)$ の値の収束条件を $O(10^{-15})$ のオーダーの精度で行う。
元の関数のTaylor展開の高次項まで考慮するため、近似解の精度の大幅な向上が見込めるが、1ステップごとに膨大な計算を行うため、高い計算コストがかかることが予想される。
結果
調和振動子ポテンシャル
各図のグラフでは横軸を”経過時間” $t$、縦軸を”運動量,座標変位”$p,q$ としている。
グラフの凡例の見方の一例として "Runge p" を取り上げる。
"Runge" の部分においてそれぞれ、"Euler" は Euler法、"Runge" は Runge-Kutta法、"LAI" は LAI法で時間発展させた結果示す。
同様に "p" の部分においてそれぞれ、"p" は運動量、"q" は座標変位を表す。
また、時間刻み幅の比較やLAIでの $a_{j}$ の次数の表記で、"$dt=0.025$" は時間刻み幅 $dt$ が $dt=0.025$ を表し、"$j=5$" はLAIでの $a_{j}$ の次数が $j=5$ であることを示す。
"real p", "real q" は共に調和振動子ポテンシャルでの解析解であり、
"real p"は運動量 $p$ の解析解 $\sin(t)$、"real q"は座標 $q$ の解析解 $\cos(t)$ を表す。
調和振動子ポテンシャル条件において、全ての試行で時間刻み幅 $dt$ を $dt=0.025$ で固定している。
図1はRunge-Kutta法、図2はLAI法の座標・運動量変化を表す曲線をそれぞれEuler法の結果と重ね合わせたグラフである。図1,図2を見るに、どちらもEuler法と解析解とのずれが肉眼で確認できる。一方、Runge-Kutta法, LAI法の結果は解析解との一致が良く、精度を詳細に調べるには、データの解析が必要である。一般化座標 $q(t)$ の解析解と近似解との差をdatファイルで調べ、全離散点に対するRMS{二乗平均平方根 : Root Mean Square}を表1にまとめた。表1より、Euler法、Runge-Kutta法、LAI法でのそれぞれの誤差は、$O(10^{-2})$, $O(10^{-8})$, $O(10^{-16})$ のオーダーであることが分かった。これらの結果から、LAI法 $\gg$ Runge-Kutta法 $\gg$ Euler法の順で計算精度が向上することを確認できた。
また、図3の曲線を見るに、$j=2$ と $j=5$ のどちらの近似解も解析解との一致が良い。表1より、$j=2$ の結果は$O(10^{-8})$, $j=5$ の結果は$O(10^{-16})$ のオーダーの精度であった。つまり、$j=5$ の方が 解析解との一致が良く、$a_{j}$ の次数が高いほど計算精度が向上することを確認できた。
LAI $\gg$ Runge-Kutta $\gg$ Euler の順で精度が高いことが分かる。
Morseポテンシャル
各図のグラフでは横軸を"経過時間" $t$、縦軸を"Hamiltonian" Energy としている。
表記は3.1節と同様に従う。
図4はEuler法、図5はRunge-Kutta法、図6-7はLAI法のMorse振動子の全エネルギーを表す曲線グラフである。図5と図6はそれぞれEuler法の結果と重ね合わせている。
図4−6より、Runge-Kutta法、LAI法の Hamiltonian がEuler法の結果より収束している。さらに、表2の結果から LAI法 $\ll$ Runge-Kutta法 $\ll$ Euler法 の順で、Hamiltonian の真値との誤差が小さくなることを確認できた。さらに図4−6の中で図4が顕著に示しているが、方法の選択だけでなく、時間刻み幅を小さくするほど Hamiltonian の収束精度が向上することも確認できた。
これらの結果から、3つの方法に共通して時間刻み幅が小さいほうが Hamiltonian が良く保存されること、LAI法 $\gg$ Runge-Kutta法 $\gg$ Euler法の順で計算精度が向上することが確認できた。
図7の結果から $j=2$ と $j=5$ のどちらもEuler法, Runge-Kutta法より収束が早いことが分かるが、 $j=5$ の方が収束精度が高いことが分かる。さらに表2から、
$j=2$ が $O(10^{-11})$、$j=5$ が $O(10^{-16})$ のオーダーの誤差を確認できた。
よって、3.1節の調和振動子と同様、$a_{j}$ の次数が高いほど計算精度が向上していた。
考察
積分の概念からも分かるように、3.2節での図4-6 の結果から、全ての時間発展法において1ステップの時間刻み幅を小さくするほど近似解とエネルギー保存の精度が向上していることが示されている。
2.1節の式(8)から分かるように、Euler法は各ステップあたりTaylor展開の打ち切り誤差として、座標に対して$O(dt^2)$ のオーダーの誤差が蓄積するため、精度が低くなる。実際に3.1節から $O(dt^2)$ もの大きな誤差が存在することが分かる。したがって、分子動力学での時間発展でEuler法は用いられるべきでないことが分かる。
Runge-Kutta法は 参考文献[2]より、$O(dt^4)$ のオーダーの精度を持つとされており、今回の試行でも3.1節及び3.2節の結果から、この精度を持つことが示されている。さらに、LAIよりも計算コストが低くなることが今回の試行から分かる。精度はLAIほど高くはないが、試行条件に応じて時間刻み幅を十分小さくすることで精度を高めての使用もできる。
最後に、LAIの収束精度はRunge-Kutta法よりも遥かに高い計算精度を持つことが。LAI法の数値積分において、$a_{j}$ の次数が高いほど計算精度が向上することが3章の結果から分かる。具体的には3.1節では解析解との一致の精度が $j=2$, $j=5$ の場合にそれぞれ$O(10^{-8})$, $O(10^{-16})$ のオーダーであった。 そして、3.2節ではLAIでのエネルギー保存の精度が $j=2$, $j=5$ の場合にそれぞれ、$O(10^{-8})$, $O(10^{-15})$ のオーダーである。これから、さらに、$O(10^{-15})$ のオーダーの精度は $a_{j}$ の収束上限でもあり、 $j=5$ という比較的低い次数でも高い収束精度をもつことが分かる。この高い収束精度から、特にトラジェクトリー(軌跡)の局所的な微小部分でも解析解との良く一致するだろう。ただし、subroutine DGESV を使うと、時間の1ステップごとに各 $a_{j}$ が収束するように膨大な計算が行われるので、$a_{j}$ の次数が低くても他の2つの手法に比べて圧倒的に計算コストが高くなってしまう。この計算コストの高さから、計算精度のさらなる向上が計算機の性能によって大きな制限を受けてしまうことが考えられる。
結論
今回用いた全ての時間発展法に共通して、時間刻み幅を小さくすることで精度が向上することが確認できた。
手法の選択では、LAI法 $\gg$ Runge-Kutta法 $\gg$ Euler法の順に計算精度が高いことが示された。しかし、計算速度はEuler法 $\gg$ Runge-Kutta法 $\gg$ LAI法の順であった。よって LAI法 を用いた数値積分プログラムの実行には膨大な計算コストがかかり、計算機の性能による大きな制限がかかることも判明した。以上から、計算精度と計算コストがトレードオフの関係であることを考慮し、(勿論Euler法は含まれないが)異なる時間発展法を、各手法の特性、試行条件、目標に合わせて臨機応変に用いることが、計算コストを効果的に緩和し、研究を円滑に進める上で重要だと分かった。