熱浴かサーモスタットか

はじめに

分子動力学法は、運動方程式を数値積分することで系を時間発展させ、種々の物理量を調べる手法です。普通にやるとエネルギーが保存するため、得られる統計集団はミクロカノニカル(NVEアンサンブル)となります。それでは不便なので運動方程式をいじって温度制御をする手法が多数考案されました。それらの温度制御法は「熱浴」もしくは「サーモスタット」と呼ばれます。ただ、どれを熱浴と呼ぶのか、サーモスタットと呼ぶのかは人によって異なり、統一されていないように見えます。また、同じ手法であっても熱浴と呼ばれたり、サーモスタットと呼ばれたりします。本稿では、この熱浴/サーモスタット問題について個人的な意見をまとめておこうと思います。

温度制御法については分子動力学法における温度制御の基礎も参考にしてください。

熱浴とサーモスタット

まず、熱浴(heatbath)の定義から入りましょう。興味のある系に比べて、十分に熱容量が大きく、温度がわかっている系を用意し、これを熱浴と呼びます。この熱浴と興味ある系を接触させ、エネルギーだけをやりとりできるようにします。すると、そのうち両者の温度が等しくなります。熱浴は十分に大きいので、温度は変わりません。従って、興味ある系は熱浴と同じ温度になります。これが熱浴のイメージです。

次にサーモスタット(thermostat)です。有名なのがバイメタル型のサーモスタットです。熱膨張率の異なる金属を貼り合わせておき、ある程度以上温度が高くなると曲がり、そうでない場合はまっすぐになるように調整しておきます。この金属と電気回路を組み合わせ、サーモスタットがまっすぐになるような温度であれば電気を流し続け、温度がある程度高くなったらサーモスタットが曲がって電気が切れることで温度を一定に保つ機構を実現します。こたつの中で「カチッ」という音が聞こえたら、おそらくバイメタル型のサーモスタットの音です。

雑にまとめると、

• 興味ある系より圧倒的に大きな内部自由度を持っており、接触することで勝手に温度が一定になるのが熱浴
• 興味ある系の温度を測定し、目標温度にむけてフィードバック制御をするのがサーモスタット

ということになります。

温度制御法の種類

分子動力学法は、普通にニュートンないしハミルトンの運動方程式を解くとエネルギーが保存します。そこで、運動方程式を修正して温度制御をします。温度制御の方法には、決定論的(deterministic)にする方法と確率論的(stochastic)にする方法の二種類があります。

決定論的の温度制御法は、速度スケーリング法、ガウシアンサーモスタット法、ベレンゼン法、能勢フーバー法などがあります。確率論的な温度制御法は、ランジュバン法、DPDなどがあります。

速度スケーリング法

速度スケーリング法は、温度制御法のなかでも最初期に提案されたものです。ある温度で平衡状態にある粒子系の速度はマクスウェル分布に従います。マクスウェル分布とは、要するに温度が高ければ高いほど分散が大きいガウス分布です。そこで、速度の分散が指定の温度のマクスウェル分布の分散になるように、無理やり速度を$\vec{v}_i \rightarrow \alpha \vec{v}_i$とスケーリングします。$\alpha$は、スケーリング後の速度の分散が、指定の温度のマクスウェル分布の分散に一致するように決めます。

こうして、速度の方を無理やり指定温度にしていれば、そのうちポテンシャルエネルギーの法も指定温度の平衡状態になるでしょ？というのが速度スケーリング法です。

ガウシアンサーモスタット法

速度スケーリング法は非常にナイーブな方法で、適当なステップごとに速度を無理やりスケーリングします。すると、運動方程式が閉じなくなります。運動方程式が閉じるとは、変数の時間微分がその変数自身の関数で書けている状態をいいます。これを自励的(autonomous)といいます。速度スケーリング法は自励的ではありません。

そこで、運動方程式に速度の分散が指定の温度のマクスウェル分布の分散になる、という制約を追加し、ラグランジュの未定定数法により、その制約を満たしながら運動するように運動方程式を変形したのがガウシアンサーモスタット法です。運動方程式に制約条件として「速度の分散」が常に指定の値になっているので、そのまま平衡状態になってくれれば、ポテンシャルエネルギーも指定温度の平衡状態になるでしょ？というのがガウシアンサーモスタット法です。

ベレンゼン法

速度スケーリング法も、ガウシアンサーモスタット法も、速度を無理やり指定の温度に合わせようとします。これは系が指定温度から離れた温度状態にあった場合、かなり強烈な外力を加えることになり、たとえばショックウェーブなどが発生してあまりよろしくありません。そこで、系の温度を測定し、その温度と指定温度の差によってゆっくり制御しましょう、というのがベレンゼン(Berendsen)法です。

ベレンゼン法は「どれくらいゆっくり温度制御をするか」を指定できるため、状況に応じて適切に制御できます。また、運動方程式は自励的です。しかし、平衡状態になった時、実現する統計集団が厳密なカノニカル分布にならない、という欠点があります。これは、正しいカノニカル分布に従う場合と、エネルギーの揺らぎ方が異なる、という意味で、エネルギーの期待値は正しい値に収束します。

また、欠点と言えるかどうかはわかりませんが、運動方程式が時間反転対称性を満たしません。

能勢フーバー法

能勢フーバー法(Nosé–Hoover)法は、これまでに述べた決定論的な温度制御法の欠点が全て克服された手法になっています。すなわち以下の性質を持っています。

