易辛模型 正則系統 漢米爾頓( 即能量)

實驗五

 目的

 原理

1.模型與正則系統

2.Metropolis Algorithm

3.臨界現象

4.熱力學極限

 實驗內容

目的

本實驗在判紹相轉變與臨界的最重要模型-易辛模型。我們的重點在學習

如何以Monte Carlo 求和的方法來計算熱力學的平均量，並觀察相轉變的性

質．

原理

1.模型與正則系統

易辛模型在統計物理中一個非常重要的模型，無論是相變理論的建立或

是對實用系統的模擬，它都佔有非常重要的地位。首先，讓我們先來看看它

的定義。

假想我們有一個具有N 個晶格點的晶格，在每一格點上都有一個自旋變

數，i=1,2,3,…,N， 只能有兩個數值。此晶格系統的狀態即由此N

個定義系統的微觀狀態，而一有N 個自旋變數的易辛模型其微觀

狀態( microstate ) 的數目顯然有個。我們有興趣的物理或熱力學量

通常都是微觀伏態的函數，A=A(S) ，例如最重要的就是能量E 及序參量

( order parameter ) M。如果我們知道任一微觀狀態S 出現的機會，則經由權重

求和即可得到這些物理量的平均值。

在統計理中，一個系統若與一溫度為T 的熱庫保持熱接觸，則此系統達

熱平時，其微觀態出現的機會就由Boltzmann 機率數來決定，因此物理量A

在熱平時的大小即為加下的平均

(5-1)

其中即為Boltzmann 函數

， (5-2)

而Z 則由機率守恆， 來決定

(5-3)

也就是說溫度T 與系統在所有狀態s 的能量結構， ，決定了該系

統的熱力學性質。這也就是為什麼我們要以模型的能量結構來作為模型

的定義。對易辛模型而言，其漢米爾頓( 即能量)為

(5-4)

其中J 代表相鄰自旋間的交互作用強度，B 代表外表磁場，因此第二

項即為Boltzmannn 能量。第一項的n. n. 表示最近鄰(nearest neighbor)，易

辛模型的序參量即為。若J>0 系統的最低能量狀態是所有的

( 或-1)，也就是(或-N)，此乃所謂鐵磁性，反之若J<0，則系統

的最低能量狀態為相鄰的自旋反向，此即反鐵磁性。

2.Metropolis Algorithm

簡單的說，統計力學要計算的就是權重求如(weighted sum)

(5-5)

要對(5-5)精確求如是非常困難的，僅有少數的模型我們知道它的精確

解。事實上，如何模型的對稱性上著手決定其是否可解(Yang-Baxter

equation) 以及如何經由特別的Bethe ansatz 求出其解已成為統計物理

中非常數學化的領域，我們不可能在此討論。我們的目的是要避免這些

繁複的數學過程，而由計算機來幫助我們求和，並從這些模擬結果來理

解相變理論。

同學們應該留意到(5-5)式的求和項數是所有的微觀態，對易辛模型

而言，則有

個狀態。統計力學是處理熱力學極限下的平

衡性質，因此，

是一個可怕的天文數字，我們不指望當今的計算機能在我們有生之年完

成這個求和計算。另一方面，我們也留意到(5-5)式正是我們討論過的高

維度的積分(或求和) ，因此Monte Carlo 法會是一個無可避免的選

擇。也正如我們討論過的，我們可以在樣本空間中均勻的抽取m 個樣

本，並且對抽出的樣本計算之和即：

(5-6)

或者，我們可以在樣本空間抽出一筆狀態使其滿足個特定習分布

，則(5-5)式可近似為

(5-7)

選取若即為Boltzmann 分布函數(5-2) 式

(5-8)

(5-9)式即表示的算數平均。因此，我們最有興趣的是如何產生一組狀態

集合使其滿足(5-8) 之分布。

通常要產生一系列滿足Boltzmann 分布的狀態是採用所謂的Markov

process 來完成的。關於下面我們要介紹的Metropolis 與Markov process 之

間的關係，我們在本實驗的附錄中詳細說明，為不影響同學們閱讀的連貫性，

我們將直接介紹Metropolis 所提出的一種產生Boltmann 分布的程序。

Metropolis algorithm ：

(1) 選擇一初始的微觀狀態。

(2) 自初始狀態得一新的嘗試狀態;例如我們可將初始狀態的某一個自旋反向

而得一新的嘗試狀態。

(3) 計覦舊狀態( old ) 與嘗試狀態( trial ) 之能量差。

(4) 若就接受嘗試狀態為新的狀態，略過下面的步驟，直接執行步驟

8。

(5) 若則計算。

(6) 產生一個在(0，1)之間的均勻亂。

(7) 若則接受嘗試狀態為新的狀態;否則舊的狀態保留為下一步的新狀

態。

(8) 對新的狀態( 無論是由(4)或(7) 所決定) 求取物理量。

(9) 重覆(2) –(8) ，直到產生足夠的狀態。

