A 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.