1. Calculate the dependence of the average magnetization on the temper- ature T (in units of J/kB), using the Monte Carlo simulation code as the below script by Python. For one of the temperatures with non-zero magnetization show the dependence of the instantaneous magnetization on the Monte Carlo step, and indicate the initial equilibration period, and the part of the trajectory used to calculate the average.

2. Using the results from problem #1, evaluate the critical temperature Tc (in units of J/kB).

3. How does your result compare to the analytic solution? What are the sources of error?

4. Does your result for problem #2 depend on the size of the lattice? Show results to support your answer.

5. Show that the magnetization at T < Tc depends on the initial configu- ration of the system.

import matplotlib import matplotlib.pyplot as plt import numpy as np import math import sys

maglist = [] N = 20 T = 2 inUT=1/T lattice = np.zeros((N, N), dtype = int) for i in range(N): for j in range(N): lattice[i,j] = np.sign(2*np.random.rand()-1)

frozenlat = lattice Mn = sum(sum(frozenlat)) print(Mn)

print(Mn) print(lattice) energy = 0 for i in range(N): for j in range(N): WF=lattice[(i+1)%N,j]+lattice[(i-1)%N,i]+lattice[i,(j+1)%N]+lattice[i,(j-1)%N] energy = energy + lattice[i,j]*WF energy += lattice[i,j]*WF energy *= 0.5

nIter = 10000 log = open('log.dat','w') for itt in range(0,nIter): i = int(N*np.random.rand()) j = int(N*np.random.rand()) S = frozenlat[i,j] WF=frozenlat[(i+1)%N,j]+frozenlat[(i-1)%N,i]+frozenlat[i,(j+1)%N]+frozenlat[i,(j-1)%N] deltaE = 2 * S * WF if math.exp(-float(deltaE)*inUT) > np.random.rand(): lattice[i,j] = -S energy += deltaE Mn -= 2*S if itt%10 ==0: log.write("%d %f %f " %(itt,float(energy/(N*N)),float(Mn/(N*N))))

data = np.loadtxt('log.dat') plt.figure(1) plt.subplot(211) plt.title('PBC Monte Carlo Ising Model') plt.ylabel('Energy') plt.plot(data[:,0],data[:,1]) plt.subplot(212) plt.ylabel('Magnetization') plt.xlabel('Monte Carlo Step') plt.plot(data[:,0],data[:,2]) plt.show()