徐 琳，陳雨澤，劉家昊

(北京師范大學(xué) 物理學(xué)系，北京 100875)

Monte-Carlo算法是本科計(jì)算物理基礎(chǔ)課程中的內(nèi)容之一，其本質(zhì)是一種將概率問(wèn)題轉(zhuǎn)化為統(tǒng)計(jì)問(wèn)題求解的隨機(jī)模擬算法，在模擬物理問(wèn)題及數(shù)學(xué)建模等方面應(yīng)用十分廣泛.但是由于課程安排有限，這部分內(nèi)容一帶而過(guò)且未做考試要求，因此大部分同學(xué)對(duì)于Monte-Carlo算法的了解并不多，甚至并不知道如何進(jìn)行利用，可以說(shuō)比較遺憾.

本文進(jìn)行了本科課程之外的嘗試，分別利用三種不同的算法模擬了二維Ising模型，從運(yùn)算效率、模擬效果、能否消除臨界慢化現(xiàn)象等方面進(jìn)行了對(duì)比，旨在深入探討Monte-Carlo算法中Metropolis、Swendsen-Wang和Wolff算法的原理和適用范圍差異以及在二維Ising模型模擬中的優(yōu)劣問(wèn)題，希望幫助同學(xué)加強(qiáng)對(duì)Monte-Carlo算法的理解.舉一反三，為學(xué)習(xí)針對(duì)其他物理問(wèn)題的模擬和探索所使用的其他計(jì)算物理相關(guān)研究方法打下了一定的基礎(chǔ)，同時(shí)希望有利于拓展老師教學(xué)以及學(xué)生自主探究問(wèn)題的思路，傳遞自主學(xué)習(xí)和探索精神.

1 Ising模型

1.1 Ising模型簡(jiǎn)介

物質(zhì)的磁性源于原子的磁性，外層電子軌道磁矩、電子自旋磁矩、原子核磁矩在宏觀上顯示出物質(zhì)的磁性.磁鐵在加熱到某個(gè)臨界溫度以上時(shí)，磁性就會(huì)消失；降溫至臨界溫度以下，磁性又會(huì)恢復(fù).這個(gè)臨界溫度叫做居里溫度，這種從有磁性到無(wú)磁性的連續(xù)相變，叫做二級(jí)相變. 由恩斯特·伊辛(Ernst Ising，德國(guó))創(chuàng)建的磁自旋線性模型——Ising模型，起初就是為了更好地研究鐵磁性并解釋“鐵磁相變”問(wèn)題.

將鐵磁性物質(zhì)假想成規(guī)則排列的箭頭，每個(gè)箭頭擁有上、下兩個(gè)方向，對(duì)應(yīng)自旋方向. 相鄰箭頭之間存在相互作用且受外界環(huán)境的影響，能夠發(fā)生隨機(jī)轉(zhuǎn)變，即箭頭上、下方向的轉(zhuǎn)變. 溫度的大小最能影響漲落，溫度越高則箭頭發(fā)生的隨機(jī)無(wú)序轉(zhuǎn)變?cè)絼×遥?dāng)達(dá)到居里溫度以上時(shí)上、下方向的箭頭磁性相互抵消，宏觀表現(xiàn)為磁性消失；反之，溫度低時(shí)系統(tǒng)宏觀表現(xiàn)出磁性.總而言之，“熱”創(chuàng)造無(wú)序，而“磁”抵抗無(wú)序，這種競(jìng)爭(zhēng)解釋了相變過(guò)程, 如圖1所示.

圖 1 一維Ising模型

除了描述晶體磁性以外，非合金中的有序-無(wú)序相變、玻璃相物質(zhì)特性等都可以用Ising模型進(jìn)行描述. 基于Ising模型的簡(jiǎn)單和普適，將箭頭向上或向下的狀態(tài)比喻成其他意義，便可以被普遍應(yīng)用于其他領(lǐng)域模擬各類臨界現(xiàn)象.

1.2 平均場(chǎng)近似

Ising模型中每一個(gè)箭頭方向只取兩個(gè)值中的一個(gè)(即磁自旋中的1或-1，記為σ)，且僅與其最近鄰的箭頭相互作用.使用哈密頓量計(jì)算系統(tǒng)能量

(1)

其中J為相互作用能，i、j代表相鄰的兩個(gè)箭頭， 表示近鄰相互作用，μ為玻爾磁子，B表示外界磁場(chǎng)強(qiáng)度.

(2)

(3)

假設(shè)外部磁場(chǎng)強(qiáng)度為0，根據(jù)式(1)可以得到系統(tǒng)的能量E、磁化強(qiáng)度M、平均比熱C分別為:

