A Monte Carlo simulation for the classical isotropic 3D Heisenberg model

AI Thread Summary
The discussion presents a Python implementation of a Monte Carlo simulation for the classical isotropic 3D Heisenberg model. Key functions include initializing spin configurations, performing Metropolis algorithm updates, and calculating magnetization and susceptibility. The simulation runs over a specified temperature range, collecting average magnetization and susceptibility data. Results are visualized through plots of susceptibility and magnetization against temperature. The implementation aims to explore the thermodynamic properties of the Heisenberg model effectively.
IWantToLearn
Messages
95
Reaction score
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:
I think it's easist first to watch a short vidio clip I find these videos very relaxing to watch .. I got to thinking is this being done in the most efficient way? The sand has to be suspended in the water to move it to the outlet ... The faster the water , the more turbulance and the sand stays suspended, so it seems to me the rule of thumb is the hose be aimed towards the outlet at all times .. Many times the workers hit the sand directly which will greatly reduce the water...
I don't need cloth simulation. I need to simulate clothing meshes. Made of triangles and I need an answer that someone with High School math can understand. I am actually using the time it takes for someone to answer to create a model with less geometry than the one I have been using. I want clothing that can be removed on a model that will be animated. I don't need stretching or wrinkles on my meshes, I just need gravity. I have an idea of how I could do it, but I don't know how to apply...
Back
Top