• 運動方程式が自励的 (速度スケーリング法は自励的ではない)
• 運動方程式が時間反転対称性を持つ (ベレンゼン法はそうではない)
• 温度制御が急激にかからないように制御できる(速度スケーリング法やガウシアンサーモスタット法はそうできない)
• 実現する定常状態が厳密なカノニカル分布となる(ベレンゼン法はそうならない)

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

$$\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$という変数を導入し、目標温度$T_0$と現在の温度$T$との差を見ます。現在の温度が目標温度より高ければ($T>T_0$)、$\dot{\zeta}>0$なので、そのうち$\zeta>0$になります。すると、$\zeta$は摩擦として働き、系が冷えて温度が下がります。現在の温度が目標温度より低ければ逆のことが起きて、$\zeta$は負の摩擦、すなわち粒子は加速され、系が温められて温度が上がります。

ランジュバン法

ベレンゼン法や能勢フーバー法が決定論的な運動方程式に従っていたのに対して、確率的な運動方程式に従う方法で温度制御する方法がランジュバン法です。運動方程式は以下で与えられます。

$$\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$は摩擦力、$\vec{R}$はランダムな力で、例えば$\vec{R}$の$x$成分$R_x(t)$が

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

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

決定論的な温度制御法が全ての粒子を一度に制御していたのに対して、確率論的な温度制御法はそれぞれの粒子(自由度)に対して個別に温度制御をかけます。

この方程式は、液体の中に浮かぶコロイドの運動を模したもので、ランジュバン方程式と呼ばれます。コロイドの運動は液体による摩擦をうけてブレーキがかかります($-\gamma \vec{p}_i$の項)。一方、コロイドが静止していても、まわりの液体分子の衝突により突き動かされます($\vec{R}$の項)。このブレーキとアクセルが釣り合ったところで温度が決まります。

ランジュバン法も、定常状態が厳密なカノニカル分布になります。ただし、粒子個別に力積を与えるため、運動量が保存しません。運動量が保存しないため、流体運動を記述することができません。

DPD

DPD (Dissipative particle dynamics)法は、散逸粒子動力学法とも呼ばれ、ランジュバン法に対して運動量が保存するように改良した方法です。ランジュバン法では粒子一つに摩擦力と揺動力をかけていたのを、粒子のペアごとにかけることで、作用反作用の法則が成り立つようにして、全体の運動量が保存するようにしています。

熱浴か？サーモスタットか？

では、これらの温度制御法は熱浴でしょうか？それともサーモスタットでしょうか？

まず、熱浴のイメージとして「対象となる系に対して十分大きな自由度があり、それと接触することでそのうち同じ温度になる」というものがあります。数値計算において、対象となる系に比べて十分大きな系を用意するわけにはいかないため、何かしら近似をすることになります。

では、能勢フーバー法は熱浴でしょうか？サーモスタットでしょうか？

能勢フーバー法は、系に一つだけ仮想的な粒子を追加し、その自由度と全粒子を相互作用させています。その仮想的な粒子の質量は十分に重いですが、これを「熱浴」と呼んで良いかは微妙です。一方、運動方程式から、現在の温度と指定の温度の差を見て、温度が低ければ上がるように、高ければ下がるように制御をかけている、という解釈ができます。そういう意味で「サーモスタット」と呼ぶのに抵抗はありません。

逆に、ランジュバン法はどうでしょうか？もともとランジュバン方程式が液体の中にあるコロイド粒子の運動を記述するものであったことを考えると、「なにか莫大な自由度を持つ系の中に、興味ある系を沈めた」という感じがあります。そういう意味で文字通り「お風呂につかる」イメージで、ランジュバン法を熱浴と呼ぶのは良さそうです。一方、別に系の温度を測定してフィードバック制御をかけているわけではないので、サーモスタットと呼ぶのはどうかな、という気がします。

まとめ

温度制御法について概観し、熱浴なのかサーモスタットなのか考えてみました。あくまで個人的な意見として、能勢フーバー法はサーモスタット、ランジュバン法は熱浴と呼ぶのがしっくりくる気がしています。

しかし、ランジュバン法は、多くの人が「Langevin thermostat」という名前で呼んでいるため、論文とかに書くならthermostatとしたほうが良いでしょう。また、能勢フーバー法の提案者であるW. G. Hooverさんは「能勢フーバー法はheatbathである」と断言していました。提案者がそう言ってるならそれに従った方が良い気もします。

というわけで、別に何が正解というわけでもないのですが、学生さんから「この温度制御法は熱浴ですか？それともサーモスタットですか？」と聞かれたのをきっかけに、考えをまとめておきました。

なにかの参考になれば幸いです。

参考文献

• 速度スケーリング法: L. V. Woodcock, Chem. Phys. Lett., vol. 10, pp. 257 (1971).
• ガウシアンサーモスタット法： W. G. Hoover, A. J. C. Ladd, and B. Moran, Phys. Rev. Lett., vol. 48, pp. 1818 (1982); D. J. Evans, J. Chem. Phys, vol. 78, pp. 3297 (1983).
• ベレンゼン法： H. J. C. Berendsen et al., J. Chem. Phys., vol. 81, pp. 3684 (1984).
• 能勢フーバー法： W. G. Hoover, Phys. Rev. A, vol. 31, pp. 1695 (1985).