Monte Carlo Simulation of 2D Ising Moldel（code）

來自專欄 Monte Carlo Method3 人贊了文章

2D伊辛模型的理論講解網上很多了，除了wiki, 我再推薦一個講解的很不錯的網站

ISING模型 - 集智百科?

wiki.swarma.net

我第一次實現Ising model用的是逐點更新的方法（此後還會用團簇更新重新算一遍），即隨機選取一個點來決定是否翻轉這個點的磁矩，下面是我自己模擬這個模型得到的一些結果：

M^2-T Graph For L=4, 8, 16, 32

In(M^2)-In(L) Graph For T=2, 2.269, 2.6

註：

T=2.269時，縱軸的四個數據縱坐標分別為： -0.267（+-0.0049），-0.434（+-0.0078），-0.602（+-0.0039），-0.804（+-0.034），之間的差值分別為 -0.167（+-0.0127），-0.168（+-0.0117），-0.202（+- 0.0379）[-0.164~-0.240,]，說明在誤差允許範圍內，In(M^2) 和 In(L) 呈線性相關關係（括弧內為errorbar）

把T=2.269的四個數據做線性擬合，如圖所示，斜率為0.257，和理論值0.25（2倍的臨界指數）非常接近。關於ising model 臨界指數的介紹見下網頁：

https://en.wikipedia.org/wiki/Ising_critical_exponents?

en.wikipedia.org

代碼（python 3）：

主函數塊

#2D Ising Model Monte Carlo Simulationfrom numpy import randomimport numpy as npfrom numba import jit#some parametersiterations=80000000dimension = 32# Energy of the latticek=dimensionextract=np.arange(1,k+1)extract[k-1]=0def total_energy(matrix): S=[matrix[m,n] for m in range(k) for n in range(k)] S1=[matrix[i,j] for i in range(k) for j in extract] S2=[matrix[i, j] for i in extract for j in range(k)] E=[S[n]*(S1[n]+S2[n]) for n in range(len(S))] return -np.sum(E)@jit(int32(int32[:,:],int32,int32))def local_energy(matrix, m, n): list = [m + 1, n + 1, m - 1, n - 1] if list[0] > k - 1: list[0] = list[0] - k if list[1] > k - 1: list[1] = list[1] - k if list[2] < 0: list[2] = list[2] + k if list[3] < 0: list[3] = list[3] + k return -matrix[m, n] * (matrix[list[0], n] + matrix[m, list[1]] + matrix[list[2], n] + matrix[m, list[3]])#Markov Chaindef Markov_chain(iteration,Temperature): # create 2d lattice x_array = random.randint(0, 2, size=(dimension, dimension)) x = x_array * 2 - 1 #some useful matrix,parameters Energy_chain = [] M_chain=[] C_chain=[] X_chain=[] MN_value=random.randint(0, dimension, size=(2, iterations)) E_total=total_energy(x) for i in range(iteration): m=MN_value[0,i] n=MN_value[1,i] E_0= local_energy(x,m,n) E_change=-2*E_0 x[m, n] = -x[m, n] E_total = E_total+E_change if E_change >0: probability = np.exp(-E_change /Temperature) if random.rand()>probability: x[m, n] = -x[m, n] E_total = E_total - E_change M_chain.append((np.sum(x)*(np.sum(x)))/(dimension*dimension)) Energy_chain.append(E_total) return [Energy_chain,M_chain]#Sampling and find the average energydef E_M_C_X_value(T): cut=round(iterations/10) list=Markov_chain(iterations,T) E=np.mean(list[0][cut:])/(dimension*dimension) M=np.mean(list[1][cut:])/(dimension*dimension) C_cut=[np.var(list[0][i:i+cut])/((dimension*dimension)*T*T) for i in range(0,iterations,cut) ] X_cut=[np.var(list[1][i:i+cut])/((dimension*dimension)*T) for i in range(0,iterations,cut)] C=np.mean(C_cut[1:]) X=np.mean(X_cut[1:]) E_error=np.sqrt(np.var(list[0][cut:])/(len(list[0][cut:]))) M_error=np.sqrt(np.var(list[1][cut:])/(len(list[1][cut:]))) C_error=np.sqrt(np.var(C_cut[1:])/(len(C_cut[1:]))) X_error=np.sqrt(np.var(X_cut[1:])/(len(X_cut[1:]))) return [E,M,C,X,E_error,M_error,C_error,X_error]

2，調用函數畫圖

（1）E, M, C, X