(10)重覆(1) –(9) 數次，並計算物量的平均值的漲落，以決定誤差是否在可

允許範圍之內。

若以易辛模型為例，我們可將上述步驟逐一說明首考慮一維N=4 的易辛

模型，假設初始狀態為

接著我們可以利用一亂數選擇將4 個自旋中之一自旋反向，則可得一新的狀

態，例如

或

前者將反向後者則為。以上即步驟( 1 ) 、( 2 ) ，其次計算

能差，對量而言， 而對

而言，則，假設J>0。根據步驟( 4 ) ，若嘗試狀態

為

則直接計算此狀態的理量如，

並接受此嘗試狀態為新的狀態。同理，若嘗試狀態為則

(步驟( 5 ) ) ，我們產生一個( 0，1 ) 的均勻亂數，若r> 則

舊狀態仍保留為新的狀態( 步驟( 7 ) ) ，故

， 。我們可能對上述步驟執行1000 次平均重覆10 次( 共做1000

次的抽樣)，對此10 個E 與M 做統計分析而得其誤差。

對於易辛模型，我們有興趣的物除了E 與M 之外，還有比熱C，

(5-10)

磁化率，

(5-11)

自旋關聯函數，

3.臨界現象

當溫度為零時，易辛模型將停留在能量最低的狀態，即所有自旋同向，如

則(假設循環邊界條件) ， 。當溫度升高時，

有些自旋因熱擾動而有一些機會翻成向下，如

則其M 值就減小為。而當溫度足夠高的時

候，則每一自旋向上與向下機會楮當，即有N/2 為，N/2 為，則

M=0。因此的典型行為如下

當時，M 的行為呈現power-law ，

(5-13)

而稱為臨界溫度( critical temperature ) ，即當M 變為0 的溫度。

同樣的，比熱與磁化率也有相對的power-law

(5-14)

(5-15)

例如比熱的圖致如下

對於不同的模型的值可能不同，但

(5-16)

卻對任何模型似手都是對的， (5-16) 即著名的Rushbrook scaling law

4.熱力學極限

同學需認識相轉變是一個熱力學極限下的現象， 才可能發生。對

有限的系統(尺度為L )，則會觀察到很大的起伏現象。例如

與

是完全等價的兩個狀態，它們有相同的能量卻有完全不同的序參量

M。若L 太小，則Monte Carlo 模擬的過程中可能會在兩個狀態附近

徘徊相同的時間，則即使溫度很低，卻可能得限的系統中，

也僅以一極大值來現。因此在使用Monte Carlo 模擬來在研究統計模

型時，我們必須對不同的晶格L 做模擬，藉由的過程來外

差其熱力學極限下的行為。基於這點，我們必須再強調Monte Carlo

的重要性，因為一旦尺度增加，(5-5)式的狀態數亦隨之指數增加，而

Monte Carlo 方法的精確度則仍以的方式收斂，這使得我們

仍有希望對大尺度的晶格進行計算。

最後在我們開始進行實驗內之前，讓我們陳列一個二維易辛模型

的Monte Carlo 模擬程式，同學們必須仔細閱讀每一令對前一節所列

舉例子的密切關係。

program Ising

implicit real*8(a-h,o-z)

real ran2,rr,rq,rtr,mag

parameter(nl=16,ns=nl**2,nnb=4,nbo=ns*nnb/2)

parameter nbin=10, nmcs=1000)

dimension ic(ns),ibond(nbo,2)

dimension isl(ns),isr(ns),isu(ns),isd(ns)

dimension iscon(ns,nnb)

dimension wt(-2*nnb:2*nnb),x(nbin)

dimension bneng(nbin),bnmag(nbin)

c

c

 c --------Lattice Connection--------

c

do i=1,ns

isl(i)=i-1

enddo

do j=1,nl

jj=nl*(j-1)+1

isl(jj)=nl*j

enddo

do i=1,ns

ll=isl(i)

isr(ll)=i

enddo

do j=1,ns

isu(j)=j+nl

if(j.gt.ns-nl) isu(j)=isu(j)-ns

enddo

c

do i=1,ns

ll=isu(i)

isd(ll)=i

enddo

c

do i=1,ns

iscon(i,1)=isl(i)

iscon(i,2)=isr(i)

iscon(i,3)=isu(i)

iscon(i,4)=isd(i)

enddo

c

c

c-------Basics Parameters for MCS--------------

c

nwormup=100

rns=1 .d0/ns

rns2=1 .d0/ns/ns

rmcs=1 .d0/nmcs

c----------initial config. ic ?----------

c

do i=1,ns

ic(i)=1 !ordered state

enddo

c

energy=0.

mag=0.

do i=1.ns

mag=mag+ic(i)*ic(i)

do j=1 ,nnb

ij=iscon(i,j)

energy=energy+ic(i)"ic( ij)

enddo

enddo

c

c

