シミュレーション工学
196分子動力学法(2) 温度制御と圧力制御シミュレーション工学慶應義塾大学大学院理工学研究科基礎理工学専攻物理情報専修渡辺宙志
View Slide
296• 分子動力学法における温度や圧力の定義• 温度や圧力制御アルゴリズム• 温度や圧力といった物理量の定義が、必ずしも自明ではないということを知る本講義の目的物理量の定義と制御• 分子動力学法における温度や圧力の定義は非自明• 温度とは？圧力とは？
396肌寒い お風呂がぬるい アイロンが熱い気体の温度 液体の温度 固体の温度これら全てに共通する「温度」とはなんだろう？どうやって数値化しているのだろう？
496ガラス温度計液体の膨張による体積変化温湿度計(バイメタル方式)熱膨張率の異なる金属を貼り合わせたものデジタル体温計(サーミスタ方式)温度により抵抗が変わる物質を利用
596熱電対二種類の金属の熱電能の差を起電力に変える(ゼーベック効果)赤外線温度計物体から放出される赤外線のエネルギーを測定これら全ては、現実物質の温度変化に対する性質の変化を利用基準となる温度を使って較正する必要があるそもそも温度とは？
696変数(Variable)我々がa prioriに認める物理量支配方程式に直接でてくる量観測量(Observable)支配方程式には含まれない量変数から導かれる量何が変数で何が観測量かは支配方程式により異なる
796熱伝導方程式𝜕𝑇𝜕𝑡= 𝜅𝜕2𝑇𝜕𝑥2位置𝑥、時刻𝑡、温度𝑇(𝑥, 𝑡)は全て変数(variable)この支配方程式を認めた時点で「この世界には温度というものがあり、この方程式に従う」と宣言
896ナビエ・ストークス方程式(非圧縮)𝜕 Ԧ𝑣𝜕𝑡+ Ԧ𝑣 ⋅ ∇ = −1𝜌∇𝑃 + 𝜈Δ Ԧ𝑣 + Ԧ𝐹速度場 Ԧ𝑣と圧力場𝑃が変数(variable)この支配方程式を認めた時点で、圧力という量の存在を認めている。「圧力とは何か？」という問いが意味をもたない
996ハミルトンの運動方程式𝑑𝑝𝑖𝑑𝑡= −𝜕𝐻𝜕𝑞𝑖𝑑𝑞𝑖𝑑𝑡=𝜕𝐻𝜕𝑝𝑖エネルギー𝐻、座標𝑞𝑖、運動量𝑝𝑖が変数(variable)温度𝑇や圧力𝑃は全て観測量(observable)この世界では、温度や圧力は定義しなければならない
1096示量変数(extensive variable)全く同じ系を二つつけた時、二倍になる量体積𝑉、エントロピー𝑆、物質量𝑁示強変数(intensive variable)全く同じ系を二つつけた時、値が変わらない量圧力𝑃、温度𝑇、化学ポテンシャル𝜇
1196示量変数 示強変数体積𝑉 圧力𝑃かけたらエネルギーになる量をお互いに共役と呼ぶ温度𝑇エントロピー𝑆物質量𝑁 化学ポテンシャル𝜇
1296𝑈(𝑆, 𝑉, 𝑁)𝐻(𝑆, 𝑃, 𝑁) 𝐹(𝑇, 𝑉, 𝑁)𝐺(𝑇, 𝑃, 𝑁) Ω(𝑇, 𝑉, 𝜇)𝐹 = 𝑈 − 𝑇𝑆𝐻 = 𝑈 + 𝑃𝑉𝐺 = 𝐻 − 𝑇𝑆𝐺 = 𝐹 + 𝑃𝑉Ω = 𝐹 − 𝜇𝑁共役な量＝ルジャンドル変換で互いに変換する量全て示量変数
1396示量変数 示強変数体積𝑉 圧力𝑃かけたらエネルギーになる量をお互いに共役と呼ぶ温度𝑇エントロピー𝑆物質量𝑁 化学ポテンシャル𝜇示量変数と示強変数、どちらが「基本的な量」だろうか？
1496a prioriな量とは、我々が未定義で使う量のこと熱力学では、共役な量のどちらかがa prioriな量数値計算では、支配方程式に含まれる変数それ以外の物理量は全てa prioriな量から定義される分子動力学法において、温度や圧力を定義しよう
1596温度𝑇は運動エネルギー𝐾に比例𝑇 =2𝐾3𝑁𝑘𝑁：粒子数𝑘：ボルツマン定数運動エネルギー𝐾は運動量の二乗和𝐾 = 𝑖𝑝𝑖22𝑚以上から𝑇 =23𝑁𝑘𝑖𝑝𝑖22𝑚
1696𝑇 =23𝑁𝑘𝑖𝑝𝑖22𝑚こちらは全部知っているのでこれが計算できる𝐾 = 𝑖𝑝𝑖22𝑚← これは単なる定義𝑇 =2𝐾3𝑁𝑘← これはどこからきたのか？
1796(𝑝, 𝑞)空間の分布関数𝑓(𝑝, 𝑞)を考えるある点(𝑝, 𝑞)におけるエネルギー𝐻(𝑝, 𝑞)を定義する規格化条件 න 𝑓𝑑𝑝𝑑𝑞 = 1内部エネルギー න 𝐻𝑓𝑑𝑝𝑑𝑞 = 𝑈
1896この系のエントロピーを以下のように定義する𝑆 = −𝑘 න 𝑓 log 𝑓 𝑑𝑝𝑑𝑞න 𝑓𝑑𝑝𝑑𝑞 = 1න 𝐻𝑓𝑑𝑝𝑑𝑞 = 𝑈以下の条件を満たしつつ、エントロピーを最大化する分布関数を求めたい規格化条件エネルギーの期待値が𝑈
1996ラグランジュの未定定数法より(※)𝐼 = 𝛼 න 𝑓𝑑𝑝𝑑𝑞 + 𝛽 න 𝐻𝑓𝑑𝑝𝑑𝑞 + න 𝑓 log 𝑓 𝑑𝑝𝑑𝑞※ あとの便利のために符号をかえたりボルツマン定数を吸収させたりしている汎関数微分から𝛿𝐼 = 0 𝛼 + 𝛽𝐻 + 1 + log 𝑓 = 0以上から𝑓 = 𝑍−1exp(−𝛽𝐻) 𝑍 = exp(𝛼 + 1)カノニカル分布 分配関数(規格化条件から決まる)
2096𝑆 = −𝑘 න 𝑓 log 𝑓 𝑑𝑝𝑑𝑞 エントロピーの定義𝑓 = 𝑍−1exp(−𝛽𝐻) log 𝑓 = −𝛽𝐻 − log 𝑍𝑆 = 𝛽𝑘 න 𝐻𝑓𝑑𝑝𝑑𝑞 + 𝐶= 𝑈カノニカル分布を代入𝑑𝑆𝑑𝑈= 𝛽𝑘
2196𝑑𝑆𝑑𝑈= 𝛽𝑘𝑑𝑆𝑑𝑈=1𝑇さきほど求めた関係式熱力学関係式逆温度と温度が結びついた𝛽 =1𝑘𝑇
2296全ての物理量は位相空間の座標(𝑝, 𝑞)の関数物理量𝐴(𝑝, 𝑞)の期待値𝐴 ≡ න 𝐴𝑓𝑑𝑝𝑑𝑞位相空間の各点(𝑝, 𝑞)で定義される物理量𝐴(𝑝, 𝑞)に、その点における分布関数𝑓(𝑝, 𝑞)をかけて全位相空間において平均せよ
2396分布がカノニカル分布である時 という𝑝𝜕𝐻𝜕𝑝𝑝𝜕𝐻𝜕𝑝= 𝑍−1 න 𝑝𝜕𝐻𝜕𝑝exp −𝛽𝐻 𝑑𝑝𝑑𝑞𝐴 ≡ න 𝐴𝑓𝑑𝑝𝑑𝑞𝑓 = 𝑍−1exp(−𝛽𝐻)物理量の期待値を考える物理量𝐴の期待値𝐴 = 𝑝𝜕𝑝𝐻期待値を求めたい物理量 カノニカル分布
2496𝜕𝜕𝑝exp −𝛽𝐻 = −𝛽𝜕𝐻𝜕𝑝exp(−𝛽𝐻) であるから𝑝𝜕𝐻𝜕𝑝= 𝑍−1 න 𝑝𝜕𝐻𝜕𝑝exp −𝛽𝐻 𝑑𝑝𝑑𝑞= 𝑍−1 න −𝑝𝛽𝜕𝜕𝑝exp(−𝛽𝐻) 𝑑𝑝𝑑𝑞=1𝛽𝑍−1 න exp(−𝛽𝐻) 𝑑𝑝𝑑𝑞部分積分= 𝑍= 1/𝛽
2596𝑝𝜕𝐻𝜕𝑝=1𝛽= 𝑘𝑇ハミルトニアンが𝐻 = 𝑝22𝑚+ 𝑉(𝑞)の形をしている場合𝑝𝜕𝐻𝜕𝑝=𝑝2𝑚= 𝑘𝑇 𝑝22𝑚=12𝑘𝑇 等分配則3次元多粒子系の場合𝑖𝑝𝑖22𝑚=32𝑁𝑘𝑇
2696𝑇 =2𝐾3𝑁𝑘𝑖𝑝𝑖22𝑚=32𝑁𝑘𝑇カノニカル分布の部分積分から𝑖𝑝𝑖22𝑚= 𝐾運動エネルギーの定義分子動力学法における(解析力学における)温度の定義
2796𝑝𝜕𝐻𝜕𝑝=1𝛽= 𝑘𝑇以下の量が温度になったのは部分積分したから全く同じ導出で以下の等式も成り立つ𝑞𝜕𝐻𝜕𝑞=1𝛽= 𝑘𝑇運動温度状態温度平衡状態では両者は同じ温度を与えるが、非平衡状態では一般に一致しない
2896運動温度時間状態温度両者が一致していなければ緩和不足(非平衡状態)平衡状態非平衡状態
2996• 温度が運動エネルギーに比例するのは等分配則のため• 逆温度はエントロピーを最大化する際のラグランジュの未定定数であり、熱力学関係を要請することで温度と結びつく• 等分配則はカノニカル分布をする系において𝑝𝜕𝑝𝐻の形の物理量の期待値が温度に比例することによる(運動温度)• 同様に𝑞𝜕𝑞𝐻も温度に比例する(状態温度)• 運動温度と状態温度は、平衡状態では一致するが、非平衡状態では必ずしも一致しない
3096圧力：粒子が壁を押す力壁がない系(周期的境界条件)で圧力は定義できるか？負の圧力はあり得るか？
3196𝐺 = 𝑖𝑝𝑖𝑞𝑖この系で以下の量を考える3次元N粒子系{𝑝𝑖, 𝑞𝑖}を考えるx,y,z座標をまとめてiのインデックスに押し込める(𝑖 = 1,2, ⋯ , 3𝑁)この量をビリアル(virial)と呼び、クラウジウスにより導入されたこの量から圧力を定義する
3296𝐺 = 𝑖𝑝𝑖𝑞𝑖ሶ𝐺 = 𝑖(𝑝𝑖ሶ𝑞𝑖+ ሶ𝑝𝑖𝑞𝑖)𝑝𝑖ሶ𝑞𝑖=𝑝𝑖2𝑚= 𝑘𝑇時間微分エネルギーの等分配則𝑖𝑝𝑖ሶ𝑞𝑖= 3𝑁𝑘𝑇 理想気体からの寄与
3396ሶ𝐺 = 𝑖(𝑝𝑖ሶ𝑞𝑖+ ሶ𝑝𝑖𝑞𝑖)ሶ𝑝𝑖= 𝑓𝑖+ 𝐹𝑖運動量の時間変化＝粒子にかかる力外力(壁から粒子iに働く力)粒子間に働く力(他粒子から粒子iに働く力の合計)𝑓𝑖𝐹𝑖
3496𝑖𝑞𝑖𝑓𝑖= 𝑖<𝑗𝑞𝑖𝑗𝑓𝑖𝑗粒子間に働く力を作用反作用を使って整理𝑞𝑖𝑗= 𝑞𝑖− 𝑞𝑗𝑓𝑖𝑗= 𝑓𝑖− 𝑓𝑗外力をガウスの定理を使って整理(※)𝑖𝑞𝑖𝐹𝑖= − 𝑃 න𝜕𝑉Ԧ𝑟 ⋅ 𝑛𝑑𝐴 = −𝑃 න𝑉∇ ⋅ 𝑛𝑑𝑉 = −3𝑃𝑉まとめると 𝑖𝑞𝑖ሶ𝑝𝑖= 𝑖<𝑗𝑞𝑖𝑗𝑓𝑖𝑗− 3𝑃𝑉※ 計算の詳細は今は気にしなくて良いです
3596ሶ𝐺 = 𝑖(𝑝𝑖ሶ𝑞𝑖+ ሶ𝑝𝑖𝑞𝑖) = 3𝑁𝑘𝑇 + 𝑖<𝑗𝑞𝑖𝑗𝑓𝑖𝑗− 3𝑃𝑉ビリアルの時間平均がゼロであることを仮定𝑃𝑉 = 3𝑁𝑘𝑇 +13𝑖<𝑗𝑞𝑖𝑗𝑓𝑖𝑗分子動力学法における圧力の定義理想気体の寄与 相互作用からの寄与
3696• 粒子に働く力を粒子間力と外力にわけたが、周期的境界条件など、壁がない系では圧力は定義できないのか？• ビリアルの時間微分の平均がゼロであるという仮定はどういう意味を持つのか？• そもそもビリアルってなに？熱力学関係式とハミルトニアンから圧力を導出
3796これから計算がごちゃごちゃするが全部追う必要はない分配関数から必要な物理量が全て求められる、という感覚だけ身に着ける
3896𝑃 = −𝜕𝐹𝜕𝑉𝑇以下の熱力学関係式から出発するヘルムホルツ自由エネルギーと分配関数𝐹 = −𝑘𝑇 log 𝑍𝑃 = 𝑘𝑇𝜕 log 𝑍𝜕𝑉=𝑘𝑇𝑍𝜕𝑍𝜕𝑉
3996𝜕𝑍𝜕𝑉この量を計算したい𝑍 = න exp −𝛽𝐻 𝑑Γ 分配関数は陽に体積に依存しない体積変化はハミルトニアンにのみ影響体積が変化したときのハミルトニアンの応答を調べる𝜕𝐻𝜕𝑉𝐿 𝛼𝐿サイズ𝛼倍体積𝛼3倍
4096𝑍 = න exp −𝛽𝐻 𝑑Γ𝑑Γ = ෑ𝑖𝑑𝑝𝑖𝑑𝑞𝑖分配関数とハミルトニアン𝜕𝑍𝜕𝑉= න −𝛽𝜕𝐻𝜕𝑉exp −𝛽𝐻 𝑑Γ𝑃 =𝑘𝑇𝑍𝜕𝑍𝜕𝑉Vで微分に代入𝑃 =𝑘𝑇𝑍න −𝛽𝜕𝐻𝜕𝑉e−𝛽𝐻 𝑑Γ = −𝜕𝐻𝜕𝑉
4196𝑃 = −𝜕𝐻𝜕𝑉を計算するために空間を一様に𝛼倍にスケールする𝑞𝑖→ 𝛼𝑞𝑖𝑝𝑖→ 𝑝𝑖/𝛼ሶ𝑞𝑖→ 𝛼 ሶ𝑞𝑖𝑝𝑖=𝜕𝐿𝜕 ሶ𝑞𝑖座標とスケーリングが異なることに注意𝑉 → 𝛼3𝑉 𝑑𝑉𝑑𝛼→ 3𝛼2𝑉座標と運動量の変化体積変化
4296もともとのハミルトニアン𝐻 = 𝑖𝑝𝑖22𝑚+ 𝑖<𝑗Φ(𝑞𝑖𝑗) = 𝐾 𝑝𝑖+ 𝑉({𝑝𝑖})スケールされたハミルトニアン𝐻(𝛼) = 𝑖𝑝𝑖22𝛼2𝑚+ 𝑖<𝑗Φ(𝛼𝑞𝑖𝑗)𝜕𝐻𝜕𝑉= lim𝛼→1𝜕𝐻𝜕𝛼𝑑𝛼𝑑𝑉=13𝑉lim𝛼→1𝜕𝐻𝜕𝛼ハミルトニアンの体積微分運動項 ポテンシャル項𝑑𝑉𝑑𝛼→ 3𝛼2𝑉
4396lim𝛼→1𝜕𝜕𝛼𝑖𝑝𝑖22𝛼2𝑚= −2𝐾 = −3𝑁𝑘𝑇運動項ポテンシャル項lim𝛼→1𝜕Φ 𝛼𝑞𝑖𝑗𝜕𝛼= Φ′ 𝑞𝑖𝑗𝑞𝑖𝑗= −𝑓𝑖𝑗𝑞𝑖𝑗まとめると𝜕𝐻𝜕𝑉=13𝑉lim𝛼→1𝜕𝐻𝜕𝛼=13𝑉−3𝑁𝑘𝑇 + 𝑖<𝑗𝑞𝑖𝑗𝑓𝑖𝑗
4496𝑃 = −𝜕𝐻𝜕𝑉 であり𝜕𝐻𝜕𝑉=13𝑉−3𝑁𝑘𝑇 + 𝑖<𝑗𝑞𝑖𝑗𝑓𝑖𝑗であるから𝑃𝑉 = 3𝑁𝑘𝑇 +13𝑖<𝑗𝑞𝑖𝑗𝑓𝑖𝑗ビリアル定理と同じ公式が得られた
4596• 分子動力学法の世界には、一般化運動量、一般化座標、そしてハミルトニアンしか存在しない→圧力は与えられるものではなく定義するもの• 熱力学関係式から出発し、分配関数から圧力の表式をもとめた• 外力や内力、ビリアルといった概念を使わずに、系のスケーリングへの応答だけから圧力を定義→周期的境界条件でも適用可能分配関数から必要な物理量をなんでも求めることができる
4696ハミルトンの運動方程式𝑑𝑝𝑖𝑑𝑡= −𝜕𝐻𝜕𝑞𝑖𝑑𝑞𝑖𝑑𝑡=𝜕𝐻𝜕𝑝𝑖エネルギーが変化しない通常は体積も変化しない熱の出入りもないミクロカノニカルアンサンブルになる(NVE)時間発展について
4796ミクロカノニカルアンサンブル粒子数N、体積V、エネルギーEが一定NVEアンサンブル熱浴(heatbath)カノニカルアンサンブル粒子数N、体積V、温度Tが一定NVTアンサンブルエネルギーが揺らぐ
4896𝐸𝑓(𝐸)エネルギーが完全に固定で揺らがないNVEアンサンブル𝐸𝑓(𝐸)エネルギーはある幅をもって揺らぐNVTアンサンブル
4996実際の物質はNが非常に大きいため、NVTアンサンブルでも事実上エネルギーは揺らがない𝐸𝑓(𝐸)NVEアンサンブル𝐸𝑓(𝐸)NVTアンサンブル𝑁 ≫ 1の時
5096我々が知りたい量は温度依存性であることが多い300Kにおける水の圧力は？300Kにおけるポリエチレンの弾性率は？ミクロカノニカル分布とカノニカル分布は本質的に変わらないので、興味ある温度𝑇に対応するエネルギー𝑈(𝑇)を持った系でNVEシミュレーションをすれば良い
5196温度𝑇における内部エネルギー𝑈(𝑇)は物質により異なる金の比熱：130 [J/kg ℃]水の比熱：4182 [J/kg ℃]同じ温度でも、温まり易い物質はエネルギーが低く温まり難い物質はエネルギーが高い※ いずれも20℃における値
5296ある温度𝑇におけるエネルギー𝑈(𝑇)を知りたい比熱が高いほど同じ温度でのエネルギーは高い温度変化に対するエネルギーの変化率が比熱𝐶 =𝜕𝑈𝜕𝑇𝑈 𝑇 = න0𝑇𝐶𝑑𝑇それを積分したものが全エネルギーこの量を知りたい
5396𝑈 𝑇 ≡ 𝐻 = 𝑍−1 න 𝐻 exp −𝛽𝐻 𝑑Γ内部エネルギーの定義𝑍 = න exp −𝛽𝐻 𝑑Γ分配関数の定義𝜕𝑍𝜕𝛽= −𝛽 න 𝐻 exp −𝛽𝐻 𝑑Γ = −𝛽𝑍 𝐻𝑈 𝑇 = −𝜕 log 𝑍𝜕𝛽𝛽で微分
5496𝑈 𝑇 = −𝜕 log 𝑍𝜕𝛽内部エネルギーの温度依存性がわかる分配関数がわかる=問題が厳密に解けている厳密に解けてない問題を解きたいのだから一般に𝑈(𝑇)はわからないフィードバック制御による温度調整
55961. 全エネルギー一定の計算を行う2. 運動エネルギーから温度を計算する3. 目標温度との差を見てエネルギーを調整4. 1.～3.を繰り返す𝑇 =2𝐾3𝑁𝑘
5696圧縮率：圧力が増えた時、どれだけ体積が減るか𝜅 = −1𝑉𝜕𝑉𝜕𝑃熱力学的に安定な系は、圧縮率が正でなくてはならない平衡状態でピストンを少し押す圧縮率が正なら押し返される(安定)圧縮率が負ならさらに引き込まれる(不安定)圧縮率が正体積を増やせば圧力が下がる体積を減らせば圧力が上がる
57961. 全エネルギー一定の計算を行う2. ビリアル定理から圧力を計算する3. 目標圧力との差を見て体積を調整4. 1.～3.を繰り返す数値シミュレーションにおける体積変化とは？周期的境界条件ではどうするか？
5896スケールされたハミルトニアン𝐻(𝛼) = 𝑖𝑝𝑖22𝛼2𝑚+ 𝑖<𝑗Φ(𝛼𝑞𝑖𝑗)スケール因子𝛼を一般化座標に含めたハミルトニアン𝐻 = 𝑖𝑝𝑖22𝛼2𝑚+ 𝑖<𝑗Φ(𝛼𝑞𝑖𝑗) +𝜋22𝑀+ 𝑃0𝛼3このハミルトニアンに従う系は、目標圧力𝑃0に制御される(Andersenのハミルトニアン)
5996𝐻 = 𝑖𝑝𝑖22𝛼𝑚+ 𝑖<𝑗Φ(𝛼𝑞𝑖𝑗) +𝜋22𝑀+ 𝑃0𝛼3𝜋と𝛼が共役な変数であると考える𝑑𝛼𝑑𝑡=𝜕𝐻𝜕𝜋=𝜋𝑀𝛼が座標、𝜋が運動量、𝑀が質量となる仮想粒子𝑑𝜋𝑑𝑡= −𝜕𝐻𝜕𝛼=1𝛼3𝑖𝑝𝑖2𝑚− 𝑖<𝑗Φ′ 𝛼𝑞𝑖𝑗𝑞𝑖𝑗− 3𝛼2𝑃0仮想粒子に働く力
6096スケールされた運動量と座標を導入して書き直す𝑝𝑖= 𝑝𝑖/𝛼 𝑞𝑖= 𝛼𝑞𝑑𝜋𝑑𝑡=1𝛼𝑖𝑝𝑖2𝑚−1𝛼𝑖<𝑗Φ′ 𝑞𝑖𝑗𝑞𝑖𝑗− 3𝛼2𝑃0= 2𝐾= − ሚ𝑓𝑖𝑗= 3𝑉/𝛼=1𝛼2𝐾 + 𝑖<𝑗ሚ𝑓𝑖𝑗𝑞𝑖𝑗− 3𝑃0𝑉= 3𝑃𝑉=3𝑉𝛼𝑃 − 𝑃0
6196𝑑𝛼𝑑𝑡=𝜕𝐻𝜕𝜋=𝜋𝑀𝑑𝜋𝑑𝑡=3𝑉𝛼𝑃 − 𝑃0𝐻 = 𝑖𝑝𝑖22𝛼𝑚+ 𝑖<𝑗Φ(𝛼𝑞𝑖𝑗) +𝜋22𝑀+ 𝑃0𝛼3追加された仮想粒子現在の圧力が設定圧力より高ければ加速そうでなければ減速仮想粒子の運動量が正なら膨張負なら収縮圧縮率が正体積を増やせば圧力が下がる体積を減らせば圧力が上がる最終的に圧力が目標圧力に収束する
6296𝐻 = 𝑖𝑝𝑖22𝛼𝑚+ 𝑖<𝑗Φ(𝛼𝑞𝑖𝑗) +𝜋22𝑀+ 𝑃0𝛼3運動を支配するハミルトニアン𝑝𝑖= 𝑝𝑖/𝛼 𝑞𝑖= 𝛼𝑞スケールされた運動量と座標圧力が制御されるのは、スケールされた世界「運動方程式に従う座標と運動量」と「我々が観測する座標と運動量」が異なる
6396• 熱力学的に安定な系では、体積が増えると圧力は下がり、体積が減ると圧力が上がる• ハミルトニアンにスケールをつかさどる仮想粒子を追加したハミルトニアンを考える(Andersenのハミルトニアン)• Andersenのハミルトニアンは、現在の圧力が目標圧力より高いと系が膨張し、低いと系が収縮するダイナミクスを記述する• スケールされた座標と運動量を観測すると、圧力が制御されているように見える• 「運動に従う変数」と「観測する変数」が異なる
6496𝐻 = 𝑖𝑝𝑖22𝛼𝑚+ 𝑖<𝑗Φ(𝛼𝑞𝑖𝑗) +𝜋22𝑀+ 𝑃0𝛼3Andersenのハミルトニアン：仮想粒子による圧力制御能勢のハミルトニアン：仮想粒子による温度制御𝐻 = 𝑖𝑝𝑖22𝑠2𝑚+ 𝑖<𝑗Φ(𝑞𝑖𝑗) +𝜋22𝑄+ 3𝑁𝑘𝑇0log 𝑠追加された仮想粒子スケールするのは運動量のみ
6596𝐻 = 𝑖𝑝𝑖22𝑠2𝑚+ 𝑖<𝑗Φ(𝑞𝑖𝑗) +𝜋22𝑄+ 3𝑁𝑘𝑇0log 𝑠𝑑𝜋𝑑𝑡= −𝜕𝐻𝜕𝑠=2𝑠𝑖𝑝𝑖22𝑠2𝑚−3𝑁𝑘𝑇0𝑠𝑝𝑖= 𝑝𝑖/𝑠スケールされた運動量を導入=1𝑠𝑖𝑝𝑖2𝑚−3𝑁𝑘𝑇0𝑠=3𝑁𝑘𝑠𝑇 − 𝑇0= 2𝐾 = 3𝑁𝑘𝑇
6696𝐻 = 𝑖𝑝𝑖22𝑠2𝑚+ 𝑖<𝑗Φ(𝑞𝑖𝑗) +𝜋22𝑄+ 3𝑁𝑘𝑇0log 𝑠𝑑𝑠𝑑𝑡=𝜕𝐻𝜕𝜋=𝜋𝑄𝑑𝜋𝑑𝑡=3𝑁𝑘𝑠𝑇 − 𝑇0現在の温度が設定温度より高ければ減速そうでなければ加速仮想粒子の運動量が正なら減速負なら加速現在の温度が目標温度より高い→ ሶ𝜋 > 0なので、いつか𝜋 > 0になる→ 𝑠が増える→ 𝑝𝑖が小さくなる→温度が下がる最終的に温度が目標温度に収束する
6796𝑝𝑖の世界時間𝑝𝑖の世界時間能勢のハミルトニアンは仮想時間をスケールすることで温度を調整している𝑠いちいち仮想時間を考えるのは面倒くさいNosé-Hoover法
6896温度制御前のハミルトニアン𝐻 𝑝, 𝑞 =𝑝22𝑚+ 𝑉(𝑞) ※簡単のために一自由度系を考える能勢のハミルトニアン𝐻𝑁= 𝐻 𝑝/𝑠, 𝑞 +𝜋22𝑄+ 𝑘𝑇 log 𝑠運動方程式𝑑𝑝𝑑𝑡= −𝜕𝐻𝑁𝜕𝑞𝑑𝑞𝑑𝑡=𝜕𝐻𝑁𝜕𝑝
6996𝑝 = 𝑝/𝑠スケールされた運動量と時間を導入ǁ𝑡 = 𝑡/𝑠𝑑𝑞𝑑𝑡=𝜕𝐻𝑁𝜕𝑝=𝜕𝐻(𝑝/𝑠, 𝑞)𝜕𝑝=𝜕𝐻( 𝑝, 𝑞)𝜕𝑝=𝜕𝐻𝜕 𝑝𝑑 𝑝𝑑𝑝=1𝑠𝜕𝐻𝜕 𝑝𝑑 ǁ𝑡 = 𝑑𝑡/𝑠であるから𝑑𝑞𝑑 ǁ𝑡=𝜕𝐻𝜕 𝑝スケール因子𝑠が消えた
7096𝑝 = 𝑝/𝑠スケールされた運動量の時間微分を考える𝑑 𝑝𝑑𝑡=1𝑠𝑑𝑝𝑑𝑡−𝑝𝑠2𝑑𝑠𝑑𝑡= −1𝑠𝜕𝐻𝜕𝑞−𝑝𝑠2𝜋𝑄≡ 𝜁𝑑 𝑝𝑑 ǁ𝑡= −𝜕𝐻𝜕𝑞− 𝑝𝜁𝑑 ǁ𝑡 = 𝑑𝑡/𝑠であるからスケール因子𝑠が消えた
7196仮想粒子の運動量の時間微分を考える 𝜁 ≡𝜋𝑄𝑑𝜁𝑑𝑡≡1𝑄𝑠𝑝𝜕𝐻𝜕𝑞− 𝑘𝑇0※計算省略𝑑 ǁ𝑡 = 𝑑𝑡/𝑠であるから𝑑𝜁𝑑 ǁ𝑡≡1𝑄𝑝𝜕𝐻𝜕𝑞− 𝑘𝑇0 スケール因子𝑠が消えた
7296スケールされた物理量 𝑝, ǁ𝑡をあらためて𝑝, 𝑡と書くと𝑑𝑞𝑑𝑡=𝜕𝐻𝜕𝑝𝑑𝑝𝑑𝑡= −𝜕𝐻𝜕𝑞− 𝑝𝜁𝑑𝜁𝑑𝑡=1𝑄𝑝𝜕𝐻𝜕𝑝− 𝑘𝑇この運動方程式による温度制御法をNosé-Hoover法と呼ぶ時間スケール因子𝑠が消え、運動空間と観測空間が一致した
7396
7496𝑝𝜕𝐻𝜕𝑝= 𝑘𝑇 であるから、摩擦係数のダイナミクスは𝑑𝑝𝑑𝑡= −𝜕𝐻𝜕𝑞− 𝑝𝜁𝑑𝜁𝑑𝑡=𝑘𝑄𝑇 − 𝑇0現在の温度が目標温度より高い→ 摩擦係数𝜁が大きくなる→ 運動量が小さくなる→ 温度が下がる摩擦項現在の温度が目標温度より低い場合は摩擦係数が負になる→温度が制御される
7596• ハミルトンの運動方程式に仮想粒子を追加し、空間スケールを制御することで圧力を制御するのがAndersenの圧力制御• ハミルトンの運動方程式に仮想粒子を追加し、時間スケールを制御することで温度を制御するのが能勢の温度制御• 変数変換にスケール因子を消去し、時間発展系と観測系を一致させた手法がNosé-Hoover法
7696Nosé-Hoover法は温度が制御できるだけでなく定常分布が厳密なカノニカル分布になる𝐸𝑓(𝐸)NVEアンサンブル𝐸𝑓(𝐸)NVTアンサンブルこの揺らぎも含めて正しい分布になる
7796𝑝𝑞O点のダイナミクス𝑝𝑞O分布のダイナミクス𝑑𝑑𝑡𝑝𝑞 =−𝜕𝑞𝐻𝜕𝑝𝐻𝜕𝑓𝜕𝑡= −∇ ⋅ (𝜌 Ԧ𝑣)
7896運動方程式：位相空間に速度場を定義初期条件：位相空間にトレーサーを置くこと方程式の解：トレーサーの軌跡𝑞𝑝O
7996𝑝𝑞O速度場に一つトレーサーを置くと軌跡を追える速度場に多数のトレーサーを置くと分布を追える𝑝𝑞O点が動く 分布が動く
8096速度場密度場速度が変化するところで密度が変化する
8196𝑥 𝑥 + ℎ𝑣(𝑥)𝜌(𝑥) 𝑣(𝑥 + ℎ)𝜌(𝑥 + ℎ)𝜌(𝑥)𝑣(𝑥 + ℎ)速度場密度場𝑣(𝑥)𝜌(𝑥 + ℎ)注目領域
8296𝑥 𝑥 + ℎ𝑣(𝑥)𝜌(𝑥) 𝑣(𝑥 + ℎ)𝜌(𝑥 + ℎ)注目領域入ってくる量 出ていく量注目領域の量の変化＝入ってくる量ー出ていく量Δ𝜌 = 𝑣 𝑥 𝜌(𝑥) − 𝑣 𝑥 + ℎ 𝜌 𝑥 + ℎ𝜕𝜌𝜕𝑡= −∇ ⋅ (𝜌 Ԧ𝑣)連続極限連続の式
8396分布関数𝑓(𝑝, 𝑞)：時刻𝑡、場所(𝑝, 𝑞)における系の存在確率運動方程式は位相空間に速度場 Ԧ𝑣(𝑝, 𝑞)を作るマクロな条件(エネルギーや体積といった熱力学変数)が等しい多数の系(統計集団)を用意したとき、そのミクロな状態が(𝑝, 𝑞)であるような確率密度確率の保存から、以下の連続の式が成り立つ𝜕𝑓𝜕𝑡= −∇ ⋅ (𝑓 Ԧ𝑣)
8496ハミルトンダイナミクスが作る速度場Ԧ𝑣 =−𝜕𝑞𝐻𝜕𝑝𝐻∇=𝜕𝑝𝜕𝑞𝜕𝑓𝜕𝑡= −∇ ⋅ (𝑓 Ԧ𝑣)= −𝜕𝜕𝑝−𝑓𝜕𝑞𝐻 +𝜕𝜕𝑞(𝑓𝜕𝑝𝐻)ナブラ演算子= − −𝜕𝑓𝜕𝑝𝜕𝐻𝜕𝑞− 𝑓𝜕2𝐻𝜕𝑝𝜕𝑞+𝜕𝑓𝜕𝑝𝜕𝐻𝜕𝑝+ 𝑓𝜕2𝐻𝜕𝑝𝜕𝑞キャンセル
8596𝜕𝑓𝜕𝑡=𝜕𝑓𝜕𝑝𝜕𝐻𝜕𝑞−𝜕𝑓𝜕𝑞𝜕𝐻𝜕𝑝分布関数𝑓が𝐻のみの関数なら𝜕𝑓𝜕𝑝=𝑑𝑓𝑑𝐻𝜕𝐻𝜕𝑝𝜕𝑓𝜕𝑞=𝑑𝑓𝑑𝐻𝜕𝐻𝜕𝑞代入すると𝜕𝑓𝜕𝑡= 0ハミルトンダイナミクスは分布関数を変化させない→ミクロカノニカルアンサンブル
8696Nosé-Hooverの運動方程式が作る速度場ሶ𝑞 =𝜕𝐻𝜕𝑝ሶ𝑝 = −𝜕𝐻𝜕𝑞− 𝑝𝜁ሶ𝜁 =1𝑄𝑝𝜕𝐻𝜕𝑝− 𝑘𝑇Ԧ𝑣 =ሶ𝑝ሶ𝑞𝜁 ሶ分布関数𝑓(𝑝, 𝑞, 𝜁)について連続の式が成り立つ𝜕𝑓𝜕𝑡= −∇ ⋅ (𝑓 Ԧ𝑣)
8796定常分布𝑓eqを考える𝜕𝑓eq𝜕𝑡= 0 ∇ ⋅ 𝑓eqԦ𝑣 = 0∇ ⋅ 𝑓eqԦ𝑣 =𝜕𝜕𝑝𝑓eqሶ𝑝 +𝜕𝜕𝑞𝑓eqሶ𝑞 +𝜕𝜕𝜁𝑓eqሶ𝜁= 0定常分布𝑓eqが存在するなら、この偏微分方程式を満たす
8896𝜕𝜕𝑝𝑓eqሶ𝑝 +𝜕𝜕𝑞𝑓eqሶ𝑞 +𝜕𝜕𝜁𝑓eqሶ𝜁 = 0ሶ𝑞 =𝜕𝐻𝜕𝑝ሶ𝑝 = −𝜕𝐻𝜕𝑞− 𝑝𝜁 ሶ𝜁 =1𝑄𝑝𝜕𝐻𝜕𝑝− 𝑘𝑇定常分布が満たすべき式Nosé-Hoover法の運動方程式𝑓eq∝ exp −𝛽𝐻 −𝛽2𝜁2
8996我々が興味ある分布は 𝑓0𝑝, 𝑞 なので追加された自由度𝜁を消去する𝑓0𝑝, 𝑞 ≡ න 𝑓eq(𝑝, 𝑞, 𝜁)𝑑𝜁𝑓eq∝ exp −𝛽𝐻 −𝛽2𝜁2 であったから𝑓0∝ exp −𝛽𝐻もとの世界で見ると、カノニカル分布が実現する
9096𝑝𝑞O𝑝𝜁Oハミルトンダイナミクス由来の流れ 熱浴由来の流れሶ𝑞 = 𝜕𝑝𝐻ሶ𝑝 = −𝜕𝑞𝐻 − 𝑝𝜁ሶ𝜁 = 𝑄−1 𝑝𝜕𝑝𝐻 − 𝑘𝑇
9196• Nosé-Hoover法は、ハミルトンダイナミクスに追加自由度を追加して温度制御をする手法• 追加自由度を消去し、もとの世界の位相空間に射影すると、分布関数が厳密にカノニカル分布となる• ハミルトンダイナミクスが位相空間に作る流れが非圧縮流になるのに対して、 Nosé-Hoover法が作る流れは圧縮流になるNosé-Hoover法は、単に温度制御をする手法にとどまらず、拡張ハミルトンダイナミクスという分野を生み出した
9296𝑑𝑝𝑑𝑡= −𝜕𝐻𝜕𝑞− 𝑝𝜁𝑑𝜁𝑑𝑡=𝑘𝑄𝑇 − 𝑇0Nosé-Hoover法は「全体の温度」しかチェックしていない𝑇 =23𝑁𝑘𝑖𝑝𝑖22𝑚温度は速度の分散で求めるFlying icecube問題高速ですれ違う氷二つの系は「温度が高い」といえるか？
9396系が相分離していると、温度が一様にならないB原子の温度A原子の温度全体の温度
9496温度の時間発展 スペクトル温度制御なし温度制御あり 温度制御あり温度制御なし非自明なピーク時間𝑝𝜁O熱浴由来の流れの「回転」を反映し、系に熱浴由来の振動が加わる
9596• Nosé-Hoover法は、系がエルゴード的かつ定常状態であれば厳密にカノニカル分布を実現する→条件によっては正しく温度制御ができない• 相分離するような系では、それぞれの部分系が異なる温度に落ち着くことがある• 熱浴由来のエネルギー振動が入るため、ダイナミクスが信頼できないことがある• 別の手法であるLangevin熱浴ではこの問題はおきない(ただし温度収束が遅い)「手法」は前提条件を理解して使わないと危険
9696• 物理量の定義は難しい• 温度や圧力の定義は自明ではない• 数値計算では、支配方程式によりa prioriな変数と観測量が異なる• ハミルトンの運動方程式はNVEアンサンブルを実現する• 運動方程式を拡張することで、温度や圧力を制御できる• 手法の性質をよく知らずに使うのは危険
kaityo256
kfujita
matleenalaakso
codeforeveryone
fullfull
utokyoissr2360
tarugoconf
asahide
dwchiang
asial_corp
motherhoodinadhd
eileencodes
revolveconf
maggiecrowley
marktimemedia
cassininazir
brianwarren
chrisshort
mongodb
marcelosomers
smashingmag
destraynor
morganepeng