Monte Carlo simulation for the classical isotropic 3D Heisenberg model

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()