- #1

IWantToLearn

- 95

- 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

Mentor's Note: Code tags added.

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: