Other Methods for sampling Ising models

Contents

Other Methods for sampling Ising models#

Wolf Algorithm (cluster methods):

Near the critical point $T_{c}$ where the system develops a magnetization, any single-spin-flip dynamics becomes very slow (the correlation time diverges). Wolff came up with a clever method to flip whole clusters of spins!

Pick a spin at random, remember its direction $s = \pm 1$ , and flip it.
For each of the four neighboring spins, if it is in the direction $s$ , flip it with probability $p = 1 - e^{β 2 J}$ .
For each of the new flipped spins, recursively flip their neighbors as in step 2

from numba import njit
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

@njit
def wolff_ising(spins, J=1.0, T=2.0, nsteps=5000):
    N = spins.shape[0]
    beta = 1.0 / T
    p_add = 1.0 - np.exp(-2 * beta * J)

    max_cluster_size = N * N
    cluster_mask = np.zeros((N, N), dtype=np.uint8)

    S_out = []
    M_out = []

    for step in range(nsteps):
        cluster_mask[:, :] = 0

        # Seed
        i, j = np.random.randint(N), np.random.randint(N)
        spin0 = spins[i, j]
        cluster_mask[i, j] = 1

        # Cluster grows from stack
        stack = [(i, j)]
        while stack:
            x, y = stack.pop()
            for dx, dy in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
                xn, yn = (x + dx) % N, (y + dy) % N
                if cluster_mask[xn, yn] == 0 and spins[xn, yn] == spin0:
                    if np.random.rand() < p_add:
                        cluster_mask[xn, yn] = 1
                        stack.append((xn, yn))

        # Flip the cluster
        for x in range(N):
            for y in range(N):
                if cluster_mask[x, y]:
                    spins[x, y] *= -1

        if step % 10 == 0:
            S_out.append(spins.copy())
            M_out.append(np.mean(spins))

    return S_out, M_out

# Run simulation
N = 20
spins_init = np.random.choice([-1, 1], size=(N, N))
S_out, M_out = wolff_ising(spins_init, J=1.0, T=2.2, nsteps=500) # at critical temp

Visualize Simulation#

# Create animation
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
im = ax[0].imshow(S_out[0], cmap="coolwarm", vmin=-1, vmax=1)
line, = ax[1].plot([], [], lw=2)
ax[1].set_xlim(0, len(M_out))
ax[1].set_ylim(-1.05, 1.05)
ax[1].set_xlabel("Step")
ax[1].set_ylabel("Magnetization")
ax[0].set_title("Spin Configuration")
ax[1].set_title("Magnetization")

def update(frame):
    im.set_data(S_out[frame])
    line.set_data(np.arange(frame + 1), M_out[:frame + 1])
    return im, line

ani = FuncAnimation(fig, update, frames=len(S_out), interval=100, blit=True)
plt.close()

HTML(ani.to_jshtml())