Commit 6dd1d73f by engelerp

### Solutions exercise 6

parent d8a0b8f0
 import numpy as np from timeit import default_timer as timer import random import matplotlib.pyplot as plt #does the next level of binning def binning_step(data): newdata = [] for i in range(int(len(data)/2)): newdata.append((data[2*i]+data[2*i+1])/2.) stdev = np.std(newdata)/np.sqrt(len(newdata)-1) return [stdev,newdata] #perform a binning analysis of the data, up to the level where at least minsamples samples remain def binning(data, minsamples=8): levels = [] stdevs = [] stdevs.append(np.std(data)/np.sqrt(len(data)-1)) levels.append(0) level = 1 while(len(data) >= 2*minsamples): stdev,data = binning_step(data) stdevs.append(stdev) levels.append(level) level += 1 return [levels, stdevs] #perform a Wolff cluster update step on sys of size H*L #absolute value of return is the cluster size, twice the return is the change in magnetisation def Wolff_step(sys,H,L,probability): flipped = 1 #choose random spin spin_v = random.randint(0,H-1) spin_h = random.randint(0,L-1) #flip the spin sys[spin_v][spin_h] *= -1 #keep track of what we do to the magnetisation flip_sign = sys[spin_v][spin_h] #check neighbours recursively spins_to_check = [[spin_v,spin_h]] while(len(spins_to_check) > 0): indices = spins_to_check.pop() spin_v = indices[0] spin_h = indices[1] current_spin = sys[spin_v][spin_h] #add neighbours to cluster and schedule them with correct probability if(current_spin * sys[(spin_v+1)%H][spin_h] == -1 and random.random() < probability): spins_to_check.append([(spin_v+1)%H,spin_h]) sys[(spin_v+1)%H][spin_h] *= -1 flipped += 1 if(current_spin * sys[(spin_v-1)%H][spin_h] == -1 and random.random() < probability): spins_to_check.append([(spin_v-1)%H,spin_h]) sys[(spin_v-1)%H][spin_h] *= -1 flipped += 1 if(current_spin * sys[spin_v][(spin_h+1)%L] == -1 and random.random() < probability): spins_to_check.append([spin_v,(spin_h+1)%L]) sys[spin_v][(spin_h+1)%L] *= -1 flipped += 1 if(current_spin * sys[spin_v][(spin_h-1)%L] == -1 and random.random() < probability): spins_to_check.append([spin_v,(spin_h-1)%L]) sys[spin_v][(spin_h-1)%L] *= -1 flipped += 1 return flipped*flip_sign #perform one single spin update on sys of extent H*L. #returns change in magnetisation def MCMC_step(sys,H,L,exptable): #pick random spin spin_v = random.randint(0,H-1) spin_h = random.randint(0,L-1) #calculate local field h = sys[spin_v][(spin_h-1)%L] + sys[spin_v][(spin_h+1)%L] + sys[(spin_v+1)%H][spin_h] + sys[(spin_v-1)%H][spin_h] #flip with metropolis probability if(random.random() < exptable[sys[spin_v][spin_h]*h]): if(sys[spin_v][spin_h] == 1): sys[spin_v][spin_h] = -1 return -2 else: sys[spin_v][spin_h] = 1 return 2 else: return 0 #returns the magnetisation of sys def magnetisation(sys): magn = 0 L = len(sys[0]) H = len(sys) for a in sys: for s in a: magn += s return float(magn)/float(H*L) #run a simulation using cluster updates def run_cluster_simulation(dT,Tmin=1.8,Tmax=3.2,meas=8192*2,equil=1000*2): ####################### ### Initialisation ####################### #arrays to store info about correlations statistics_magnetisation = [] statistics_magnetisationq = [] #arrays to store the expectation values and their standard deviations for all temperatures magnetisations = [] magnetisationsq = [] magnetisations_std = [] magnetisationsq_std = [] temperatures = [] #define system parameters L = 16 H = 16 J = 1. #initialise the system sys = [] for i in range(H): a = [] for j in range(L): a.append(1) sys.append(a) #initialise the magnetisation total_magnetisation = L*H T = Tmin ####################### ### Simulation ####################### #simulate for the specified temperatures while(T <= Tmax): #arrays to store all samples magdata = [] #magnetisation magqdata = [] #magnetisation squared T += dT beta = 1./T probability = 1. - np.exp(-2.*beta) #probability of accepting the addition of a site to the cluster #equilibrate for i in range(int(equil)): total_magnetisation += 2*Wolff_step(sys,H,L,probability) #measure for i in range(int(meas)): wolff_return = Wolff_step(sys,H,L,probability) m = np.abs(wolff_return)/float(L*H) total_magnetisation += 2*wolff_return magdata.append(m) magqdata.append((total_magnetisation/float(L*H))**2) #extract means and their standard deviations magnetisations.append(np.mean(magdata)) magnetisationsq.append(np.mean(magqdata)) magnetisations_std.append(np.std(magdata)) magnetisationsq_std.append(np.std(magqdata)) #save temperature temperatures.append(T) #do a binning analysis of the data and store it statistics_magnetisation.append(binning(magdata)) statistics_magnetisationq.append(binning(magqdata)) #calculate correlation times for all temperatures from the binning results. binning level 5 seems reasonable for most temperatures taus = [] tausq = [] for i in range(len(temperatures)): taus.append(0.5*(statistics_magnetisation[i][1][5]**2/magnetisations_std[i]**2*(meas-1.)-1.)) tausq.append(0.5*(statistics_magnetisationq[i][1][5]**2/magnetisationsq_std[i]**2*(meas-1.)-1.)) ####################### ### Plotting ####################### #correlation time magnetisation plt.plot(temperatures,taus) plt.title("Correlation time Wolff, magnetization") plt.xlabel("Temperature") plt.ylabel("Correlation time") #plt.show() plt.savefig("wolff_corr_times_m.pdf") plt.clf() #correlation time magnetisation squared plt.plot(temperatures,taus) plt.title("Correlation time Wolff, magnetization\$^2\$") plt.xlabel("Temperature") plt.ylabel("Correlation time") #plt.show() plt.savefig("wolff_corr_times_mq.pdf") plt.clf() #binning analysis results for i in range(len(statistics_magnetisation)): plt.plot(statistics_magnetisation[i][0],statistics_magnetisation[i][1]) leg = [] for i in range(len(temperatures)): leg.append("T=" + str(round(temperatures[i],3))) plt.legend(leg) plt.title("Binning analysis Wolff, magnetisation") plt.xlabel("Binning Level") plt.ylabel("Error") #plt.show() plt.savefig("wolff_binning_m.pdf") plt.clf() for i in range(len(statistics_magnetisationq)): plt.plot(statistics_magnetisationq[i][0],statistics_magnetisationq[i][1]) leg = [] for i in range(len(temperatures)): leg.append("T=" + str(round(temperatures[i],3))) plt.legend(leg) plt.title("Binning analysis Wolff, magnetisation\$^2\$") plt.xlabel("Binning Level") plt.ylabel("Error") #plt.show() plt.savefig("wolff_binning_mq.pdf") plt.clf() #calculate errorbars #binning level 5 seems reasonable for most temperatures errorbars_m = [] errorbars_mq = [] for i in range(len(temperatures)): errorbars_m.append(statistics_magnetisation[i][1][5]) errorbars_mq.append(statistics_magnetisationq[i][1][5]) #plot magnetisation with errorbars plt.errorbar(temperatures,magnetisations,yerr=errorbars_m,fmt='o',ms=1) plt.title("Magnetisation, Wolff") plt.xlabel("Temperature") plt.ylabel("Magnetisation") #plt.show() plt.savefig("wolff_m.pdf") plt.clf() #plot magnetisation squared with errorbars plt.errorbar(temperatures,magnetisationsq,yerr=errorbars_mq,fmt='o',ms=1) plt.title("Magnetisation\$^2\$, Wolff") plt.xlabel("Temperature") plt.ylabel("Magnetisation^2") #plt.show() plt.savefig("wolff_mq.pdf") plt.clf() #run a simulation using single spin updates def run_metropolis_simulation(dT,Tmin=1.8,Tmax=3.2,meas=131072*4,equil=8192*2): ####################### ### Initialisation ####################### #arrays to store info about correlations statistics_magnetisation = [] statistics_magnetisationq = [] #arrays to store the expectation values and their standard deviations for all temperatures magnetisations = [] magnetisationsq = [] magnetisations_std = [] magnetisationsq_std = [] temperatures = [] #define system parameters L = 16 H = 16 J = 1. #initialise the system sys = [] for i in range(H): a = [] for j in range(L): a.append(1) sys.append(a) #initialise the magnetisation total_magnetisation = L*H T = Tmin ####################### ### Simulation ####################### #simulate for the specified temperatures while(T <= Tmax): #arrays to store all samples magdata = [] #magnetisation magqdata = [] #magnetisation squared T += dT beta = 1./T exptable = [1.,-1.,np.exp(-4.*J*beta),-1.,np.exp(-8.*J*beta),1.,-1.,1.,-1.] #lookuptable for the exponents we may need #equilibrate for i in range(int(equil)): total_magnetisation += MCMC_step(sys,H,L,exptable) #measure for i in range(int(meas)): total_magnetisation += MCMC_step(sys,H,L,exptable) m = np.abs(total_magnetisation)/float(L*H) magdata.append(m) magqdata.append((m)**2) #extract means and their standard deviations magnetisations.append(np.mean(magdata)) magnetisationsq.append(np.mean(magqdata)) magnetisations_std.append(np.std(magdata)) magnetisationsq_std.append(np.std(magqdata)) #save temperature temperatures.append(T) #do a binning analysis of the data and store it statistics_magnetisation.append(binning(magdata)) statistics_magnetisationq.append(binning(magqdata)) #calculate correlation times for all temperatures from the binning results. binning level 14 seems okay, but still not good at certain temperatures (CSD) taus = [] tausq = [] for i in range(len(temperatures)): taus.append(0.5*(statistics_magnetisation[i][1][14]**2/magnetisations_std[i]**2*(meas-1.)-1.)) tausq.append(0.5*(statistics_magnetisationq[i][1][14]**2/magnetisationsq_std[i]**2*(meas-1.)-1.)) ####################### ### Plotting ####################### #correlation time magnetisation plt.plot(temperatures,taus) plt.title("Correlation time Metropolis, magnetization") plt.xlabel("Temperature") plt.ylabel("Correlation time") #plt.show() plt.savefig("metropolis_corr_times_m.pdf") plt.clf() #correlation time magnetisation squared plt.plot(temperatures,taus) plt.title("Correlation time Metropolis, magnetization\$^2\$") plt.xlabel("Temperature") plt.ylabel("Correlation time") #plt.show() plt.savefig("metropolis_corr_times_mq.pdf") plt.clf() #binning analysis results for i in range(len(statistics_magnetisation)): plt.plot(statistics_magnetisation[i][0],statistics_magnetisation[i][1]) leg = [] for i in range(len(temperatures)): leg.append("T=" + str(round(temperatures[i],3))) plt.legend(leg) plt.title("Binning analysis Metropolis, magnetisation") plt.xlabel("Binning Level") plt.ylabel("Error") #plt.show() plt.savefig("metropolis_binning_m.pdf") plt.clf() for i in range(len(statistics_magnetisationq)): plt.plot(statistics_magnetisationq[i][0],statistics_magnetisationq[i][1]) leg = [] for i in range(len(temperatures)): leg.append("T=" + str(round(temperatures[i],3))) plt.legend(leg) plt.title("Binning analysis Metropolis, magnetisation\$^2\$") plt.xlabel("Binning Level") plt.ylabel("Error") #plt.show() plt.savefig("metropolis_binning_mq.pdf") plt.clf() #calculate errorbars #binning level 14 seems reasonable for most temperatures errorbars_m = [] errorbars_mq = [] for i in range(len(temperatures)): errorbars_m.append(statistics_magnetisation[i][1][14]) errorbars_mq.append(statistics_magnetisationq[i][1][14]) #plot magnetisation with errorbars plt.errorbar(temperatures,magnetisations,yerr=errorbars_m,fmt='o',ms=1) plt.title("Magnetisation, Metropolis") plt.xlabel("Temperature") plt.ylabel("Magnetisation") #plt.show() plt.savefig("metropolis_m.pdf") plt.clf() #plot magnetisation squared with errorbars plt.errorbar(temperatures,magnetisationsq,yerr=errorbars_mq,fmt='o',ms=1) plt.title("Magnetisation\$^2\$, Metropolis") plt.xlabel("Temperature") plt.ylabel("Magnetisation^2") #plt.show() plt.savefig("metropolis_mq.pdf") plt.clf() #run run_metropolis_simulation(0.1) run_cluster_simulation(0.1)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment