Ising and Potts updates.

I build one local Monte Carlo update before repeating it. The square lattice is periodic. In the Ising model a spin has value -1 or +1; the dimensionless interaction strength is one and temperature is measured in units of the interaction energy divided by Boltzmann’s constant. For Potts grain growth the integers identify grains, not compositions. Neighbor-copy proposals describe the supplied grain-growth kinetics and should not be presented as unbiased equilibrium sampling.

Python notebook ↓ · Julia notebook ↓ · Complete collection ↓

Run notebook cells in order. Extract the complete collection without changing its folder structure: the phase-field lesson imports the original kernel from its parent folder. Use Python with NumPy, Matplotlib and Jupyter, or Julia with IJulia. Julia Fourier examples also require FFTW. The download includes installation commands and the actual check report.

Step 1: Inspect periodic neighbors

I first test the boundary indexing on a grid of known values. Opposite edges are adjacent. Rows and columns identify sites; they are dimensionless indices.

Python

import numpy as np
rng = np.random.default_rng(7)
def neighbor_values(state, row, column):
    rows, columns = state.shape
    return [state[(row + 1) % rows, column], state[(row - 1) % rows, column],
            state[row, (column + 1) % columns], state[row, (column - 1) % columns]]

test_grid = np.arange(9).reshape(3, 3)
assert neighbor_values(test_grid, 0, 0) == [3, 6, 1, 2]
print("Periodic neighbor indexing passed.")

using Random
rng = MersenneTwister(7)
function neighbor_values(state, row, column)
    rows, columns = size(state)
    return [state[mod1(row+1,rows),column], state[mod1(row-1,rows),column],
            state[row,mod1(column+1,columns)], state[row,mod1(column-1,columns)]]
end
test_grid = [0 1 2; 3 4 5; 6 7 8]
@assert neighbor_values(test_grid, 1, 1) == [3, 6, 1, 2]
println("Periodic neighbor indexing passed.")

Step 2: Calculate the Ising energy change

Flipping a spin changes each of its four bond energies. In a uniform positive-spin neighborhood the change is +8. I check this before using a probability rule.

Python

def ising_energy_change(old_spin, neighbors):
    return 2 * old_spin * sum(neighbors)

assert ising_energy_change(1, [1, 1, 1, 1]) == 8
assert ising_energy_change(-1, [1, 1, 1, 1]) == -8
print("Ising energy changes passed.")

function ising_energy_change(old_spin, neighbors)
    return 2 * old_spin * sum(neighbors)
end
@assert ising_energy_change(1, [1,1,1,1]) == 8
@assert ising_energy_change(-1, [1,1,1,1]) == -8
println("Ising energy changes passed.")

Step 3: Accept or reject a trial

I always accept a decrease in energy. For an increase I accept with probability exp(-ΔE/T). The dimensionless temperature T must be positive in this implementation. This rule permits occasional uphill moves; energy need not decrease at every update.

Python

def accept_trial(energy_change, temperature, rng):
    if temperature <= 0:
        raise ValueError("Use positive temperature in this example.")
    if energy_change <= 0:
        return True
    return rng.random() < np.exp(-energy_change / temperature)

assert accept_trial(-8, 1.5, rng)
assert accept_trial(0, 1.5, rng)
assert 0 < np.exp(-8 / 1.5) < 1
print("Acceptance checks passed.")

function accept_trial(energy_change, temperature, rng)
    temperature > 0 || error("Use positive temperature in this example.")
    if energy_change <= 0
        return true
    end
    return rand(rng) < exp(-energy_change / temperature)
end
@assert accept_trial(-8, 1.5, rng)
@assert accept_trial(0, 1.5, rng)
@assert 0 < exp(-8 / 1.5) < 1
println("Acceptance checks passed.")

Step 4: A single Ising trial and repeated trials

I choose a random site with replacement. One sweep consists of as many attempted updates as there are sites, so some sites can be chosen more than once. A spin flip does not conserve total spin.