(4)

(5)

(6)

其中N為系統(tǒng)的規(guī)格，在二維情況下可以設(shè)定為n×m的網(wǎng)格結(jié)構(gòu).

平均磁化率可寫為

(7)

2 Monte-Carlo方法模擬二維Ising模型

2.1 主要步驟

利用Monte-Carlo法模擬二維Ising模型的步驟主要分為以下幾步：1) 設(shè)置規(guī)格N、溫度; 2) 隨機(jī)產(chǎn)生初始狀態(tài)，選擇一個(gè)初始格點(diǎn)i作為翻轉(zhuǎn)算法的起點(diǎn)；3) 計(jì)算翻轉(zhuǎn)產(chǎn)生的能量差δE; 4) 計(jì)算翻轉(zhuǎn)的轉(zhuǎn)移概率；5) 產(chǎn)生(0,1)區(qū)間內(nèi)的隨機(jī)數(shù)r，若r

對(duì)于二維Ising模型來(lái)說(shuō)，與初始格點(diǎn)i近鄰的有四個(gè)格點(diǎn)1、2、3、4.因此翻轉(zhuǎn)一個(gè)自旋產(chǎn)生的能量差為

δEi=2Jσi(σ1+σ2+σ3+σ4)

(8)

顯然，能量差的可能值可以被直接應(yīng)用以節(jié)省運(yùn)算時(shí)間:

δEi=-8J,-4J,0,4J,8J

(9)

2.2 遍歷性與細(xì)致平衡條件

利用Monte-Caro法進(jìn)行計(jì)算時(shí)需要進(jìn)行概率抽樣，對(duì)于二維Ising模型來(lái)說(shuō)，就是進(jìn)行自旋的翻轉(zhuǎn)，而這樣的“翻轉(zhuǎn)”還滿足兩個(gè)條件.

首先對(duì)于任意的一個(gè)初始狀態(tài)，經(jīng)過(guò)n次翻轉(zhuǎn)算法后，最終都會(huì)達(dá)到一個(gè)穩(wěn)態(tài).無(wú)論從哪個(gè)狀態(tài)出發(fā)，達(dá)到另一種狀態(tài)的所有路徑中總有一條路徑的概率不為0.這就是遍歷性.

其次，從狀態(tài)ρ到狀態(tài)τ的概率流與從狀態(tài)τ到狀態(tài)ρ的概率流相等，也就是

Jρτ=Jτρ

(10)

將狀態(tài)ρ與狀態(tài)τ的穩(wěn)態(tài)概率分別寫作Pρ與Pτ，前者與后者之間的轉(zhuǎn)移概率分別寫作Pρτ與Pτρ.則有

Jρτ=Jτρ=PρPρτ=PτPτρ

(11)

當(dāng)整個(gè)過(guò)程達(dá)到穩(wěn)態(tài)時(shí)，狀態(tài)ρ與狀態(tài)τ之間的相互轉(zhuǎn)化也達(dá)到平衡，這就是細(xì)致平衡.

2.3 利用Metropolis算法模擬

Metropolis算法的特點(diǎn)是一次翻轉(zhuǎn)一個(gè)自旋，如圖2所示.隨機(jī)選取初始格點(diǎn)i后，假設(shè)此格點(diǎn)的自旋被翻轉(zhuǎn)，計(jì)算此轉(zhuǎn)變的δE.設(shè)定一個(gè)接受概率，決定是否要進(jìn)行此次翻轉(zhuǎn),記為Pacc.顯然，只要令Pacc≠0就滿足遍歷性.

圖 2 (a)翻轉(zhuǎn)一次后到達(dá)(b)的狀態(tài)

(12)

根據(jù)式(11)，在Boltzmann系統(tǒng)中有

(13)

結(jié)合式(12)與式(13)，得到

(14)

Nicholas Metropolis給出了符合式(14)的接受概率的解的形式，正是Metropolis接受準(zhǔn)則

(15)

利用Metropolis接受準(zhǔn)則能達(dá)到更高的效率，且更不容易以局部穩(wěn)態(tài)結(jié)局.

完成上述理論基礎(chǔ)后，按照本文2.1的步驟編寫程序進(jìn)行仿真模擬，可以得到非常接近理論情況的結(jié)果.

以能量隨溫度變化圖像為例，如圖3所示.在溫度很低時(shí)平均能量最低，此時(shí)自旋平行排列；隨著溫度升高經(jīng)過(guò)臨界點(diǎn)，由于比熱對(duì)數(shù)發(fā)散的緣故，能量隨溫度變化很快.當(dāng)溫度超過(guò)臨界溫度，平均自旋磁矩為0，但是由于近鄰自旋仍存在相互作用，平均能量并不趨于0.因此，我們的模擬結(jié)果與物理圖像符合的很好.

