import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation %matplotlib inline

# set some global parameters for the model l = 20 # number of spins per dimension n = l*l # total number of spins h = 0. # external field J = 1. # strength of the spin-spin coupling # let's use units in which k_B = 1. Tc = 2.269 # critical temperature for 2D Ising with J=1, h=0, k_B=1. T = 2.0 # initial temperature beta = 1/T # initial value of beta

# functions to initialize the configuration def init_lattice(l): return 2*np.random.randint(2, size=(l,l))-1 def draw_lattice(lattice_conf): plt.imshow(lattice_conf)

# sum over the neighbors of spin i # note that we use "periodic boundary conditions" so a spin at the right edge is # "neighboring" the a spin on the left edge. This ensures that there are not boundary effects. def compute_neighsum(lattice_conf, i, j): neighsum = lattice_conf[(i+1)%l,j] neighsum += lattice_conf[(i-1)%l,j] neighsum += lattice_conf[i,(j+1)%l] neighsum += lattice_conf[i,(j-1)%l] return neighsum # compute the total energy for the lattice model def compute_energy(lattice_conf): energy = 0. for i in range(l): for j in range(l): energy += -h * lattice_conf[i,j] energy += -0.5 * J * compute_neighsum(lattice_conf,i,j) * lattice_conf[i,j] # factor of 2 because the pair i,j is visited twice in the loop return energy

def compute_delta_e(lattice_conf, i, j): # implement this function! # return the energy difference associated with flipping # the spin indexed by i,j # Note: you do not need to compute the total energy! energy_difference = (-h * -lattice_conf[i,j]) - (-h * lattice_conf[i,j]) energy_difference += (-0.5 * J * compute_neighsum(lattice_conf,i,j) * -lattice_conf[i,j]) - (-0.5 * J * compute_neighsum(lattice_conf,i,j) * lattice_conf[i,j]) return energy_difference

def monte_carlo_step(lattice_conf, beta): # select a random site to flip i,j = np.random.randint(l, size=2) # compute the energy difference of flipping spin i,j delta_e = compute_delta_e(lattice_conf,i,j) # TODO! accept or reject the move using the Metropolis criterion p_acc = np.minimum(1, np.e**(-1 * beta * delta_e)) if np.random.random(1)[0] <= p_acc: lattice_conf[i,j] *= -1 return lattice_conf else: return lattice_conf

def run_mc_sweep(lattice_conf, beta): for attempt in range(n): lattice_conf = monte_carlo_step(lattice_conf, beta) return lattice_conf

def run_trajectory(lattice_conf, n_sweeps, beta): traj = np.zeros([n_sweeps, l, l], dtype=np.int) for sweep in range(n_sweeps): lattice_conf = run_mc_sweep(lattice_conf, beta) traj[sweep] = lattice_conf return traj

# example plotting the initial lattice conf = init_lattice(l) draw_lattice(conf)

# run a trajectory with 1000 mc sweeps n_sweeps = 1000 traj = run_trajectory(conf, n_sweeps, beta)

# You can run this function to visualize the trajectory fig = plt.figure() ax = plt.gca() im = ax.imshow(traj[0]) # set initial display dimensions def animate_func(i): im.set_array(traj[i]) return [im] ani = FuncAnimation(fig, animate_func)

def average_m(lattice_conf): mag = 0. for i in range(l): for j in range(l): mag += lattice_conf[i,j] return mag/n

def compute_running_m(traj, burn_in): step_averages = np.zeros(traj.shape[0]-burn_in) for i in range(traj.shape[0]-burn_in): step_averages[i] = average_m(traj[i+burn_in]) running_average= np.zeros(traj.shape[0]-burn_in) for i in range(traj.shape[0]-burn_in): j= i while (j >= 0): running_average[i] += step_averages[j] j += -1 running_average[i] = running_average[i]/(i+1) return running_average def compute_m(traj, burn_in): m = compute_running_m(traj,burn_in)[-1] return m def compute_avg_e(traj, burn_in=100): step_energies = np.zeros(traj.shape[0]-burn_in) for i in range(traj.shape[0]-burn_in): step_energies[i] = compute_energy(traj[i+burn_in]) average_energy = np.mean(step_energies) return average_energy def compute_var_e(traj, burn_in): step_energies = np.zeros(traj.shape[0]-burn_in) for i in range(traj.shape[0]-burn_in): step_energies[i] = compute_energy(traj[i+burn_in]) squared_energies = np.square(step_energies) average_squared_energies = np.mean(squared_energies) average_energy = np.mean(step_energies) squared_average_energy = average_energy**2 Cv = (T**(-2)) * ( average_squared_energies - squared_average_energy) return Cv

# TODO, plot the running average <M> as a function of the number of samples y= compute_running_m(traj,0) x=np.arange(y.size) plt.plot(x,y) plt.title('Average Magnitization as a function of the number to samples') plt.ylabel('Average Magnitization') plt.xlabel('Number of Samples') plt.show() u= compute_running_m(traj,100) x=np.arange(u.size) plt.plot(x,u) plt.title('Average Magnitization as a function of the number to samples with a burn in of 100') plt.ylabel('Average Magnitization') plt.xlabel('Number of Samples') plt.show()

# Example, compute the average of <M> for h in the range(-0.5, 0.5) T=2.0 beta = 1/T hs = np.linspace(-0.5, 0.5, 10) ms = np.zeros(10) i=0 for h_field in hs: h = h_field conf = init_lattice(l) # important to reinitialize traj = run_trajectory(conf, 500, beta) ms[i] = compute_m(traj, burn_in=100) i+=1

# TODO plot <M> vs. h plt.plot(hs,ms) plt.title('Average Magnitization as a function of h') plt.ylabel('Average Magnitization') plt.xlabel('h') plt.show()

# make sure to set h=0! #Since beta = 1/T Cv = 1/T^2(<E^2>-<E>^2) h=0. Ts = np.linspace(1.0, 4.0, 10) es = np.zeros(10) var_es = np.zeros(10) i=0 for Temp in Ts: T = Temp conf = init_lattice(l) traj = run_trajectory(conf, 500, 1/T) es[i] = compute_avg_e(traj, 100) var_es[i] = compute_var_e(traj, 100) i+=1

# TODO plot <E> vs T plt.plot(Ts,es) plt.plot(Ts,np.square(es)) plt.title('Average Energy as a Function of Temperature') plt.ylabel('Average Energy') plt.xlabel('Temperature') plt.show()

# TODO using the relation between var(E) and C_v, plot C_v vs. T plt.plot(Ts,var_es) plt.title('Heat Capacity as a function of Temperature') plt.ylabel('Heat Capacity') plt.xlabel('Temperature') plt.show()