Python

def ising_trial(state, temperature, rng):
    row = rng.integers(state.shape[0])
    column = rng.integers(state.shape[1])
    old = state[row, column]
    change = ising_energy_change(old, neighbor_values(state, row, column))
    if accept_trial(change, temperature, rng):
        state[row, column] = -old

spins = rng.choice([-1, 1], (16, 16))
for trial in range(spins.size * 10):
    ising_trial(spins, 1.5, rng)
assert set(np.unique(spins)) <= {-1, 1}
print("Ising state values:", np.unique(spins))

function ising_trial!(state, temperature, rng)
    row = rand(rng, axes(state,1))
    column = rand(rng, axes(state,2))
    old = state[row,column]
    change = ising_energy_change(old, neighbor_values(state,row,column))
    if accept_trial(change, temperature, rng)
        state[row,column] = -old
    end
end
spins = rand(rng, [-1,1], 16,16)
for trial in 1:length(spins)*10
    ising_trial!(spins, 1.5, rng)
end
@assert all(value in (-1,1) for value in spins)
println("Ising state values: ", unique(spins))

Step 5: Change the local rule for Potts grain growth

For a Potts trial I propose a neighbor’s label and count unlike neighbors before and after the change. The energy change is the difference in these counts, in units of the grain-boundary interaction. No new grain label can be created by copying a neighbor. I check both the local change and that property.

Python

def potts_energy_change(old, candidate, neighbors):
    before = sum(value != old for value in neighbors)
    after = sum(value != candidate for value in neighbors)
    return after - before

assert potts_energy_change(1, 2, [2, 2, 2, 1]) == -2

def potts_trial(state, temperature, rng):
    row = rng.integers(state.shape[0])
    column = rng.integers(state.shape[1])
    neighbors = neighbor_values(state, row, column)
    candidate = neighbors[rng.integers(4)]
    change = potts_energy_change(state[row, column], candidate, neighbors)
    if accept_trial(change, temperature, rng):
        state[row, column] = candidate

grains = rng.permutation(np.arange(1, 257)).reshape(16, 16)
initial_labels = set(grains.ravel())
for trial in range(grains.size * 10):
    potts_trial(grains, 0.3, rng)
assert set(grains.ravel()) <= initial_labels
print("Grains initially:", len(initial_labels), "remaining labels:", len(np.unique(grains)))

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(9, 4))
axes[0].imshow(spins, vmin=-1, vmax=1, cmap="coolwarm")
axes[0].set_title("Ising: spins -1 and +1")
axes[1].imshow(grains, cmap="tab20")
axes[1].set_title("Potts: grain labels")
for ax in axes:
    ax.set_xlabel("Column index")
    ax.set_ylabel("Row index")
plt.show()

function potts_energy_change(old, candidate, neighbors)
    return count(!=(candidate), neighbors) - count(!=(old), neighbors)
end
@assert potts_energy_change(1,2,[2,2,2,1]) == -2
function potts_trial!(state, temperature, rng)
    row = rand(rng, axes(state,1))
    column = rand(rng, axes(state,2))
    neighbors = neighbor_values(state,row,column)
    candidate = rand(rng, neighbors)
    change = potts_energy_change(state[row,column],candidate,neighbors)
    if accept_trial(change,temperature,rng)
        state[row,column] = candidate
    end
end
grains = reshape(randperm(rng,256),16,16)
initial_labels = Set(grains)
for trial in 1:length(grains)*10
    potts_trial!(grains,0.3,rng)
end
@assert issubset(Set(grains),initial_labels)
println("Grains initially: ", length(initial_labels), "; remaining labels: ", length(unique(grains)))

What I checked

I executed all Python notebook cells in a clean kernel and all Julia script stages. The checks apply to the short examples shown here. They do not establish grid convergence or reproduce every long simulation. Julia notebook-kernel execution is not included in these checks; the notebook code cells are identical to the tested Julia script.

Return to the code collection →