圖3 平均能量-溫度關(guān)系曲線

然而多次改變參數(shù)條件和循環(huán)次數(shù)后就會(huì)發(fā)現(xiàn)，這種一次只翻轉(zhuǎn)一個(gè)自旋的方式效率很低，要經(jīng)過(guò)很高的迭代次數(shù)后才能達(dá)到穩(wěn)態(tài).如圖4所示的結(jié)果是循環(huán)2 000次后得到的.

此外，溫度較高時(shí)Metropolis算法往往能夠得到符合度較好的結(jié)果，但是在溫度接近相變溫度時(shí)得到的結(jié)果很不穩(wěn)定，收斂很慢，且重復(fù)性很差.甚至達(dá)到錯(cuò)誤的穩(wěn)定狀態(tài).

圖4 (a)狀態(tài)經(jīng)過(guò)多次翻轉(zhuǎn)達(dá)到(b)狀態(tài)后穩(wěn)定，仍可能因接受一次翻轉(zhuǎn)回到不穩(wěn)定狀態(tài)(c)

如圖4所示，一次只能翻轉(zhuǎn)一個(gè)自旋的局限性，導(dǎo)致了上述不穩(wěn)定、收斂慢的問(wèn)題.因此一個(gè)一次翻轉(zhuǎn)多個(gè)自旋的方法似乎更加適用，Wolff算法就是如此.

2.4 利用Wolff算法模擬

Wolff算法的特點(diǎn)是一次翻轉(zhuǎn)多個(gè)自旋，且被翻轉(zhuǎn)的自旋的初始值是相同的.如圖5所示，將(a)中的部分自旋進(jìn)行翻轉(zhuǎn)，得到(b)狀態(tài).

圖5 一次翻轉(zhuǎn)多個(gè)自旋示意圖

圖6演示了Wolff算法的實(shí)現(xiàn)過(guò)程.假設(shè)翻轉(zhuǎn)的多個(gè)自旋稱為一個(gè)自選群. (a)為初始狀態(tài)，選定一個(gè)初始格點(diǎn)i開始進(jìn)行翻轉(zhuǎn)，以概率Pi依次判斷與i近鄰的4個(gè)格點(diǎn)是否被接受加入自旋群；一次翻轉(zhuǎn)后，原來(lái)的初始格點(diǎn)被標(biāo)記，選取新加入自旋群的格點(diǎn)作為新的初始格點(diǎn)；重復(fù)上述步驟，直到自選群中所有格點(diǎn)都被標(biāo)記；最后，以接受概率Pacc決定是否翻轉(zhuǎn)這個(gè)自旋群.

圖6 Wolff算法示意圖(√表示接受，×表示不接受)

(16)

又由式(8)得

δE=Eτ-Eρ=2J(g-h)

(17)

因此結(jié)合式(16)、(17)后得到

(18)

接受概率為1，意味著只要構(gòu)造出一個(gè)待翻轉(zhuǎn)的自旋群，那么它就一定會(huì)被翻轉(zhuǎn).

完成以上參數(shù)的確定后，同樣按照本文2.1的步驟進(jìn)行編程模擬，得到圖7的結(jié)果.值得注意的是在求磁化強(qiáng)度時(shí)，由于在低溫時(shí)磁化強(qiáng)度會(huì)在±1之間震蕩，求和結(jié)果趨近于0，因此采用絕對(duì)值更加方便.

圖7 利用Wollf算法模擬得到的(a)平均能量(b)磁化強(qiáng)度(c)比熱(d)平均磁化率隨溫度變化圖像

Wolff算法的模擬結(jié)果與Metropolis算法的結(jié)果可以相統(tǒng)一，具體可以將圖7和圖3中對(duì)于平均能量的模擬進(jìn)行對(duì)比.此外可以發(fā)現(xiàn)，在平均能量隨溫度變化最快的位置對(duì)應(yīng)著比熱的峰值，這與比熱為能量對(duì)溫度的導(dǎo)數(shù)的物理圖像也符合.

明顯的區(qū)別在于，在用Wolff算法進(jìn)行模擬的過(guò)程中，在任何給定的溫度區(qū)域都不會(huì)出現(xiàn)明顯的慢化現(xiàn)象或者不收斂的情況.尤其是在相變溫度附近，Wolff算法的運(yùn)行速度明顯快于Metropolis算法.

2.5 利用Swendsen-Wang算法模擬

