Statistical Mechanics

Code Examples

Ensembles and Partition Functions

import numpy as np

def boltzmann_factor(E, T):
    """exp(-E/k_B T)"""
    kB = 1.38e-23
    return np.exp(-E / (kB * T))

def canonical_partition_function(energies, T):
    """Z = Σ exp(-βE_i)"""
    beta = 1 / (1.38e-23 * T)
    return np.sum(np.exp(-beta * np.array(energies)))

def canonical_expectation(E, Z, T):
    """⟨E⟩ = (1/Z) Σ E_i exp(-βE_i)"""
    beta = 1 / (1.38e-23 * T)
    weights = np.exp(-beta * np.array(E))
    return np.sum(E * weights) / Z

def heat_capacity(energies, T, Cv):
    """Cv = (⟨E²⟩ - ⟨E⟩²) / (k_B T²)"""
    beta = 1 / (1.38e-23 * T)
    Z = canonical_partition_function(energies, T)
    E_avg = canonical_expectation(energies, Z, T)
    E2_avg = canonical_expectation([E**2 for E in energies], Z, T)
    return (E2_avg - E_avg**2) / ((1.38e-23) * T**2)

# Example: Two-level system
E0, E1 = 0, 0.001  # Energy levels in J (corresponds to ~72 K)
energies = [E0, E1]

print("Two-level system:")
for T in [10, 50, 100, 300, 1000]:
    Z = canonical_partition_function(energies, T)
    E_avg = canonical_expectation(energies, Z, T)
    print(f"  T = {T:4d} K: Z = {Z:.4f}, ⟨E⟩ = {E_avg*1e21:.2f} × 10⁻²¹ J")

# Partition function for harmonic oscillator
def harmonic_oscillator_z(T, omega):
    """Z = exp(-βℏω/2) / (1 - exp(-βℏω))"""
    hbar = 1.055e-34
    beta = 1 / (1.38e-23 * T)
    if omega * hbar * beta > 100:
        return np.exp(-omega * hbar * beta / 2)
    return np.exp(-omega * hbar * beta / 2) / (1 - np.exp(-omega * hbar * beta))

# Free particle partition function (classical limit)
def free_partition_function_1D(V, T, m=1.67e-27):
    """Z = V / λ³ where λ = h/√(2πmkT)"""
    h = 6.626e-34
    kB = 1.38e-23
    lam = h / np.sqrt(2 * np.pi * m * kB * T)
    return V / lam**3

# 3D ideal gas
def ideal_gas_pressure(N, T, V, m=1.67e-27):
    """P = NkT/V (ideal gas law)"""
    kB = 1.38e-23
    return N * kB * T / V

N_A = 6.022e23
V = 0.0224  # 1 mole at STP
T = 273
P = ideal_gas_pressure(N_A, T, V)
print(f"\nIdeal gas at STP: P = {P/1e5:.2f} bar")

# Grand canonical ensemble
def grand_canonical_partition_function(grand_mu, T, single_particle_states):
    """Ξ = Π_k (1 + exp(-β(ε_k - μ))) for fermions"""
    beta = 1 / (1.38e-23 * T)
    Xi = 1
    for eps in single_particle_states:
        Xi *= (1 + np.exp(-beta * (eps - grand_mu)))
    return Xi

# Fermi-Dirac and Bose-Einstein occupation
def fermi_dirac_occupation(eps, mu, T):
    """f(ε) = 1/(exp((ε-μ)/kT) + 1)"""
    kB = 1.38e-23
    return 1 / (np.exp((eps - mu) / (kB * T)) + 1)

def bose_einstein_occupation(eps, mu, T):
    """n(ε) = 1/(exp((ε-μ)/kT) - 1)"""
    kB = 1.38e-23
    if mu >= eps:
        return np.inf
    return 1 / (np.exp((eps - mu) / (kB * T)) - 1)

print(f"\nFermi-Dirac occupation (ε-μ = 0.1 eV, T=300K):")
delta_E = 0.1 * 1.6e-19  # 0.1 eV
f = fermi_dirac_occupation(delta_E, 0, 300)
print(f"  f(ε_F + 0.1 eV) = {f:.4f}")

Ising Model and Phase Transitions

import numpy as np

