分子動力学法における温度制御の基礎
渡辺宙志
東京大 学情 報基 盤セ ンタ ー [email protected] 1 はじ め に 1.1 こ の文 章 につ いて 分子動力学法における温度制御法、特にNos´e–Hoover法による温度制御の意味、および使 用上の注意をまとめた。間違いなどがあれば連絡をしていただければ修整する。なお、本稿 は第三回分子シミュレーションスクールのテキストとして用意した文章に加筆修正したもの であ る。 1.2 温 度と は 何か 温度とはなんであろうか。夏は暑く冬は寒い、お湯は熱く氷は冷たいなど、我々は日常的 に温度を体感している。温度計や体温計で温度を測ることもできるし、温度とは至極あたり まえの概念のように思える。しかし、深く考えると温度とはそれほど自明の概念ではない。 我々は既に熱がエネルギーの一形態であり、温度が高い物体は、それだけ構成粒子が激しく運 動していることを知っている。では、氷がものすごい高スピードですれ違う系は高温の系と 解釈してよいのだろうか(図1)。また、水に浮かぶ微粒子の運動を考える。ミクロに見れば微 粒子は水分子と相互作用しているため、系が平衡状態にあれば微粒子と水分子は同じ温度で あるはずである。しかし、いま水分子を忘れ、微粒子の間のみに働く相互作用を考えたらど うなるであろうか。この系に温度は定義できるだろうか。また、その温度は水温と同じであ ろうか。本講義の主眼は分子動力学法における温度制御であるが、温度についてきちんと理 解せずに温度を制御するのは危険である。そこで、少し温度についておさらいしておこう。 二つの物体を接したときに、その間に熱のやり取りが無い場合、その二つの物体の温度は 等しいと定義する。これがいわゆる熱力学第零法則であり、この法則によって温度という概 念が定義できる。逆に、温度が異なる二つの物体が接すると、その間には熱が流れることに なる。熱はミクロに見ればエネルギーの一形態であり、温度が高い物質はそれだけ大きな分(a)
(b)
図1:左：宇宙空間で高速ですれ違う氷二つの系は、温度が高いといえるか？右：水に浮かぶ 微粒子 だけ に 注目 した系 の温度は水温と同じとするべき か？子運動エネルギーを持っている。従って、温度と分子運動にはなんらかの関係があるはずで ある 。実 際に 、温度Tの粒 子系における運動量分布は、マク スウェル分 布 f (px, py, pz)∝ exp ( −p2x+ p2y+ p2z 2mkBT ) (1) に従 うこ と が 知 られ ており、ここからいわゆるエネル ギー 等分 配則 、 * p2x 2m + = * p2y 2m + = * p2z 2m + = 1 2kBT (2) が導かれる。すなわち温度は運動量分布の二次のモーメントに比例し、その比例定数がボル ツマン定数kBである。実際に数値計算を行い、運動量のヒストグラムを取ったものが図2で ある 。フィッティングにより、この 系は温度T = 1.3を持つ ことが わ かる。 0 0.01 0.02 0.03 0.04 -5 -4 -3 -2 -1 0 1 2 3 4 5 Momentum 図2:粒子系における運動量分布。黒丸が数値計算データ、点線がガウス分布exp(−p2/(2kBT )) によるフィッティング。これにより温度T = 1.3のマクスウェル分布をしていることがわかる。 さて、図2のシミュレーションは、温度制御をしていない分子動力学法の結果である。温 度制御をしていない、つまり普通の分子動力学法は、ハミルトンの運動方程式を用いる。ハ ミルトンの運動方程式系ではエネルギーが保存量となるため、時間発展の結果、統計集団と してミクロカノニカル分布を実現する。統計力学では、全エネルギー一定の系がミクロカノ ニカル、温度一定の系がカノニカル分布なのであった。ミクロカノニカル系の温度を測定す るとはどういうことであろうか？これらを理解するためには、温度やカノニカル分布の定義 を再確 認 す る必 要 があ る。 1.3 変 分 原理 と 分布 まずミクロカノニカル、カノニカルという二つの分布の違いを簡単に見るために、それぞ れを変分原理から導いておこう。位相空間Γ = (q1, q2,· · · , p1, p2,· · ·)を考える。この空間に ∫ f dΓ = 1 (3) を満たす分布関数f (Γ)を定義する。ただし、dΓ≡ dq1dq2· · · dp1dp2· · ·であり、積分は全空間 につ い て行 う。分布関数fが定ま ると、こ の系の物理量Aの期 待 値hAiは 、 hAi ≡ ∫ Af dΓ (4)
と定 義 され る 。さ て、この系のエントロピーSを次 のよ うに 定義 しよ う。 S≡ −kBhln fi = −kB ∫ f ln f dΓ (5) さて、この系のエネルギーを与える関数H(Γ)を定義したとき、エネルギーを一定にする条件 の下でエントロピーを最大化することを考える1。条件H = Eの元でδS = 0の変分を取ると、 f = 1 Ω Ω = ∫ δ(H− E)dΓ (6) を得る。つまり、分布関数はH = Eを満たす等エネルギー面でのみ0でない値をとり、かつ その 値は ど こ で も一 定である。こ れがミクロカノニカル 分布で あ る。 次に、エネルギーの期待値が一定であるという条件でエントロピーを最大化してみよう。 つま り、 hHi = U (7) の条件の下でエントロピーを最大化するのであるから、ラグランジュの未定乗数βを用いて I = βU− S (8) を定 義し 、その 変 分δI = 0を満 たすようにfを 決めると、 f = Z−1exp(−βH) Z = ∫ exp(−βH)dΓ (9) を得る。fは全位相空間で0でない値を持ち、その値はボルツマン重みexp(−βH)に比例す る。これがカノニカル分布である。エントロピーの定義式(5)にエネルギーの表式(7)、及びf の表 式(9)を用 い る と、SとUには dS = βkBdU (10) が成 り立 つ こ と がわ かる。こ れが熱力学関係式 dS = dU T (11) と等 価 であ る とい う要請 から、ラグランジュの 未定定数βが β = 1 kBT (12) のように温度と結びつくのであった。つまり、温度はカノニカル分布で定義されるもので あって、こ の枠 組 みでは ミ クロカノニカル分布では温 度は定 義 されな い 。 次に、温度Tを持つの運動量分布がマクスウェル分布(1)となることを示そう。先に導入した Hは、古典力学におけるハミルトニアンに対応する。簡単のため、一粒子系のハミルトニアン H = p 2 x+ p2y+ p2z 2m + V (qx, qy, qz) を考 え る。こ の と き、 px ∂H ∂px = p 2 x m (13) 1なぜこのエントロピーを選び、なぜこのエントロピーを最大化しようとするのかを考えて見よ。考えても 答 え は 出 な い で あ ろ う が 、きっと 統 計 力 学 の 知 見 が 深 ま る で あ ろ う。
とい う 量を 考 えよ う。この量 の平均は px ∂H ∂px = Z−1 ∫ px ∂H ∂px exp(−βH)dΓ = Z−1 ∫ px (−β) ∂ ∂px (exp(−βH)) dΓ = 1 β ∫ ∂px ∂pxexp(−βH)dΓ = 1 β と計算できる。途中で部分積分を用いた。pyやpzについてもまったく同様に計算できて、等 分配 則(2)と等 価 な 関係式 * p2x m + = * p2y m + = * p2z m + = 1 β = kBT (14) を得 る 。 ここで注意したいのは、平均h· · ·iはアンサンブル平均であり、同じ系を複数用意して同様 な観測を行った場合にその平均が温度を与えるという意味であって、多数の粒子を含む一つ の系の、それぞれの粒子が従う分布を意味してはいないことである。系がエルゴード的であ れば、一つの粒子について長時間観測することで同様な平均も得られるが、図2はある時刻 におけるスナップショットから得たデータである。ミクロカノニカル系の、しかもあるス ナップショットから温度が得られるという事実は局所平衡(Local Equilibrium)の概念を導入す ることで初めて理解される。つまり、ミクロカノニカル系は全体としてエネルギーを保存し ているが、そのうちの一つの粒子のエネルギーは他の粒子と相互作用することで揺らいでい る。注目している粒子のエネルギーが系全体のエネルギーに比べて無視できるほど小さいと き、粒子にとって他の粒子は熱浴として働くことになる。同様に他の粒子についても残りの 粒子が熱浴として働くことで、すべての粒子が同じ温度の熱浴と相互作用していると解釈す ることができ、その結果、粒子による平均がアンサンブル平均と同一視される。このとき、 その 熱浴 の 温 度 は、系 全体 の粒子の運動エネルギーの平 均から 定 まる。 1.4 温 度を 制 御す る前に 統計力学においては、温度は分布から定まる量である。分布がボルツマン分布をしている とき、運動エネルギーと温度に簡単な関係式が成立する。したがって、一般の場合に運動量 の分布の二次のモーメントが温度を与えるかどうかは非自明である。またp2/m = kBTの 場合 と まった く 同様 にして、 q∂H ∂q = kBT (15) も成立する。運動量から定まる温度を運動温度(Kinetic Temperature)、座標から定まる温度を 状態温度(Configuration Temeprature)と呼んで区別することもある。状態温度がビリアルその ものであることにも注意をしておきたい。もちろん、この二つの温度は平衡状態では一致す るべき量であるが、一般に運動温度の方が状態温度よりも緩和時間が短いため、単に運動温 度だけを観測して温度が平衡状態に至ったと判断すると、実はまだ状態温度は緩和の途中で 計算結果がおかしなことになる、ということも起きる。また温度勾配がある系やシェアがか かった系など、非平衡状態においてはこの二つの温度は一致する保証はない。いずれにせ よ、数値計算において温度を制御する前に、自分が観測しているのはどんな量か、それをど うい う 意味 で 温度 と呼ん でいるのかをしっか り認識し ておく 必 要があ る 。
2 分子 動 力 学 法 と温度 2.1 温度 制御の 必 要 性
分子動力学法(Molecular Dynamics method, MD)は、粒子間にかかる相互作用を力として表 現し、その時間発展を追うことで系の様々な物性を調べる手法である。支配方程式としてハ ミルトンの運動方程式を使うと全エネルギーが保存量となり、時間発展の結果、統計集団と してミクロカノニカル分布が得られる。最初に与えるエネルギーを変えて物理量を計算すれ ば、物理量の全エネルギー依存性がわかることになる。しかし、実際に知りたいのは物理量 の全エネルギー依存性ではなく、温度依存性であることが多い。たとえばたんぱく質の折り たたみを計算しようとしたら、体温付近でのシミュレーションをしたくなるであろう。そこ で、な ん らか の 温 度制御法 が必要となる。 ここで、温度制御は必ずしもカノニカル分布の実現手段である必要はないことに注意した い。自由度が大きければ、カノニカル分布における全エネルギーの揺らぎはほぼ無視できる ほど小さくなり、ミクロカノニカル分布と本質的な違いはなくなる。したがって、ある温度 Tにおけるエネルギーの期待値hHiが分かっていれば、ミクロカノニカル分布においても全 エネルギーがE =hHiとなるように設定することで温度Tのカノニカル分布と区別はつかな い結 果 を得 る こと ができ る。温 度とエネルギーの関係 は、比 熱Cを用い て C = ∂E ∂T (16) と表 さ れる の であった 。し たがって 、温 度Tにおける全系の エ ネルギ ーEは積 分 E = ∫ T 0 CdT (17) で与えられる。すなわち、比熱の温度依存性C(T )が分かっていれば任意の温度に対応する エネルギーが分かるため、温度制御の必要はない。さて、エネルギーの期待値は分配関数Z を用 い て hHi = −∂ ln Z ∂β = kBT 2∂ ln Z ∂T (18) で与えられる。式(16)と式(18)から、比熱の温度依存性がわかる、ということは系の分配関 数がわかっていることと等価である。系の分配関数が分かるということは、その問題が解け ているということであり、その系の情報はすべて分かっていることになる。もちろん分配関 数をあらかじめ求めることは困難であるため、比熱の温度依存性もシミュレーションを行う 前には分からない。そこで、最初は適当なエネルギー(温度)からシミュレーションをはじめ、 最終 的に 望 む 温 度に 収束するように温度制御をかける必 要が ある 。 2.2 温 度制 御 法 分子動力学法は運動方程式を数値積分する手法であるから、分子動力学法に温度制御を入 れるということは、温度が制御されるように適切に運動方程式を修正することを意味する。 得ら れた 新 し い 運動 方程式は、以 下のような性質を持ってい る ことが 望 ましい 。 カノ ニカ ル 分 布(Canonical Distribution) 時間発展の結果、指定温度のカノニカル分布を実 現す る 。 時間反 転 対 称 性(Time Reversiblity) 修 正された運動方程式 が時間 反 転対称 性 を保持 す る。 自励 性(Autonomous) 運動方程式が自己完結していること。たとえば外力など、時間に陽に 依存す る 項 が 運動方 程式に含まれないこ と。
エル ゴ ード 性(Ergodicity) 軌道が位相空間を埋め尽くす。等エネルギー面を覆うミクロカノ ニカルの場合と異なり、カノニカル分布には保存量が無いため、軌道が全位相空間を 覆う必 要 があ ること に注意。
計算コ ス ト(Efficiency) たとえ温度が適切に制御できる方法であっても、数値計算コスト(プ ログ ラ ム作 成 コス トも 含む)が高い方法は望ましく ない 。
分子動力学法での温度制御の初出は、おそらくL. V. WoodcockによるVelocity Scaling method であろう[1]。これは適当なステップごとに平均速度を指定温度に強制的にあわせることで、 平衡状態においてポテンシャルエネルギーも指定温度になるという方法であった。この手法 は 速 度 の 変 更 が 人 為 的 に 行 わ れ る た め 、運 動 方 程 式 が 自 励 的 で な い 。つ い で 、Hooverや Evansらが運動エネルギーを一定とする変分を取ることで、温度一定の運動方程式を導いて い る(Gaussian Thermostat法)[2]。こ の 方 法 で は 運 動 方 程 式 は 自 励 的 で あ る が 、運 動 エ ネ ル ギーは一定で、ポテンシャルエネルギーのみがボルツマン重みにしたがって揺らぐ。能勢 は 新 た な 自 由 度 を 付 け 加 え る こ と で 、運 動 エ ネ ル ギ ー も ボ ル ツ マ ン 重 み で 揺 ら ぐ 方 法 、
Extended System法を提案した[3]。この方法をHooverが変数変換し、仮想時間を扱わずにすむ
ようにしたのがNos´e–Hoover法である[4]。Nos´e–Hoover法の運動方程式は自励的かつ時間反 転対称であり、厳密にカノニカル分布が不変分布となる。かつ、ハミルトンダイナミクスに 一自由度のみを付け加える簡便さと実装の容易さから広く使われている。他には、速度を適 当な時間遅れをもってスケーリングするBerendsenの方法も収束が早くかつ安定であるため に、こちらも特に生物分野で広く使われているようである[5]。ただし、この方法では定常分 布がミクロカノニカル、カノニカル、どちらの分布にもならない。最近ではハミルトンダイ ナミクスに従いつつカノニカル分布を実現するNos´e–Poincar´e法という手法も提案された[6]。 Nos´e–Poincar´e法は適当なハミルトニアンを用意することで、系の時間発展はハミルトンの運 動方程式によって行うが、変数変換により興味ある系にカノニカル分布を実現する方法であ る。温度制御をかけつつシンプレクティック積分が適用できるため、こちらも広く使われて いる よ うで あ る。 2.3 温度 制御と カ ノニ カル分 布 Nos´e–Hoover法が広く使われている理由は実装が容易であるとともに、指定温度のカノニ カル分布が厳密に定常分布となることにある。以下、カノニカル分布がNos´e–Hoover法にお いて定常分布となる様子を見てみよう。簡単のため、一自由度系のハミルトニアンH(p, q)を 考え る 。この ハ ミ ルトニ アンにNos´e–Hoover法を適用した運 動 方程式 は ˙ q = ∂H ∂p ˙ p = −∂H ∂p − pζ ˙ ζ = 1 τ2 ( p∂H ∂p − 1 β ) (19) で与えられる。ここでβは逆温度β = 1/(kBT )であるが、後の計算の便利のためにβのままに しておく。τは温度制御の時定数を決め、小さいほど強く、大きいほど弱く温度制御がかか る。この系の位相空間は(p, q, ζ)で張られる三次元空間であり、分布関数もその空間において 定 義 さ れ る 。分 布 関 数 の 定 義 の 仕 方 は い く つ か あ る が 、こ こ で は 簡 単 に 系 の 状 態 が 点 (p, q, ζ)の周囲にある微小体積要素dΓ≡ dpdqdζに存在する確率がf (p, q, ζ)dΓと表現できると
き、fをこの系の分布関数と呼ぶことにしよう2。確率は保存量であるため、確率分布関数で あるfは 以下の 連 続 の式を 満たす。 ∂f ∂t = −divJ = −div ( ˙ Γf ) = −∂( ˙pf ) ∂p − ∂( ˙qf ) ∂q − ∂( ˙ζ)f ∂ζ ただしJは位相空間における確率流であり、速度場Γ˙と密度場fの積で与えられる。ここで、 fが定 常 分 布 であ れ ば時 間による偏微分の項はゼロとなる た め、fは ( ∂ ˙p ∂p+ ∂ ˙q ∂q + ∂ ˙ζ ∂ζ ) f + ∂f ∂pp +˙ ∂f ∂qq +˙ ∂f ∂ζ ˙ ζ = 0 を満 た す必 要 があ る。運動 方程式からp˙、q˙などを代入 する と、fに 関 する偏 微 分方程 式 ζf = ( ∂H ∂p − pζ ) ∂f ∂p + ∂H ∂q ∂f ∂q + 1 τ2 ( p∂H ∂p − 1 β ) ∂f ∂ζ を得 る 。この 定 常 解は、 f = Z−1exp [ −β ( H + ζ 2 2τ2 )] である3。さて、実際に興味のある系の分布関数は(p, q)に関する分布関数f 0である。そこでf のζに 関 する 自 由 度 を積 分することで消去すると4 f0 = ∫ ∞ −∞f dζ = Z −1 0 exp(−βH) と、興味ある系の分布関数f0(p, q)が指定温度βによるカノニカル分布となることがわかる。 つ ま り、Nos´e–Hoover法 は 、ま ず 興 味 の あ る 系(p, q)に 対 し て 自 由 度 を 一 つ 追 加 し た 世 界 (p, q, ζ)において定常分布f ∝ exp(−β(H + ζ2/(2τ2)))を実現する。追加自由度の部分がガウ ス分布となるために積分で消去することができ、結果として興味のある系に対してカノニカ ル分 布を 実 現 す る手法 である。 3 実際 の 使用 に お ける 注意 前節までで温度制御の必要性とその理論的な背景を説明してきた。しかし実際に使う場合 には数値積分法や熱浴をつけることによる影響など、事前に知っておいた方が良いことが多 い。以下、温度制御付きの分子動力学法を使う際の数値計算上の注意についてまとめる。 3.1 温 度制 御 と数 値積分 法 ハミルトンの運動方程式を数値積分する手法としては、シンプレクティック積分が広く使 われている。しかし、温度制御された系ではシンプレクティック積分が適用できない5。以下 2実際には、軌道は一次元であり、空間は三次元であるから、どんなに長い軌道をとったとしても軌道が空間 に占める体積は0であるため、この定義では系の状態の存在確率はゼロである。軌道の密度と分布関数を結び つ け る に は 、も う 少 し 精 緻 な 議 論 が 必 要 と な る 。 3代 入 し て 確 か め て み よ 4こ れ は (p, q, ζ)空 間 か ら(p, q)空 間 へ の 射 影 に 対 応 す る 。 5たまに市販のパッケージソフトを使っている人で、温度制御をしているにも関わらずシンプレクティック積 分をしている人がいる。おそらくプルダウンメニューから選んでいるのだろうが、温度制御した系にシンプレ ク ティック 積 分 を 適 用 す る の は ナ ン セ ン ス で あ る 。
ではまず、シンプレクティック積分とは何かを簡単に解説し、その後温度制御された系にお ける 数 値積 分 法に ついて 解説する。 簡単のため、一自由度系のハミルトニアンH = p2/2 + V (q)の時間発展を考えよう。ただ し質 量 を1とし てい る。ハ ミルトンの運動方程式は { ˙ p = −∂V∂q ˙ q = p (20) で与 え られ る 。こ れを、リュービル演算子Lを用 いて { ˙ p = iLp ˙ q = iLq (21) と書 こう6。リュービル 演算 子の表式は iL = −∂V ∂q ∂ ∂p+ p ∂ ∂q (22) であ る 。式(21)を形 式的に積分すると、 { p(t) = eitLp(0) q(t) = eitLq(0) (23) とな る。た だし 、 eitL= 1 + itL +(itL) 2 2 +· · · (24) である。p(0)やq(0)に演算したらp(t)やq(t)を得るのだから、U (t) = eiLtは時刻をtだけ進める 演算子、すなわち時間発展演算子であることがわかる。一般に時間発展演算子を厳密に計算 することは不可能であるため、適当な時間刻み∆tをとり、U (∆t)の近似表式U (∆t)˜ を使って { p(n∆t) = U (∆t)˜ np(0) q(n∆t) = U (∆t)˜ nq(0) (25) として時間発展を計算するのが数値積分である。U (∆t)˜ は∆tを小さくすればするほど元の演 算子U (∆t)に近づき、結果として数値積分の精度があがるが、その分、計算したい時刻まで のステップ数が多くなる。また、近似的な時間発展演算子を使っているため、ステップ数が 増えれば数値誤差も増える。そこで、数値誤差を押さえつつ高速な積分を行うため、数多く の数値積分法が提案されてきた。一般の数値積分ではルンゲ・クッタ法や予測子−修正子法 などが広く使われているが、ハミルトン系の時間発展を追う場合にはより性質のよいシンプ レク ティック積 分 が用 いら れることが多い。 一般のリュービル演算子を指数の肩に乗せた形式を計算することは困難だが、リュービル 演算子を適切に分解すると、展開が有限次で打ち切られ、厳密に計算できる場合がある。い ま、リュービル演算子Lを、運動に関する項LKとポテンシャルからの寄与LUに分けよう。具 体的 には iLK = p ∂ ∂q iLU = −∂V ∂q ∂ ∂p (26) 6ここで i≡√−1をつけているのは、演算子をエルミートにするためである。エルミート性を気にしなけれ ば つ け な く て も よ い 。そ の 場 合 は 演 算 子 の 表 式 か ら もiを 取 れ ば よ い 。
であ る 。ここ で 、iLKをqに二回演算してみると、 (iLK)2q = iLK ( p ∂ ∂qq ) = iLKp = p∂p ∂q = 0 であることがわかる。つまり(iLK)2 = 0である。同様な計算から(iLU)2 = 0であることもわ かる。演算子の二乗が零ということは、演算子を指数関数に乗せた場合の展開項が一次で切 れる とい う こ と なの で、 { eitLK = 1 + itLK eitLU = 1 + itLU (27)
が厳密に成り立つ。そこで、exp(iLt) = exp(i(LK+ LU)t)をexp(iLKt)とexp(iLUt)を使って
近似することを考える。演算子が指数の肩に乗ったものを分解する方法には鈴木−トロッ ター 分 解と 呼 ばれ る一 般論があるが、ここでは低次の場 合の 結果 のみ 示す と、
exp(i∆t(LK+ LU)) = exp(i∆tLK) exp(i∆tLU) + O(∆t) 一次近 似
exp(i∆t(LK+ LU)) = exp(i ∆t 2 LK) exp(i∆tLU) exp(i ∆t 2 LK) + O(∆t 2) 二次 近 似 となる7。このようにして得られた演算子により数値積分を行う手法をシンプレクティック積 分(Symplectic Integration)と呼ぶ8。シンプレクティックの名は、系のシンプレクティック性を 保存したまま離散化できることによる。その結果、もとのハミルトニアンの値Hが時間発展 において厳密に保存する。より正確にはn次のシンプレクティック積分では全エネルギーが O(∆tn)のオーダーで揺らぐが、一方的に増えたり減ったりせず、安定して長時間積分でき る9。ルンゲクッタ法や予測子−修正子法など、シンプレクティック性を持たない積分法で は、全系のエネルギーが保存せず、長時間積分すると徐々にエネルギーが増えていったり、 逆に減っていったりしてしまう。エネルギーの厳密な保存、これがシンプレクティック積分 の大き な 特 徴で あ る。 シンプレクティック積分の具体的な実行手続きを見てみよう。一次近似の場合、時間発展 演 算 子 はU (∆t) = exp(iLK˜ ∆t) exp(iLU∆t)で あ る 。こ れ は ま ずexp(iLU∆t)を 演 算 し 、次 に exp(iLK∆t)を演 算 する という 意味である。 ei∆tLK = 1 + i∆tLK= 1 + ∆tp ∂ ∂q ei∆tLU = 1 + i∆tLU = 1− ∆t∂V ∂q ∂ ∂p (28) であ る から 、具 体 的な積分 手順は p ← p −∂V ∂q∆t q ← q + p∆t 7指数分解は二次までは簡単だが、三次以上の高次公式はとたんに式が面倒になる。展開係数も1や1/2のよ う な 簡 単 な も の に な ら ず、項 数 も(指 数 関 数 を 展 開 し て い る の だ か ら 当 た り 前 だ が)指 数 関 数 的 に 増 え る 。 8よ り 正 確 に は 、こ こ で 示 し た 積 分 法 は シ ン プ レ ク ティック 積 分 の 陽 解 法 で あ る 。 9もっと 正 確 に 言 う な ら 、離 散 化 し た 時 間 発 展 演 算 子 の ヤ コ ビ ア ン の 行 列 式 が 厳 密 に1と な る 。
で与 え られ る 。二 次近似の 場合にも同様に、 q ← q + p∆t 2 p ← p −∂V ∂q∆t q ← q + p∆t 2
であることがわかる。なお、積分の次数としてはexp(i∆t2 LK) exp(i∆tLU) exp(i∆t2 LK)として もexp(i∆t2 LU) exp(i∆tLK) exp(i∆t2 LU)と し て も 同 じ 二 次 近 似 で あ る 。し か し 、前 者 は 一 ス テップの間に力の計算が一度であるのに対し、後者は二回行わなければならない。分子動力 学法 で は力 の 計算 がもっとも重いので、一 般には前者 を用い る とよい 。 さて、シンプレクティック積分はエネルギーを保存する積分法であった。温度制御を行う と、当然全系のエネルギーは保存しない。また、温度制御がある系の時間発展はシンプレク ティック性を持たないため、そもそもシンプレクティック積分をすることはできない。しか し、シンプレクティック積分で行った指数分解の手法を用いて、同様な数値積分法を構築す るこ と がで き る。 ハミルトニアンHで記述され、かつNos´e–Hoover法で温度制御されている系を考えよう。 運動 方 程式 は ˙ p = −∂H ∂q − pζ ˙ q = ∂H ∂p ˙ ζ = 1 τ2 ( p2− 1 β ) (29) で与 え られ る 。こ の系のリュービル演算子Lは L = −i [ p ∂ ∂q − q ∂ ∂p− pζ ∂ ∂p+ 1 τ2 ( p2− 1 β ) ∂ ∂ζ ] (30) であ る 。 シンプレクティック積分においては、リュービル演算子を簡単に積分できる部分に分ける ことで数値積分公式を導出していた。その本質は、指数の肩にリュービル演算子を乗せた時 に項が有限次で切れるということにあったが、この運動方程式には−pζ ∂ ∂pという非可換な項 があ るた め に 同 様な ことはできない。し かし、こ の場 合は たま たま exp ( −pζ ∂ ∂p ) = p exp (−ζ) ∂ ∂p (31) と、指数の展開項を無限次まで足し上げた場合の厳密な表式を求めることができるため、こ れを 使って 積 分法 を構 築することができる。具 体的に、一 次 の積分 法 は以下 の 通り。 q ← q +∂H ∂p∆t (32) p ← p −∂H ∂q∆t (33) p ← p exp (−ζ∆t) (34) ζ ← ζ + 1 τ2 ( p2− 1 β ) ∆t (35)
同様 に 、二次 の 対 称分解に 対応する公式も以下のよう に求 める こと がで きる 。 q ← q +∂H ∂p ∆t 2 (36) p ← p exp (−ζ∆t/2) (37) ζ ← ζ + 1 τ2 ( p2− 1 β ) ∆t 2 (38) p ← p −∂H ∂q ∆t (39) ζ ← ζ + 1 τ2 ( p2− 1 β ) ∆t 2 (40) p ← p exp (−ζ∆t/2) (41) q ← q +∂H ∂p ∆t 2 (42) シンプレクティック積分の場合と同様に、exp(i∆t2 LU)に対応する項を中央に持ってくること で、１ ス テップ に 力 の計 算が一度でよいようにしている 。 このように時間発展演算子の指数分解を行うことで数値積分法を構築する手法をRESPA (REversible System Propagator Algorithm)と呼ぶ[7]。シンプレクティック積分と異なり、一般の
RESPAにおいてはエネルギーが保存量とはならないが、もとの運動方程式が時間反転対称で あれば、離散化した後もその対称性を厳密に保つという性質がある。その性質がシミュレー ションにどれだけうれしいかは、実は定かではない。ただ、温度制御された系に普通のシン プレクティック積分を適用しても、それは単なる修正オイラー法であって特にうれしい性質 は持ってい な い、と いう ことは知っておいた方がよい 。 3.2 温 度制 御 に起 因する 振動 モード 系の時間発展を追うことができる、ということがモンテカルロ法などに比べた分子動力学 法の利点の一つであった。しかし、温度制御を行った場合、系に内在しない時間スケールを 持ち込む場合があるため注意が必要である。以下、Nos´e–Hoover法を用いた場合に系にどの よう なダ イ ナ ミ クスが 導入されるか調べてみよう。 ハミ ル ト ニア ン が H = 1 2m ∑ i p2i + V (q1, q2,· · ·) (43) で与えられるような、ポテンシャルVに運動が支配される系を考える。この系をNos´e–Hoover 法で 温 度制 御 した 場合 の運動方程式は ˙ pi = ∂V ∂qi − piζ ˙ qi = pi m ˙ ζ = 1 τ2 ( ∑ i p2i m − NkBT ) (44) で与えられるのであった。ただし、τは緩和時間パラメータである。一般にこの微分方程式 系は非可積分であるが、τが系の特徴的な時間スケール(分子運動のスケール)に比べて十分 長い場合、系はτに比例する特徴的な振動モードを持つ。実際のシミュレーション結果を図3 に示す。３次元Lennard-Jonesポテンシャル、N = 64000、熱浴をつけた場合とつけていない 場合のエネルギーの時間発展を比べた。実時間ではその差は分からないが、Nos´e–Hoover法
0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 40 50 60 70 80 90 100 110 120 130 Kinetic Energy Time Without Thermostat With thermostat 0 0.5 1 1.5 2 2.5 20 40 60 80 100 120 140 Power Spectra Frequency Without Thermostat With thermostat 図3:温度制御の有無と時間スペクトル。３次元Lennard-Jonesポテンシャル、64000粒子系。 Nose-HooverをつけないNVEアンサンブルとNose-HooverをつけたNVTアンサンブル。設定温
度はT = 1とした。左図が運動エネルギー(温度)。温度制御が無いものはT = 1.33程度だが、 温度制御のあるものはT = 1に制御されている。左図が温度の時間発展のフーリエ変換。温 度制御の無いものはほぼ白色ノイズで、原点に鋭いピークが一つ立つのみだが、温度制御の ある も のは 熱 浴に 起因す るピークが立つ。 により温度制御された系にはフーリエスペクトルにおいて特徴的なピークが立つことが分 かる。一般に、この振動モードの周期はτに比例、振幅はτに反比例することが知られてい る。この振動モードの起源を熱浴に起因するモードを断熱近似にて導出しよう[8]。簡単の ため、一自由度の調和振動子系を考えるが、多自由系においても自然ハミルトン系(運動エ ネル ギー がpの二次 で与え られる系)なら同じである。 設定 温 度をT = 1とし 、kB= 1の単位系をとると、運 動方程 式 は ˙ p = −q − pζ ˙ q = p ˙ ζ = 1 τ2(p 2− 1) (45) で与 え られ る 。ζの 項を時 間で一階微分して整理する と、 τ2ζ = 2p ˙¨ p (46) ˙ p = q− pζであ る から 、 τ2ζ¨ = 2p(−q − pζ) = −2p2ζ− 2pq ここ でζ˙を 用 い てp2を 消去 すれば、 τ2ζ =¨ −2ζ(τ2ζ + 1)˙ − 2pq (47) 整理 して 、 τ2ζ + 2τ¨ 2ζ ˙ζ + 2ζ =−2pq (48) ここで、τの値が大きければζの変化はpqに比べて十分遅く(断熱近似)、かつpqの平均は0で ある か ら、 τ2ζ + 2τ¨ 2ζ ˙ζ + 2ζ = 0 (49)
p
q
ζ
図4: Nos´e–Hoover法の位相空間の流れ場。付け加えられた項は、p− ζ平面における回転を表 現している。これが、等エネルギー面の間の遷移を表現し、結果として温度一定の運動を表 現す る。 さて 、こ こ で ζ(t) = τ−1ζ(t/τ )¯ (50) とい うス ケ ー リ ングを 考えると、式(49)は ¨¯ ζ + 2¯ζ ˙¯ζ + 2 ¯ζ = 0 (51) と、τが消去される。すなわち、この系には振幅がτに反比例し、周期がτに比例するモード があ るこ と が 分 かる。ま た 、x = ¯ζ, y = ˙¯ζと表記すると、 ˙ x = y (52) ˙ y = −2x(y + 1) (53) という力学系に帰着される。これは変数分離できて時間不変量E = x2+ y− log(y + 1)を得 る。曲線x2+ y− log(y + 1) = Eは閉曲線であるので、この力学系は周期運動をすることがわ かる。以上から、断熱近似をするとζはτに比例する周期を持つ周期関数であることがわか る。こ のモ ー ドが 図3のフーリエスペクトルに現れてい る。 Nos´e–Hoover系では振動モードが存在し、これが系に特徴的な時間スケールを持ち込んで いる。ここで、変数変換した系が時間不変量を持つのは、温度制御に用いる変数が一自由度 であることに起因しており、Nos´e–Hoover–Chain法やKinetic-Moments法など、多変数熱浴で は式(49)が非可積分となり、特徴的な振動モードは消える。しかし、系に非自明なモードを 持ち込んでいることには変わりがなく、鋭いピークの変わりに幅の広いモードがフーリエス ペクトルに現れるであろう。系が十分多自由度であればフーリエスペクトルはデルタ関数的 であるべきで、これは、熱浴制御された系ではフーリエモードで書いた線形応答理論が成り 立たないことを示唆する。このようなモードとNos´e–Poincar´e法との関連は未だ明らかでな い。また、Gaussian-Thermostat法は、特徴的な時間スケールを系に持ち込まないので線形応 答理論に悪影響を及ぼさない可能性があるが、まだ詳しいことは明らかになっていない。 ここでは計算を簡単にするために調和振動子を用いたが、実際のシミュレーションでどう なるかを見るために3次元Lennard-Jones粒子について温度制御をした計算結果を紹介してお0.85 0.9 0.95 1 1.05 1.1 0 50 100 150 200 250 300 350 400 450 500 Tau=1 Tau=10 Tau=100 図5:温度の時間発展のτ依存性。3次元Lennard-Jones粒子、設定温度T = 0.9で、密度0.7。粒 子数N = 681472。熱浴 の緩和時間はτ = 1, 10, 100の三種類 に ついて 表 示して あ る。 こう。運動 方程 式 は 以下の 通り。ただ しmは粒 子の質 量、Tは設 定 温度で あ る。 ˙ qi = pi m ˙ pi = ∂V ∂qi − ζpi ˙ ζ = 1 τ2 ( 1 3N 3N ∑ i p2i − kBT ) (54) 式(44)と異なり、運動エネルギーの項を3Nで割ってあることに注意。式(44)では同じτの値 でも粒子数が異なると振る舞いが大きく変わるが、式(54)では粒子数依存性はほぼなくなる ので、実際には式(54)の形を用いた方が便利である。計算結果を図5に示す。τの値として 1, 10, 100の三通りの場合について表示してある。調和振動子の場合と異なり、振幅にはτ依 存性 は ほと ん どな いが、振 動数はτに比 例するこ とが 分か る10。 3.3 エ ルゴ ー ド性 調和振動子にNos´e–Hooverを適用するとエルゴード性が失われ、結果としてカノニカル分 布を再現しないことが早くから知られていた(図6参照)。そこで、エルゴード性の実現のた め にNos´e–Hoover–Chain法[9]やKinetic-Moments法[10]な ど の 多 変 数 熱 浴 法 が 提 案 さ れ た 。 Nos´e–Hoover法では自由度を一つ追加することで温度を制御するが、多変数熱浴法では追加 自由度を二つ以上加えることで、温度制御しつつエルゴード性の実現も行う。しかし、なぜ Nos’e–Hooverではエルゴード性が失われ、多変数熱浴法ではエルゴード的となるのか、そし て付け加える変数はいくつ必要なのかはわかっていなかった。ここでは温度制御とエルゴー ド性 に つい て 議論 する[11, 12]。 もともとエルゴード性とは、マクロに言えば「物理量のアンサンブル平均と長時間平均が 等しい 」 hAi = lim T→∞ 1 T ∫ T 0 A(t)dt (55) 10実際にexp(−γt) sin(ωt)型の減衰振動を仮定してフィッティングをすると、時間スケールが t/τでスケールさ れ る こ と が 確 か め ら れ る 。
-0.8 -0.4 0 0.4 0.8 -0.8 -0.4 0 0.4 0.8 q p -0.8 -0.4 0 0.4 p -0.4 0 0.4 0.8 q -0.08 0 0.08 ζ 図6:調和振動子系にNos´e–Hoover法を適用した場合の例。左図は(p, q)平面をプロットしたも の。軌道が位相空間を覆い尽くさず、エルゴード的でないことがわかる。特にエネルギーに 初期条件に依存する上限と下限が存在する。右図は(p, q, ζ)空間における軌道の構造。軌道 がト ー ラス を 構成 してお り、何ら かの保存量があること を 示唆す る 。 というのが定義である。これはミクロには「すべての等価な微視的状態がすべて同じ確率で 実現する」ということを意味する。ハミルトンダイナミクスにおいてはエネルギーが保存量 となるので、運動がエルゴード的であるとは、位相空間における長時間軌道が等エネルギー 面を 覆 い尽 く すこ とであ る。 さて、温度制御された系でのエルゴード性の定義は、ミクロカノニカルの場合と少し異な るので注意が必要である。温度制御された系では保存量が存在しないため、軌道は位相空間 上すべてを覆い尽くす時に系がエルゴード的であると定義される。ただし、軌道がある位相 空間の微小部分を通る確率はその微小部分に対応するエネルギーをEとして、ボルツマン重 みexp(−E/kBT )に比例しなくてはならない。もし系に保存量が存在した場合、それはその ままエルゴード性の破れを意味する。調和振動子系H = p2/2 + q2/2にNos´e–Hoover系を適用 した 場合 の 運 動 方程 式として式(45) ˙ p = −q − pζ ˙ q = p ˙ ζ = 1 τ2(p 2− 1) からはじめよう。ここから断熱近似により時間不変量を構築し、エルゴード性が破れている こ と を 示 す。本 質 的 に は3.2節 と 同 じ 計 算 だ が 、よ り 分 か り や す く す る た め に 変 数 変 換 p = r cos θ、q = r sin θを施す。すると、運 動方程式は ˙r = −rζ cos2θ ˙ ζ = 1 τ2(r 2cos2θ− 1) となる。ここでτが十分大きければrやζに比べてθの運動は十分に速く、cos2θはほぼ定数で ある とみ な す こ とがで きる。そ こで、0 < t < 2πの間で 平均 を取 るこ とで ˙r = −πrζ ˙ ζ = π τ2(r 2− 2)
とな る 。これ は 変 数分離形 となっており、直ちに時間 不変 量 r2 2 + τ2ζ2 2 − ln r 2 = C (56) を 得 る 。Cは 初 期 条 件 か ら 決 ま る 定 数 で あ る 。こ れ が 図6の ト ー ラ ス に 対 応 し て い る 。 τ2ζ2/2≥ 0で あ り、かつH = r2/2で あるから、全 エネルギーHは不等 式 H−1 2ln H ≤ C 0 (57) を満たす。ただしC0 ≡ C + ln 2である。式(57)の左辺は下に凸の関数であるから、この不等 式はHに最大値と最小値が存在することを意味する。もし系がカノニカル分布をするのであ れば、Hは0から無限大の値までボルツマン重みexp(−βH)に比例する確率で揺らぐはずで あった。しかし、系に内在する周期と熱浴の時間スケールが解離していると漸近的に構成さ れる保存量を通じて、エネルギーに上限と下限が生じる。その結果得られる分布はカノニカ ル分布にならない。つまり、長時間の計算をしても位相平均と時間平均が一致せず、エル ゴー ド性 を 失 う こと になる。 ここで、エルゴード性を失う本質は、式(3.3)が変数分離形となって保存量が構成されるこ とであった。し たがって系がどんなに多自由度であっても 、系 に周 期運 動が 内在 し、かつ Nos´e–Hoover法に内在する周期と時間スケールが分離していればエルゴード性が破れる可能 性 が あ る 。ま たKinetic–Moments法 やNos´e–Hoover–Chain法 な ど の 多 変 数 熱 浴 で は 、(最 低 で も)三自由度あるために、断熱近似した場合においても運動方程式が求積できず、その結果 エルゴード性が失われないことがわかる。すなわち、熱浴変数は二変数あればエルゴード性 は失 われ な い 。 3.4 そ の他 の 問題 温度制御をかけているはずなのに、温度が正しく制御されない場合がある。これは、系に 熱浴とカップルしやすい自由度とカップルしにくい自由度が存在する時に起きることが多 い。Nos´e–Hoover法は、全ての粒子の運動エネルギーの平均を制御する手法であるため、「系 の右端にある粒子が遅くなった」という事象が、系の左端にある(実際には相互作用してい ないはずの)粒子に影響を与える。このため系に「熱を受け取りやすい場所」と「熱を受け 取り に くい 場 所」があると 、 1. まず 熱を 受 け取り や すい場所が加熱され る。 2. 平均 エネ ル ギーが 上 がるためにNos´e–Hooverにより 全 体が冷 却 される 。 3. する と 、相 対 的に 熱を 受け取りにくい場所 の温度 が下 がる 。 4. 以上 を 繰り 返 し 、系に 人為的な熱流が生じ る。 ということがおきる(図7)。ここで、「熱を吸収しやすい/しづらい場所」と書いたが、実際に は系に熱の吸収のしやすさの異なる自由度が同居する場合も同様なことがおきる。この際、 ある自由度が加熱、他の自由度が冷却され、全体として等分配則が成り立たないような状態 に な り う る 。こ の よ う な 場 合 に は 、自 由 度 ご と に う ま く 熱 浴 を 設 定 す る か 、も し く は Langevin熱浴のような局所的に働く温度制御法を用いるのが良い。Langevin熱浴は粒子ごと に独立に熱浴がつくため、たとえば空間的に離れた粒子が影響を及ぼしあうことはない。た だし確率的な力を導入するために決定論的な時間発展は失われ、かつ摩擦項の存在で運動方 程式 が時 間 反 転 対称 ではなくなる。
熱を受け取りやすい場所 熱を受け取りにくい場所 熱流 図7:一様で無い系、特に系に熱を吸収しやすい場所と吸収しづらい場所が共存する系では、 うまく 温度 制 御が できな い場合がある。 4 おわ り に 分子動力学法における温度制御の意味と、その実用上の注意を述べてきた。もともとハミ ルトンの運動方程式は、極小作用の原理から導かれている。しかし、温度を制御するという ことは運動方程式を修正するということであり、それはすなわち対応する原理を失っている ということを意味する。分子動力学法は時間発展を追えるのがメリットであるが、温度制御 をかけた運動方程式の時間発展の物理的な意味は必ずしも明確ではない。さらに温度制御に よる人為的な時間スケールも導入される場合がある。そこで、実際の数値計算では温度制御 を全エネルギーの調整手段であると割り切って、温度がある程度緩和したら熱浴制御を切 り、ミク ロカ ノニ カルで 数値 計算を行うと良い。 計算機の能力の増大と共にパッケージソフトウェアも増え、自分ではコードを書かずに分 子動力学法計算を使ってできることも増えてきた。数値計算は手段であって目的ではない。 しかし、自分がいったい何をやっているかよく理解していないととんでもない誤りを犯すこ ともあるし、そうでなくとも新しい結果が数値計算のアヤか新発見か区別できないのでは困 る。数値計算の中身をブラックボックスにせず、自分が行っている計算がどんなアルゴリズ ムでど ん な計 算 を行っているのかを常に把握するよう にし てお きた い。 参考文 献
[1] L. V. Woodcock, Chem. Phys. Lett., 10, 257 (1971).
[2] W. G. Hoover, A. J. C. Ladd, and B. Moran, Phys. Rev. Lett., 48, 1818 (1982); D. J. Evans, J.
Chem. Phys, 78, 3297 (1983).
[3] S. Nos´e, Mol. Phys. 52 255 (1984).
[4] W. G. Hoover, Phys. Rev. A, 31, 1695 (1985).
[5] H. J. C. Berendsen et al., J. Chem. Phys., 81, 3684 (1984).
[6] S. D. Bond, B. J. Leimkuhler, and B. B. Lairdy, J. Comput. Phys., 151, 114 (1999). [7] M. Tuckerman, B. J. Berne, and G. J. Martyna, J. Chem. Phys., 97, 1990 (1992) [8] H. Watanabe and H. Kobayashi, Molecular Simulation, 33, 77 (2007).
[9] G. J. Martyna, M. L. Klein, and M. Tuckerman, J. Chem. Phys., 97, 2635 (1992). [10] W. G. Hoover and B. L. Holian, Phys. Lett. A, 211, 253 (1996).
[11] 渡辺 宙 志,小林 礼人,日本物理学会誌 第62巻第10号, 785 (2007) [12] H. Watanabe and H. Kobayashi, Phys. Rev. E, 75, 040102(R), (2007).