#2D Ising Model Monte Carlo Simulationimport numpy as npimport matplotlib.pyplot as pltimport monte_carlo_functionimport timestart=time.process_time()#derive the relation E(T)T=[1.0,2.2, 2.269, 2.35, 2.4, 2.45]E=[]M=[]C=[]X=[]E_error=[]M_error=[]C_error=[]X_error=[]print(Temperature list is,T)for i in range(len(T)): print(Process:,i+1,/,len(T)) list=monte_carlo_function.E_M_C_X_value(T[i]) E.append(list[0]) M.append(list[1]) C.append(list[2]) X.append(list[3]) E_error.append(list[4]) M_error.append(list[5]) C_error.append(list[6]) X_error.append(list[7])print(E is,E)print(M is,M)print(C is,C)print(X is,X)print(E_error is,E_error)print(M_error is,M_error)print(C_error is,C_error)print(X_error is,X_error)# plot the E(T) graphplt.figure()plt.errorbar(T,E,yerr=E_error,fmt=o,ecolor=r,color=b,elinewidth_=2,capsize=4)plt.xlabel(Temperature / K)plt.ylabel(Energy / J)plt.title(2D Ising Model-E(T) Graph)plt.grid(True)plt.show()#plot the M(T) graphplt.figure()plt.errorbar(T,M,yerr=M_error,fmt=o,ecolor=r,color=b,elinewidth_=2,capsize=4)plt.xlabel(Temperature / K)plt.ylabel( Magnetization )plt.title(2D Ising Model-M(T) Graph)plt.grid(True)plt.show()#plot the C(T) graphplt.figure()plt.errorbar(T,C,yerr=C_error,fmt=o,ecolor=r,color=b,elinewidth_=2,capsize=4)plt.xlabel(Temperature / K)plt.ylabel(Specific Heat Capacity)plt.title(2D Ising Model-C(T) Graph)plt.grid(True)plt.show()#plot the X(T) graphplt.figure()plt.errorbar(T,X,yerr=X_error,fmt=o,ecolor=r,color=b,elinewidth_=2,capsize=4)plt.xlabel(Temperature / K)plt.ylabel(Permeability)plt.title(2D Ising Model-X(T) Graph)plt.grid(True)plt.show()end=time.process_time()print(end-start)

（2）計算臨界指數

