Monte Carlo simulations 京都大学工学研究科材料工学科西谷滋人平成 14 年 6 月 25 日 1 Monte Carlo の原理 [Wm .G. Hoover, 計算統計力学, (森北出版，1999), pp.68–70.] 平衡統計力学の目標は，どんな観測できる量も，対応する微視的状態関数の平均として表すことである．正準集団では，状態の重みは exp(−Ei /kT ) とわかっているので，この重みに比例する頻度で微視的状態を発生させるのはやさしい．そのような状態を生成し，重みつき平均を計算すれば，その平均は次第に正準集団平均 Obs(q, p) = Obsi exp(−Ei /kT ) exp(−Ei /kT ) (1) に収束するだろう．直接この平均をとる方法は，自由度が数個以上ある大抵の興味深い系では実用的でない．なぜなら，相空間のほとんどの部分が無視できる確率しか持たないからである．修正モンテカルロ標本化法は Metropolis, Resenbluth 夫妻 Teller 夫妻によって剛体円盤の状態方程式に応用された． [N.A.Metropolis, A.W. and M.N. Rosenbluth, A.H. andE. Tellr, J. Chem. Phys. 21, (1953),1087.] この方法は，配置空間における比較的小さなエネルギー変化 ∆E = kT と exp(−∆E/kT ) に比例した状態遷移確率をもとにした，単純で使える方法， Obs = Obsi 1 (2) であり，陽な重みがいらない．モンテカルロ配置の列は，非対称な動きから生成される．ポテンシャルエネルギーを減らすような動きはすべて受け入れられる．ポテンシャルエネルギーを増やすような動きは，確率 exp(−∆E/kT ) で受け入れられる．この方法は確率密度が最終的に exp(−E/kT ) に収束することを保証する (以下で詳述)． 2 巡回セールスマン問題ここでは，巡回セールスマン問題を例にモンテカルロシミュレーションのアルゴリズムを見ておきましょう．巡回セールスマン問題とは，ある街から出発していくつかの街を次々とめぐって元の街に戻ってくる最短の経路を求める問題です．訪れる街の数が少ないときにはすべての経路を数え上げればいいのですが，数が増えるとその計算時間は指数関数的に増えてしまうと予想される問題です．このような問題はセールスマンだけでなく，コンピュータの CPU の配置や，都市ガスの配管設計などでも出会う問題です． 1 2 巡回セールスマン問題対象となる経路の長さの合計を E(a) = r[ai ] − r[ai+1 ] (3) i=1..N としておきます．ここで a はそれぞれの街の順番あるいは配置を示しています．r はそれぞれ街の座標で r で距離を求めます．するとこの関数は模式的に図 1 のようになると考えることができます．すべての配置について E(a) を求めれば，最小値が求まります．あるいは初期の配置を a = [1, 2, 3, . . . , N, 1] (4) として，一定の手順で変更 δa を加え，E(a) が下がった場合にその配置を採用するという方法をとることも出来ます．しかし，この方法は極小値に落ちてしまうことが分かるでしょう．最小値を探すには時として坂を駆け登る必要があるのです．このような問題の一つの解法としてアニーリング法 (simulated annealing) があります．これは格子欠陥を多く含んだ金属を高温へ上げて欠陥を掃き出し柔らかくする熱処理 (焼きなまし，annealing) からの類推で名付けられた手法です． E(a) a 図 1: 配置 a とその経路の総和 E(a) を示す模式図アルゴリズムは以下のとおりです． 1. 配置 a を仮定し E(a) を求める． 2. a からすこし違った配置 a + δa を作る． 3. ∆E = E(a + δa) − E(a) を求める． 4. ∆E < 0 なら新たな配置を採用する． 5. ∆E > 0 なら新たな配置を exp(−∆E/T ) の確率で受け入れる． 6. 手順 2 以下を適当な回数繰り返す．ここで T は温度から類推される制御パラメータで，十分大きい場合はすべての状態が採用されます．一方，T を下げるにつれて採用される試行が少なくなり，徐々に状態が凍結されていきます． a + δa の生成方法と T の下げ方をうまく取れば最小値に近い状態が確率的に高く出現し，最小値かそれに近い状態へ落ち着くというアルゴリズムです．N = 16 の場合の (a) 初期配置，(b) 単純に 2 3 MONTE CARLO のアルゴリズム (a) (b) E(a)=9.601317879 a=[1,2,3,4,5,6,7,8,9,10, 11,12,13,14,15,16,1] (c) E(a)=4.379272386 a=[1,5,15,9,11,7,13,2,1 4,16,3,10,8,4,12,6,1] E(a)=3.572634278 a=[1,13,7,11,9,15,10,5, 3,14,16,8,4,12,6,2,1] 図 2: ランダムに生成した 16 個の街に適用した巡回セールスマン問題の (a) 初期配置，(b) 単純に E(a) が下がった場合だけを採用した結果，および (c) simulated annealing の計算結果． E(a) が下がった場合だけを採用した結果，および (c) simulated annealing の計算結果を図 2 に示しました．以下は街の位置をランダムに生成し，表示する Maple script です． restart; > with(plots): > N_city:=16; > Path:=[seq(i,i=1..N_city),1]; > Position:=seq([evalf(rand()/10^12),evalf(rand()/10^12)],i=1..N_city): > Real_Path:=[seq(Position[Path[i]],i=1..N_city),Position[1]]: > PLOT(CURVES(Real_Path)); 街の間の距離をあらかじめ計算して，2 次元の配列に入れておき，配置 a から任意に 2 都市を取りだして入れ替える試行を δa として上述のアルゴリズムにしたがって最短経路を探すプログラムを作成してください． 3 Monte Carlo のアルゴリズム (N,V,T) が規定された正準集合 (canonical ensemble) を対象にした正準モンテカルロ法のアルゴリズムにしたがって作成したコードを以下に示す．アルゴリズムとしてはモンテカルロ法と焼きなまし法は全く同じです．しかし，熱平衡物性を求めるモンテカルロ法と，最安定状態を求める焼きなまし法とは目的が違い，温度制御の仕方が違うという点は認識しておいてください．最大のステップ幅 d をどのように選択するかが常に問題となる．d が大きすぎると，試行ステッ 3 4 マルコフ過程と遷移確率 List 3.1 NVT 一定 MC の Code のコア部分 int reject=0; double Sum_U=0,TotalU; double u_i=Energy.Total(); TotalU=u_i; // STEP_2 u_i for(int iter=1;iter<=iter_max;iter++){ x=rand_generate(&ix); int select_atom=x*nbase; // STEP_3 select atom try_move(select_atom); // STEP_4 try_move double u_j=Energy.Total(); // STEP_5 u_j double dU=u_j-u_i; if( dU <= 0 ){ // STEP_6 TotalU = u_i = u_j; } else { x=rand_generate(&ix); // STEP_7 if( exp(-dU/kT) > x){ // STEP_8 TotalU = u_i = u_j; } else { // STEP_9 reset_move(select_atom); reject++; } } Sum_U += TotalU; } プ数のわずかしか受け入れられないので確率密度関数 p(x) のサンプリングが不十分ということになる．一方，d が小さすぎると，p(x) のサンプリングが少なすぎるということになる． d のおおよその値は試行ステップ数の 3 分の 1 から 2 分の 1 といったところが妥当とされている．また，初期値は p(x) ができるだけ速やかに非対称分布に近づくように設定する．これは p(x) が最大値をとるような配置をとることで達成できる． 4 マルコフ過程と遷移確率 [Ueda p.90] サンプリングの手続きから明らかなように，ある時刻に状態 m を系がとるか否かは，その 1 ステップ前の系の状態 l にだけ依存し，それ以前に系がどんな状態を経てきたかにはよらない．このように系の状態が確率の法則に従い，次々遷移していく確率過程において，次のステップの状態が現在の状態だけで確率的に決まるとき，この確率過程をマルコフ過程 (Markoff process) と呼び，生成された状態の連鎖はマルコフの連鎖を構成するという．前項で生成された状態の出現確率が，ステップ数 n → ∞ のとき，カノニカル分布 exp (−U/kB T ) に比例することは，以下のようにして示される．いま状態 l にある系が k ステップで状態 m に遷移する確率を (k) (1) Plm , Plm ≡ Plm (5) とおくと，マルコフ過程の場合，この遷移確率に対してチャップマン- コルモゴロフの式 (Chapman- Kolmogorov equation) (k) (k−1) Plm = Plj Pjm (6) j が成立する．ここで和は k − 1 ステップで系のとりうるすべての状態についてとる．マルコフ過程の遷移確率について次の定理がある． 4 4 マルコフ過程と遷移確率もしすべての状態が同一のエルゴード状態に属するならば， (k) lim Plm = wm , m = 1, 2, · · · , s k→∞ (7) なる極限分布 wm がすべての状態 m(s は状態の総数) に対して存在し，wm は初期状態によらず， wm > 0, m = 1, 2, · · · , s (8) である．また極限分布 wm は，(6) 式で k → ∞ とした方程式 wm = wj Pjm (9) j を満足する． k > 0 となここでエルゴード状態とは，エルゴード状態に属する任意の状態 l と m に対して，Plm る有限の k が存在する．つまり状態 l から m へ有限回の遷移で到達できるような状態の集合をさす．この定理の証明は [S. Karlin: A First Course in Stochastic Processes(Academic Press, 1966), 訳：佐藤健一，佐藤由身子：確率過程講義 (産業図書，1974)] を参照されたい．そこでカノニカル分布を実現するには，wm = exp (−Um /kB T ) とおき，(9) 式を満足する遷移確率 Pjm を求めればよい．(9) 式を exp (−Um /kB T ) = exp (−Uj /kB T ) Pjm j = exp (−Um /kB T ) Pmm + exp (−Uj /kB T ) Pjm j=m と書き換えて，この左辺に 1= Pmj = Pmm + j Pmj j=m を掛けると exp (−Um /kB T ) Pmj = exp (−Uj /kB T ) Pjm j=m j=m を得る．遷移確率 Pmj はこの式を満足するように選べばよい．そのためには exp (−Um /kB T ) Pmj = exp (−Uj /kB T ) Pjm (10) が成り立てばよい．これは十分条件である．j の代わりに l を用いて Pml : Plm = 1 : exp (− (Um − Ul ) /kB T ) と書き換えると，Um > Ul ならば，状態 m から l へは確率 1 で系を遷移させ， l から m へは確率 exp (− (Um − Ul ) /kB T ) で遷移させれば良いことが分かる．さきに述べたサンプリングのアルゴリズムは，この結果に基づいている．関係 (10) 式は，詳細釣り合いの原理 (principle of detailed balance) と呼ばれる関係であって，粒子のミクロな動きによる状態 m から l への遷移の起こる頻度が，逆の l から m への頻度と等しいことを示す．カノニカル分布で表された熱平衡状態はこの原理によって実現することが保証されるのである． 5 参考文献参考文献参考文献 [Kubo] 久保亮五編，大学演習熱学・統計力学，(1961 裳華房)． [Kamiyama and Sato] 神山新一，佐藤明，モンテカルロ・シミュレーション，(1997 朝倉書店)． [Ueda] 上田顕，コンピュータシミュレーション (-マクロな系の中の原子運動-)，(1990 朝倉書店)． [Gould] Harvey Gould and Jan Tobochnik，計算物理学入門，鈴木増雄監訳，溜渕継博監翻訳，石川正勝，宮島佐介訳 (ピアソン・エデュケーション，2000) ． 6