分子動力学法における温度制御の基礎
はじめに
分子動力学法(Molecular Dynamics method, MD)における温度制御の基礎の話を書きます。簡単のため、かなり雑な議論をします。細かい式の導出などは、2019年の金沢大学における集中講義のノートなどを参照してください。
分子動力学法における温度
ただし、
運動エネルギーは以下のように書けます。
さて、この系が温度
いま、原子が
のエネルギーが分配されます。分子動力学シミュレーションでは、全原子の速度ベクトルと位置ベクトルを追いかけるので、運動エネルギー
を得ます。これを、分子動力学シミュレーションにおける温度と定義します。
さて、分子動力学法では、
となります。
となります。要するに温度とは、原子あたりの平均運動エネルギーの2/3倍のことです。
なぜ温度制御が必要なのか
分子動力学法は、与えられたハミルトニアンから導出されるハミルトンの運動方程式を数値積分することで、原子の座標や速度を更新していく手法です。
上記の式から、ただちに
つまり、普通にシミュレーションをすると、全エネルギー
しかし、我々が知りたいのはエネルギー依存性ではなく、温度依存性であることが多いです。例えば、大気圧下で水が沸騰するのは原子あたりのエネルギーがXXジュールの時である、という情報は使いづらいです。温度100度の時に沸騰する、の方が知りたい情報です。そこで、原子の数
原子数が多ければ、
従って、温度
さて、
シミュレーションでは全エネルギーが保存するため、最終的に平衡状態では、全エネルギーの何割かが運動エネルギーに、残りがポテンシャルエネルギーに分配されます。そして、運動エネルギーと温度の関係は
とわかっているため、これで温度がわかります。問題は、全エネルギー
- 適当な全エネルギー
でシミュレーションをしてみるE - 平衡状態になった時の運動エネルギー
を求めるK - 運動エネルギー
からこの系の温度K を求めるT' - 温度
と目標温度T' を比べて、全エネルギーを調整するT - 目標温度
に一致するまで1-4を繰り返すT
という作業をする必要があります。不可能ではありませんが、極めて面倒です。そこで、「温度を指定してシミュレーションしたら、そのうち指定の温度に緩和して欲しい」という動機が生まれます。これが分子動力学法における温度制御です。
全エネルギーと温度制御
分子動力学シミュレーションにおいて、制御可能なのはエネルギーです。しかし我々は、全エネルギーの温度依存性
さて、エネルギーの温度依存性
「まともな系」であれば、一般に比熱は正になります。つまり、
- 温度を上げたければ全エネルギーを上げれば良い
- 温度を下げたければ全エネルギーを下げれば良い
ということになります。全エネルギーは運動エネルギーとポテンシャルエネルギーの和ですから、そのどちらかを増やしたり減らしたりすれば良いことになります。
さて、ポテンシャルエネルギーは位置の関数ですが、例えば相互作用する原子を引き離した時にエネルギーが増えるか減るかは相互作用の詳細に依存します。もし斥力が働いていれば引き離せばエネルギーは減るでしょうし、引力が働いていればエネルギーは増えます。Lennard-Jones原子などは、途中で引力と斥力が切り替わったりするのでとても面倒です。
一方、運動エネルギーの制御は簡単です。速度が速くなれば運動エネルギーは増え、速度が遅くなれば運動エネルギーが減ります。そこで、以下のように運動エネルギーを制御することにします。何か一定の定数
現在の温度を測定し、もし目標温度が低ければ
が満たされるように
さて、ある温度
その他の温度制御法
ベレンゼン法(Berendsen Thermostat)
速度スケーリング法は、「とりあえず運動エネルギーを指定温度の平衡状態の値に無理やり設定しておけば、そのうちポテンシャルエネルギーもその温度の値になるよね」という乱暴な方法でした。この方法だと、初期温度が指定温度から外れていた場合、かなり強烈な温度制御がかかることになります。すると、例えば衝撃波が発生したりして、シミュレーション的にうれしくありません。そこで、もう少しゆっくり温度制御しましょう、という発想で生まれたのがベレンゼン法(Berendsen Thermstat)です[2]。
Berendsen目標温度を
が成立するように温度を制御します(上記が成り立つように運動エネルギーを変化させます)。
能勢フーバー法(Nosé–Hoover thermostat)
ベレンゼンの方法は簡単ですが、系が小さい場合に正しい統計集団であるカノニカル分布を再現しない、という問題があります。そこで、系が小さい場合でも厳密にカノニカル分布を再現する手法として、能勢フーバー法が提案されました。もともと仮想時間を用いる能勢の方法というものがあり、それをフーバー(W. G. Hoover)が改良して運動方程式を簡単にしたものが能勢フーバー法です[3]。なお、能勢の英語表記を「Nose」にすると「ノーズ」になってしまうため、「Nosé」と表記します。
能勢フーバー法の運動方程式は以下で与えられます。
新たに
ランジュバン法(Langevin thermostat)
ベレンゼン法や能勢フーバー法は、決定論的な運動方程式に従っていましたが、確率的な運動方程式に従う方法で温度制御する方法があります。それがランジュバン法(Langevin thermostat)です[4]。
運動方程式はこんな感じになります。
ここで、
を満たすように決められます(ホワイトノイズ)。
ベレンゼン法や能勢フーバー法は全運動エネルギーだけを見て制御していたのに対し、ランジュバン法は、全ての自由度(例えば
温度制御法の使い分け
能勢フーバー法やランジュバン法は温度
一般に、ベレンゼン法や能勢フーバー法の方がランジュバン法よりも緩和が早いため、問題がなければそちらを使えば良いでしょう。ただし、ベレンゼン法や能勢フーバー法は「全運動エネルギー」しか見ないため、例えば系が相分離するなどして一様でない場合、定常状態において温度が一様でなくなることがあります。例えば系の右側が高温、左側が低温のまま定常になっているが、平均温度が指定温度になっている、なんてことが起きたりします。一方、ランジュバン法はすべての自由度に対して作用するため、そのようなことは起きません。不安なら、ランジュバン法と能勢フーバー法の両方で時間発展させ、同じ状態に落ち着くか確認すると良いでしょう。
また、温度の「落ちつき方」も違います。ベレンゼン法やランジュバン法は指定の温度に指数関数的に緩和していくのに対して、能勢フーバー法は振動しながら緩和します。また、定常状態に落ち着いたあとも、手法由来の振動がエネルギーに乗ってしまうため、たとえばエネルギーの時間揺らぎのスペクトルを測定すると、そこに手法由来のピークが現れてしまいます。時間に関してフーリエ変換する場合は注意したほうが良いでしょう。
まとめ
以上、駆け足で分子動力学法における温度制御法をまとめてみました。雑にまとめるとこんな感じです。
- 等分配則から、平衡状態における温度と運動エネルギーの関係がわかる
-
アンサンブルもNVE アンサンブルも対して変わらないので、温度NVT と全エネルギーT の関係さえわかっていれば、その温度におけるエネルギーE でE(T) アンサンブルで計算すれば良い。しかしNVE の関数形は事前にわからない。E(T) - 温度制御が必要なのは、その温度での平衡状態において運動エネルギー
とポテンシャルエネルギーK の比が事前にわからないからV - 運動エネルギーを増減させることで温度を上下させ、最終的に指定温度に緩和するようにする手法が温度制御法
- 温度制御法には決定論的な手法(速度スケーリング法、ベレンゼン法、能勢フーバー法)と、確率的な方法(ランジュバン法)がある
- ベレンゼン法や能勢フーバー法の方がランジュバン法より緩和が早い(ことが多い)
- ベレンゼン法とランジュバン法は温度が指数関数的に緩和するが、能勢フーバー法は振動しながら緩和する
- ベレンゼン法や能勢フーバー法は、温度が一様にならない場合がある
他にも能勢フーバー法でのエルゴード性の破れとか、それに対応する能勢フーバー法チェイン法とかいろいろマニアックな話題はありますが、それは例えば以下の記事参照。
- 調和振動子に作用させた熱浴のエルゴード性 調和振動子系をNose-Hooverで制御すると温度が正しく制御されない理由と回避方法
- 分子動力学法における熱浴の保存量 Nose-Hoover法や多変数熱浴における数値積分の保存量について
-
L. V. Woodcock, Chem. Phys. Lett., vol. 10, p. 257 (1971). ↩︎
-
H. J. C. Berendsen et al., J. Chem. Phys., vol. 81,p. 3684 (1984). ↩︎
-
S. Nosé, Mol. Phys. vol. 52 p. 255 (1984)., W. G. Hoover, Phys. Rev. A, vol. 31, p. 1695 (1985). ↩︎
-
W. G. Hoover, A. J. C. Ladd, and B. Moran, Phys. Rev. Lett. vol. 48, p. 1818 (1982). ↩︎
Discussion