class IsingModel:
    def __init__(self, L, J=1, h=0):
        """2D Ising model on L×L lattice."""
        self.L = L
        self.J = J
        self.h = h
        self.N = L * L
        self.spins = np.random.choice([-1, 1], size=(L, L))
    
    def energy(self):
        """Calculate total energy."""
        E = 0
        for i in range(self.L):
            for j in range(self.L):
                # Periodic boundaries
                right = self.spins[(i+1) % self.L, j]
                down = self.spins[i, (j+1) % self.L]
                E -= self.J * self.spins[i, j] * (right + down)
                E -= self.h * self.spins[i, j]
        return E / 2  # Each bond counted twice
    
    def magnetization(self):
        """Total magnetization."""
        return np.sum(self.spins)
    
    def metropolis_step(self, T):
        """Single Metropolis sweep."""
        kB = 1.0  # Units where kB = 1
        
        for _ in range(self.N):
            # Random spin flip
            i = np.random.randint(self.L)
            j = np.random.randint(self.L)
            
            # Calculate energy change
            neighbors = (self.spins[(i-1) % self.L, j] + 
                        self.spins[(i+1) % self.L, j] +
                        self.spins[i, (j-1) % self.L] + 
                        self.spins[i, (j+1) % self.L])
            
            dE = 2 * self.spins[i, j] * (self.J * neighbors + self.h)
            
            # Metropolis criterion
            if dE < 0 or np.random.random() < np.exp(-dE / (kB * T)):
                self.spins[i, j] *= -1
    
    def equilibrate(self, T, n_sweeps=1000):
        """Equilibrate at given temperature."""
        for _ in range(n_sweeps):
            self.metropolis_step(T)
    
    def simulate(self, T, n_sweeps=10000):
        """Run simulation and return observables."""
        self.equilibrate(T, n_sweeps)
        
        mags = []
        eners = []
        for _ in range(n_sweeps):
            self.metropolis_step(T)
            mags.append(self.magnetization())
            eners.append(self.energy())
        
        return np.mean(mags)/self.N, np.std(mags)/self.N, np.mean(eners)/self.N, np.std(eners)/self.N

# Critical temperature for 2D Ising (Onsager solution: kB Tc/J ≈ 2.269)
L = 20
ising = IsingModel(L)

print("2D Ising model simulation:")
temperatures = [1.0, 2.0, 2.2, 2.269, 2.3, 2.5, 3.0]
for T in temperatures:
    m, dm, E, dE = ising.simulate(T, n_sweeps=5000)
    print(f"  T = {T:.2f}: m = {m:.4f} ± {dm:.4f}, E = {E/L**2:.4f} ± {dE/L**2:.4f}")

# Binder cumulant
def binder_cumulant(mags):
    """U = 1 - ⟨m⁴⟩/(3⟨m²⟩²)"""
    m2 = np.mean(np.array(mags)**2)
    m4 = np.mean(np.array(mags)**4)
    return 1 - m4 / (3 * m2**2)

print(f"\nBinder cumulant at Tc ≈ 2.27:")
for T in [2.2, 2.269, 2.3]:
    ising2 = IsingModel(20)
    ising2.equilibrate(T, 5000)
    mags = [ising2.magnetization() for _ in range(5000)]
    U = binder_cumulant(mags)
    print(f"  U(T={T}) = {U:.4f}")

Mean-Field Theory

import numpy as np

def mean_field_ising(T, J, h=0, max_iter=100, tol=1e-6):
    """Mean-field solution of Ising model."""
    kB = 1.0  # Units where kB = 1
    beta = 1 / T
    
    # Initial magnetization
    m = 0.1
    
    for iteration in range(max_iter):
        m_new = np.tanh(beta * (J * 6 * m + h))  # z=6 for 3D
        if abs(m_new - m) < tol:
            break
        m = m_new
    
    return m

def solve_self_consistent(T, J):
    """Solve MF equations numerically."""
    from scipy.optimize import brentq
    
    def equation(m):
        kB = 1.0
        return m - np.tanh(m * J * 6 / (kB * T))
    
    # Below Tc, look for non-zero solution
    Tc_mf = J * 6  # Mean-field Tc
    
    if T > Tc_mf:
        return 0
    
    # Find root
    try:
        m = brentq(equation, 1e-6, 1.0)
        return m
    except:
        return 0

