分子動力学法における温度制御の基礎

はじめに

分子動力学法(Molecular Dynamics method, MD)における温度制御の基礎の話を書きます。簡単のため、かなり雑な議論をします。細かい式の導出などは、2019年の金沢大学における集中講義のノートなどを参照してください。

分子動力学法における温度

$NN$原子から構成される三次元の系を考えます。原子に$11$から$NN$まで通し番号をつけ、
$ii$番目の原子を「$ii$原子」と呼ぶことにしましょう。$ii$原子の速度を$v⃗i\vec{v}_i$、位置を$r⃗i\vec{r}_i$としましょう。この系のハミルトニアン、すなわち全エネルギーが以下のように与えられています。

$$H=K({v⃗i})+V({r⃗i}) H = K(\{\vec{v}_i\}) + V(\{\vec{r}_i\})$$

ただし、$KK$は運動エネルギー、$VV$はポテンシャル(相互作用)エネルギーです。ここではポテンシャルエネルギーについての詳細は考えません。

運動エネルギーは以下のように書けます。

$$K=∑i12mv⃗i2 K = \sum_i \frac{1}{2}m \vec{v}_i^2$$

さて、この系が温度$TT$であるとしましょう。すると、等分配則により、自由度あたり$kBT/2k_B T/2$のエネルギーが分配されるのでした。ただし$kBk_B$はボルツマン定数です。

いま、原子が$NN$個あり、三次元空間なので、運動エネルギーは$3N3N$の自由度を持っています。従って、運動エネルギーには全体として

$$K=32NkBT K = \frac{3}{2} N k_B T$$

のエネルギーが分配されます。分子動力学シミュレーションでは、全原子の速度ベクトルと位置ベクトルを追いかけるので、運動エネルギー$KK$は簡単に計算することができます。そこで、先程の式を温度について解くと、

$$T=2K3NkB T= \frac{2K}{3N k_B}$$

を得ます。これを、分子動力学シミュレーションにおける温度と定義します。

さて、分子動力学法では、$kB=1k_B=1$とする単位系をとることが多いです。すると、

$$T=2K3N T= \frac{2K}{3N}$$

となります。$K/NK/N$は、原子あたりの運動エネルギーです。原子あたりの運動エネルギー(平均運動エネルギー)を$Kˉ\bar{K}$と書くと、最終的に

$$T=23Kˉ T= \frac{2}{3} \bar{K}$$

となります。要するに温度とは、原子あたりの平均運動エネルギーの2/3倍のことです。

なぜ温度制御が必要なのか

分子動力学法は、与えられたハミルトニアンから導出されるハミルトンの運動方程式を数値積分することで、原子の座標や速度を更新していく手法です。$ii$原子の運動量$p⃗i=mv⃗i\vec{p}_i = m\vec{v}_i$を用いると、ハミルトンの運動方程式は以下のように書けます。

$$dp⃗idt=−∂H∂r⃗idr⃗idt=∂H∂p⃗i \begin{aligned} \frac{d \vec{p}_i}{dt} &= -\frac{\partial H}{\partial \vec{r}_i}\\ \frac{d \vec{r}_i}{dt} &= \frac{\partial H}{\partial \vec{p}_i} \end{aligned}$$

上記の式から、ただちに$HH$の時間微分がゼロとなることがわかります。

$$dHdt=∑i(∂H∂r⃗idr⃗idt+∂H∂p⃗idp⃗idt)=∑i(∂H∂r⃗i∂H∂p⃗i−∂H∂p⃗i∂H∂r⃗i)=0 \begin{aligned} \frac{d H}{dt} &= \sum_i \left( \frac{\partial H}{\partial \vec{r}_i} \frac{d\vec{r}_i}{dt} + \frac{\partial H}{\partial \vec{p}_i} \frac{d\vec{p}_i}{dt} \right) \\ &= \sum_i \left( \frac{\partial H}{\partial \vec{r}_i} \frac{\partial H}{\partial \vec{p}_i} - \frac{\partial H}{\partial \vec{p}_i} \frac{\partial H}{\partial \vec{r}_i} \right)\\ &=0 \end{aligned}$$

つまり、普通にシミュレーションをすると、全エネルギー$EE$が保存します。通常、体積$VV$一定の条件でシミュレーションをするため、原子の数$NN$、体積$VV$、エネルギー$EE$が一定のシミュレーションをしていることになります。この結果得られる統計集団を$NVENVE$アンサンブルといいます。

