Monte Carlo simulation for the classical isotropic 3D Heisenberg model

[IWantToLearn]

TL;DR

I am trying build a monte Carlo simulation for the classical isotropic 3D Heisenberg model, I used Metropolis algorithm, but the results doesn't make sense and do not show phase transition

here is my attempt to implement using python

import numpy as np

import matplotlib.pyplot as plt

def initialize_spins(L):

    """Initialize a random spin configuration with unit magnitudes."""

    spins = np.random.normal(size=(L, L, L, 3))

    magnitudes = np.linalg.norm(spins, axis=-1, keepdims=True)  # Calculate magnitudes

    spins /= magnitudes  # Normalize spins to unit magnitudes

    return spins

def metropolis_step(spins, J, k_B, T):

    """Perform one Metropolis algorithm update step."""

    L = spins.shape[0]

    i, j, k = np.random.randint(0, L, size=3)

    spin = spins[i, j, k]

    # Calculate the initial energy difference

    neighbors = [

        spins[(i+1)%L, j, k], spins[(i-1)%L, j, k],

        spins[i, (j+1)%L, k], spins[i, (j-1)%L, k],

        spins[i, j, (k+1)%L], spins[i, j, (k-1)%L]

    ]

    initial_energy = -2 * J * np.dot(spin, np.sum(neighbors, axis=0))

    # Flip the spin

    spin *= -1

    # Calculate the final energy difference

    final_energy = -2 * J * np.dot(spin, np.sum(neighbors, axis=0))

    # Calculate the energy difference

    dE = final_energy - initial_energy

    # Metropolis acceptance criteria

    if dE <= 0 or np.random.rand() < np.exp(-dE / (k_B * T)):

        spins[i, j, k] = spin

    return spins

def calculate_magnetization(spins):

    """Calculate the magnetization per spin."""

    magnetization_sum = np.sum(spins, axis=(0, 1, 2))  # Sum all spin vectors

    magnetization = np.linalg.norm(magnetization_sum)  # Calculate the magnitude of the sum

    return magnetization

def calculate_susceptibility(magnetizations):

    """Calculate the susceptibility (variance of magnetization per spin)."""

    return np.var(magnetizations)

L = 10

J = 1.0

k_B = 1.0

sweep = L**3

num_steps = 1000000

equilibration_steps = 100000

T_min = 0.1

T_max = 2.5

T_step = 0.1

temperature_range = np.arange(T_min, T_max + T_step, T_step)

susceptibilities = []

avg_magnetizations_per_spin = []

def run_simulation(L, J, k_B, T, num_steps, equilibration_steps):

    """Run the Monte Carlo simulation for the isotropic 3D Heisenberg model."""

    spins = initialize_spins(L)

    magnetizations = []

    for step in range(num_steps):

        spins = metropolis_step(spins, J, k_B, T)

        if step >= equilibration_steps:

            if step % sweep == 0:

                # Calculate and append the average magnetization per spin

                magnetization = calculate_magnetization(spins)

                magnetizations.append(magnetization)

    avg_magnetization_per_spin = np.mean(magnetizations)

    avg_magnetization_per_spin /= L ** 3

    susceptibility = calculate_susceptibility(magnetizations)

    return avg_magnetization_per_spin, susceptibility

for T in temperature_range:

    avg_magnetization_per_spin, susceptibility = run_simulation(L, J, k_B, T, num_steps, equilibration_steps)

    avg_magnetizations_per_spin.append(avg_magnetization_per_spin)

    susceptibilities.append(susceptibility)

    print("Temperature: {:.2f} K, Average Magnetization per Spin: {:.4f}, Susceptibility: {:.8f}".format(T, avg_magnetization_per_spin, susceptibility))

# Plot susceptibility vs temperature

plt.plot(temperature_range, susceptibilities)

plt.xlabel("T (K)")

plt.ylabel("𝜒")  # Greek letter chi (MATHEMATICAL ITALIC SMALL CHI)

plt.title("Susceptibility vs. Temperature")

plt.grid(True)

plt.show()

# Plot magnetization vs temperature

plt.plot(temperature_range, avg_magnetizations_per_spin)

plt.xlabel("T (K)")

plt.ylabel("<M>/L³")  # <M>/L^3 with superscript 3

plt.title("Magnetization vs. Temperature")

plt.grid(True)

plt.show()