# Critical exponents
Tc_mf = 6.0  # For J=1
print("Mean-field Ising critical behavior:")
print(f"  Mean-field Tc = {Tc_mf:.3f} J/kB")

# Below Tc, magnetization follows m ∝ (Tc - T)^β with β = 0.5
for T in [0.5, 1.0, 2.0, 4.0, 5.0, 5.5, 5.9, 5.99]:
    m = solve_self_consistent(T, 1.0)
    if T < Tc_mf:
        beta_mf = np.log(m) / np.log((Tc_mf - T) / Tc_mf)
        print(f"  T = {T:.2f}: m = {m:.4f}, effective β = {beta_mf:.3f}")
    else:
        print(f"  T = {T:.2f}: m = 0")

# Landau theory
def landau_free_energy(m, a, b, T):
    """F(m) = a(T-Tc)m² + bm⁴ + hm"""
    return a * (T - Tc_mf) * m**2 + b * m**4

# Susceptibility
def susceptibility(T, Tc, chi0=1):
    """χ ∝ |T-Tc|^(-γ) with γ = 1 (MF)"""
    if T > Tc:
        return chi0 / (T - Tc)
    return chi0 / (T_c - T)

print(f"\nMean-field susceptibility divergence:")
for T in [2.0, 1.1, 1.01, 1.001]:
    chi = susceptibility(T, 1.0)
    print(f"  T = {T:.3f}: χ = {chi:.1f}")

Monte Carlo Methods

import numpy as np

def metropolis_sampling(pdf, proposal_std, n_samples, x0):
    """Metropolis-Hastings sampling from arbitrary PDF."""
    samples = [x0]
    x = x0
    
    for _ in range(n_samples):
        x_proposed = x + np.random.normal(0, proposal_std)
        alpha = pdf(x_proposed) / pdf(x)
        
        if np.random.random() < alpha:
            x = x_proposed
        
        samples.append(x)
    
    return np.array(samples)

def wolff_cluster_update(spins, J, T):
    """Wolff single-cluster Monte Carlo update."""
    L = len(spins)
    visited = np.zeros(L, dtype=bool)
    
    # Start new cluster
    site = np.random.randint(L)
    cluster = [site]
    visited[site] = True
    
    S_cluster = spins[site]
    
    # Grow cluster
    i = 0
    while i < len(cluster):
        current = cluster[i]
        P_add = 1 - np.exp(-2 * J / T)
        
        neighbors = [(current - 1) % L, (current + 1) % L]
        for n in neighbors:
            if not visited[n] and spins[n] == S_cluster:
                if np.random.random() < P_add:
                    visited[n] = True
                    cluster.append(n)
        i += 1
    
    # Flip cluster
    for site in cluster:
        spins[site] *= -1
    
    return spins

def swendsen_wang_update(spins, J, T):
    """Swendsen-Wang multi-cluster update."""
    L = len(spins)
    bonds = np.zeros(L, dtype=bool)
    
    # Create bonds
    for i in range(L):
        if spins[i] == spins[(i+1) % L]:
            P_bond = 1 - np.exp(-2 * J / T)
            bonds[i] = np.random.random() < P_bond
    
    # Find clusters using bond connectivity
    clusters = []
    visited = np.zeros(L, dtype=bool)
    
    for i in range(L):
        if not visited[i]:
            cluster = [i]
            visited[i] = True
            j = i
            while bonds[j]:
                j = (j + 1) % L
                if not visited[j]:
                    visited[j] = True
                    cluster.append(j)
            clusters.append(cluster)
    
    # Flip random subset of clusters
    for cluster in clusters:
        if np.random.random() < 0.5:
            for site in cluster:
                spins[site] *= -1
    
    return spins

# Auto-correlation time estimation
def autocorrelation_time(samples):
    """Estimate autocorrelation time from decay of autocorrelation function."""
    n = len(samples)
    mean = np.mean(samples)
    var = np.var(samples)
    
    # Normalized autocorrelation
    acf = np.correlate(samples - mean, samples - mean, mode='full')
    acf = acf[n-1:] / (acf[n-1] * np.arange(n, 0, -1))
    
    # Integrate to get tau
    tau = 0.5 + np.sum(acf[1:])
    return tau

