4x4 Ising 模型/ exact and Monte Carlo

Exact: all spin states all generated for Boltzmann sum

# ex_Ising_4x4.py;
# exec time, dt= 31.1 (PC216, Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz)
# RASP4B executing time=1m50.195s

import time
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

def SQL(L):
    N=L*L; nbo=N*2; bond=[[0,0] for i in range(nbo)]
    nn=0
    for j in range(L):
        for i in range(L):
            I=L*j+i; J=I+1; K=I+L
            bond[nn][0]=I; bond[nn][1]=J
            bond[N+nn][0]=I; bond[N+nn][1]=K
            if(i==L-1): J=L*j; bond[nn][0]=I; bond[nn][1]=J
            if(j==L-1): bond[N+nn][1]=bond[N+nn][1]-N
            #print('nn=',nn,I,J,bond[nn][:])
            nn+=1
    nnb=[[0,0,0,0] for i in range(N)];
    for i in range(N): nnb[i][0]=bond[i][1]; nnb[i][1]=bond[N+i][1]
    for i in range(N): i1=nnb[i][0]; i2=nnb[i][1]; nnb[i1][2]=i; nnb[i2][3]=i;
    #for i in range(N): print('i=',i,nnb[i])
    return bond,nnb

def ex_Ising(N,m):
    TS=[]
    for n in range(m):
        S=[]; #print('\n',n,'S=',end='')
        for i in range(N):
            ni=int(n/2**i)%2
            #print(ni,end=' ')
            nii=ni*2-1
            S.append(nii)
        TS.append(S)
    print('N=',N,m,' len(TS)=',len(TS))
    TE=[]; TM=[]; TM2=[]
    TT=np.linspace(0.2,8,40)
    #TT=[0.1,1]
    for T in TT:
        Z=0; E=0; mag=0; mag2=0;
        for S in TS:
            eng=0; M=0; M2=0
            for i in range(nbo):
                i0=bond[i][0]; i1=bond[i][1]; eng+=-S[i0]*S[i1]
            for i in range(N): M+=S[i]
            M2=M**2; pr=np.exp(-eng/T);
            Z+=pr; E+=eng*pr; mag+=M*pr; mag2+=M2*pr;
            #print('%6.2f %12.5e %12.5e %12.5e %12.5e' %(T,eng,pr,E,Z))
        Ea=E/Z/N; maga=mag/Z; mag2a=mag2/Z/N**2
        TE.append(Ea); TM.append(maga); TM2.append(mag2a)
        print('%6.2f %12.5f %12.5f %12.5f' %(T,Ea,maga,mag2a))
    return TT,TE,TM,TM2

L=4; N=L*L; m=2**N; nbo=N*2; print('L,N,nbo,m=',L,N,nbo,m)
bond,nnb=SQL(L); #input('bbbb')
t0=time.time()
TT,TE,TM,TM2=ex_Ising(N,m)
t1=time.time(); dt=round(t1-t0,3)
print('exec time, dt=',dt)

plt.subplot(121)
plt.title("2D Ising 4x4(Exact)")
plt.xlabel('T'); plt.ylabel('E')
plt.plot(TT,TE,'r--o',label='N=4x4')
plt.legend()
plt.subplot(122)
plt.title("2D Ising 4x4(Exact)")
plt.xlabel('T'); plt.ylabel('M2')
plt.plot(TT,TM2,'r--o',label='N=4x4')
plt.legend()
plt.savefig('ex_4x4_Ising.png')

Monte Carlo for 4x4 Ising

#--- MCMC.py

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import time