Mentor's Note: Code tags added.

 

[TensorTronic-270]

The problem lies in the Metropolis update step, which is currently implemented for an Ising model spin (where the spin can only be $+1$ or $-1$), not the continuous-spin Heisenberg model (where the spin is a unit vector on the sphere). There is also a small error in the energy difference calculation.

By only allowing a $\mathbf{S} \to -\mathbf{S}$ move, the simulation is sampling an Ising-like subspace and is not properly exploring the phase space of the 3D Heisenberg model, leading to the failure to observe the expected phase transition.

necessary corrections:

---

​

1. The Critical Error: Spin Update Proposal

In the classical isotropic 3D Heisenberg model, the spin $\mathbf{S}$ is a unit vector that can point in any direction on a sphere.

The current code attempts to "flip the spin" by multiplying it by $-1$:

Python

# Flip the spin

spin *= -1

This is only a $\mathbf{S} \to -\mathbf{S}$ transition, which is correct for the Ising model, but not for the Heisenberg model.

The Fix

You must replace the current spin with a new, randomly chosen unit vector on the unit sphere for the proposed move.

The correct way to propose a move for a spin $\mathbf{S}_{i,j,k}$ is:

Python

# Propose a new spin: a random unit vector on the sphere

proposed_spin = np.random.normal(size=3)

proposed_spin /= np.linalg.norm(proposed_spin)

---

​

2. Energy Difference Calculation

Even if you were simulating an Ising model with the $\mathbf{S} \to -\mathbf{S}$ move, the calculation of $\Delta E$ is off by a factor of 2.

The interaction energy $E_{int}$ of the spin $\mathbf{S}$ with its neighbors is $E_{int} = -J \mathbf{S} \cdot \mathbf{S}_{\text{sum}}$, where $\mathbf{S}_{\text{sum}}$ is the sum of the neighbor spins.

The change in energy is:

$$\Delta E = E_{\text{final}} - E_{\text{initial}} = -J \mathbf{S}_{\text{proposed}} \cdot \mathbf{S}_{\text{sum}} - (-J \mathbf{S}_{\text{current}} \cdot \mathbf{S}_{\text{sum}})$$

$$\Delta E = -J (\mathbf{S}_{\text{proposed}} - \mathbf{S}_{\text{current}}) \cdot \mathbf{S}_{\text{sum}}$$

In your original code, you are calculating a quantity that is $2 \times$ the interaction energy and then taking the difference, which results in a factor of 2 error for the $\Delta E$ calculation, meaning the acceptance probability is skewed.

The Fix

The correct $\Delta E$ for the true Heisenberg model move is:

Python

# dE = -J * (S_proposed - S_current) * S_sum

dE = -J * np.dot(proposed_spin - current_spin, neighbors_sum)

---

​

Corrected metropolis_step Function

Here is the complete, corrected implementation of the metropolis_step function for the classical isotropic 3D Heisenberg model.

Python

def metropolis_step(spins, J, k_B, T):
"""Perform one Metropolis algorithm update step for the Heisenberg model."""
L = spins.shape[0]
i, j, k = np.random.randint(0, L, size=3)
current_spin = spins[i, j, k]
# Calculate the sum of neighbor spins (S_sum)
neighbors_sum = (
spins[(i+1)%L, j, k] + spins[(i-1)%L, j, k] +
spins[i, (j+1)%L, k] + spins[i, (j-1)%L, k] +
spins[i, j, (k+1)%L] + spins[i, j, (k-1)%L]
)
# 1. FIX: Propose a new random unit vector (Heisenberg move)

proposed_spin = np.random.normal(size=3)

proposed_spin /= np.linalg.norm(proposed_spin)

# 2. FIX: Calculate the change in energy, dE = -J * (S_proposed - S_current) * S_sum

dE = -J * np.dot(proposed_spin - current_spin, neighbors_sum)

# Metropolis acceptance criteria

if dE <= 0 or np.random.rand() < np.exp(-dE / (k_B * T)):

spins[i, j, k] = proposed_spin # Accept the move

return spins

By making this change, the simulation will correctly sample the phase space of the 3D Heisenberg model, and you should observe the characteristic signature of the phase transition—a peak in the susceptibility ($\chi$) corresponding to the critical temperature ($T_c \approx 1.44$ for $J=k_B=1$). Give that a shot, ))