しかし、我々が知りたいのはエネルギー依存性ではなく、温度依存性であることが多いです。例えば、大気圧下で水が沸騰するのは原子あたりのエネルギーがXXジュールの時である、という情報は使いづらいです。温度100度の時に沸騰する、の方が知りたい情報です。そこで、原子の数$NN$、体積$VV$、温度$TT$が一定のとなる統計集団を考えます。これを$NVTNVT$アンサンブルといいます。

原子数が多ければ、$NVENVE$アンサンブルと$NVTNVT$アンサンブルも本質的には同じものになります。従って、シミュレーションは$NVENVE$アンサンブルで行って、その時の運動エネルギーから温度$TT$を求めれば、そのシミュレーションで得られた観測量は、温度$TT$における値であるとみなすことができます。

従って、温度$TT$における平衡状態の運動エネルギー$KK$さえわかれば、$NVENVE$アンサンブルから$NVTNVT$アンサンブルの情報を引き出すことができます。

さて、$NVENVE$アンサンブルにおいて、我々が調整できるパラメータは原子数$NN$、体積$VV$、そして全エネルギー$EE$です。全エネルギーは、運動エネルギーとポテンシャルエネルギーの和でした。

$$E=K+V E = K + V$$

シミュレーションでは全エネルギーが保存するため、最終的に平衡状態では、全エネルギーの何割かが運動エネルギーに、残りがポテンシャルエネルギーに分配されます。そして、運動エネルギーと温度の関係は

$$K=32NkBT K = \frac{3}{2} N k_B T$$

とわかっているため、これで温度がわかります。問題は、全エネルギー$EE$が$KK$と$VV$にどのように分配されるかはシミュレーション前にはわからないということです。全エネルギー$EE$が与えられた時、それがどのくらい運動エネルギー$KK$に分配され、どのくらいがポテンシャルエネルギー$VV$に分配されるかは系の詳細に依存します。従って、温度$TT$における物理量(例えば圧力)の値を知りたければ

1. 適当な全エネルギー$EE$でシミュレーションをしてみる
2. 平衡状態になった時の運動エネルギー$KK$を求める
3. 運動エネルギー$KK$からこの系の温度$T′T'$を求める
4. 温度$T′T'$と目標温度$TT$を比べて、全エネルギーを調整する
5. 目標温度$TT$に一致するまで1-4を繰り返す

という作業をする必要があります。不可能ではありませんが、極めて面倒です。そこで、「温度を指定してシミュレーションしたら、そのうち指定の温度に緩和して欲しい」という動機が生まれます。これが分子動力学法における温度制御です。

全エネルギーと温度制御

分子動力学シミュレーションにおいて、制御可能なのはエネルギーです。しかし我々は、全エネルギーの温度依存性$E(T)E(T)$を知りません。その状態で目標温度$TT$になるようにエネルギー$EE$を増やしたり減らしたりすることになります。

さて、エネルギーの温度依存性$E(T)E(T)$を、温度$TT$で微分すると、この系の比熱が得られます。例えば体積一定の条件で微分すると定積比熱が得られます。

$$CV=(∂E∂T)V C_V = \left( \frac{\partial E}{\partial T}\right)_V$$

「まともな系」であれば、一般に比熱は正になります。つまり、$E(T)E(T)$は$TT$に関して増加関数です。従って、

• 温度を上げたければ全エネルギーを上げれば良い
• 温度を下げたければ全エネルギーを下げれば良い

ということになります。全エネルギーは運動エネルギーとポテンシャルエネルギーの和ですから、そのどちらかを増やしたり減らしたりすれば良いことになります。

さて、ポテンシャルエネルギーは位置の関数ですが、例えば相互作用する原子を引き離した時にエネルギーが増えるか減るかは相互作用の詳細に依存します。もし斥力が働いていれば引き離せばエネルギーは減るでしょうし、引力が働いていればエネルギーは増えます。Lennard-Jones原子などは、途中で引力と斥力が切り替わったりするのでとても面倒です。

一方、運動エネルギーの制御は簡単です。速度が速くなれば運動エネルギーは増え、速度が遅くなれば運動エネルギーが減ります。そこで、以下のように運動エネルギーを制御することにします。何か一定の定数$α\alpha$を考え、速度をスケールします。

$$v⃗i→αv⃗i \vec{v}_i \rightarrow \alpha \vec{v}_i$$

現在の温度を測定し、もし目標温度が低ければ$α>1\alpha > 1$に、温度が高ければ$α<1\alpha < 1$に設定すればそのうち温度が目標温度に近づくことになります。$α\alpha$の決め方ですが、「目標温度における運動エネルギーになるように決めてしまう」のが一番楽です。運動エネルギー$KK$は$α2K\alpha^2K$に変換されるため、目標温度を$TT$、現在の運動エネルギーを$KK$として、

