Research entry
The Ising Model — Metropolis Monte Carlo
2018 · Academic Archive
Monte Carlo simulation of the 2D and 3D Ising model using the Metropolis algorithm — computing energy, magnetisation, heat capacity, and susceptibility across a temperature sweep to find the ferromagnetic phase transition and Curie temperature.
Open 3D Visualiser → · GitHub → · Full Research Writeup →
Interactive Simulation
Hamiltonian: H = −J Σ σᵢσⱼ − h Σ σᵢ
Cv = β²(⟨E²⟩ − ⟨E⟩²)N² · χ = β(⟨M²⟩ − ⟨M⟩²)N²
Square lattice: Tc ≈ 2.27 J/k_B (Onsager exact) · Triangular: Tc ≈ 3.64 J/k_B
Overview
A computational physics assignment from JS TP (Junior Sophister Theoretical Physics) at Trinity College Dublin, 2018. The project implemented the Metropolis algorithm to simulate the Ising model across a range of geometries (2D square, 2D triangular, 3D cubic) and investigate the ferromagnetic phase transition.
The Ising model assigns a spin σᵢ ∈ {−1, +1} to each site of a lattice. The energy of a configuration is given by the Hamiltonian:
H = −J Σ_{⟨i,j⟩} σᵢσⱼ − h Σᵢ σᵢ- J > 0: ferromagnetic — neighbouring spins prefer to align
- J < 0: antiferromagnetic — neighbouring spins prefer to anti-align
- h: external magnetic field
The Metropolis Algorithm
- Initialise a random lattice of spins
s₀ = {σ₁, ..., σₙ} - For each spin σᵢ, compute the energy change from flipping it:
ΔE = 2σᵢ ( J Σⱼ σⱼ + h ) - Accept the flip if ΔE < 0. Otherwise accept with probability e^{−βΔE}
- Repeat until the lattice converges (
|⟨M⟩| → 1for a ferromagnet) - After warm-up, collect averages of E, M, E², M² over many iterations
Observables
All four observables are computed from a single run (no re-running needed):
| Observable | Formula |
|---|---|
| Average energy per spin | ⟨E⟩ = (1/N²) Σ (−J σᵢ Σⱼ σⱼ − h σᵢ) |
| Average magnetisation per spin | ⟨M⟩ = (1/N²) Σ σᵢ |
| Specific heat capacity | Cv = β²(⟨E²⟩ − ⟨E⟩²) N² |
| Magnetic susceptibility | χ = β(⟨M²⟩ − ⟨M⟩²) N² |
Key Results
Phase Transition and Curie Temperature
Below the Curie temperature Tc, the lattice is ferromagnetic (all spins aligned, |⟨M⟩| ≈ 1). Above Tc, thermal fluctuations disorder the lattice (paramagnetic). At Tc itself, both Cv and χ show a sharp spike — the signature of a second-order phase transition.
| Geometry | Neighbours | Ground state energy | Tc (measured) | Tc (theory) |
|---|---|---|---|---|
| 2D Square | 4 | −2J | 2.7 ± 0.35 J/k_B | 2.269 J/k_B (Onsager) |
| 2D Triangular | 6 | −3J | 3.7 ± 0.35 J/k_B | 3.641 J/k_B |
| 3D Cubic | 6 | −3J | 4.1 ± 0.35 J/k_B | ~4.511 J/k_B |
The ground state energies matched theory exactly. The triangular and 3D cubic lattices share 6 neighbours each, giving the same ground state energy but slightly different Curie temperatures.
Boundary Conditions
Replacing periodic boundary conditions with hard walls on a 10×10 lattice reduced Tc from 2.6 to 2.2 J/k_B. Edge spins only have 2–3 neighbours, weakening the cooperative effect. The relative impact scales as (4N−4) / (N−2)² for boundary vs. interior spins.
External Magnetic Field
With h ≠ 0, the average energy shifts to −2|J| − |h|. A larger field raises the effective Curie temperature — more thermal energy is needed to break the combined field + exchange alignment. Magnetic susceptibility decreases with larger h (the lattice is already saturated).
Hysteresis Loop
Sweeping h from positive to negative and back produces a hysteresis loop — the lattice retains its magnetisation longer than naively expected. Larger J widens the loop. Negative J (antiferromagnetic) collapses the loop to near zero unless |h| is large enough to force full alignment.
Coupling Constant J
The relationship between J and Tc is approximately linear — stronger spin-spin interactions require higher temperatures to destroy order. Ground state energies follow −2|J| regardless of the sign of J (antiferromagnets anti-align perfectly, cancelling the sign).
Original Python Code (2018)
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d
class IsingSolver:
def __init__(self, N, J=1.0, h=0.0, geometry='square'):
self.N = N
self.J = J
self.h = h
self.geometry = geometry
self.grid = np.random.choice([-1, 1], size=(N, N))
def get_neighbours(self, i, j):
N = self.N
neighbours = [
((i-1) % N, j), ((i+1) % N, j),
(i, (j-1) % N), (i, (j+1) % N)
]
if self.geometry == 'triangular':
if i % 2 == 0:
neighbours += [((i+1) % N, (j-1) % N), ((i-1) % N, (j-1) % N)]
else:
neighbours += [((i+1) % N, (j+1) % N), ((i-1) % N, (j+1) % N)]
return neighbours
def metropolis_step(self, T):
beta = 1.0 / T
for _ in range(self.N**2):
i, j = np.random.randint(0, self.N, size=2)
spin = self.grid[i, j]
neighbour_sum = sum(self.grid[ni, nj] for ni, nj in self.get_neighbours(i, j))
dE = 2 * spin * (self.J * neighbour_sum + self.h)
if dE < 0 or np.random.rand() < np.exp(-beta * dE):
self.grid[i, j] *= -1
def avg_magnetisation(self):
return np.mean(self.grid)
def run(self, T_range, warmup=100, samples=50):
results = {'T': [], 'E': [], 'M': [], 'Cv': [], 'chi': []}
for T in T_range:
for _ in range(warmup):
self.metropolis_step(T)
E_list, M_list = [], []
for _ in range(samples):
self.metropolis_step(T)
E = self._total_energy() / self.N**2
M = abs(self.avg_magnetisation())
E_list.append(E); M_list.append(M)
avgE = np.mean(E_list); avgM = np.mean(M_list)
Cv = (1/T**2) * (np.mean(np.array(E_list)**2) - avgE**2) * self.N**2
chi = (1/T) * (np.mean(np.array(M_list)**2) - avgM**2) * self.N**2
results['T'].append(T); results['E'].append(avgE)
results['M'].append(avgM); results['Cv'].append(Cv)
results['chi'].append(chi)
return results
def _total_energy(self):
E = 0
for i in range(self.N):
for j in range(self.N):
spin = self.grid[i, j]
neighbour_sum = sum(self.grid[ni, nj] for ni, nj in self.get_neighbours(i, j))
E += -self.J * spin * neighbour_sum / 2 - self.h * spin
return E
def find_curie_temperature(self, results):
M_smooth = gaussian_filter1d(results['M'], sigma=2)
dM = np.gradient(M_smooth, results['T'])
return results['T'][np.argmin(dM)]
# Example usage
solver = IsingSolver(N=10, J=1.0)
T_range = np.arange(0.5, 5.0, 0.1)
results = solver.run(T_range)
Tc = solver.find_curie_temperature(results)
print(f"Curie temperature: Tc ≈ {Tc:.2f} J/k_B")
fig, axes = plt.subplots(2, 2, figsize=(10, 8))
axes[0,0].plot(results['T'], results['E'])
axes[0,0].set_title('⟨E⟩ vs T'); axes[0,0].set_xlabel('T'); axes[0,0].set_ylabel('⟨E⟩')
axes[0,1].plot(results['T'], results['M'], color='green')
axes[0,1].axvline(Tc, color='red', linestyle='--', label=f'Tc={Tc:.2f}')
axes[0,1].set_title('|⟨M⟩| vs T'); axes[0,1].legend()
axes[1,0].plot(results['T'], results['Cv'], color='orange')
axes[1,0].set_title('Cv vs T')
axes[1,1].plot(results['T'], results['chi'], color='purple')
axes[1,1].set_title('χ vs T')
plt.tight_layout()
plt.show()What the Interactive Simulation Shows
Use the controls above to explore the physics:
- Run “Live Simulation” and watch the lattice spontaneously order at low T (all blue or all red), and disorder at high T (mixed)
- Run “Temp Sweep” to generate all four observable plots and read off Tc from the |⟨M⟩| inflection point and the Cv/χ spikes
- Switch to Triangular geometry to see the higher Tc (~3.7 vs ~2.7) from 6 neighbours vs 4
- Set J < 0 for antiferromagnetic behaviour — the lattice forms a checkerboard pattern
- Set h ≠ 0 to see the external field bias the magnetisation and shift the effective Tc