Swendsen-Wang算法和Wolff算法原理類似，也是一種對(duì)“團(tuán)簇”進(jìn)行翻轉(zhuǎn)的算法.區(qū)別在于它會(huì)先按照一定規(guī)則將系統(tǒng)分成多個(gè)團(tuán)簇，每個(gè)團(tuán)簇隨機(jī)取值，然后嘗試翻轉(zhuǎn)；而Wolff算法僅僅構(gòu)造一個(gè)團(tuán)簇并翻轉(zhuǎn)它.因此我們可以簡(jiǎn)單理解為，Swendsen-Wang算法兼具了Wolff算法和Metropolis算法的特點(diǎn).那么，理論上也應(yīng)當(dāng)兼具兩者的優(yōu)勢(shì).

以平均能量隨溫度變化圖像為例進(jìn)行對(duì)比，我們發(fā)現(xiàn)在比較理想的情況下，不改變模擬條件能夠達(dá)到和Wolff算法及Metropolis算法完全相同的結(jié)果.在臨界溫度附近，Swendsen-Wang算法并沒有出現(xiàn)臨界慢化現(xiàn)象，在一定的重復(fù)次數(shù)下也不會(huì)出現(xiàn)不收斂的情況.只是與Wolff算法比較，它的運(yùn)行效率要更慢.圖8的結(jié)果是循環(huán)了3 000次得到的.

圖8 平均能量-溫度關(guān)系曲線的SW算法模擬

原因可想而知，由于翻轉(zhuǎn)“團(tuán)簇”，整個(gè)系統(tǒng)的關(guān)聯(lián)區(qū)域并不大，每一次翻轉(zhuǎn)被接受的概率很大，但不是1.因此Swendsen-Wang算法可以較好地克服臨界慢化現(xiàn)象，但是并不能像Wolff算法一樣高效.在達(dá)到相同的模擬效果的情況下，Swendsen-Wang算法要弛豫更多次的Monte-Carlo步.

3 對(duì)比總結(jié)

Monte-Carlo法是研究Ising模型這種自旋系統(tǒng)相變的一個(gè)重要方法，本文采用了經(jīng)典的“重要性抽樣”.

簡(jiǎn)單來(lái)講，Metropolis算法和Wolff算法最主要的區(qū)別在于前者每步嘗試翻轉(zhuǎn)“一個(gè)”自旋，而后者每步嘗試翻轉(zhuǎn)“一團(tuán)”自旋.概率選擇的原則是要使最終系統(tǒng)各個(gè)狀態(tài)出現(xiàn)的概率符合正則分布.對(duì)于溫度很高的系統(tǒng)兩種算法無(wú)明顯差別，但由于臨界系統(tǒng)中存在很大的關(guān)聯(lián)區(qū)域，單個(gè)自旋的翻轉(zhuǎn)很難被接受，因此此時(shí)Metropolis算法會(huì)變得異常緩慢，這也是模擬過(guò)程中在相變溫度附近結(jié)果無(wú)法收斂的原因.

Wolff算法的優(yōu)勢(shì)顯而易見.首先，Wolff算法的接受概率為1，這使得系統(tǒng)自旋轉(zhuǎn)變的效率更高、收斂更快，很大程度上消除了相變溫度附近的臨界慢化現(xiàn)象.此外，一次翻轉(zhuǎn)多個(gè)自旋使系統(tǒng)構(gòu)型變化更快，在同等計(jì)算步長(zhǎng)的情況下，Wolff算法的信噪比更高，效率更高.與Swendsen-Wang算法相比，主要差別在于Wolff算法模擬過(guò)程中會(huì)以概率1翻轉(zhuǎn)初始格點(diǎn)所在的自旋群，而對(duì)于Swendsen-Wang算法而言需要以小于1的概率翻轉(zhuǎn)每一個(gè)自旋群.因此，Wolff算法的效率要比Swendsen-Wang算法更高.

4 結(jié)論

通過(guò)對(duì)二維Ising模型的模擬，我們發(fā)現(xiàn)與Metropolis算法和Swendsen-Wang算法相比，Wolff算法不僅模擬效率高、效果好，還能夠更好地克服臨界慢化現(xiàn)象.本文通過(guò)分析對(duì)比的方式，能夠更深刻地理解Wolff算法被稱是最好的聚類翻轉(zhuǎn)Monte-Carlo算法的原因，明白Metropolis算法和Swendsen-Wang算法的劣勢(shì)和不同的適用范圍.本文是在本科計(jì)算物理課程之外自主探索的結(jié)果，希望能對(duì)學(xué)生自主學(xué)習(xí)、拓展學(xué)習(xí)計(jì)算物理的思路起到啟發(fā)作用.