分子動力学法(2) 温度制御や圧力制御 / Simulation 06

シミュレーション工学

Transcript

1. 1 96 分子動力学法(2) 温度制御と圧力制御 シミュレーション工学 慶應義塾大学大学院理工学研究科基礎理工学専攻物理情報専修 渡辺宙志
2. 2 96 • 分子動力学法における温度や圧力の定義 • 温度や圧力制御アルゴリズム • 温度や圧力といった物理量の定義が、必ずし も自明ではないということを知る 本講義の目的 物理量の定義と制御 • 分子動力学法における温度や圧力の定義は 非自明 • 温度とは？圧力とは？
3. 3 96 肌寒い お風呂がぬるい アイロンが熱い 気体の温度 液体の温度 固体の温度 これら全てに共通する「温度」とはなんだろう？ どうやって数値化しているのだろう？
4. 4 96 ガラス温度計 液体の膨張による体積変化 温湿度計(バイメタル方式) 熱膨張率の異なる金属を貼り合わせた もの デジタル体温計(サーミスタ方式) 温度により抵抗が変わる物質を利用
5. 5 96 熱電対 二種類の金属の熱電能の差を起電 力に変える(ゼーベック効果) 赤外線温度計 物体から放出される赤外線のエネ ルギーを測定 これら全ては、現実物質の温度変化に対する性質の変化を利用 基準となる温度を使って較正する必要がある そもそも温度とは？
6. 6 96 変数(Variable) 我々がa prioriに認める物理量 支配方程式に直接でてくる量 観測量(Observable) 支配方程式には含まれない量 変数から導かれる量 何が変数で何が観測量かは支配方程式により異なる
7. 7 96 熱伝導方程式 𝜕𝑇 𝜕𝑡 = 𝜅 𝜕2𝑇 𝜕𝑥2 位置𝑥、時刻𝑡、温度𝑇(𝑥, 𝑡)は全て変数(variable) この支配方程式を認めた時点で「この世界には温度 というものがあり、この方程式に従う」と宣言
8. 8 96 ナビエ・ストークス方程式(非圧縮) 𝜕 Ԧ 𝑣 𝜕𝑡 + Ԧ 𝑣 ⋅ ∇ = − 1 𝜌 ∇𝑃 + 𝜈Δ Ԧ 𝑣 + Ԧ 𝐹 速度場 Ԧ 𝑣と圧力場𝑃が変数(variable) この支配方程式を認めた時点で、圧力という量 の存在を認めている。 「圧力とは何か？」という問いが意味をもたない
9. 9 96 ハミルトンの運動方程式 𝑑𝑝𝑖 𝑑𝑡 = − 𝜕𝐻 𝜕𝑞𝑖 𝑑𝑞𝑖 𝑑𝑡 = 𝜕𝐻 𝜕𝑝𝑖 エネルギー𝐻、座標𝑞𝑖 、運動量𝑝𝑖 が変数(variable) 温度𝑇や圧力𝑃は全て観測量(observable) この世界では、温度や圧力は定義しなければならない
10. 10 96 示量変数(extensive variable) 全く同じ系を二つつけた時、二倍になる量 体積𝑉、エントロピー𝑆、物質量𝑁 示強変数(intensive variable) 全く同じ系を二つつけた時、値が変わらない量 圧力𝑃、温度𝑇、化学ポテンシャル𝜇
11. 11 96 示量変数 示強変数 体積𝑉 圧力𝑃 かけたらエネルギーになる量をお互いに共役と呼ぶ 温度𝑇 エントロピー𝑆 物質量𝑁 化学ポテンシャル𝜇
12. 12 96 𝑈(𝑆, 𝑉, 𝑁) 𝐻(𝑆, 𝑃, 𝑁) 𝐹(𝑇, 𝑉, 𝑁) 𝐺(𝑇, 𝑃, 𝑁) Ω(𝑇, 𝑉, 𝜇) 𝐹 = 𝑈 − 𝑇𝑆 𝐻 = 𝑈 + 𝑃𝑉 𝐺 = 𝐻 − 𝑇𝑆 𝐺 = 𝐹 + 𝑃𝑉 Ω = 𝐹 − 𝜇𝑁 共役な量＝ルジャンドル変換で互いに変換する量 全て示量変数
13. 13 96 示量変数 示強変数 体積𝑉 圧力𝑃 かけたらエネルギーになる量をお互いに共役と呼ぶ 温度𝑇 エントロピー𝑆 物質量𝑁 化学ポテンシャル𝜇 示量変数と示強変数、どちらが「基本的な量」だろうか？
14. 14 96 a prioriな量とは、我々が未定義で使う量のこと 熱力学では、共役な量のどちらかがa prioriな量 数値計算では、支配方程式に含まれる変数 それ以外の物理量は全てa prioriな量から定義される 分子動力学法において、温度や圧力を定義しよう
15. 15 96 温度𝑇は運動エネルギー𝐾に比例 𝑇 = 2𝐾 3𝑁𝑘 𝑁：粒子数 𝑘：ボルツマン定数 運動エネルギー𝐾は運動量の二乗和 𝐾 = ෍ 𝑖 𝑝𝑖 2 2𝑚 以上から 𝑇 = 2 3𝑁𝑘 ෍ 𝑖 𝑝𝑖 2 2𝑚
16. 16 96 𝑇 = 2 3𝑁𝑘 ෍ 𝑖 𝑝𝑖 2 2𝑚 こちらは全部知っているので これが計算できる 𝐾 = ෍ 𝑖 𝑝𝑖 2 2𝑚 ← これは単なる定義 𝑇 = 2𝐾 3𝑁𝑘 ← これはどこからきたのか？
17. 17 96 (𝑝, 𝑞)空間の分布関数𝑓(𝑝, 𝑞)を考える ある点(𝑝, 𝑞)におけるエネルギー𝐻(𝑝, 𝑞)を定義する 規格化条件 න 𝑓𝑑𝑝𝑑𝑞 = 1 内部エネルギー න 𝐻𝑓𝑑𝑝𝑑𝑞 = 𝑈
18. 18 96 この系のエントロピーを以下のように定義する 𝑆 = −𝑘 න 𝑓 log 𝑓 𝑑𝑝𝑑𝑞 න 𝑓𝑑𝑝𝑑𝑞 = 1 න 𝐻𝑓𝑑𝑝𝑑𝑞 = 𝑈 以下の条件を満たしつつ、エントロピーを最大化する分布 関数を求めたい 規格化条件 エネルギーの期待値が𝑈
19. 19 96 ラグランジュの未定定数法より(※) 𝐼 = 𝛼 න 𝑓𝑑𝑝𝑑𝑞 + 𝛽 න 𝐻𝑓𝑑𝑝𝑑𝑞 + න 𝑓 log 𝑓 𝑑𝑝𝑑𝑞 ※ あとの便利のために符号をかえたりボルツマン定数を吸収させたりしている 汎関数微分から 𝛿𝐼 = 0 𝛼 + 𝛽𝐻 + 1 + log 𝑓 = 0 以上から 𝑓 = 𝑍−1exp(−𝛽𝐻) 𝑍 = exp(𝛼 + 1) カノニカル分布 分配関数(規格化条件から決まる)
20. 20 96 𝑆 = −𝑘 න 𝑓 log 𝑓 𝑑𝑝𝑑𝑞 エントロピーの定義 𝑓 = 𝑍−1exp(−𝛽𝐻) log 𝑓 = −𝛽𝐻 − log 𝑍 𝑆 = 𝛽𝑘 න 𝐻𝑓𝑑𝑝𝑑𝑞 + 𝐶 = 𝑈 カノニカル分布を代入 𝑑𝑆 𝑑𝑈 = 𝛽𝑘
21. 21 96 𝑑𝑆 𝑑𝑈 = 𝛽𝑘 𝑑𝑆 𝑑𝑈 = 1 𝑇 さきほど求めた関係式 熱力学関係式 逆温度と温度が結びついた 𝛽 = 1 𝑘𝑇
22. 22 96 全ての物理量は位相空間の座標(𝑝, 𝑞)の関数 物理量𝐴(𝑝, 𝑞)の期待値 𝐴 ≡ න 𝐴𝑓𝑑𝑝𝑑𝑞 位相空間の各点(𝑝, 𝑞)で定義される物理量𝐴(𝑝, 𝑞) に、その点における分布関数𝑓(𝑝, 𝑞)をかけて全位 相空間において平均せよ
23. 23 96 分布がカノニカル分布である時 という 𝑝 𝜕𝐻 𝜕𝑝 𝑝 𝜕𝐻 𝜕𝑝 = 𝑍−1 න 𝑝 𝜕𝐻 𝜕𝑝 exp −𝛽𝐻 𝑑𝑝𝑑𝑞 𝐴 ≡ න 𝐴𝑓𝑑𝑝𝑑𝑞 𝑓 = 𝑍−1exp(−𝛽𝐻) 物理量の期待値を考える 物理量𝐴の期待値 𝐴 = 𝑝𝜕𝑝 𝐻 期待値を求めたい物理量 カノニカル分布
24. 24 96 𝜕 𝜕𝑝 exp −𝛽𝐻 = −𝛽 𝜕𝐻 𝜕𝑝 exp(−𝛽𝐻) であるから 𝑝 𝜕𝐻 𝜕𝑝 = 𝑍−1 න 𝑝 𝜕𝐻 𝜕𝑝 exp −𝛽𝐻 𝑑𝑝𝑑𝑞 = 𝑍−1 න − 𝑝 𝛽 𝜕 𝜕𝑝 exp(−𝛽𝐻) 𝑑𝑝𝑑𝑞 = 1 𝛽 𝑍−1 න exp(−𝛽𝐻) 𝑑𝑝𝑑𝑞 部分積分 = 𝑍 = 1/𝛽
25. 25 96 𝑝 𝜕𝐻 𝜕𝑝 = 1 𝛽 = 𝑘𝑇 ハミルトニアンが𝐻 = 𝑝2 2𝑚 + 𝑉(𝑞)の形をしている場合 𝑝 𝜕𝐻 𝜕𝑝 = 𝑝2 𝑚 = 𝑘𝑇 𝑝2 2𝑚 = 1 2 𝑘𝑇 等分配則 3次元多粒子系の場合 ෍ 𝑖 𝑝𝑖 2 2𝑚 = 3 2 𝑁𝑘𝑇
26. 26 96 𝑇 = 2𝐾 3𝑁𝑘 ෍ 𝑖 𝑝𝑖 2 2𝑚 = 3 2 𝑁𝑘𝑇 カノニカル分布の部分積分から ෍ 𝑖 𝑝𝑖 2 2𝑚 = 𝐾 運動エネルギーの定義 分子動力学法における(解析力学における) 温度の定義
27. 27 96 𝑝 𝜕𝐻 𝜕𝑝 = 1 𝛽 = 𝑘𝑇 以下の量が温度になったのは部分積分したから 全く同じ導出で以下の等式も成り立つ 𝑞 𝜕𝐻 𝜕𝑞 = 1 𝛽 = 𝑘𝑇 運動温度 状態温度 平衡状態では両者は同じ温度を与えるが、 非平衡状態では一般に一致しない
28. 28 96 運動温度 時間 状態温度 両者が一致していなければ緩和不足(非平衡状態) 平衡状態 非平衡状態
29. 29 96 • 温度が運動エネルギーに比例するのは等分配則 のため • 逆温度はエントロピーを最大化する際のラグラ ンジュの未定定数であり、熱力学関係を要請す ることで温度と結びつく • 等分配則はカノニカル分布をする系において 𝑝𝜕𝑝 𝐻の形の物理量の期待値が温度に比例するこ とによる(運動温度) • 同様に𝑞𝜕𝑞 𝐻も温度に比例する(状態温度) • 運動温度と状態温度は、平衡状態では一致する が、非平衡状態では必ずしも一致しない
30. 30 96 圧力：粒子が壁を押す力 壁がない系(周期的境界条件)で圧力は定義できるか？ 負の圧力はあり得るか？
31. 31 96 𝐺 = ෍ 𝑖 𝑝𝑖 𝑞𝑖 この系で以下の量を考える 3次元N粒子系{𝑝𝑖 , 𝑞𝑖 }を考える x,y,z座標をまとめてiのインデックスに押し込める (𝑖 = 1,2, ⋯ , 3𝑁) この量をビリアル(virial)と呼び、クラウジウスにより導入された この量から圧力を定義する
32. 32 96 𝐺 = ෍ 𝑖 𝑝𝑖 𝑞𝑖 ሶ 𝐺 = ෍ 𝑖 (𝑝𝑖 ሶ 𝑞𝑖 + ሶ 𝑝𝑖 𝑞𝑖 ) 𝑝𝑖 ሶ 𝑞𝑖 = 𝑝𝑖 2 𝑚 = 𝑘𝑇 時間微分 エネルギーの等分配則 ෍ 𝑖 𝑝𝑖 ሶ 𝑞𝑖 = 3𝑁𝑘𝑇 理想気体からの寄与
33. 33 96 ሶ 𝐺 = ෍ 𝑖 (𝑝𝑖 ሶ 𝑞𝑖 + ሶ 𝑝𝑖 𝑞𝑖 ) ሶ 𝑝𝑖 = 𝑓𝑖 + 𝐹𝑖 運動量の時間変化＝粒子にかかる力 外力 (壁から粒子iに働く力) 粒子間に働く力 (他粒子から粒子iに働く力の合計) 𝑓𝑖 𝐹𝑖
34. 34 96 ෍ 𝑖 𝑞𝑖 𝑓𝑖 = ෍ 𝑖<𝑗 𝑞𝑖𝑗 𝑓𝑖𝑗 粒子間に働く力を作用反作用を使って整理 𝑞𝑖𝑗 = 𝑞𝑖 − 𝑞𝑗 𝑓𝑖𝑗 = 𝑓𝑖 − 𝑓𝑗 外力をガウスの定理を使って整理(※) ෍ 𝑖 𝑞𝑖 𝐹𝑖 = − 𝑃 න 𝜕𝑉 Ԧ 𝑟 ⋅ 𝑛𝑑𝐴 = −𝑃 න 𝑉 ∇ ⋅ 𝑛𝑑𝑉 = −3𝑃𝑉 まとめると ෍ 𝑖 𝑞𝑖 ሶ 𝑝𝑖 = ෍ 𝑖<𝑗 𝑞𝑖𝑗 𝑓𝑖𝑗 − 3𝑃𝑉 ※ 計算の詳細は今は気にしなくて良いです
35. 35 96 ሶ 𝐺 = ෍ 𝑖 (𝑝𝑖 ሶ 𝑞𝑖 + ሶ 𝑝𝑖 𝑞𝑖 ) = 3𝑁𝑘𝑇 + ෍ 𝑖<𝑗 𝑞𝑖𝑗 𝑓𝑖𝑗 − 3𝑃𝑉 ビリアルの時間平均がゼロであることを仮定 𝑃𝑉 = 3𝑁𝑘𝑇 + 1 3 ෍ 𝑖<𝑗 𝑞𝑖𝑗 𝑓𝑖𝑗 分子動力学法における圧力の定義 理想気体の寄与 相互作用からの寄与
36. 36 96 • 粒子に働く力を粒子間力と外力にわけたが、 周期的境界条件など、壁がない系では圧力 は定義できないのか？ • ビリアルの時間微分の平均がゼロであると いう仮定はどういう意味を持つのか？ • そもそもビリアルってなに？ 熱力学関係式とハミルトニアンから圧力を導出
37. 37 96 これから計算がごちゃごちゃするが全部追う必要はない 分配関数から必要な物理量が全て求められる、 という感覚だけ身に着ける
38. 38 96 𝑃 = − 𝜕𝐹 𝜕𝑉 𝑇 以下の熱力学関係式から出発する ヘルムホルツ自由エネルギーと分配関数 𝐹 = −𝑘𝑇 log 𝑍 𝑃 = 𝑘𝑇 𝜕 log 𝑍 𝜕𝑉 = 𝑘𝑇 𝑍 𝜕𝑍 𝜕𝑉
39. 39 96 𝜕𝑍 𝜕𝑉 この量を計算したい 𝑍 = න exp −𝛽𝐻 𝑑Γ 分配関数は陽に体積に依存しない 体積変化はハミルトニアンにのみ影響 体積が変化したときのハミルトニアンの応答を調べる 𝜕𝐻 𝜕𝑉 𝐿 𝛼𝐿 サイズ𝛼倍 体積𝛼3倍
40. 40 96 𝑍 = න exp −𝛽𝐻 𝑑Γ 𝑑Γ = ෑ 𝑖 𝑑𝑝𝑖 𝑑𝑞𝑖 分配関数とハミルトニアン 𝜕𝑍 𝜕𝑉 = න −𝛽 𝜕𝐻 𝜕𝑉 exp −𝛽𝐻 𝑑Γ 𝑃 = 𝑘𝑇 𝑍 𝜕𝑍 𝜕𝑉 Vで微分 に代入 𝑃 = 𝑘𝑇 𝑍 න −𝛽 𝜕𝐻 𝜕𝑉 e−𝛽𝐻 𝑑Γ = − 𝜕𝐻 𝜕𝑉
41. 41 96 𝑃 = − 𝜕𝐻 𝜕𝑉 を計算するために空間を一様に𝛼倍にスケールする 𝑞𝑖 → 𝛼𝑞𝑖 𝑝𝑖 → 𝑝𝑖 /𝛼 ሶ 𝑞𝑖 → 𝛼 ሶ 𝑞𝑖 𝑝𝑖 = 𝜕𝐿 𝜕 ሶ 𝑞𝑖 座標とスケーリングが 異なることに注意 𝑉 → 𝛼3𝑉 𝑑𝑉 𝑑𝛼 → 3𝛼2𝑉 座標と運動量の変化 体積変化
42. 42 96 もともとのハミルトニアン 𝐻 = ෍ 𝑖 𝑝𝑖 2 2𝑚 + ෍ 𝑖<𝑗 Φ(𝑞𝑖𝑗 ) = 𝐾 𝑝𝑖 + 𝑉({𝑝𝑖 }) スケールされたハミルトニアン 𝐻(𝛼) = ෍ 𝑖 𝑝𝑖 2 2𝛼2𝑚 + ෍ 𝑖<𝑗 Φ(𝛼𝑞𝑖𝑗 ) 𝜕𝐻 𝜕𝑉 = lim 𝛼→1 𝜕𝐻 𝜕𝛼 𝑑𝛼 𝑑𝑉 = 1 3𝑉 lim 𝛼→1 𝜕𝐻 𝜕𝛼 ハミルトニアンの体積微分 運動項 ポテンシャル項 𝑑𝑉 𝑑𝛼 → 3𝛼2𝑉
43. 43 96 lim 𝛼→1 𝜕 𝜕𝛼 ෍ 𝑖 𝑝𝑖 2 2𝛼2𝑚 = −2𝐾 = −3𝑁𝑘𝑇 運動項 ポテンシャル項 lim 𝛼→1 𝜕Φ 𝛼𝑞𝑖𝑗 𝜕𝛼 = Φ′ 𝑞𝑖𝑗 𝑞𝑖𝑗 = −𝑓𝑖𝑗 𝑞𝑖𝑗 まとめると 𝜕𝐻 𝜕𝑉 = 1 3𝑉 lim 𝛼→1 𝜕𝐻 𝜕𝛼 = 1 3𝑉 −3𝑁𝑘𝑇 + ෍ 𝑖<𝑗 𝑞𝑖𝑗 𝑓𝑖𝑗
44. 44 96 𝑃 = − 𝜕𝐻 𝜕𝑉 であり 𝜕𝐻 𝜕𝑉 = 1 3𝑉 −3𝑁𝑘𝑇 + ෍ 𝑖<𝑗 𝑞𝑖𝑗 𝑓𝑖𝑗 であるから 𝑃𝑉 = 3𝑁𝑘𝑇 + 1 3 ෍ 𝑖<𝑗 𝑞𝑖𝑗 𝑓𝑖𝑗 ビリアル定理と同じ公式が得られた
45. 45 96 • 分子動力学法の世界には、一般化運動量、一般 化座標、そしてハミルトニアンしか存在しない →圧力は与えられるものではなく定義するもの • 熱力学関係式から出発し、分配関数から圧力の 表式をもとめた • 外力や内力、ビリアルといった概念を使わずに、 系のスケーリングへの応答だけから圧力を定義 →周期的境界条件でも適用可能 分配関数から必要な物理量をなんでも求めることができる
46. 46 96 ハミルトンの運動方程式 𝑑𝑝𝑖 𝑑𝑡 = − 𝜕𝐻 𝜕𝑞𝑖 𝑑𝑞𝑖 𝑑𝑡 = 𝜕𝐻 𝜕𝑝𝑖 エネルギーが変化しない 通常は体積も変化しない 熱の出入りもない ミクロカノニカルアンサンブルになる (NVE) 時間発展について
47. 47 96 ミクロカノニカルアンサンブル 粒子数N、体積V、エネルギーEが一定 NVEアンサンブル 熱浴(heatbath) カノニカルアンサンブル 粒子数N、体積V、温度Tが一定 NVTアンサンブル エネルギーが揺らぐ
48. 48 96 𝐸 𝑓(𝐸) エネルギーが完全に固定で 揺らがない NVEアンサンブル 𝐸 𝑓(𝐸) エネルギーはある幅を もって揺らぐ NVTアンサンブル
49. 49 96 実際の物質はNが非常に大きいため、NVTアン サンブルでも事実上エネルギーは揺らがない 𝐸 𝑓(𝐸) NVEアンサンブル 𝐸 𝑓(𝐸) NVTアンサンブル 𝑁 ≫ 1の時
50. 50 96 我々が知りたい量は温度依存性であることが多い 300Kにおける水の圧力は？ 300Kにおけるポリエチレンの弾性率は？ ミクロカノニカル分布とカノニカル分布は本質的に 変わらないので、興味ある温度𝑇に対応するエネル ギー𝑈(𝑇)を持った系でNVEシミュレーションをすれ ば良い
51. 51 96 温度𝑇における内部エネルギー𝑈(𝑇)は 物質により異なる 金の比熱：130 [J/kg ℃] 水の比熱：4182 [J/kg ℃] 同じ温度でも、温まり易い物質はエネルギーが低く 温まり難い物質はエネルギーが高い ※ いずれも20℃における値
52. 52 96 ある温度𝑇におけるエネルギー𝑈(𝑇)を知りたい 比熱が高いほど同じ温度でのエネルギーは高い 温度変化に対するエネルギーの変化率が比熱 𝐶 = 𝜕𝑈 𝜕𝑇 𝑈 𝑇 = න 0 𝑇 𝐶𝑑𝑇 それを積分したものが全エネルギー この量を知りたい
53. 53 96 𝑈 𝑇 ≡ 𝐻 = 𝑍−1 න 𝐻 exp −𝛽𝐻 𝑑Γ 内部エネルギーの定義 𝑍 = න exp −𝛽𝐻 𝑑Γ 分配関数の定義 𝜕𝑍 𝜕𝛽 = −𝛽 න 𝐻 exp −𝛽𝐻 𝑑Γ = −𝛽𝑍 𝐻 𝑈 𝑇 = − 𝜕 log 𝑍 𝜕𝛽 𝛽で微分
54. 54 96 𝑈 𝑇 = − 𝜕 log 𝑍 𝜕𝛽 内部エネルギーの温度依存性がわかる 分配関数がわかる=問題が厳密に解けている 厳密に解けてない問題を解きたいのだから 一般に𝑈(𝑇)はわからない フィードバック制御による温度調整
55. 55 96 1. 全エネルギー一定の計算を行う 2. 運動エネルギーから温度を計算する 3. 目標温度との差を見てエネルギーを調整 4. 1.～3.を繰り返す 𝑇 = 2𝐾 3𝑁𝑘
56. 56 96 圧縮率：圧力が増えた時、どれだけ体積が減るか 𝜅 = − 1 𝑉 𝜕𝑉 𝜕𝑃 熱力学的に安定な系は、圧縮率が正でなくてはならない 平衡状態でピストンを少し押す 圧縮率が正なら 押し返される(安定) 圧縮率が負なら さらに引き込まれる (不安定) 圧縮率が正 体積を増やせば圧力が下がる 体積を減らせば圧力が上がる
57. 57 96 1. 全エネルギー一定の計算を行う 2. ビリアル定理から圧力を計算する 3. 目標圧力との差を見て体積を調整 4. 1.～3.を繰り返す 数値シミュレーションにおける体積変化とは？ 周期的境界条件ではどうするか？
58. 58 96 スケールされたハミルトニアン 𝐻(𝛼) = ෍ 𝑖 𝑝𝑖 2 2𝛼2𝑚 + ෍ 𝑖<𝑗 Φ(𝛼𝑞𝑖𝑗 ) スケール因子𝛼を一般化座標に含めたハミルトニアン 𝐻 = ෍ 𝑖 𝑝𝑖 2 2𝛼2𝑚 + ෍ 𝑖<𝑗 Φ(𝛼𝑞𝑖𝑗 ) + 𝜋2 2𝑀 + 𝑃0 𝛼3 このハミルトニアンに従う系は、目標圧力𝑃0 に制御される (Andersenのハミルトニアン)
59. 59 96 𝐻 = ෍ 𝑖 𝑝𝑖 2 2𝛼𝑚 + ෍ 𝑖<𝑗 Φ(𝛼𝑞𝑖𝑗 ) + 𝜋2 2𝑀 + 𝑃0 𝛼3 𝜋と𝛼が共役な変数であると考える 𝑑𝛼 𝑑𝑡 = 𝜕𝐻 𝜕𝜋 = 𝜋 𝑀 𝛼が座標、𝜋が運動量、𝑀が質量となる 仮想粒子 𝑑𝜋 𝑑𝑡 = − 𝜕𝐻 𝜕𝛼 = 1 𝛼3 ෍ 𝑖 𝑝𝑖 2 𝑚 − ෍ 𝑖<𝑗 Φ′ 𝛼𝑞𝑖𝑗 𝑞𝑖𝑗 − 3𝛼2𝑃0 仮想粒子に働く力
60. 60 96 スケールされた運動量と座標を導入して書き直す ෤ 𝑝𝑖 = 𝑝𝑖 /𝛼 ෤ 𝑞𝑖 = 𝛼𝑞 𝑑𝜋 𝑑𝑡 = 1 𝛼 ෍ 𝑖 ෤ 𝑝𝑖 2 𝑚 − 1 𝛼 ෍ 𝑖<𝑗 Φ′ ෤ 𝑞𝑖𝑗 ෤ 𝑞𝑖𝑗 − 3𝛼2𝑃0 = 2𝐾 = − ሚ 𝑓𝑖𝑗 = 3𝑉/𝛼 = 1 𝛼 2𝐾 + ෍ 𝑖<𝑗 ሚ 𝑓𝑖𝑗 ෤ 𝑞𝑖𝑗 − 3𝑃0 𝑉 = 3𝑃𝑉 = 3𝑉 𝛼 𝑃 − 𝑃0
61. 61 96 𝑑𝛼 𝑑𝑡 = 𝜕𝐻 𝜕𝜋 = 𝜋 𝑀 𝑑𝜋 𝑑𝑡 = 3𝑉 𝛼 𝑃 − 𝑃0 𝐻 = ෍ 𝑖 𝑝𝑖 2 2𝛼𝑚 + ෍ 𝑖<𝑗 Φ(𝛼𝑞𝑖𝑗 ) + 𝜋2 2𝑀 + 𝑃0 𝛼3 追加された仮想粒子 現在の圧力が設定圧力より高ければ加速 そうでなければ減速 仮想粒子の運動量が正なら膨張 負なら収縮 圧縮率が正 体積を増やせば圧力が下がる 体積を減らせば圧力が上がる 最終的に圧力が目標圧力に収束する
62. 62 96 𝐻 = ෍ 𝑖 𝑝𝑖 2 2𝛼𝑚 + ෍ 𝑖<𝑗 Φ(𝛼𝑞𝑖𝑗 ) + 𝜋2 2𝑀 + 𝑃0 𝛼3 運動を支配するハミルトニアン ෤ 𝑝𝑖 = 𝑝𝑖 /𝛼 ෤ 𝑞𝑖 = 𝛼𝑞 スケールされた運動量と座標 圧力が制御されるのは、スケールされた世界 「運動方程式に従う座標と運動量」と 「我々が観測する座標と運動量」が異なる
63. 63 96 • 熱力学的に安定な系では、体積が増えると圧力は下がり、 体積が減ると圧力が上がる • ハミルトニアンにスケールをつかさどる仮想粒子を追加し たハミルトニアンを考える(Andersenのハミルトニアン) • Andersenのハミルトニアンは、現在の圧力が目標圧力より 高いと系が膨張し、低いと系が収縮するダイナミクスを記 述する • スケールされた座標と運動量を観測すると、圧力が制御さ れているように見える • 「運動に従う変数」と「観測する変数」が異なる
64. 64 96 𝐻 = ෍ 𝑖 𝑝𝑖 2 2𝛼𝑚 + ෍ 𝑖<𝑗 Φ(𝛼𝑞𝑖𝑗 ) + 𝜋2 2𝑀 + 𝑃0 𝛼3 Andersenのハミルトニアン：仮想粒子による圧力制御 能勢のハミルトニアン：仮想粒子による温度制御 𝐻 = ෍ 𝑖 𝑝𝑖 2 2𝑠2𝑚 + ෍ 𝑖<𝑗 Φ(𝑞𝑖𝑗 ) + 𝜋2 2𝑄 + 3𝑁𝑘𝑇0 log 𝑠 追加された仮想粒子 スケールするのは 運動量のみ
65. 65 96 𝐻 = ෍ 𝑖 𝑝𝑖 2 2𝑠2𝑚 + ෍ 𝑖<𝑗 Φ(𝑞𝑖𝑗 ) + 𝜋2 2𝑄 + 3𝑁𝑘𝑇0 log 𝑠 𝑑𝜋 𝑑𝑡 = − 𝜕𝐻 𝜕𝑠 = 2 𝑠 ෍ 𝑖 𝑝𝑖 2 2𝑠2𝑚 − 3𝑁𝑘𝑇0 𝑠 ෤ 𝑝𝑖 = 𝑝𝑖 /𝑠 スケールされた運動量を導入 = 1 𝑠 ෍ 𝑖 ෤ 𝑝𝑖 2 𝑚 − 3𝑁𝑘𝑇0 𝑠 = 3𝑁𝑘 𝑠 𝑇 − 𝑇0 = 2𝐾 = 3𝑁𝑘𝑇
66. 66 96 𝐻 = ෍ 𝑖 𝑝𝑖 2 2𝑠2𝑚 + ෍ 𝑖<𝑗 Φ(𝑞𝑖𝑗 ) + 𝜋2 2𝑄 + 3𝑁𝑘𝑇0 log 𝑠 𝑑𝑠 𝑑𝑡 = 𝜕𝐻 𝜕𝜋 = 𝜋 𝑄 𝑑𝜋 𝑑𝑡 = 3𝑁𝑘 𝑠 𝑇 − 𝑇0 現在の温度が設定温度より高ければ減速 そうでなければ加速 仮想粒子の運動量が正なら減速 負なら加速 現在の温度が目標温度より高い → ሶ 𝜋 > 0なので、いつか𝜋 > 0になる → 𝑠が増える → ෤ 𝑝𝑖 が小さくなる→温度が下がる 最終的に温度が目標温度に収束する
67. 67 96 𝑝𝑖 の世界 時間 ෤ 𝑝𝑖 の世界 時間 能勢のハミルトニアンは仮想時間を スケールすることで温度を調整している 𝑠 いちいち仮想時間を考えるのは面倒くさい Nosé-Hoover法
68. 68 96 温度制御前のハミルトニアン 𝐻 𝑝, 𝑞 = 𝑝2 2𝑚 + 𝑉(𝑞) ※簡単のために一自由度系を考える 能勢のハミルトニアン 𝐻𝑁 = 𝐻 𝑝/𝑠, 𝑞 + 𝜋2 2𝑄 + 𝑘𝑇 log 𝑠 運動方程式 𝑑𝑝 𝑑𝑡 = − 𝜕𝐻𝑁 𝜕𝑞 𝑑𝑞 𝑑𝑡 = 𝜕𝐻𝑁 𝜕𝑝
69. 69 96 ෤ 𝑝 = 𝑝/𝑠 スケールされた運動量と時間を導入 ǁ 𝑡 = 𝑡/𝑠 𝑑𝑞 𝑑𝑡 = 𝜕𝐻𝑁 𝜕𝑝 = 𝜕𝐻(𝑝/𝑠, 𝑞) 𝜕𝑝 = 𝜕𝐻( ෤ 𝑝, 𝑞) 𝜕𝑝 = 𝜕𝐻 𝜕 ෤ 𝑝 𝑑 ෤ 𝑝 𝑑𝑝 = 1 𝑠 𝜕𝐻 𝜕 ෤ 𝑝 𝑑 ǁ 𝑡 = 𝑑𝑡/𝑠であるから 𝑑𝑞 𝑑 ǁ 𝑡 = 𝜕𝐻 𝜕 ෤ 𝑝 スケール因子𝑠が消えた
70. 70 96 ෤ 𝑝 = 𝑝/𝑠 スケールされた運動量の時間微分を考える 𝑑 ෤ 𝑝 𝑑𝑡 = 1 𝑠 𝑑𝑝 𝑑𝑡 − 𝑝 𝑠2 𝑑𝑠 𝑑𝑡 = − 1 𝑠 𝜕𝐻 𝜕𝑞 − 𝑝 𝑠2 𝜋 𝑄 ≡ 𝜁 𝑑 ෤ 𝑝 𝑑 ǁ 𝑡 = − 𝜕𝐻 𝜕𝑞 − ෤ 𝑝𝜁 𝑑 ǁ 𝑡 = 𝑑𝑡/𝑠であるから スケール因子𝑠が消えた
71. 71 96 仮想粒子の運動量の時間微分を考える 𝜁 ≡ 𝜋 𝑄 𝑑𝜁 𝑑𝑡 ≡ 1 𝑄𝑠 ෤ 𝑝 𝜕𝐻 𝜕𝑞 − 𝑘𝑇0 ※計算省略 𝑑 ǁ 𝑡 = 𝑑𝑡/𝑠であるから 𝑑𝜁 𝑑 ǁ 𝑡 ≡ 1 𝑄 ෤ 𝑝 𝜕𝐻 𝜕𝑞 − 𝑘𝑇0 スケール因子𝑠が消えた
72. 72 96 スケールされた物理量 ෤ 𝑝, ǁ 𝑡をあらためて𝑝, 𝑡と書くと 𝑑𝑞 𝑑𝑡 = 𝜕𝐻 𝜕𝑝 𝑑𝑝 𝑑𝑡 = − 𝜕𝐻 𝜕𝑞 − 𝑝𝜁 𝑑𝜁 𝑑𝑡 = 1 𝑄 𝑝 𝜕𝐻 𝜕𝑝 − 𝑘𝑇 この運動方程式による温度制御法をNosé-Hoover法と呼ぶ 時間スケール因子𝑠が消え、運動空間と観測空間が一致した
73. 73 96
74. 74 96 𝑝 𝜕𝐻 𝜕𝑝 = 𝑘𝑇 であるから、摩擦係数のダイナミクスは 𝑑𝑝 𝑑𝑡 = − 𝜕𝐻 𝜕𝑞 − 𝑝𝜁 𝑑𝜁 𝑑𝑡 = 𝑘 𝑄 𝑇 − 𝑇0 現在の温度が目標温度より高い → 摩擦係数𝜁が大きくなる → 運動量が小さくなる → 温度が下がる 摩擦項 現在の温度が目標温度より低い場合は摩擦係数が負になる →温度が制御される
75. 75 96 • ハミルトンの運動方程式に仮想粒子を追加し、 空間スケールを制御することで圧力を制御す るのがAndersenの圧力制御 • ハミルトンの運動方程式に仮想粒子を追加し、 時間スケールを制御することで温度を制御す るのが能勢の温度制御 • 変数変換にスケール因子を消去し、時間発展 系と観測系を一致させた手法がNosé-Hoover 法
76. 76 96 Nosé-Hoover法は温度が制御できるだけでなく 定常分布が厳密なカノニカル分布になる 𝐸 𝑓(𝐸) NVEアンサンブル 𝐸 𝑓(𝐸) NVTアンサンブル この揺らぎも含めて正しい分布になる
77. 77 96 𝑝 𝑞 O 点のダイナミクス 𝑝 𝑞 O 分布のダイナミクス 𝑑 𝑑𝑡 𝑝 𝑞 = −𝜕𝑞 𝐻 𝜕𝑝 𝐻 𝜕𝑓 𝜕𝑡 = −∇ ⋅ (𝜌 Ԧ 𝑣)
78. 78 96 運動方程式：位相空間に速度場を定義 初期条件：位相空間にトレーサーを置くこと 方程式の解：トレーサーの軌跡 𝑞 𝑝 O
79. 79 96 𝑝 𝑞 O 速度場に一つトレーサーを 置くと軌跡を追える 速度場に多数のトレーサーを 置くと分布を追える 𝑝 𝑞 O 点が動く 分布が動く
80. 80 96 速度場 密度場 速度が変化するところで密度が変化する
81. 81 96 𝑥 𝑥 + ℎ 𝑣(𝑥)𝜌(𝑥) 𝑣(𝑥 + ℎ)𝜌(𝑥 + ℎ) 𝜌(𝑥) 𝑣(𝑥 + ℎ) 速度場 密度場 𝑣(𝑥) 𝜌(𝑥 + ℎ) 注目領域
82. 82 96 𝑥 𝑥 + ℎ 𝑣(𝑥)𝜌(𝑥) 𝑣(𝑥 + ℎ)𝜌(𝑥 + ℎ) 注目領域 入ってくる量 出ていく量 注目領域の量の変化＝入ってくる量ー出ていく量 Δ𝜌 = 𝑣 𝑥 𝜌(𝑥) − 𝑣 𝑥 + ℎ 𝜌 𝑥 + ℎ 𝜕𝜌 𝜕𝑡 = −∇ ⋅ (𝜌 Ԧ 𝑣) 連続極限 連続の式
83. 83 96 分布関数𝑓(𝑝, 𝑞)：時刻𝑡、場所(𝑝, 𝑞)における系の存在確率 運動方程式は位相空間に速度場 Ԧ 𝑣(𝑝, 𝑞)を作る マクロな条件(エネルギーや体積といった熱力学変数) が等しい多数の系(統計集団)を用意したとき、そのミ クロな状態が(𝑝, 𝑞)であるような確率密度 確率の保存から、以下の連続の式が成り立つ 𝜕𝑓 𝜕𝑡 = −∇ ⋅ (𝑓 Ԧ 𝑣)
84. 84 96 ハミルトンダイナミクスが作る速度場 Ԧ 𝑣 = −𝜕𝑞 𝐻 𝜕𝑝 𝐻 ∇= 𝜕𝑝 𝜕𝑞 𝜕𝑓 𝜕𝑡 = −∇ ⋅ (𝑓 Ԧ 𝑣) = − 𝜕 𝜕𝑝 −𝑓𝜕𝑞 𝐻 + 𝜕 𝜕𝑞 (𝑓𝜕𝑝 𝐻) ナブラ演算子 = − − 𝜕𝑓 𝜕𝑝 𝜕𝐻 𝜕𝑞 − 𝑓 𝜕2𝐻 𝜕𝑝𝜕𝑞 + 𝜕𝑓 𝜕𝑝 𝜕𝐻 𝜕𝑝 + 𝑓 𝜕2𝐻 𝜕𝑝𝜕𝑞 キャンセル
85. 85 96 𝜕𝑓 𝜕𝑡 = 𝜕𝑓 𝜕𝑝 𝜕𝐻 𝜕𝑞 − 𝜕𝑓 𝜕𝑞 𝜕𝐻 𝜕𝑝 分布関数𝑓が𝐻のみの関数なら 𝜕𝑓 𝜕𝑝 = 𝑑𝑓 𝑑𝐻 𝜕𝐻 𝜕𝑝 𝜕𝑓 𝜕𝑞 = 𝑑𝑓 𝑑𝐻 𝜕𝐻 𝜕𝑞 代入すると 𝜕𝑓 𝜕𝑡 = 0 ハミルトンダイナミクスは分布関数を変化させない →ミクロカノニカルアンサンブル
86. 86 96 Nosé-Hooverの運動方程式が作る速度場 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝑝𝜁 ሶ 𝜁 = 1 𝑄 𝑝 𝜕𝐻 𝜕𝑝 − 𝑘𝑇 Ԧ 𝑣 = ሶ 𝑝 ሶ 𝑞 𝜁 ሶ 分布関数𝑓(𝑝, 𝑞, 𝜁)について連続の式が成り立つ 𝜕𝑓 𝜕𝑡 = −∇ ⋅ (𝑓 Ԧ 𝑣)
87. 87 96 定常分布𝑓eq を考える 𝜕𝑓eq 𝜕𝑡 = 0 ∇ ⋅ 𝑓eq Ԧ 𝑣 = 0 ∇ ⋅ 𝑓eq Ԧ 𝑣 = 𝜕 𝜕𝑝 𝑓eq ሶ 𝑝 + 𝜕 𝜕𝑞 𝑓eq ሶ 𝑞 + 𝜕 𝜕𝜁 𝑓eq ሶ 𝜁 = 0 定常分布𝑓eq が存在するなら、この偏微分方程式を満たす
88. 88 96 𝜕 𝜕𝑝 𝑓eq ሶ 𝑝 + 𝜕 𝜕𝑞 𝑓eq ሶ 𝑞 + 𝜕 𝜕𝜁 𝑓eq ሶ 𝜁 = 0 ሶ 𝑞 = 𝜕𝐻 𝜕𝑝 ሶ 𝑝 = − 𝜕𝐻 𝜕𝑞 − 𝑝𝜁 ሶ 𝜁 = 1 𝑄 𝑝 𝜕𝐻 𝜕𝑝 − 𝑘𝑇 定常分布が満たすべき式 Nosé-Hoover法の運動方程式 𝑓eq ∝ exp −𝛽𝐻 − 𝛽 2 𝜁2
89. 89 96 我々が興味ある分布は 𝑓0 𝑝, 𝑞 なので 追加された自由度𝜁を消去する 𝑓0 𝑝, 𝑞 ≡ න 𝑓eq (𝑝, 𝑞, 𝜁)𝑑𝜁 𝑓eq ∝ exp −𝛽𝐻 − 𝛽 2 𝜁2 であったから 𝑓0 ∝ exp −𝛽𝐻 もとの世界で見ると、カノニカル分布が実現する
90. 90 96 𝑝 𝑞 O 𝑝 𝜁 O ハミルトンダイナミクス由来の流れ 熱浴由来の流れ ሶ 𝑞 = 𝜕𝑝 𝐻 ሶ 𝑝 = −𝜕𝑞 𝐻 − 𝑝𝜁 ሶ 𝜁 = 𝑄−1 𝑝𝜕𝑝 𝐻 − 𝑘𝑇
91. 91 96 • Nosé-Hoover法は、ハミルトンダイナミクス に追加自由度を追加して温度制御をする手法 • 追加自由度を消去し、もとの世界の位相空間 に射影すると、分布関数が厳密にカノニカル 分布となる • ハミルトンダイナミクスが位相空間に作る流 れが非圧縮流になるのに対して、 Nosé- Hoover法が作る流れは圧縮流になる Nosé-Hoover法は、単に温度制御をする手法にとどまらず、 拡張ハミルトンダイナミクスという分野を生み出した
92. 92 96 𝑑𝑝 𝑑𝑡 = − 𝜕𝐻 𝜕𝑞 − 𝑝𝜁 𝑑𝜁 𝑑𝑡 = 𝑘 𝑄 𝑇 − 𝑇0 Nosé-Hoover法は「全体の温度」 しかチェックしていない 𝑇 = 2 3𝑁𝑘 ෍ 𝑖 𝑝𝑖 2 2𝑚 温度は速度の分散で求める Flying icecube問題 高速ですれ違う氷二つの系は「温度が高い」といえるか？
93. 93 96 系が相分離していると、温度が一様にならない B原子の温度 A原子の温度 全体の温度
94. 94 96 温度の時間発展 スペクトル 温度制御なし 温度制御あり 温度制御あり 温度制御なし 非自明なピーク 時間 𝑝 𝜁 O 熱浴由来の流れの「回転」を反映し、 系に熱浴由来の振動が加わる
95. 95 96 • Nosé-Hoover法は、系がエルゴード的かつ定 常状態であれば厳密にカノニカル分布を実現 する→条件によっては正しく温度制御ができ ない • 相分離するような系では、それぞれの部分系 が異なる温度に落ち着くことがある • 熱浴由来のエネルギー振動が入るため、ダイ ナミクスが信頼できないことがある • 別の手法であるLangevin熱浴ではこの問題は おきない(ただし温度収束が遅い) 「手法」は前提条件を理解して使わないと危険
96. 96 96 • 物理量の定義は難しい • 温度や圧力の定義は自明ではない • 数値計算では、支配方程式によりa priori な変数と観測量が異なる • ハミルトンの運動方程式はNVEアンサン ブルを実現する • 運動方程式を拡張することで、温度や圧 力を制御できる • 手法の性質をよく知らずに使うのは危険