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

# 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

Problem 2a

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! return 2*h*lattice_conf[i,j] + 2*J*compute_neighsum(lattice_conf,i,j)*lattice_conf[i,j] #This code below ran way too slow... took over an hour to finish #Kept it here just for my reference #new_lattice = copy.deepcopy(lattice_conf) #new_lattice[i,j] = -lattice_conf[i, j] #initial_energy = compute_energy(lattice_conf) #final_energy = compute_energy(new_lattice)

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 prob_acc = min([np.exp(-beta*delta_e), 1]) acc_param = np.random.uniform(0,1) if acc_param < prob_acc: lattice_conf[i,j] *= -1 return lattice_conf else: return lattice_conf

Problem 2b

After performing 'N' sweeps, we will have a lattice in a new/different configuration. This new configuration can represent how the lattice changed over some period of time. If we defined a 'new' lattice state after only one sweep, then the new lattice would differ from the original by only a single spin. It makes more sense to define a 'new' state after performing 'N' sweeps, because it is likely that multiple spins within the lattice will be changing at a time.

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 compute_m(traj, burn_in): return np.mean(traj[burn_in:]) def compute_running_m(traj, burn_in): vals = [] for i in range(burn_in, len(traj)): M_avg = compute_m(traj[burn_in:i], burn_in) vals.append(M_avg) return vals def compute_avg_e(traj, burn_in=100): es = [] for config in traj[burn_in:]: es.append(compute_energy(config)) return np.mean(es) def compute_var_e(traj, burn_in=100): es = [] for config in traj[burn_in:]: es.append(compute_energy(config)) return np.var(es)

#just for reference #plotting <M> (not running <M>) vs number of sweeps #y_vals = compute_running_m(traj, burn_in=100) #x_vals = np.arange(0, len(y_vals)) #plt.plot(x_vals, y_vals) #plt.ylabel('<M>') #plt.xlabel('N_sweeps') #<M> seems to tend towards zero, but running <M> plot does not

Problem 2c

The plot after tossing the first 100 sweeps looks much better. Perhaps the lattice configurations in the first 100 sweeps are too similar to the original lattice state. This may be a problem if the original randomly generated lattice was actually in a configuration that was energetically unfavorable or otherwise physically unreasonable. It would take some time for the lattice to revert to a more physically reasonable configuration.

# TODO, plot the running average <M> as a function of the number of samples burn_in_val = 100 toss_sweeps = compute_running_m(traj, burn_in_val) a_vals = np.arange(burn_in_val, len(toss_sweeps) + burn_in_val) b_vals = toss_sweeps plt.plot(a_vals, b_vals) plt.ylabel('Running <M>') plt.xlabel('N_sweeps') no_toss_sweeps = compute_running_m(traj, burn_in=0) c_vals = np.arange(0, len(no_toss_sweeps)) d_vals = no_toss_sweeps plt.plot(c_vals, d_vals)

/shared-libs/python3.7/py/lib/python3.7/site-packages/numpy/core/fromnumeric.py:3373: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
/shared-libs/python3.7/py/lib/python3.7/site-packages/numpy/core/_methods.py:170: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
/shared-libs/python3.7/py/lib/python3.7/site-packages/numpy/core/fromnumeric.py:3373: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
/shared-libs/python3.7/py/lib/python3.7/site-packages/numpy/core/_methods.py:170: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)

Problem 2d

By re-initializing for every point on our plot, we can the relationship between <M> and h as independent of the initial lattice configuration.

# 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, 100) i+=1

# TODO plot <M> vs. h plt.plot(hs, ms) plt.ylabel('<M>') plt.xlabel('h')

Problem 2e

# make sure to set h=0! 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.xlabel('T') plt.ylabel('<E>')

Problem 2f

The function is maximized at around T = 2.3. If we look back to the parameters we set at the beginning of this problem, this is the critical temperature! This is what we'd expect to see from a plot of C_v vs T.

# TODO using the relation between var(E) and C_v, plot C_v vs. T #C_v = 1/(kT**2) * var(E) plt.plot(Ts, var_es/(T**2)) plt.xlabel('T') plt.ylabel('C_v')