Lab 12: MCMC Tutorial#
In this lab, we will study the 2D Ising model using Markov chain Monte Carlo. We demonstrate how to implement the Metropolis-Hastings algorithm to sample the Ising model at different temperatures and study the phase transition.
2D Ising Model#
The set of lattice sites is \(\Lambda\), usually a square grid. At each site, where \(\sigma_k \in \{ -1, +1\}\)is the spin at each site \(k\). For any two adjacent sites \(\langle i j \rangle\) (with no double counting) there is a spin interaction \(J\). The energy of the system is given by
\[E(\sigma) = -J\sum_{\langle ij\rangle} \sigma_i\sigma_j\]
The probability of a given configuration \(\sigma\) is given by the Boltzmann distribution
\[p(\sigma) = \frac{1}{Z}\exp \left(-\frac{E(\sigma)}{k_\mathrm{B}T}\right),\]
where \(Z = \sum_{\sigma} \exp \left(-\frac{E(\sigma)}{k_\mathrm{B}T}\right)\) is the partition function, \(T\) is the temperature, and \(k_\mathrm{B}\) is the Boltzmann constant. We assume periodic boundary conditions.
import numpy as np
import matplotlib.pyplot as plt
import tqdm
Energy of a spin configuration#
We will do a warmup problem by computing the energy of the Ising model spin configuration shown in lecture. Blue squares represent spins pointing up, and red squares represent spins pointing down. First, we initialize the lattice with the following spin configuration:
N = 10
J = 1
lattice_spins = np.zeros((N, N))
for i in range(N):
for j in range(N):
if i < 5:
lattice_spins[i, j] = -1
else:
lattice_spins[i, j] = 1
lattice_spins[1, 3] = 1
plt.figure()
plt.imshow(lattice_spins, cmap="RdYlBu")
# show gridlines
for i in range(N):
plt.axhline(i + 0.5, color="black", lw=0.1)
plt.axvline(i + 0.5, color="black", lw=0.1)
plt.show()
energy = 0
for i in range(N):
for j in range(N):
for k, l in [(-1, 0), (1, 0), (0, 1), (0, -1)]: #steps to take when navigating through the map, to obtain the sites of neighbors
i_neigh = i + k if i + k < N else 0
j_neigh = j + l if j + l < N else 0
energy += -J * lattice_spins[i, j] * lattice_spins[i_neigh, j_neigh]
energy = energy / 2 #to avoid double counting
print(f"{energy=}")
energy=np.float64(-152.0)
MCMC Sampling#
As discussed in lecture, we will use the Metropolis-Hastings algorithm to sample the Ising model at a given temperature.
As discussed in class the critical temperature for the 2D Ising model is \(T_\mathrm{c}=\frac{2}{\ln(1+\sqrt{2})}\approx 2.269\) in units where \(J=k_\mathrm{B}=1\). For \(T<T_\mathrm{c}\), the system is in the ordered phase (spins tend to align), and for \(T>T_\mathrm{c}\), the system is in the disordered phase. (The order parameter of this phase transition is M)
Let’s start with \(T=2\), a \(100\times100\) lattice, and 1 million MCMC steps.
To apply the MCMC method to the Ising model, we design a Markov process using the Metropolis algorithm as follows
On step \(k\), randomly choose one of the spins \(\sigma_k\) and consider flipping it \(\sigma_k \to -\sigma_k\).
Calculate the change in energy that would result from flipping spin \(k\), i.e. the quantity: $\(\Delta E = - \left[J \sum_{\langle ij\rangle} \sigma_j\right] \Delta \sigma_i,\)\( where \)\Delta \sigma_i = -2\sigma_i$ is the change due to the spin flip
If \(\Delta E \leq 0\), accept the spin flip.
If \(\Delta E > 0\), accept the spin flip with probability \(\exp(-\Delta E/k_\mathrm{B}T)\).
Update the moving average of whatever quantity we are interested in.
Repeat.
STEPS = 1_000_000
N = 100
J = 1
KB = 1
T = 2
# initialize random lattice
# lattice_spins = 2 * (np.random.randint(2, size=(N, N)) - 0.5)
# initialize ordered lattice with all spins up
lattice_spins = np.ones((N, N))
num_accept = 0
m_values = []
for t in tqdm.tqdm(range(STEPS)):
i, j = np.random.randint(N), np.random.randint(N)
# we only need to consider the neighbors of
# (i, j) to calculate the change in energy
delta_energy = 0
for k, l in [(-1, 0), (1, 0), (0, 1), (0, -1)]:
i_neigh = i + k if i + k < N else 0
j_neigh = j + l if j + l < N else 0
delta_energy += -J * -2 * lattice_spins[i, j] * lattice_spins[i_neigh, j_neigh]
if delta_energy <= 0:
lattice_spins[i, j] *= -1
num_accept += 1
elif delta_energy > 0:
prob = np.exp(-delta_energy / (KB * T))
if np.random.random() < prob:
lattice_spins[i, j] *= -1
num_accept += 1
m_values.append(np.mean(lattice_spins))
print(f"acceptance_rate = {num_accept/STEPS}")
0%| | 0/1000000 [00:00<?, ?it/s]
1%| | 5592/1000000 [00:00<00:17, 55913.12it/s]
1%| | 11184/1000000 [00:00<00:17, 54941.75it/s]
2%|▏ | 16680/1000000 [00:00<00:17, 54877.60it/s]
2%|▏ | 22553/1000000 [00:00<00:17, 56386.91it/s]
3%|▎ | 28194/1000000 [00:00<00:17, 56170.22it/s]
3%|▎ | 33988/1000000 [00:00<00:17, 56764.98it/s]
4%|▍ | 39666/1000000 [00:00<00:16, 56631.54it/s]
5%|▍ | 45366/1000000 [00:00<00:16, 56745.23it/s]
5%|▌ | 51249/1000000 [00:00<00:16, 57393.71it/s]
6%|▌ | 57117/1000000 [00:01<00:16, 57787.94it/s]
6%|▋ | 62897/1000000 [00:01<00:16, 56626.18it/s]
7%|▋ | 68770/1000000 [00:01<00:16, 57255.92it/s]
7%|▋ | 74501/1000000 [00:01<00:16, 56497.62it/s]
8%|▊ | 80420/1000000 [00:01<00:16, 57293.75it/s]
9%|▊ | 86155/1000000 [00:01<00:16, 57042.44it/s]
9%|▉ | 91970/1000000 [00:01<00:15, 57370.83it/s]
10%|▉ | 97710/1000000 [00:01<00:16, 55781.74it/s]
10%|█ | 103339/1000000 [00:01<00:16, 55927.11it/s]
11%|█ | 109260/1000000 [00:01<00:15, 56892.41it/s]
11%|█▏ | 114957/1000000 [00:02<00:15, 56020.54it/s]
12%|█▏ | 120727/1000000 [00:02<00:15, 56512.68it/s]
13%|█▎ | 126640/1000000 [00:02<00:15, 57285.70it/s]
13%|█▎ | 132374/1000000 [00:02<00:15, 57258.13it/s]
14%|█▍ | 138104/1000000 [00:02<00:15, 57264.19it/s]
14%|█▍ | 143909/1000000 [00:02<00:14, 57493.37it/s]
15%|█▍ | 149804/1000000 [00:02<00:14, 57926.61it/s]
16%|█▌ | 155599/1000000 [00:02<00:14, 56499.15it/s]
16%|█▌ | 161321/1000000 [00:02<00:14, 56708.55it/s]
17%|█▋ | 166999/1000000 [00:02<00:15, 54470.54it/s]
17%|█▋ | 172904/1000000 [00:03<00:14, 55791.98it/s]
18%|█▊ | 178671/1000000 [00:03<00:14, 56338.94it/s]
18%|█▊ | 184498/1000000 [00:03<00:14, 56905.42it/s]
19%|█▉ | 190201/1000000 [00:03<00:14, 56349.19it/s]
20%|█▉ | 196125/1000000 [00:03<00:14, 57199.44it/s]
20%|██ | 201853/1000000 [00:03<00:14, 56478.92it/s]
21%|██ | 207758/1000000 [00:03<00:13, 57235.37it/s]
21%|██▏ | 213488/1000000 [00:03<00:13, 56941.82it/s]
22%|██▏ | 219447/1000000 [00:03<00:13, 57724.86it/s]
23%|██▎ | 225396/1000000 [00:03<00:13, 58247.16it/s]
23%|██▎ | 231354/1000000 [00:04<00:13, 58643.19it/s]
24%|██▎ | 237315/1000000 [00:04<00:12, 58929.31it/s]
24%|██▍ | 243210/1000000 [00:04<00:13, 57989.38it/s]
25%|██▍ | 249014/1000000 [00:04<00:12, 58002.91it/s]
25%|██▌ | 254944/1000000 [00:04<00:12, 58385.64it/s]
26%|██▌ | 260786/1000000 [00:04<00:13, 53371.36it/s]
27%|██▋ | 266712/1000000 [00:04<00:13, 55018.18it/s]
27%|██▋ | 272671/1000000 [00:04<00:12, 56322.57it/s]
28%|██▊ | 278387/1000000 [00:04<00:12, 56564.14it/s]
28%|██▊ | 284084/1000000 [00:05<00:12, 56374.54it/s]
29%|██▉ | 289750/1000000 [00:05<00:12, 55270.77it/s]
30%|██▉ | 295648/1000000 [00:05<00:12, 56348.76it/s]
30%|███ | 301486/1000000 [00:05<00:12, 56943.77it/s]
31%|███ | 307433/1000000 [00:05<00:12, 57687.90it/s]
31%|███▏ | 313214/1000000 [00:05<00:12, 55891.35it/s]
32%|███▏ | 318869/1000000 [00:05<00:12, 56080.59it/s]
32%|███▏ | 324775/1000000 [00:05<00:11, 56954.55it/s]
33%|███▎ | 330482/1000000 [00:05<00:11, 56031.12it/s]
34%|███▎ | 336160/1000000 [00:05<00:11, 56247.96it/s]
34%|███▍ | 342085/1000000 [00:06<00:11, 57131.18it/s]
35%|███▍ | 347865/1000000 [00:06<00:11, 57327.48it/s]
35%|███▌ | 353714/1000000 [00:06<00:11, 57670.02it/s]
36%|███▌ | 359485/1000000 [00:06<00:11, 57491.97it/s]
37%|███▋ | 365427/1000000 [00:06<00:10, 58065.80it/s]
37%|███▋ | 371322/1000000 [00:06<00:10, 58328.43it/s]
38%|███▊ | 377157/1000000 [00:06<00:10, 58013.39it/s]
38%|███▊ | 382993/1000000 [00:06<00:10, 58115.73it/s]
39%|███▉ | 388886/1000000 [00:06<00:10, 58355.88it/s]
39%|███▉ | 394723/1000000 [00:06<00:10, 57126.42it/s]
40%|████ | 400534/1000000 [00:07<00:10, 57413.64it/s]
41%|████ | 406477/1000000 [00:07<00:10, 58009.60it/s]
41%|████ | 412368/1000000 [00:07<00:10, 58276.42it/s]
42%|████▏ | 418243/1000000 [00:07<00:09, 58417.05it/s]
42%|████▏ | 424087/1000000 [00:07<00:09, 58390.19it/s]
43%|████▎ | 429928/1000000 [00:07<00:09, 57091.16it/s]
44%|████▎ | 435645/1000000 [00:07<00:09, 56621.16it/s]
44%|████▍ | 441520/1000000 [00:07<00:09, 57244.34it/s]
45%|████▍ | 447450/1000000 [00:07<00:09, 57850.98it/s]
45%|████▌ | 453240/1000000 [00:07<00:09, 57856.13it/s]
46%|████▌ | 459153/1000000 [00:08<00:09, 58233.87it/s]
47%|████▋ | 465055/1000000 [00:08<00:09, 58464.95it/s]
47%|████▋ | 470904/1000000 [00:08<00:09, 56985.60it/s]
48%|████▊ | 476701/1000000 [00:08<00:09, 57272.75it/s]
48%|████▊ | 482436/1000000 [00:08<00:09, 55584.57it/s]
49%|████▉ | 488009/1000000 [00:08<00:09, 54767.42it/s]
49%|████▉ | 493830/1000000 [00:08<00:09, 55765.57it/s]
50%|████▉ | 499417/1000000 [00:08<00:09, 54895.85it/s]
51%|█████ | 505098/1000000 [00:08<00:08, 55450.68it/s]
51%|█████ | 510962/1000000 [00:08<00:08, 56385.94it/s]
52%|█████▏ | 516608/1000000 [00:09<00:08, 55837.36it/s]
52%|█████▏ | 522198/1000000 [00:09<00:08, 54052.94it/s]
53%|█████▎ | 527618/1000000 [00:09<00:08, 53794.69it/s]
53%|█████▎ | 533007/1000000 [00:09<00:08, 53492.61it/s]
54%|█████▍ | 538879/1000000 [00:09<00:08, 55020.71it/s]
54%|█████▍ | 544674/1000000 [00:09<00:08, 55880.24it/s]
55%|█████▌ | 550537/1000000 [00:09<00:07, 56693.16it/s]
56%|█████▌ | 556305/1000000 [00:09<00:07, 56985.62it/s]
56%|█████▌ | 562008/1000000 [00:09<00:07, 56908.04it/s]
57%|█████▋ | 567885/1000000 [00:10<00:07, 57462.94it/s]
57%|█████▋ | 573814/1000000 [00:10<00:07, 58006.14it/s]
58%|█████▊ | 579705/1000000 [00:10<00:07, 58274.60it/s]
59%|█████▊ | 585635/1000000 [00:10<00:07, 58579.44it/s]
59%|█████▉ | 591494/1000000 [00:10<00:07, 58123.07it/s]
60%|█████▉ | 597413/1000000 [00:10<00:06, 58439.07it/s]
60%|██████ | 603349/1000000 [00:10<00:06, 58711.52it/s]
61%|██████ | 609255/1000000 [00:10<00:06, 58812.05it/s]
62%|██████▏ | 615187/1000000 [00:10<00:06, 58961.83it/s]
62%|██████▏ | 621084/1000000 [00:10<00:06, 58171.15it/s]
63%|██████▎ | 626944/1000000 [00:11<00:06, 58297.42it/s]
63%|██████▎ | 632776/1000000 [00:11<00:06, 56906.39it/s]
64%|██████▍ | 638475/1000000 [00:11<00:06, 55842.87it/s]
64%|██████▍ | 644227/1000000 [00:11<00:06, 56330.55it/s]
65%|██████▍ | 649906/1000000 [00:11<00:06, 56463.14it/s]
66%|██████▌ | 655715/1000000 [00:11<00:06, 56940.81it/s]
66%|██████▌ | 661552/1000000 [00:11<00:05, 57363.68it/s]
67%|██████▋ | 667292/1000000 [00:11<00:05, 56950.43it/s]
67%|██████▋ | 673088/1000000 [00:11<00:05, 57247.96it/s]
68%|██████▊ | 678816/1000000 [00:11<00:05, 56942.77it/s]
68%|██████▊ | 684615/1000000 [00:12<00:05, 57249.90it/s]
69%|██████▉ | 690532/1000000 [00:12<00:05, 57819.49it/s]
70%|██████▉ | 696316/1000000 [00:12<00:05, 57673.56it/s]
70%|███████ | 702085/1000000 [00:12<00:05, 54700.34it/s]
71%|███████ | 707696/1000000 [00:12<00:05, 55103.11it/s]
71%|███████▏ | 713245/1000000 [00:12<00:05, 55212.86it/s]
72%|███████▏ | 719058/1000000 [00:12<00:05, 56067.38it/s]
72%|███████▏ | 724698/1000000 [00:12<00:04, 56162.04it/s]
73%|███████▎ | 730462/1000000 [00:12<00:04, 56597.23it/s]
74%|███████▎ | 736301/1000000 [00:12<00:04, 57129.84it/s]
74%|███████▍ | 742074/1000000 [00:13<00:04, 57307.98it/s]
75%|███████▍ | 747890/1000000 [00:13<00:04, 57560.62it/s]
75%|███████▌ | 753649/1000000 [00:13<00:04, 57347.21it/s]
76%|███████▌ | 759386/1000000 [00:13<00:04, 56080.24it/s]
77%|███████▋ | 765002/1000000 [00:13<00:04, 55609.85it/s]
77%|███████▋ | 770866/1000000 [00:13<00:04, 56497.30it/s]
78%|███████▊ | 776795/1000000 [00:13<00:03, 57322.64it/s]
78%|███████▊ | 782533/1000000 [00:13<00:03, 57284.67it/s]
79%|███████▉ | 788360/1000000 [00:13<00:03, 57574.70it/s]
79%|███████▉ | 794121/1000000 [00:13<00:03, 57427.25it/s]
80%|███████▉ | 799866/1000000 [00:14<00:03, 57050.59it/s]
81%|████████ | 805822/1000000 [00:14<00:03, 57794.13it/s]
81%|████████ | 811711/1000000 [00:14<00:03, 58118.19it/s]
82%|████████▏ | 817525/1000000 [00:14<00:03, 56878.05it/s]
82%|████████▏ | 823320/1000000 [00:14<00:03, 57192.43it/s]
83%|████████▎ | 829200/1000000 [00:14<00:02, 57666.37it/s]
83%|████████▎ | 834971/1000000 [00:14<00:02, 55826.84it/s]
84%|████████▍ | 840589/1000000 [00:14<00:02, 55928.46it/s]
85%|████████▍ | 846437/1000000 [00:14<00:02, 56677.34it/s]
85%|████████▌ | 852114/1000000 [00:14<00:02, 55108.31it/s]
86%|████████▌ | 857639/1000000 [00:15<00:02, 54827.63it/s]
86%|████████▋ | 863342/1000000 [00:15<00:02, 55469.23it/s]
87%|████████▋ | 869208/1000000 [00:15<00:02, 56407.64it/s]
87%|████████▋ | 874857/1000000 [00:15<00:02, 54882.50it/s]
88%|████████▊ | 880441/1000000 [00:15<00:02, 55160.51it/s]
89%|████████▊ | 886317/1000000 [00:15<00:02, 56216.58it/s]
89%|████████▉ | 891948/1000000 [00:15<00:02, 52580.46it/s]
90%|████████▉ | 897755/1000000 [00:15<00:01, 54132.39it/s]
90%|█████████ | 903415/1000000 [00:15<00:01, 54840.13it/s]
91%|█████████ | 908934/1000000 [00:16<00:01, 54119.11it/s]
91%|█████████▏| 914463/1000000 [00:16<00:01, 54456.61it/s]
92%|█████████▏| 920320/1000000 [00:16<00:01, 55661.54it/s]
93%|█████████▎| 925902/1000000 [00:16<00:01, 54763.98it/s]
93%|█████████▎| 931608/1000000 [00:16<00:01, 55433.28it/s]
94%|█████████▎| 937482/1000000 [00:16<00:01, 56408.79it/s]
94%|█████████▍| 943132/1000000 [00:16<00:01, 55137.68it/s]
95%|█████████▍| 948830/1000000 [00:16<00:00, 55674.78it/s]
95%|█████████▌| 954407/1000000 [00:16<00:00, 54948.14it/s]
96%|█████████▌| 960090/1000000 [00:16<00:00, 55496.75it/s]
97%|█████████▋| 965923/1000000 [00:17<00:00, 56332.50it/s]
97%|█████████▋| 971563/1000000 [00:17<00:00, 56200.35it/s]
98%|█████████▊| 977484/1000000 [00:17<00:00, 57092.18it/s]
98%|█████████▊| 983398/1000000 [00:17<00:00, 57699.08it/s]
99%|█████████▉| 989171/1000000 [00:17<00:00, 54916.54it/s]
99%|█████████▉| 994692/1000000 [00:17<00:00, 53980.70it/s]
100%|██████████| 1000000/1000000 [00:17<00:00, 56610.83it/s]
acceptance_rate = 0.076903
The final distribution of spins looks like this
plt.figure()
plt.imshow(lattice_spins, cmap="RdYlBu")
# show gridlines
for i in range(N):
plt.axhline(i + 0.5, color="black", lw=0.1)
plt.axvline(i + 0.5, color="black", lw=0.1)
plt.show()
Magnetization#
We can measure the magnetization, defined as the average spin over the lattice
\[M = \frac{1}{N}\sum_{i=1}^N \sigma_i,\]
where \(N\) is the number of spins.
Note we use a burn-in period of 100,000 steps to allow the system to reach equilibrium.
BURNIN = 100_000
m_mean = np.mean(m_values[BURNIN:])
m_std = np.std(m_values[BURNIN:])
plt.figure()
plt.plot(range(STEPS)[:BURNIN], m_values[:BURNIN], label="Burn in")
plt.plot(range(STEPS)[BURNIN:], m_values[BURNIN:], label="Sampling")
# show mean as dashed line
plt.plot(range(STEPS), m_mean * np.ones((STEPS)), "--", color="black")
# show variance as filled box
plt.fill_between(range(STEPS), m_mean - m_std, m_mean + m_std, color="gray", alpha=0.3)
plt.xlabel("Steps")
plt.ylabel("Magnetization")
plt.legend(title=f"T={T}")
plt.show()
print(f"magnetization mean = {m_mean}")
print(f"magnetization std = {m_std}")
magnetization mean = 0.9139394173333335
magnetization std = 0.007050481607383647
np.random.random()
0.24046307940991352