$$T=2α2K3N T= \frac{2\alpha^2 K}{3N}$$

が満たされるように$α\alpha$を決めてやります。すると、常に運動エネルギーは「指定の温度における等分配則が成立した場合の運動エネルギー」となります。

さて、ある温度$TT$において、全エネルギー$E(T)E(T)$は、平衡状態において$K(T)K(T)$と$V(T)V(T)$に分配されます。逆に、今、全エネルギー$EE$とポテンシャルエネルギー$VV$が、温度$TT$における平衡状態における値でない時、上記の方法で「常に無理やり運動エネルギーだけ温度$TT$の平衡状態の値$K(T)K(T)$」を維持してやると、そのうち$VV$や$EE$も温度$TT$における平衡状態での値になるでしょ、というのが温度制御の基本です。この方法を速度スケーリング法(Velocity Scaling method)と呼びます^[1]。速度スケーリング法は、もっともプリミティブな温度制御法です。

その他の温度制御法

ベレンゼン法(Berendsen Thermostat)

速度スケーリング法は、「とりあえず運動エネルギーを指定温度の平衡状態の値に無理やり設定しておけば、そのうちポテンシャルエネルギーもその温度の値になるよね」という乱暴な方法でした。この方法だと、初期温度が指定温度から外れていた場合、かなり強烈な温度制御がかかることになります。すると、例えば衝撃波が発生したりして、シミュレーション的にうれしくありません。そこで、もう少しゆっくり温度制御しましょう、という発想で生まれたのがベレンゼン法(Berendsen Thermstat)です^[2]。

Berendsen目標温度を$T0T_0$、現在の温度を$TT$とした時、緩和の時定数$τ\tau$を使って

$$dTdt=T0−Tτ \frac{dT}{dt} = \frac{T_0 - T}{\tau}$$

が成立するように温度を制御します(上記が成り立つように運動エネルギーを変化させます)。$τ\tau$を大きく取ればとるほど、ゆっくり温度を制御することになります。ベレンゼン法は非常に簡単であり、安定しているために広く使われています。

能勢フーバー法(Nosé–Hoover thermostat)

ベレンゼンの方法は簡単ですが、系が小さい場合に正しい統計集団であるカノニカル分布を再現しない、という問題があります。そこで、系が小さい場合でも厳密にカノニカル分布を再現する手法として、能勢フーバー法が提案されました。もともと仮想時間を用いる能勢の方法というものがあり、それをフーバー(W. G. Hoover)が改良して運動方程式を簡単にしたものが能勢フーバー法です^[3]。なお、能勢の英語表記を「Nose」にすると「ノーズ」になってしまうため、「Nosé」と表記します。

能勢フーバー法の運動方程式は以下で与えられます。

$$dp⃗idt=−∂H∂r⃗i−p⃗iζdr⃗idt=∂H∂p⃗idζdt=1Q(T−T0) \begin{aligned} \frac{d \vec{p}_i}{dt} &= -\frac{\partial H}{\partial \vec{r}_i} - \vec{p}_i \zeta\\ \frac{d \vec{r}_i}{dt} &= \frac{\partial H}{\partial \vec{p}_i} \\ \frac{d \zeta}{dt} &= \frac{1}{Q}\left(T - T_0\right) \end{aligned}$$

新たに$ζ\zeta$という変数が導入されています。$QQ$は熱浴の質量と呼ばれ、緩和の時定数を決めます。ここでは証明しませんが、この運動方程式は目標温度$T0T_0$におけるカノニカル分布が厳密な定常状態となります。

ランジュバン法(Langevin thermostat)

ベレンゼン法や能勢フーバー法は、決定論的な運動方程式に従っていましたが、確率的な運動方程式に従う方法で温度制御する方法があります。それがランジュバン法(Langevin thermostat)です^[4]。

運動方程式はこんな感じになります。

$$dp⃗idt=−∂H∂r⃗i−γp⃗i+R⃗dr⃗idt=∂H∂p⃗i \begin{aligned} \frac{d \vec{p}_i}{dt} &= -\frac{\partial H}{\partial \vec{r}_i} - \gamma \vec{p}_i + \vec{R}\\ \frac{d \vec{r}_i}{dt} &= \frac{\partial H}{\partial \vec{p}_i} \\ \end{aligned}$$