open(7,nle='Il 6.datl)

write(7,") 'Q=',nqs,' SQ lattice, L=',nl,' MCS=',nmcs

c

c

do 6000 nt=0,30

c

c

c

temp=0. 1 +0. 10d0*nt

g=1. d0/temp !Ferromagnet:IC

c g=-1 .d0/temp !AF

c

do i=-2*nnb,2*nnb

wt(i)=dexp(1 .d0*g*i)

enddo

c

c

c##############################################

c

c

do 710 kmcs=1 ,nbin

c

c

teng=0.d0

tmag=0.d0

c

c

c============================================

c

c

do 720 jmcs=1 ,nmcs

c

c

do 601 is=1,nmcs

c

c

nsi=0

do i=1,nnb

js=iscon(is,i)

nsi=nsi+ic(js)

enddo

nde=-2*ic(is)*nsi

trpr=wt(nde)

c

c

if(trpr.it. 1 .d0) then

rtr=ran2(iseed)

if(rtr.gt.trpr) cycle

endif

c

c-----ACCEPT the new state----

c

ic(is)=-ic(is)

energyzenergy+1 ."nde

mag=mag+2."ic(is)

c

601 continue

c-----------------------------------

c

c stop

teng=teng+energy

tmag=tmag+mag*mag

c

720 continue

c=============================================

c

aeng=teng*rmcs*rns/nnb

amag=tmag*rmcs*rns*rns

bneng(kmcs)=aeng

bnmag(kmcs)=amag

710 continue

c#########################################

c

call binavg(bneng,nbin,avgeng,sde)

call binavg(bnmag,nbin,avgmag,sdm)

write(7,752) temp,avgeng,sde,avgnlag,sdm

write(*,752) temp,avgeng,sde,avgmag,sdm

752 format( f.3,2x,3(f9.5,' ',f9.5,2x))

6000 continue

7788 continue

end

c

c**********BIN AVERAGE******************

c

subroutine binavg(x, nbin,xav, dx)

implicit real*8(a-h,o-z)

dimension x(1 0)

sum=0.d0

do i=1,nbin

sum=sum+x(i) '

enddo

xav=sum/nbin

sum=0.d0

do i=1,nbin

sum=sum+(x( i)-xav)""2

enddo

dx=dsqn(sumfnbin')

end

c

c

c

c****************************************

c

c

c

Function ran2(idum)

integer

idum,IMI ,IM2,IMMI ,IAI ,k2,IQI ,IQ2,IRI ,!R2,NTAB,NDIV

real ran2,AM,EPS, RNMX

parameter (IM1 =2147483563, IM2=2147483399,AM=l .0/IM1 ,

& IMMI =IMI -1 , IAI =40014, IA2=40692, IQI =53668, IQ2=52774,

& IRI =12211 , IR2=3791 ,NT/kB=32,NDIV=1 +IMMI/NTAB,

& EPS=1 .2e-7,RNMX=1 .0-EPS)

C Long period (>2el 8) random number generator. Initialize

C by setting idum to a negative integer; thereafter, do not

C alter idum between successive deviates in a sequence.

C Returns a random deviate contained in (0.0,1 .0).

C Source: Numerical Recipes in Fortran, 2nd Edition.

integer idum2,j,k,iv(NTAB),iy

save iv,iy,idum2

data idum2/23456789/. ivfNTAB*0/, iy/0/

If (idum .le. 0) then

idum = max(-idum,1 )

idum2 = idum

Do j = NTAB+8, 1, -1 .

k = idumfIQI

idum = IA1*(idum-k*IQ1 )-k*IRI

if (idum.0It.0) idum = idum+IM1

if tj.le.NTAB) iv(j) = idum

End Do

iy = iv(1)

End If

k = idum/Q1

idum = IA1*(idum-k*IQ1 )-k*IR1

if (idum.It.0) idum=idum+IM1

k = idum2/IQ2

idum2 = IA2*(idum2-k*IQ2)-k*IR2

if (idum2.It.0) idum2=idum2+IM2

j = 1+iy/NDIV

iy = iv(j)-idum2

iv(i)=idum

if (iy.lt.1 ) iy=iy+IMM1

ran2 = min(AM*iy,RNMX)

return

End

實驗內容

1．請寫一個程式可以精確計算一個易辛模型的E，M、C、

等物理量在不同溫度的值。也就是說，此程式必須能產生所

有的狀態， ，並且正確地將以及

算出。

2．利用program Ising 計算上一題的晶格的物理量，並與1

的結果比較得其true error , 再與以平均值的統計誤差比較，以

了解我們如何掌握Monte Carlo 模擬所得數據的誤差。

3．program Ising 中尚未考慮加入磁場，請您將磁場造成的Zeeman

Energy 加入程式的計算中。利用修改過的程式，計算一晶

格的序參量M 在不同溫度下，隨著磁場

的變化曲線。

4．追蹤、與晶格之序參量M 模擬過程中，在

時，是否能觀察到M 有正到負的變遷。對不同晶格尺度，

此變遷所需的時間有何不同。

5．對二維方形晶格模擬以得到其及。

6．重覆5 但為三維簡單立方晶格。

7.重覆5 但為二維三角形晶格。