Theoretical Physics Project 2-Ising model Solution

Project 2: Ising model
1. 1D Ising model 
Write a program that implements the Metropolis algorithm for single spin flips in a 1D Ising model of N spins producing a new spin configuration Xk+1 from the present spin configuration Xk. Use the canonical ensemble for a heat bath of temperature T. Set up periodic boundary conditions (effectively bending the chain into circle, such that the first and the last spin are adjacent to each other). Choose units such that J = 1. The thermal energy kBT is given in units of J.
a) Set N = 20, and kBT = 1 and simulate L = 500 individual trial spin flips. Begin from a “cold” initial state, where all spins are pointing to the same direction. Set the external magnetic field to H = 0. Visualise how the spin configuration changes with the number of trial spin flips (i.e. with time). It might be advantageous for the visualisation to show the configuration only after every 5th trial or so. What can you observe? (4 points)
b) Repeat the same simulation for three different thermal energies kBT = 0.1, 1 and 10 and keep track of the energy E at each trial spin flip. For a given configuration X = {s1,...,sN}, it is given by
N
E = −J Xsisi+1,
i=1
First run the simulation for at least L = 1000 trial spin flips. Plot the time evolution of the energy for each temperature. Your result will be strongly fluctuating. Improve it by repeating each simulation M = 100 times, each time starting with a different random seed. Then plot the time evolution of the mean energy ⟨E⟩ and its
Monte Carlo error estimate p(⟨E2⟩ − ⟨E⟩2)/M.
Discuss your result. When is equilibrium reached? What can you observe after the system has reached equilibrium? (4 points)
c) Make a plot of the mean energy per particle N1 ⟨E⟩t (averaged over simulation time t,
e.g. over L = 1000 trial spin flips) versus the thermal energy kBT after equilibrium has been reached (i.e. ignoring the first L ∼ 1000 trial spin flips). Choose kBT = 1...10 and average over M = 100 independent simulations to obtain a smooth result. Compare to the analytical result for the thermodynamic limit:
,
and discuss your results. (4 points)
d) Make a plot of the specific heat per particle at constant volume cV = N1 CV versus the thermal energy kBT after equilibrium has been reached. Use the same parameters as in c). As discussed in the lecture, you can calculate CV as follows:

Compare to the analytical result for the thermodynamic limit:
,
and discuss your results. (4 points)
e) Finally, make a plot of the magnetisation per particle m = N1 ⟨M⟩ versus the thermal energy kBT after equilibrium has been reached. Use the same parameters as in c), but repeat your measurement for different external magnetic fields H = 0, 0.1, 1 and 10. For a given configuration X, the magnetisation M is given by
N
M = Xsi.
i=1
Compare to the analytical result for the thermodynamic limit:
,
and discuss your results. (4 points)
f) (Bonus) A domain is a region of d adjacent spins that all point into the same direction. Make a plot of the average number of domains ⟨ndomains⟩ versus the thermal energy kBT, after equilibrium has been reached. Use the same parameters as in c), Discuss your result. (4 bonus points)
2. 2D Ising model (8 points)
While the 1D Ising model does not yet feature a phase transition, the 2D Ising model already does. This can be seen e.g. in the analytical result for the absolute magnetisation per particle for H = 0:
,
where Tc is the Curie temperature at which the phase transition occurs.√
Implement a 2D Ising model on a square lattice with n = N = 30 spins per side, again using periodic boundary conditions, sn+1,j = s1,j and si,n+1 = si,1, connecting opposite sides with each other. Calculate the absolute magnetisation of the model for H = 0 in dependence of the thermal energy kBT with the same techniques used for the 1D model. Note that more time might be needed to reach equilibrium. Discuss your results. Can you reproduce the predicted phase transition?
3. Bonus: Wang–Landau Sampling
Download the original paper on the Wang-Landau algorithm from Stud.IP ▶ Ubung¨ : Methods of Computational Physics ▶ Files ▶ Additional Material. Read it and answer the following questions:
a) (Bonus) Does the WLS-Algorithm satisfy Detailed Balance? Comment on the size of possible violations. (2 bonus points)
b) (Bonus) To obtain the absolute density of states the numerically determined relative density needs to be normalised. Quote two independent normalisation constraints that can be used for that purpose. (2 bonus points)
c) (Bonus) Quote a practical solution to quantify the ‘flatness’ of the energy density histogram on the fly. (2 bonus points)