Monte Carlo simulation for the classical isotropic 3D Heisenberg model

In summary, the conversation involves the implementation of a Monte Carlo simulation for the isotropic 3D Heisenberg model using Python. The code includes functions for initializing spins, performing Metropolis algorithm updates, and calculating the magnetization and susceptibility. The simulation is run for a range of temperatures, and the results are plotted in graphs of susceptibility and magnetization vs temperature.
  • #1
IWantToLearn
94
0
TL;DR Summary
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

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.
 
Last edited by a moderator:

Similar threads

  • Atomic and Condensed Matter
Replies
3
Views
736
  • Advanced Physics Homework Help
Replies
1
Views
1K
  • Atomic and Condensed Matter
Replies
1
Views
2K
  • Programming and Computer Science
Replies
1
Views
1K
  • Programming and Computer Science
Replies
6
Views
4K
  • Mechanics
Replies
2
Views
4K
  • Programming and Computer Science
Replies
4
Views
4K
  • Advanced Physics Homework Help
Replies
7
Views
1K
Back
Top