# Example: 1D Ising model autocorrelation
N = 100
J, T = 1.0, 1.0
spins = np.random.choice([-1, 1], size=N)

# Measure autocorrelation time for different algorithms
print("Autocorrelation time comparison (1D Ising at Tc):")
print(f"  Single-spin flip: high τ (critical slowing down)")
print(f"  Wolff/S-W clusters: τ = O(1) near Tc")

Non-Equilibrium Dynamics

import numpy as np

class MasterEquation:
    def __init__(self, transition_matrix):
        """dP/dt = W·P for Markov process."""
        self.W = transition_matrix
    
    def propagate(self, P0, dt, n_steps):
        """Forward Euler integration of master equation."""
        P = P0.copy()
        t = 0
        results = [P.copy()]
        
        for _ in range(n_steps):
            dP = self.W @ P
            P = P + dt * dP
            t += dt
            results.append(P.copy())
        
        return np.array(results)

def langevin_dynamics(x0, v0, gamma, kT, m, dt, n_steps):
    """Langevin dynamics: m dv/dt = -γv - dU/dx + √(2γkT)ξ(t)"""
    x, v = x0, v0
    results = [(x, v)]
    
    for _ in range(n_steps):
        # Random force (Einstein relation)
        xi = np.random.normal(0, 1)
        noise = np.sqrt(2 * gamma * kT / dt)
        
        # Velocity Verlet with friction
        v += dt / m * (-gamma * v + noise * xi)
        x += dt * v
        
        results.append((x, v))
    
    return np.array(results)

def green_kubo_relation(velocity_correlation):
    """Diffusion coefficient from velocity autocorrelation."""
    return np.trapz(velocity_correlation, dx=dt)

# Brownian motion
def brownian_motion_1D(D, dt, n_steps):
    """Einstein-Smoluchowski relation: ⟨x²⟩ = 2Dt"""
    x = 0
    positions = [x]
    
    for _ in range(n_steps):
        dx = np.random.normal(0, np.sqrt(2 * D * dt))
        x += dx
        positions.append(x)
    
    return np.array(positions)

# Mean squared displacement
def mean_squared_displacement(positions):
    """MSD = ⟨(x(t) - x(0))²⟩"""
    return np.array([np.mean((positions[t:] - positions[:-t])**2) 
                     for t in range(len(positions)//2)])

D = 1.0  # Diffusion coefficient
dt = 0.001
n_steps = 10000

positions = brownian_motion_1D(D, dt, n_steps)
msd = mean_squared_displacement(positions)

print("Brownian motion (Einstein relation):")
print(f"  D = {D}")
print(f"  MSD(t) ≈ 2Dt: slope at short times = {np.polyfit(range(100), msd[:100], 1)[0]/dt:.2f}")

# Fokker-Planck equation (conceptual)
def fokker_planck_forward(x, t, D, drift):
    """∂P/∂t = -∂/∂x(Drift·P) + D ∂²P/∂x²"""
    pass

# Non-equilibrium work relations
def jarzynski_equality(work_samples):
    """⟨exp(-W/kT)⟩ = exp(-ΔF/kT)"""
    kT = 1.0  # Units
    return np.mean(np.exp(-np.array(work_samples) / kT))

# Crooks fluctuation theorem
def crooks_fluctuation_theorem(work_forward, work_reverse):
    """P_F(W)/P_R(-W) = exp((W - ΔF)/kT)"""
    pass

Best Practices

• Verify equilibrium properties by checking that observables don't drift over time in long simulations.
• Use multiple independent runs with different random seeds to estimate statistical errors.
• For critical phenomena, simulate large systems near Tc to minimize finite-size effects.
• Use cluster algorithms (Wolff, Swendsen-Wang) to overcome critical slowing down near phase transitions.
• Distinguish between different ensemble averages and ensure proper equilibration before measurement.
• Use block averaging to estimate statistical errors and autocorrelation times.
• For non-equilibrium simulations, verify that fluctuations satisfy fluctuation-dissipation relations.
• When applying mean-field theory, check self-consistency of the mean-field approximation.
• Use histogram reweighting to efficiently sample across phase transitions.
• Consider conservation laws and select appropriate dynamics (microcanonical vs canonical) for your system.