# Ising Monte Carlo simulation # Version 2.12.2021 # - latest: optional numba # Vesa Apaja # Disclaimer: only for educational purposes, very ineffective near Tc import numpy as np import sys import matplotlib.pyplot as plt from time import process_time as timer try: from numba import njit except: print('Warning: You have no numba module installed, so the program will run slower.') print(' Consider installing it. (Python 3.10 may have issues in installing numba).') # dummy decorator def njit(*args, **kwargs): def decorator(func): return func return decorator # local: from analyze import analyze # # Examples: # # ideal paramagnet # python3 ising.py 0 50 100 30 0.1 0.1 40 0.1 0.1 40 1 0 3 # # Ferromagnet, B=0, no lattice plots # python3 ising.py 1.0 100 1000 50 0.1 0.1 50 0 0 1 0 0 2 # # Fast from high-T to low T # python3 ising.py 1.0 100 10 20 5.0 -0.1 50 0 0 1 1 0 2 # # Ferromagnet, B=1 # python3 ising.py 1.0 50 1000 30 0.1 0.1 50 1 0 1 1 0 2 # # Ferromagnet, B=3 # python3 ising.py 1.0 50 1000 30 0.1 0.1 50 3 0 1 1 0 2 # # Antiferromagnet # python3 ising.py -1.0 50 100 30 0.1 0.1 100 0 0.1 1 1 0 3 # # Antiferromagnet from high-T to low-T # python3 ising.py -1.0 50 100 30 5.0 -0.1 100 0 0.1 1 2 0 3 # # Hysteresis # ---------- # far below Tc # increase B: # python3 ising.py 1.0 50 500 100 1.0 0.0 1 -5.0 0.1 100 0 0 3 # decrease B # python3 ising.py 1.0 50 500 100 1.0 0.0 1 5.0 -0.1 100 0 0 2 # just below Tc # increase B: # python3 ising.py 1.0 50 500 100 2.2 0.0 1 -5.0 0.1 100 0 0 3 # decrease B # python3 ising.py 1.0 50 500 100 2.2 0.0 1 5.0 -0.1 100 0 0 2 # just above Tc # increase B: # python3 ising.py 1.0 50 500 100 2.4 0.0 1 -5.0 0.1 100 0 0 3 # decrease B # python3 ising.py 1.0 50 500 100 2.4 0.0 1 5.0 -0.1 100 0 0 2 # boiling # increase temperature: # python3 ising.py 1 100 1000 0 0.8 0.01 100 -1 0 1 2 0 2 # # decrease pressure: # python3 ising.py 1 100 500 0 1 0 1 -0.8 -0.01 100 2 0 2 # arguments value in the example meaning # --------- -------------------- ------- # J 1.0 coupling constant # L 20 lattice size (L x L) # MCSTEPS 100 # of MC steps # THERM 30 # of thermalization steps # T 0.1 first temperature # DT 0.1 temperature step # Tnum 50 # of temperature values # B 0 first magnetic field value # DB 0 magnetic field step # Bnum 1 # of magnetic field values # plotfreq 1 plotting frequency # pauses 0 pause between plots # init 1 lattice initialization # # If you have problems, ask Vesa Apaja, vesa.apaja@jyu.fi (Nanoscience Center 2nd floor) # # fix random seed for reproducibility #np.random.seed(7) def print_ETA(time): if time>60: print(' Ready in {:d} minutes'.format(int(time/60.))) else: if int(time)>1: print(' Ready in {:d} seconds'.format(int(time))) else: print(' Ready in no time') @njit() def measure(J,L,lattice,xleft,yup,B,E,E2,M,M2,ABS_M): Sx = 0 Sy = 0 for i1 in range(L): for i2 in range(L): Sy += lattice[i2,yup[i1]]*lattice[i2,i1] Sx += lattice[xleft[i2],i1]*lattice[i2,i1] S = np.sum(lattice) energy = -J*Sx-J*Sy-B*S E += energy E2 += energy**2 M += S M2 += S**2 ABS_M += abs(S) return E,E2,M,M2,ABS_M @njit() def sweep(L,lattice,xleft,xright,yup,ydown,w): for kk in range(L*L): # sweep over the lattice # pick random lattice point i = np.random.randint(L) j = np.random.randint(L) spinX=lattice[xleft[i],j]+lattice[xright[i],j] spinY=lattice[i,yup[j]]+lattice[i,ydown[j]] ijspin=lattice[i,j] spin= int((ijspin+3)/2)-1 if spinX==-2: if spinY==-2: nnsumma=0 elif spinY==0: nnsumma=1 elif spinY==2: nnsumma=2 elif spinX==0: if spinY==-2: nnsumma=3 elif spinY==0: nnsumma=4 elif spinY==2: nnsumma=5 elif spinX==2: if spinY==-2: nnsumma=6 elif spinY==0: nnsumma=7 elif spinY==2: nnsumma=8 if w[nnsumma,spin] > np.random.random(): # accept split flip lattice[i,j] = -ijspin """ Wolff algorithm for a single cluster: 1. Choose random site i. 2. Find neighbouring sites j If S[j]=S[i], join j to cluster with probability p = 1-exp(-2.0/T) 3. If site j was joined to te cluster, mark i=j and repeat step 2 until cluster stops growing 4. Invert spins in the cluster """ @njit() def wolff_sweep(L,p,S,nbr): for rep in range(10): k = np.random.randint(0, L*L - 1) Pocket, Cluster = [k], [k] while len(Pocket)>0: j = Pocket[np.random.randint(len(Pocket))] for l in nbr[j]: if np.random.uniform(0.0, 1.0) < p: if S[l] == S[j]: if l not in Cluster: Pocket.append(l) Cluster.append(l) Pocket.remove(j) for j in Cluster: S[j] *= -1 def ising(J,L,MCSTEPS,THERM,T,DT,Tnum,B,DB,Bnum,plotfreq,pauses,init): print('You gave parameters:') print(' J = ',J) print(' L = ',L) print(' MCSTEPS = ',MCSTEPS) print(' THERM = ',THERM) print(' T = ',T) print(' DT = ',DT) print(' Tnum = ',Tnum) print(' B = ',B) print(' DB = ',DB) print(' Bnum = ',Bnum) print('plotfreq = ',plotfreq) print(' pauses = ',pauses) print(' init = ',init) print('') print('Info:') print(' - This program will make plots according to the parameter plotfreq') print(' - Some plots may be utterly useless') print(' - Not all bad parameter choices are trapped in the beginning, the run may crash') print(' - All results except the lattice spins will be saved to an ascii data file results.dat,') print(' which you can use in any plotting program you wish.') print(' - You may want to backup the the old results.dat file now, it will be overwritten') ideal = False if abs(J)<0.0001: ideal = True # # spin values in the simulation # ispin = np.array([-1,1]) mu = 2 # sanity tests # if not init in [0,1,-1,2,-2,3,4,5]: print('bad init lattice, choose 0, 1, -1, 2, -2, 3, 4 or 5') return if not plotfreq in [0,1,2]: print('bad value for plotfreq, choose 0,1, or 2') return if not pauses in [0,1]: print('bad value for pauses (not 0 or 1); assuming no pauses') pauses = 0; print('') if ideal: print('J=0 : This is an IDEAL PARAMAGNET') print('Program uses parameter x=mu B/(kB T) and keeps T=1 fixed') Tnum = 1 T = 1 dT = 0 #B = 0 #DB = 0.2*mu; #Bnum = 16; if Bnum<=1: print('WARNING: Just one x value (taken from B table) will not show much.') if J>0: print('coupling constant J = ',J); print('J>0: This is a FERROMAGNETIC system') if J<0: print('coupling constant J = ',J); print('J<0: This is an ANTIFERROMAGNETIC system') print('Lattice dimensions ',L,'x',L) # # tabulate T and B values # Bs = np.array([B+i*DB for i in range(Bnum)]) # B's Ts = np.array([T+i*DT for i in range(Tnum) if T+i*DT >0.0 ]) # T's if ideal: print('Table of x = values') #multiply B's with mu Bs = mu*Bs print(', '.join(f'{B:.2f}' for B in Bs)) else: print('Table of magnetic field values') print(Bs) print('Table of temperature values') print(', '.join(f'{T:.2f}' for T in Ts)) if init==0: print('init=0: setting spins to random values +/- 1 between MC runs') if init==1: print('init=1: resetting spins up between MC runs') if init==-1: print('init=-1: resetting spins down between MC runs') if init==2: print('init=2: keeping previous spin configuration between MC runs (the first one is spins up)') if init==-2: print('init=-2: keeping previous spin configuration between MC runs (the first one is spins down)') if init==3: print('init=3: keeping previous spin configuration between MC runs (the first one is random spins)') if init==5: print('init=5, keeping previous spin configuration between MC runs (the first one is half up, half down)') if plotfreq == 0: print('plotfreq = 0, not making plots during MC runs') elif plotfreq==1: print('making plots after every MC cycle: this slows down the computation but may be useful') elif plotfreq==2: print('making plots after every MC sweep: this slows down the computation but may be useful') else: print('Bad plotfreq (not 0,1)') return #input("Press Enter to continue...") print('========== END OF INPUT PARAMETERS =======================================') # # Nearest neighbors as a dictionary # --------------------------------- # lattice is # 0 1 2 ... L-1 # L L+1 L+2 ... 2L-1 # ... # (L-1)*L ... L*L-1 # # with periodic boundaries, each point has 4 neightbors, # so the neighbors are (linear indicing) N = L*L nbr_dict = {i : ((i // L) * L + (i + 1) % L, (i + L) % N, (i // L) * L + (i - 1) % L, (i - L) % N) for i in range(N)} # convert to list nbr =[] for i in range(N): nbr.append(nbr_dict[i]) nbr = np.array(nbr) # checking if(False): print (nbr) print('nbr vs lattice checking') exit(0) ####################################################### # translation vectors for periodic boundary conditions ####################################################### xright = np.zeros((L),dtype=int) xleft = np.zeros((L),dtype=int) for a1 in range(L): xright[a1] = a1+1 xleft[a1] = a1-1 xright[L-1] = 0 xleft[0] = L-1 yup = np.zeros((L),dtype=int) ydown = np.zeros((L),dtype=int) for a1 in range(L): yup[a1] = a1+1 ydown[a1] = a1-1 yup[L-1] = 0 ydown[0] = L-1 ############################################################# nnsum = np.array([-2,0,2,-2,0,2]) MCSTEP= range(MCSTEPS) print('========== Monte Carlo begins ============================================') # # start lattice # if init==1 or init==2: lattice=np.ones((L,L),dtype=int) print('first lattice is spins up') elif init==0 or init==3: lattice=2*np.random.randint(2,size=(L,L))-1; print('first lattice has random spins') elif init==-1 or init==-2: lattice=-np.ones((L,L),dtype=int) print('first lattice is spins down') elif init==4: print('first lattice from file') print(' no done') exit(1) elif init==5: print('first lattice is half up half down') lattice = np.ones((L,L),dtype=int) lattice[L/2:L,:] = -lattice[L/2:L,:] w = np.zeros((3*3,2)) dEperkT = np.zeros((3*3,2)) ENERGY_Y = np.zeros((L*L)) ENERGY_X = np.zeros((L*L)) res = [] # initialize plot # --------------- if plotfreq>0: latticefig=plt.figure(1) ax = latticefig.gca() latticeplot = ax.imshow(lattice,vmin=-1.0, vmax=1.0) ax.set_axis_off() T_out = "{:10.3f}".format(T) B_out = "{:10.3f}".format(B) tit = 'T = '+str(T_out)+' B = '+str(B_out) ax.set_title(tit) latticeplot.set_data(lattice) plt.show(block=False) #plt.draw() plt.pause(1e-3) for T in Ts: # T loop Ttic = timer() for B in Bs: # B loop Btic = timer() print(50*'-') if abs(J)<0.0001: print('mu B/(kB T) = {:.3f}'.format(B/mu)) else: print('T = {:.3f} B = {:.3f}'.format(T,B)) BperkT = B/T JxperkT = J/T JyperkT = J/T # energy from nearest neighbour spin sums (nnsum) q_apu = 0 for q1 in range(3): # x direction for q2 in range(3): # y direction for q3 in range(2): # spin values up and down w[q_apu,q3] = 1 dEperkT[q_apu,q3]=2*JxperkT*nnsum[q1]*ispin[q3]+BperkT*ispin[q3]+2*JyperkT*nnsum[q2]*ispin[q3] ; if dEperkT[q_apu,q3] > 0: w[q_apu,q3] = np.exp(-dEperkT[q_apu,q3]); q_apu +=1 # # initial lattice: random, all spins up (1) or keep previous # if init==0: lattice=2*np.random.randint(2,size=(L,L))-1 print(' init lattice to random spins') if init==1: lattice=np.ones((L,L)) print(' init lattice to spin ups') if abs(init)>=2: pass # init : do nothing here M = 0.0 E = 0.0 E2 = 0.0 M2 = 0.0 ABS_M = 0.0 p = 1.0 - np.exp(-2.0*J/T) # Wolff algorithm probability # MC loop # ======= for mcs in range(-THERM,MCSTEPS): # Loop over MC steps if True: sweep(L,lattice,xleft,xright,yup,ydown,w) else: # OPTIONAL: Wolff algorithm S = lattice.reshape(L*L) wolff_sweep(L,p,S,nbr) lattice = S.reshape(L,L) # still thermalizing ? if mcs < 0: continue # measure and add to sums E,E2,M,M2,ABS_M = measure(J,L,lattice,xleft,yup,B,E,E2,M,M2,ABS_M) # snapshot between MC steps if plotfreq==2: latticeplot.set_data(lattice) T_out = "{:10.3f}".format(T) B_out = "{:10.3f}".format(B) tit = 'T = '+str(T_out)+' B = '+str(B_out) ax.set_title(tit) plt.show(block=False) #plt.draw() plt.pause(1e-3) # ---------------------- # MC loop ends # norm = L*L*MCSTEPS E /= norm M /= norm ABS_M /= norm norm = (L*L)**2*MCSTEPS E2 /= norm M2 /= norm print(' M = {:.4f} E = {:.4f}'.format(M,E)) # collect results for later analysis res.append([T,B,M,M2,ABS_M,E,E2]) # snapshot after MC steps if plotfreq==1: latticeplot.set_data(lattice) T_out = "{:10.3f}".format(T) B_out = "{:10.3f}".format(B) tit = 'T = '+str(T_out)+' B = '+str(B_out) ax.set_title(tit) plt.draw(),plt.pause(1e-3) # timing info # toc = timer() if len(Bs)>1: time = toc-Btic time = time*(Bs[-1]-B)/DB elif len(Ts)>1: time = toc-Ttic time = time*(Ts[-1]-T)/DT print_ETA(time) #------------- # B loop ends # ----------- # T loop ends print('========== Monte Carlo ends ==============================================') print('analysis and plotting the data follows') analyze(ideal,L,MCSTEPS,J,Ts,Bs,res) if __name__ == '__main__': args = sys.argv[1:] #print("args=",args) try: J,L,MCSTEPS,THERM,T,DT,Tnum,B,DB,Bnum,plotfreq,pauses,init = args J = float(J) L = int(L) MCSTEPS = int(MCSTEPS) THERM = int(THERM) T = float(T) DT = float(DT) Tnum = int(Tnum) B = float(B) DB = float(DB) Bnum = int(Bnum) plotfreq = int(plotfreq) pauses = int(pauses) init = int(init) except: if len(sys.argv)!= 13: print ('Wrong number of arguments: got ', len(sys.argv), ' expecting 13') print ('Argument List:', str(sys.argv[1:])) print('usage:') print('python3 sys.argv[0] J L MCSTEPS THERM T DT Tnum B DB Bnum plotfreq pauses init') exit(1) ising(J,L,MCSTEPS,THERM,T,DT,Tnum,B,DB,Bnum,plotfreq,pauses,init)