a='0.20 -2.00000 0.00000 1.00000 0.40 -2.00000 0.00000 1.00000 0.60 -1.99999 0.00000 0.99999 0.80 -1.99963 0.00000 0.99983 1.00 -1.99716 0.00000 0.99865 1.20 -1.98857 0.00000 0.99448 1.40 -1.96786 0.00000 0.98413 1.60 -1.92742 0.00000 0.96317 1.80 -1.85872 0.00000 0.92622 2.00 -1.75538 0.00000 0.86892 2.20 -1.61874 0.00000 0.79163 2.40 -1.46066 0.00000 0.70140 2.60 -1.29881 0.00000 0.60907 2.80 -1.14828 0.00000 0.52398 3.00 -1.01707 0.00000 0.45099 3.20 -0.90680 0.00000 0.39101 3.40 -0.81558 -0.00000 0.34274 3.60 -0.74030 -0.00000 0.30414 3.80 -0.67782 0.00000 0.27320 4.00 -0.62549 0.00000 0.24822 4.20 -0.58116 -0.00000 0.22784 4.40 -0.54318 -0.00000 0.21104 4.60 -0.51029 -0.00000 0.19704 4.80 -0.48153 0.00000 0.18524 5.00 -0.45614 0.00000 0.17520 5.20 -0.43354 0.00000 0.16658 5.40 -0.41328 0.00000 0.15911 5.60 -0.39500 -0.00000 0.15258 5.80 -0.37840 0.00000 0.14684 6.00 -0.36326 -0.00000 0.14176 6.20 -0.34938 -0.00000 0.13724 6.40 -0.33659 -0.00000 0.13318 6.60 -0.32478 0.00000 0.12953 6.80 -0.31382 -0.00000 0.12623 7.00 -0.30362 0.00000 0.12322 7.20 -0.29411 0.00000 0.12049 7.40 -0.28520 -0.00000 0.11798 7.60 -0.27685 -0.00000 0.11567 7.80 -0.26900 -0.00000 0.11355 8.00 -0.26160 0.00000 0.11158'

X=a.split()
Tx=[float(X[i]) for i in range(0,160,4)]
Ex=[float(X[i]) for i in range(1,160,4)]
print(len(X))
print('Tx=',Tx)
print('Ex=',Ex)

def SQL(L):
    N=L*L; nbo=N*2; bond=[[0,0] for i in range(nbo)]
    nn=0
    for j in range(L):
        for i in range(L):
            I=L*j+i; J=I+1; K=I+L
            bond[nn][0]=I; bond[nn][1]=J
            bond[N+nn][0]=I; bond[N+nn][1]=K
            if(i==L-1): J=L*j; bond[nn][0]=I; bond[nn][1]=J
            if(j==L-1): bond[N+nn][1]=bond[N+nn][1]-N
            #print('nn=',nn,I,J,bond[nn][:])
            nn+=1
    nnb=[[0,0,0,0] for i in range(N)];
    for i in range(N): nnb[i][0]=bond[i][1]; nnb[i][1]=bond[N+i][1]
    for i in range(N): i1=nnb[i][0]; i2=nnb[i][1]; nnb[i1][2]=i; nnb[i2][3]=i;
    return bond,nnb
            
            
def EM(S):
    eng=0; 
    for i in range(NBO): 
        i1=bond[i][0]; i2=bond[i][1]
        eng+=-S[i1]*S[i2]; 
    M=0;
    for i in range(N): M+=S[i]
    M2=M**2
    return eng,M,M2

def MCMC(S):
    for i in range(N):
        r=int(np.random.rand()*N)
        DE=0
        for n in range(4):
            DE+=S[nnb[r][n]]
        DE=DE*S[r]*2
        tr=np.exp(-DE/T)
        if(tr>=1): S[r]=-S[r]
        else: 
            rn=np.random.rand()
            if(tr>rn): S[r]=-S[r]
    return S

L=4; N=L*L; NBO=N*2; m=2**N; NMCS=10000;
S=[1 for i in range(N)]
bond,nnb=SQL(L)
eng,M,M2=EM(S)
print('initial S:',S,eng,M,M2)

TES=[]; TM2S=[]; t00=time.time()
for jN in range(1):
    TE=[]; TM=[]; TM2=[]
    TT=np.linspace(0.1,4,40)
    TT=Tx.copy(); jT=0
    for T in TT:
        Z=0; E=0; mag=0; mag2=0; t0=time.time()
        for j in range(NMCS):
            S=MCMC(S)
            eng,M,M2=EM(S)
            E+=eng; mag+=M; mag2=M2;
        Ea=E/N/NMCS; maga=mag/N/NMCS; mag2a=mag2/N/NMCS
        TE.append(Ea); TM.append(maga); TM2.append(mag2a)
        t1=time.time(); dt=round(t1-t0,4)
        sr='%6.2f %12.5f %12.5f %12.5f %12.5f %8.3f'
        print(sr %(T,Ea,Ex[jT],maga,mag2a,dt))
        jT+=1
    TES.append(TE); TM2S.append(TM2)
t1=time.time(); dta=round(t1-t00,4)
print('MC_YC run time=',dta)
#plt.xlim(0,5)
#plt.ylim(-1,-0.4)
plt.plot(TT,TES[0],'bo', label='MC')
plt.plot(TT,Ex,'r--', label='Exact')
plt.legend()
plt.savefig("Ising_4x4_YCMC.png")