Ising models and Metropolis algorithm

Ising models and Metropolis algorithm#

What you will learn

Ising models are minimal models used to phase transitions.
The Ising model is defined by a Hamiltonian that encodes both the interaction between neighboring spins $i$ and $j$ (for nearest neighbors, $|i - j| = 1$), and the interaction of spins with an external field $B$:
Markov Chain Monte Carlo (MCMC) presents an straightforward and affordable way of simulating ising models.
The particularly simple and efficient formulation of MCMC for ising models is given by Metropolis selection criteria.

Ising model#

Ising Model

\[{H([s]) = - \sum_{<ij>} J s_i s_j - h \sum_i s_i}\]

\[{Z = \sum_{[s]} e^{-H([s])/k_B T}}\]

Each spin variable takes values $s_i = \pm 1$
$[s] = [s_1, s_2, \dots]$ denotes a microstate defined by the configuration of all spins
Parameters of the Ising model: $h$ is the external magnetic field, and $J$ is the spin–spin coupling strength
The summation $<ij>$ is over nearest-neighbor spin pairs

Microstates of Ising model#

import matplotlib.pyplot as plt
import numpy as np

spins = np.random.choice([-1,1],size=(8,8)) # every time you sammple a new microstate

print(spins)

plt.imshow(spins)

[[ 1 -1  1 -1 -1  1 -1  1]
 [-1  1  1 -1 -1  1 -1 -1]
 [ 1 -1 -1  1 -1 -1 -1  1]
 [ 1  1  1  1 -1 -1 -1 -1]
 [ 1 -1  1  1  1 -1 -1  1]
 [-1 -1 -1  1  1  1 -1 -1]
 [ 1 -1  1  1  1  1 -1 -1]
 [ 1  1  1  1  1 -1  1  1]]

<matplotlib.image.AxesImage at 0x7efc1c6fbbb0>

../_images/68f2b86e6b59847dc1e883cafc219c7d9e2563b74a399d78d3f8de50b3d6c791.png

Sampling microstates of Ising Model#

Monte Carlo appraoch to computing observables is simple: all we have to do is sample from the exponential distribution.

\[P(E_i)= \frac{e^{-\beta E_i}}{\sum_i e^{-\beta E_i}}\]

\[\langle A\rangle = \sum_i A_iP(E_i) \approx\frac{1}{N}\sum_{i=1}^N \underbrace{A_i}_{{\rm chosen}\ {\rm acc.}\ {\rm to}\ P(E_i)}\]

Problems with regular MC
- There is an astronomical number of states, $i = 1... N^2$!
- Most of these microstates make exponentially small contribution to probability distribution! $P(E_i)$
- Brute force MC will not be efficient for sampling important portions of $P(E_i)$.

Extracting thermodynamics from Ising models#

Mangetization#

Magnetization $M$ or magnetization per spin measures the net alignment of spins.

\[M([s]) = \sum_i s_i\]

\[m([s]) = \frac{M}{N}\]

Free Energy#

We can calculate probability of macrostate defined to include all spin configurations with a fixed value of magnetization $M$.

\[P(M) = \frac{Z(M)}{Z} = \frac{\sum_{[s]} e^{-\beta H([s])} \, \delta(M - M([s]))}{Z}\]

From this, we can compute the free energy as a function of magnetization $M$

\[F(M) = -\beta^{-1} \log Z(M) + const\]

This magnetization-dependent free energy tells us about the relative stability and likelihood of different magnetization states in the system.
If there are multiple coexisting phases we will expect to see multiple minima!

Response functions#

To quantify fluctuations in energy and magnetization, we compute the following response functions:

Heat capacity, $C_V$, measures energy fluctuations:

\[C_V(T) = \left( \frac{\partial \langle E \rangle}{\partial T} \right)_V = \frac{1}{k_B T^2} \left( \langle E^2 \rangle - \langle E \rangle^2 \right)\]

Magnetic susceptibility, $\chi_T$, measures magnetization fluctuations in response to an external field $B$:

\[\chi_T = \left( \frac{\partial \langle M \rangle}{\partial B} \right)_T = \frac{1}{k_B T} \left( \langle M^2 \rangle - \langle M \rangle^2 \right)\]