ここで、$γ\gamma$は摩擦力、$R⃗\vec{R}$はランダムな力で、例えば$R⃗\vec{R}$の$xx$成分$Rx(t)R_x(t)$が

$$<Rx(t)Rx(t′)>=2Dδ(t−t′) \left< R_x(t) R_x(t')\right> = 2 D \delta(t-t')$$

を満たすように決められます(ホワイトノイズ)。$DD$は揺動力の強さであり、$DD$と$γ\gamma$の比が目標温度になるように決められます。

ベレンゼン法や能勢フーバー法は全運動エネルギーだけを見て制御していたのに対し、ランジュバン法は、全ての自由度(例えば$ii$原子の速度の$x,y,zx,y,z$成分)に対して個別に温度制御をするのが特徴です。こちらも、温度$TT$におけるカノニカル分布が厳密な定常状態となります。

温度制御法の使い分け

能勢フーバー法やランジュバン法は温度$TT$におけるカノニカル分布が厳密な定常状態になる一方、ベレンゼン法はそうなりません。しかし、その差異は極めて小さいため、あまり気にする必要はありません。

一般に、ベレンゼン法や能勢フーバー法の方がランジュバン法よりも緩和が早いため、問題がなければそちらを使えば良いでしょう。ただし、ベレンゼン法や能勢フーバー法は「全運動エネルギー」しか見ないため、例えば系が相分離するなどして一様でない場合、定常状態において温度が一様でなくなることがあります。例えば系の右側が高温、左側が低温のまま定常になっているが、平均温度が指定温度になっている、なんてことが起きたりします。一方、ランジュバン法はすべての自由度に対して作用するため、そのようなことは起きません。不安なら、ランジュバン法と能勢フーバー法の両方で時間発展させ、同じ状態に落ち着くか確認すると良いでしょう。

また、温度の「落ちつき方」も違います。ベレンゼン法やランジュバン法は指定の温度に指数関数的に緩和していくのに対して、能勢フーバー法は振動しながら緩和します。また、定常状態に落ち着いたあとも、手法由来の振動がエネルギーに乗ってしまうため、たとえばエネルギーの時間揺らぎのスペクトルを測定すると、そこに手法由来のピークが現れてしまいます。時間に関してフーリエ変換する場合は注意したほうが良いでしょう。

まとめ

以上、駆け足で分子動力学法における温度制御法をまとめてみました。雑にまとめるとこんな感じです。

• 等分配則から、平衡状態における温度と運動エネルギーの関係がわかる
• $NVENVE$アンサンブルも$NVTNVT$アンサンブルも対して変わらないので、温度$TT$と全エネルギー$EE$の関係さえわかっていれば、その温度におけるエネルギー$E(T)E(T)$で$NVENVE$アンサンブルで計算すれば良い。しかし$E(T)E(T)$の関数形は事前にわからない。
• 温度制御が必要なのは、その温度での平衡状態において運動エネルギー$KK$とポテンシャルエネルギー$VV$の比が事前にわからないから
• 運動エネルギーを増減させることで温度を上下させ、最終的に指定温度に緩和するようにする手法が温度制御法
• 温度制御法には決定論的な手法(速度スケーリング法、ベレンゼン法、能勢フーバー法)と、確率的な方法(ランジュバン法)がある
• ベレンゼン法や能勢フーバー法の方がランジュバン法より緩和が早い(ことが多い)
• ベレンゼン法とランジュバン法は温度が指数関数的に緩和するが、能勢フーバー法は振動しながら緩和する
• ベレンゼン法や能勢フーバー法は、温度が一様にならない場合がある

他にも能勢フーバー法でのエルゴード性の破れとか、それに対応する能勢フーバー法チェイン法とかいろいろマニアックな話題はありますが、それは例えば以下の記事参照。

• 調和振動子に作用させた熱浴のエルゴード性　調和振動子系をNose-Hooverで制御すると温度が正しく制御されない理由と回避方法
• 分子動力学法における熱浴の保存量 Nose-Hoover法や多変数熱浴における数値積分の保存量について

脚注

1. L. V. Woodcock, Chem. Phys. Lett., vol. 10, p. 257 (1971).

2. H. J. C. Berendsen et al., J. Chem. Phys., vol. 81,p. 3684 (1984).

3. S. Nosé, Mol. Phys. vol. 52 p. 255 (1984)., W. G. Hoover, Phys. Rev. A, vol. 31, p. 1695 (1985).

4. W. G. Hoover, A. J. C. Ladd, and B. Moran, Phys. Rev. Lett. vol. 48, p. 1818 (1982).