import numpy as npimport matplotlib.pyplot as plt#data storageT=[1.0, 1.4, 1.7, 2.0, 2.05, 2.1, 2.15, 2.2, 2.269, 2.35, 2.4, 2.45, 2.6, 2.8, 3.2, 3.6, 4]M1=[0.9984284895833333, 0.9839847222222222, 0.9470226388888889, 0.8675694791666667, 0.8513590451388889, 0.8301091840277778, 0.8113871354166666, 0.789850625, 0.7658821527777778, 0.7218853993055555, 0.7043654340277777, 0.6803491666666667, 0.6109082465277778, 0.5249506076388889, 0.39011116319444444, 0.3008007465277778, 0.24716826388888888]M2=[0.9986219631076388, 0.9828937730034722, 0.9431147625868056, 0.8446861345486111, 0.8168956751302083, 0.7800217745225695, 0.7483616121961806, 0.7093512170138889,0.6479291383463541, 0.5689366056857639, 0.5210536182725695, 0.477038625, 0.35005108289930553, 0.24247147352430556, 0.12705987565104165, 0.08772311284722223, 0.06686821245659723]M3=[0.9986008398301867, 0.9829366279771593, 0.9416635219658746, 0.8333712326388889, 0.800862243197971, 0.757239956753201, 0.7157412747395834, 0.649229851453993, 0.5478898898767542, 0.4113860016140408, 0.32151242104085287, 0.2407911782023112, 0.12640475245496963, 0.07073438117133246, 0.033179971503363714, 0.02246400141737196, 0.017040189276801215]M4=[0.9985473033322229, 0.9828611836227841, 0.9409495011992984, 0.8311556779265934, 0.793744361501058, 0.7456932663826412, 0.7047716192192501, 0.6135119233641094, 0.447242976656278, 0.22510332686551413, 0.15203926261234282, 0.08819298803551992, 0.03659119166946411, 0.016650539838578966, 0.008184549450556437, 0.005154495224210951, 0.00424778634707133]M1_error=[0.000328525148509109, 0.001160639141834489, 0.002268002174078354, 0.003788017679697471, 0.003982838449476648, 0.0042855034286385, 0.0044602113171121754, 0.004707674996794537, 0.004886777734098844, 0.005227656562248408, 0.005339283044349196, 0.005460619496380141, 0.005683464311277928, 0.0058040574223762334, 0.005432797630985786, 0.0048638796251160715, 0.004362337546551876]M2_error=[0.0002883358847623713, 0.0011469476962931806, 0.0023421270468071023, 0.004751437769646956, 0.0053314077705554136, 0.006057186964279236, 0.006648656651882294, 0.0055427413839558804, 0.007788211921715816, 0.008390379891119435, 0.008547705642675598, 0.00857381960654556, 0.008145504045975525, 0.006910511470924287, 0.004509820080195631, 0.0033295942767259975, 0.0026302868516527656]M3_error=[0.0001469163484009741, 0.0016077506262109095, 0.003437961623707859, 0.007834226155112338, 0.009887610137224065, 0.011373103390126461, 0.012691878703365909, 0.010792466191136125, 0.0039002144811601234, 0.004194260803423116, 0.004268640721794869, 0.004295472609456829, 0.01226026031266311, 0.00742401475509187, 0.003853021456862876, 0.0026369100639426174, 0.002014098505653076]M4_error=[0.00042507497304176804, 0.002287257777618821, 0.004964842887854574, 0.011590453666701823, 0.016316665617351074, 0.02057847723724343, 0.02124988973538772, 0.023202377308610075, 0.034277706844300794, 0.030884789006276992, 0.024295588262924557, 0.016587076950623718, 0.011541722908068179, 0.005363282107165037, 0.002809385670743479, 0.0017687361482280927, 0.0014402802280875499]InL=[np.log(4),np.log(8),np.log(16),np.log(32)]#T=2.0InM1=[np.log(M1[3]),np.log(M2[3]),np.log(M3[3]),np.log(M4[3])]M1_error_cut=[M1_error[3],M2_error[3],M3_error[3],M4_error[3]]#T=2.269InM2=[np.log(M1[8]),np.log(M2[8]),np.log(M3[8]),np.log(M4[8])]M2_error_cut=[M1_error[8],M2_error[8],M3_error[8],M4_error[8]]#T=2.45InM3=[np.log(M1[11]),np.log(M2[11]),np.log(M3[11]),np.log(M4[11])]M3_error_cut=[M1_error[11],M2_error[11],M3_error[11],M4_error[11]]print(InM2)print(M2_error_cut)#plot the M-T graphplt.figure()plt.errorbar(T,M1,M1_error,fmt=o,ecolor=gray,color=g,elinewidth_=2,capsize=4,label=L=4)plt.errorbar(T,M2,M2_error,fmt=o,ecolor=gray,color=b,elinewidth_=2,capsize=4,label=L=8)plt.errorbar(T,M3,M3_error,fmt=o,ecolor=gray,color=y,elinewidth_=2,capsize=4,label=L=16)plt.errorbar(T,M4,M4_error,fmt=o,ecolor=gray,color=r,elinewidth_=2,capsize=4,label=L=32)plt.ylabel(<m^2>)plt.xlabel(T/J)plt.grid(True)plt.legend()plt.show()#In(m^2)-In(L) Graphplt.figure()plt.errorbar(InL,InM1,M1_error_cut,fmt=o,ecolor=r,color=g,elinewidth_=2,capsize=4,label=T/J=2)plt.errorbar(InL,InM2,M2_error_cut,fmt=o,ecolor=r,color=b,elinewidth_=2,capsize=4,label=T/J=2.269)plt.errorbar(InL,InM3,M3_error_cut,fmt=o,ecolor=r,color=gold,elinewidth_=2,capsize=4,label=T/J=2.45)plt.plot(InL,InM1,c=g)plt.plot(InL,InM2,c=b)plt.plot(InL,InM3,c=gold)plt.xlabel(In(L))plt.ylabel(In(<m^2>))plt.grid(True)plt.legend()plt.show()#plot the M(T) graph for L=32plt.figure()plt.errorbar(T,M4,yerr=M4_error,fmt=o,ecolor=r,color=b,elinewidth_=2,capsize=4)plt.xlabel(Temperature / K)plt.ylabel( Magnetization )plt.title(2D Ising Model-M(T) Graph)plt.grid(True)plt.show()#線性擬合z=np.polyfit(InL,InM2,1)def y(x): return -0.25*x+0.08y_list=[]for i in range(len(InL)): y_list.append(y(InL[i]))plt.figure()plt.plot(InL,y_list)plt.errorbar(InL,InM2,M2_error_cut,fmt=o,ecolor=r,color=b,elinewidth_=2,capsize=4,label=T/J=2.269)plt.show()