These response functions become particularly important near critical points, where fluctuations are enhanced, and $C_V$ and $\chi_T$ may diverge in the thermodynamic limit.

Enforcing periodic boundary conditions#

In the Ising model with periodic boundary conditions, the modulo operation is used to wrap around the edges of the lattice—so that spins on one edge interact with those on the opposite edge.
This effectively turns a 2D grid into a torus (doughnut shape), avoiding edge effects and mimicking an infinite system.
Consider a 2D lattice of size $ L \times L $. For a site at position $ (i, j) $, its right neighbor is at $ (i, j+1) $. But if $ j = L - 1 $, then $ j+1 $ exceeds the lattice. To wrap around, we use:

\[ j_{\text{right}} = (j + 1) \mod L \]

This brings the index back to $ 0 $ when it exceeds $ L - 1 $.

L = 10  # Lattice size
spins = np.random.choice([-1,1], size=(L, L))

i, j = 3, 9  # Current spin

# Neighbors with periodic boundary conditions
right = spins[i][(j + 1) % L]
left  = spins[i][(j - 1) % L]
up    = spins[(i - 1) % L][j]
down  = spins[(i + 1) % L][j]

../_images/d64189c003eb37b4690f06fba214b637d97637163ef8b16ddcf7acf512fbe84f.png

def get_dE(spins, i, j, J=1, B=0):
    '''Compute change in energy of 2D spin lattice 
    after flipping a spin at a location (i,j)''' 
    
    N = len(spins)
    
    z  = spins[(i+1)%N, j] + spins[(i-1)%N, j] + \
             spins[i, (j+1)%N] + spins[i, (j-1)%N] 
             
    dE = 2*spins[i,j]*(J*z + B)
    
    return dE

Generating chains! Master Equation and Detailed Balance#

Random chain

In simple MC each sampled state $X_i$ at step $i$ is independet of the next or previous steps.
Thus, MC sampling generates totally uncorrelated samples which is good for rapid convergence according to Central Limit Theorem. But as we remarked above the samples are most likely not covering important areas!

\[P(X_N, ..., X_2, X_1) = P_N(X_N)\cdot ... P_2(X_2) \cdot P_1(X_1)\]

Markov chain:

In simple MCMC each sampled state $X_i$ at steo $i$ is generated from $i-1$.
This introduces correlations between samples which means slower covnergence to the mean.
On the other hand MCMC find the important areas for sampling much faster so in the hand it wins big compared to MC.

\[P(X_N, ..., X_2, X_1) = P_N(X_N)\cdot ... P_3(X_3| X_2) \cdot P_2(X_2| X_1)... P_1(X_1)\]

The probability of picking state $i+1$ given taht we started from $i$ is referred to as transition probability $T_{ij}$

\[T_{21} = T(X_2 \leftarrow X_1) = P(X_2|X_1)\]

As a probability it should be normalized:

\[\sum_{X} T_{XX'} = 1\]

Probability of being at $X$ at $t+\Delta t$ given a prior state at $(X',t)$ can be written in terms of $T_{XX'}$ as:

\[P(X, t+\Delta t) = \sum_{X'} T_{XX'} P(X', t)\]

Subtracting P(X,t) from both states we can obtain equation of motion for Markov chain.

\[P(X, t+\Delta t) -P(X,t)= \sum_{X'} T_{XX'} P(X', t) - P(X, t)\sum_{X'} T_{X'X}\]

Defining transition rates as the limit of $w_{X'X} = lim_{\Delta t \rightarrow 0}\frac{T_{X'X} (\Delta t)}{\Delta t}$ we arrive at the master equation”

Master equation: A continuity equation in probability space.

\[\frac{\partial P(X,t)}{\partial t} = \sum_{X'} w_{X X'} P(X', t) - \sum_{X'} w_{X' X} P(X, t)\]

From maset equation we see that in order to reach equilibrium probability must stop depending on time. There are two possibilities:

\[\frac{\partial P(X,t)}{\partial t} = \sum_{X'} P(X', t) w_{XX'} - \sum_{X'} P(X, t) w_{X' X}=0\]

When the sum is equl zero becasue influx and outflux of probabiity currents cancel each out. This corresponds to steady state condition and is not true equilbroum for there are macroscopic fluxes in the system.
When every single pairs is identically zero, this is known as detailed balance condition which is equivalent to thermodynamic equilibrium.

Detailed Balance $\equiv$ Condition for equilibrium

\[{ w_{ij}\cdot p_j = w_{ji} \cdot p_i }\]

Master equation for a two state dynamics

\[\frac{dP(x_1,t)}{dt} = w_{12} p_2 - w_{21} p_1\]

\[\frac{dP(x_2,t)}{dt} = w_{21} p_1 - w_{12} p_2\]

In the equilibrium we see that equilibrium is established when probability ratio matches the ratio of tranistions $\frac{p_1}{p_2} = \frac{w_{12}}{w_{21}}$
In the $NVT$ ensemble we have an explicit requirment for transition rates $\frac{p_1}{p_2} = e^{-\beta (E_1-E_2)} = \frac{w_{12}}{w_{21}}$

How to pick the moves for Markov chain?#

Now it is time to consider practical aspects of conducting MCMC. How to pick states and transitions?

\[[transition] = [accept/reject] \cdot [move]\]

\[T(X' \leftarrow X) = \alpha_{X' X}\cdot A_{X' X} \]

We have great freedom in picking states and moves as long as we satisfy the detailed balance! E.g as long as the ratio of transition probabilities matches ration of Botlzman factors!

The simplest case for move is to pick one spins at random per iteration:

\[\alpha_{X' X} = \frac{1}{L^2}\]

For transitions we adopt criteria that favors our chain to explore low energy (high probability) configurations:
If $p(X') < p(X),\,\,\,\,$ $A_{X'X}=\frac{p(X')}{p(X)}$
If $p(X') \geq p(X),\,\,\,$ $A_{X'X}=1$

\[\frac{T_{X'X}}{T_{XX'}} = \frac{1/L^2}{1/L^2} \cdot \frac{p(X')/p(X)}{1}\]

Metropolis algorithm

${ i. \ Initialization.}$
Generate an initial spin configuration $[s_0] = (s_1, \dots, s_N)$, either by assigning all spins the same value or choosing each spin randomly as $ +1 $ or $ -1 $.

${ ii. \ Attempt\ spin\ flip.}$
Randomly select a spin and flip it (i.e., multiply by $-1$). This produces a trial configuration $[s_1]$.

${iii. \ Accept/Reject\ move.}$
Compute the energy difference between the trial and original configurations:

\[\Delta E = E_1 - E_0\]

Determine the acceptance probability:

\[w = \frac{P[s_1]}{P[s_0]} = \exp(-\beta \Delta E)\]

Generate a random number $ r \in [0, 1] $ and apply the Metropolis criterion:

(a) If $ r \leq w $, accept the move: set $ [s_0] \leftarrow [s_1] $
(b) If $ r > w $, reject the move: retain $ [s_0] $

Return to step ii to continue sampling.

Code for running 2D ising MCMC Simulations#

If numba is installed, uncomment to benfit from jit acceleration!
To speed benchmarks to assess how long simulation would take

Show code cell source Hide code cell source

from numba import njit
import numpy as np
import ipywidgets as widgets

@njit
def compute_total_energy(spins, J, B):
    N = spins.shape[0]
    E = 0.0
    for i in range(N):
        for j in range(N):
            S = spins[i, j]
            neighbors = spins[(i+1)%N, j] + spins[i, (j+1)%N]
            E -= J * S * neighbors
            E -= B * S
    return E

@njit
def run_ising2d(N, T, J, B, n_steps, out_freq):
    """
    Blazing fast 2D Ising Metropolis MCMC.

    Parameters:
    N        : Lattice size (N x N)
    T, J, B  : Temperature, coupling, field
    n_steps  : Number of Monte Carlo steps
    out_freq : Frequency of measurements

    Returns:
    S_out    : List of sampled spin states (sparse)
    E_out    : Per-spin energy list
    M_out    : Per-spin magnetization list
    """

    spins = 2 * (np.random.rand(N, N) < 0.5) - 1
    E_t = compute_total_energy(spins, J, B)
    M_t = np.sum(spins)

    # Preallocate memory
    n_out = n_steps // out_freq
    S_out = np.empty((n_out, N, N), dtype=np.int8)
    E_out = np.empty(n_out)
    M_out = np.empty(n_out)

    k = 0  # Output index

    for step in range(n_steps):
        
        i = np.random.randint(N)
        j = np.random.randint(N)

        S = spins[i, j]
        neighbors = spins[(i+1)%N, j] + spins[(i-1)%N, j] + spins[i, (j+1)%N] + spins[i, (j-1)%N]
        dE = 2 * S * (J * neighbors + B)

        if dE <= 0 or np.random.rand() < np.exp(-dE / T):
            spins[i, j] = -S
            E_t += dE
            M_t += -2 * S

        if step % out_freq == 0:
            S_out[k] = spins
            E_out[k] = E_t / N**2
            M_out[k] = M_t / N**2
            k += 1

    return S_out, E_out, M_out

def plot_ising(i):
    
    fig, ax  = plt.subplots(ncols=3, figsize=(10,4))
    
    ax[0].pcolor(S[i])
    ax[1].plot(E)
    ax[2].plot(M)
    
    ax[0].set(ylabel='$i$', xlabel='j')
    ax[1].set(ylabel='$E$', xlabel='steps')
    ax[2].set(ylabel='$M$', xlabel='steps')
    fig.tight_layout()
    plt.show()

#Simulation Parameters
params = {'N':40,
          'J':1,
          'B':0, 
          'T': 4,
          'n_steps': 100000, 
          'out_freq': 10}

# Simulation
S, E, M = run_ising2d(**params)

np.array(S).shape, np.array(E).shape, np.array(M).shape

((10000, 40, 40), (10000,), (10000,))

widgets.interact(plot_ising, i =(0, params['n_steps']//params['out_freq']-1))

<function __main__.plot_ising(i)>

Carry out explorative simulations#

How do we know if the simulation has done enough sampling? How do we assess convergence and errors?
Test the dependence of observables on system size.
Test the dependence of observables on initital conditions.
Vary intensive parameters, e.g., temperature and field strength. Investigate changes in observables such as magnetization, energy, susceptibility, heat capacity.
Think about alternative ways of accelerating and enhancing the sampling.

Parameter sweeps#

#Simulation Parameters
params = {'N':40,
          'J':1,
          'B':0, 
          'T': 4,
          'n_steps': 1000000, 
          'out_freq': 10}

Ts = np.linspace(1, 4, 40) 
Es, Ms, Cs, Xs = [], [], [], []  

for T in Ts:
    
    params['T']=T
    
    S, E, M = run_ising2d(**params)

    # Save last 90% percent of data
    idx = int(len(E) * 0.1) 
    
    E = E[idx:]
    M = M[idx:]
    
    Es.append(np.mean(E))
    Ms.append(np.mean(M))
    
    Cs.append(np.var(E)/T**2) 
    Xs.append(np.var(M)/T)

fig, ax  = plt.subplots(ncols=2, nrows=2, figsize=(8,6))

ax[0,0].scatter(Ts, Es)
ax[0,0].set(ylabel='$E(T)$')

ax[0,1].scatter(Ts, Ms)
ax[0,1].set(ylabel='$M(T)$')

ax[1,0].scatter(Ts, Cs)
ax[1,0].set(ylabel='$C_v(T)$')

ax[1,1].scatter(Ts, Xs)
ax[1,1].set(ylabel='$\Xi(T)$')
fig.tight_layout()

../_images/64a343fd05c181d6bc9f0e1ddfbe3905b60aa28b0758644d6694f3308baf5607.png

Autocorrelation functions#

For any observable $ A(t) $, you can compute:

\[ C_A(n) = \langle A(t) A(t+n) \rangle - \langle A(t) \rangle^2 \]

or the normalized version of autocorrelation funtion:

\[ \rho_A(n) = \frac{\langle A(t) A(t+n) \rangle - \langle A(t) \rangle^2}{\langle A(t)^2 \rangle - \langle A(t) \rangle^2} \]

Observable	Type	Purpose
$ M(t) $, $ E(t) $	Global observables	System-wide relaxation
$ s_i s_{i+r} $	Spatial correlation	Correlation length / criticality
$ C(n) $	General ACF	Basis for estimating correlation time

Integrated and exponential autocorrelation time#

Once you compute $ C(n) $, you can also estimate:
Integrated autocorrelation time $ \tau_{\text{int}} $: $$ \tau_{\text{int}} = \frac{1}{2} + \sum_{n=1}^{n_{\text{cut}}} \frac{C(n)}{C(0)} $$
Exponential autocorrelation time by fitting $ C(n) \sim \exp(-n / \tau_{\text{exp}}) $

#Simulation Parameters
params = {'N':40,
          'J':1,
          'B':0, 
          'T': 4,
          'n_steps': 100000, 
          'out_freq': 10}

# Simulation
S, E, M = run_ising2d(**params)

# Example usage
#M = np.random.randn(1000)  # Replace with your magnetization time series
acorr = autocorr(M)

rho_t = acorr / acorr[0]

# Plotting
plt.plot(rho_t[:10000])
plt.xlabel('Steps, $n$')
plt.ylabel(r'$\langle M(t) M(t+n) \rangle$')
plt.title('Autocorrelation of Magnetization')
plt.grid(True)
plt.tight_layout()
plt.show()

../_images/ee574e2c85e451a192e830287bd9353cb9d89651776a0367be545ea2c39f3d22.png

# Example usage:
x = np.random.randn(1000).cumsum()  # Strongly correlated data
tau_exp, N_eff = estimate_exp_autocorr_time(x)

print(f"Exponential autocorrelation time: {tau_exp:.2f}")
print(f"Effective sample size: {N_eff:.1f} (out of {len(x)} total samples)")

../_images/584fe0f7f9807ef3235cb3ee1dd39741af363767ba839b507bddacbc6224477c.png

Exponential autocorrelation time: 180.09
Effective sample size: 2.8 (out of 1000 total samples)

Problems#

Problem-1#

Revisit the example MCMC simulation for determining $\pi$ value. Vary the size of the displacement to determine the optimal size that generates quickest convergence to the value of $\pi$

Problem-2#

Carry out MC simulation of 2D ising spin model for various lattice sizes $N= 16,32, 64$ at temperatures above and below critical e.g $T<T_c$ and $T>T_c$.

How long does it take to equilibrate system as a function of size and as a function of T?
Plot some observables as a function of number of samples states to show that the system is carrying out some sort of random walk in the configurational space.
How do profiles of Energy vs T, Magnetization vs T and heat capacity vs T, and susceptibility vs T change as a function of size of our lattice.
Does $J>0$ and $J<0$ change the nature of phase transition?

Problem-3#

Compute correlation functions of spin variable, that is how correlated are spins as a function of distance on a lattice, $L$. $C(L)=\langle s_i s_{i+L}\rangle -\langle s_i\rangle \langle s_{i+L}\rangle $ Make sure to account for the periodic boundary conditions!

Note that you can pick a special tagged spin and calculate correlation function of taged spin ($s_13$ for instance) with any other as a function of lattice spearation by averaging over produced MC configurations. Or you can take advantage of the fact that there are no priviledged spins and average over many spins and average over MC configruations e.g $s_1, s_2, ...$. E.g you can pick a horizontal line of spins and run a summation for each fixed r_ab distance.

Problem-4#

Take a 20 by 20 lattice and equilibriate the system with a value of extneral field B equal to +1. Now slowly change h to −1 in discrete steps during each of these steps, use the previously equilibriated configuration as an input to the system to undergo equilibriation again.

Caluclate average and variance quantities (e.g E, M, C etc). Notice anything interesing :)

Observable	Type	Purpose
\( M(t) \), \( E(t) \)	Global observables	System-wide relaxation
\( s_i s_{i+r} \)	Spatial correlation	Correlation length / criticality
\( C(n) \)	General ACF	Basis for estimating correlation time