FYS-STK4155/doc/Programs/HeisenbergIsingModels/Heisenberg.ipynb

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import os
import sys

In [2]:

# Define parameters
N = 10  # Size of the cubic lattice (N x N x N)
J = 1.0  # Interaction strength
kB = 1.0  # Boltzmann constant

T_values = [2.5, 2.4, 2.3, 2.2, 2.1, 2.0, 1.9, 1.8, 1.7, 1.6, 1.54, 1.52, 1.50, 1.49, 1.48, 1.47, 1.46,
             1.45, 1.44, 1.43, 1.42, 1.41, 1.40, 1.39, 1.38, 1.37, 1.36, 1.34, 1.32, 1.30, 1.25, 1.20,
             1.1, 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1]

total_steps = 10**6
equilibration_steps = total_steps/10

start_save_config = total_steps - 1000

progress = 0

In [3]:

# Function to initialize the lattice with random unit vectors
def initialize_lattice(N):
    lattice = np.random.rand(N, N, N, 3) * 2 - 1.0  # Random values in [-1.0, 1.0]
    lattice /= np.linalg.norm(lattice, axis=-1, keepdims=True)
    return lattice

In [4]:

# Function to get the neighbors of a spin
def get_neighbors(lattice, i, j, k):
    neighbors = []
    for di, dj, dk in [(-1, 0, 0), (1, 0, 0), (0, -1, 0), (0, 1, 0), (0, 0, -1), (0, 0, 1)]:
        ni, nj, nk = (i + di) % N, (j + dj) % N, (k + dk) % N  # Apply periodic boundary conditions
        neighbors.append(lattice[ni, nj, nk])
    return neighbors

In [5]:

# Function to calculate the energy of a spin and its neighbors
def calculate_energy(spin, neighbors):
    neighbor_sum = np.sum(neighbors, axis=0)
    energy = -J * np.dot(spin, neighbor_sum)
    return energy

In [6]:

# Function to calculate the magnetization magnitude of the lattice
def calculate_magnetization(lattice):
    magnetization = [0.0, 0.0, 0.0]
    for i in range(N):
        for j in range(N):
            for k in range(N):
                # Access the lattice element at position (i, j, k)
                magnetization[0] += lattice[i, j, k, 0]
                magnetization[1] += lattice[i, j, k, 1]
                magnetization[2] += lattice[i, j, k, 2]

    magnitude = magnetization[0]*magnetization[0] + magnetization[1]*magnetization[1] + magnetization[2]*magnetization[2]
    magnitude = np.sqrt(magnitude)
    magnitude = magnitude/(N**3)
    return magnitude

In [7]:

# Function to perform a Metropolis update
def metropolis_update(lattice, temperature):
    global acc_rate

    i, j, k = np.random.randint(0, N, size=3)
    spin = lattice[i, j, k]
    neighbors = get_neighbors(lattice, i, j, k)
    
    # Calculate energy before the update
    energy_before = calculate_energy(spin, neighbors)
    
    # Propose a new spin configuration
    new_spin = (np.random.rand(3) * 2 - 1.0)  # Random values between -1.0 and 1.0 for each component
    new_spin /= np.linalg.norm(new_spin)

    # Calculate energy after the update
    energy_after = calculate_energy(new_spin, neighbors)
    
    # Calculate energy difference
    delta_energy = energy_after - energy_before
    
    # Metropolis acceptance criteria
    if delta_energy <= 0 or np.random.rand() < np.exp(-delta_energy / (kB * temperature)):
        lattice[i, j, k] = new_spin
        acc_rate += 1

In [8]:

# Function to save Configurations as an .npy file (Only the X Components)
def save_configuration_X(config, temperature, step, lattice_size):
    folder_name = f"{lattice_size}X{lattice_size}X{lattice_size}_Steps{total_steps}/{temperature:.2f}/X_Comp"
    file_name = f"X_ConfigFile_Size{lattice_size}x{lattice_size}x{lattice_size}_M0.00_T{temperature:.2f}_StepNum{step}.npy"
    folder_path = os.path.join(os.getcwd(), folder_name)
    file_path = os.path.join(folder_path, file_name)
    if not os.path.exists(folder_path):
        os.makedirs(folder_path)

    np.save(file_path, config)

In [9]:

# Function to save Configurations as an .npy file (Only the Y Components)
def save_configuration_Y(config, temperature, step, lattice_size):
    folder_name = f"{lattice_size}X{lattice_size}X{lattice_size}_Steps{total_steps}/{temperature:.2f}/Y_Comp"
    file_name = f"Y_ConfigFile_Size{lattice_size}x{lattice_size}x{lattice_size}_M0.00_T{temperature:.2f}_StepNum{step}.npy"
    folder_path = os.path.join(os.getcwd(), folder_name)
    file_path = os.path.join(folder_path, file_name)
    if not os.path.exists(folder_path):
        os.makedirs(folder_path)

    np.save(file_path, config)

In [10]:

# Function to save Configurations as an .npy file (Only the Z Components)
def save_configuration_Z(config, temperature, step, lattice_size):
    folder_name = f"{lattice_size}X{lattice_size}X{lattice_size}_Steps{total_steps}/{temperature:.2f}/Z_Comp"
    file_name = f"Z_ConfigFile_Size{lattice_size}x{lattice_size}x{lattice_size}_M0.00_T{temperature:.2f}_StepNum{step}.npy"
    folder_path = os.path.join(os.getcwd(), folder_name)
    file_path = os.path.join(folder_path, file_name)
    if not os.path.exists(folder_path):
        os.makedirs(folder_path)

    np.save(file_path, config)

In [11]:

# Function to save the lattice as an .npy file
def save_lattice(lattice, temperature, step, lattice_size):
    folder_name = f"{lattice_size}X{lattice_size}X{lattice_size}_Steps{total_steps}/{temperature:.2f}/Lattice"
    file_name = f"LatticeFile_Size{lattice_size}x{lattice_size}x{lattice_size}_T{temperature:.2f}_StepNum{step}.npy"
    folder_path = os.path.join(os.getcwd(), folder_name)
    file_path = os.path.join(folder_path, file_name)
    if not os.path.exists(folder_path):
        os.makedirs(folder_path)

    np.save(file_path, lattice)

In [12]:

# Function to perform the Monte Carlo simulation
def monte_carlo_simulation(lattice, temperature, equilibration_steps, total_steps):
    magnetization_values = []

    for step in range(total_steps):
        metropolis_update(lattice, temperature)
        
        if step >= equilibration_steps and (step - equilibration_steps) % (total_steps // 100) == 0:
            # Calculate magnetization and save measurements
            magnetization = calculate_magnetization(lattice)
            magnetization_values.append(magnetization)
        
        if step >= start_save_config:
            x_components = lattice[:, :, :, 0]  # Extract x components (index 0)
            y_components = lattice[:, :, :, 1]  # Extract y components (index 1)
            z_components = lattice[:, :, :, 2]  # Extract z components (index 2)
            save_configuration_X(x_components,temperature,step,N)
            save_configuration_Y(y_components,temperature,step,N)
            save_configuration_Z(z_components,temperature,step,N)

            # Save the lattice as well
            save_lattice(lattice, temperature, step, N)

    # Calculate average magnetization and error bar
    avg_magnetization = np.mean(magnetization_values)
    error_bar = np.std(magnetization_values) / np.sqrt(len(magnetization_values))
    # Calculating magnetic susceptibility
    Chi = (np.var(magnetization_values))/temperature

    # Return average magnetization, error bar
    return avg_magnetization, error_bar, Chi

In [13]:

# Main simulation loop
global lattice
lattice = initialize_lattice(N)

avg_magnetizations = []
error_bars = []
susceptibility_values = []

file_name_ = f"simulation_results_{N}X{N}X{N}_Steps{total_steps}.txt" 
file = open(file_name_ , 'w')

for T in T_values:
    acc_rate = 0
    avg_mag, error, Chi = monte_carlo_simulation(lattice, T, equilibration_steps, total_steps)
    avg_magnetizations.append(avg_mag)
    error_bars.append(error)
    susceptibility_values.append(Chi)

    progress = progress + 1
    num_of_temp = len(T_values)
    print(f"{progress:02d}/{num_of_temp} - Temperature: {T:.2f}, Magnetization: {avg_mag:.8f}, Error: {error:.8f}, Susceptibility: {Chi:.8f}, Acc_Rate: {acc_rate}")


    line = f"Temperature: {T:.2f}, Magnetization: {avg_mag:.8f}, Error: {error:.8f}, Susceptibility: {Chi:.8f}, Acc_Rate: {acc_rate}\n"
    sys.stdout.flush()  # Force flushing the output buffer

    file.write(line)
    file.flush()  # Force flushing the file buffer

file.close()

01/43 - Temperature: 2.50, Magnetization: 0.05067800, Error: 0.00231309, Susceptibility: 0.00019261, Acc_Rate: 694993
02/43 - Temperature: 2.40, Magnetization: 0.05701614, Error: 0.00233589, Susceptibility: 0.00020461, Acc_Rate: 682533
03/43 - Temperature: 2.30, Magnetization: 0.05745391, Error: 0.00228808, Susceptibility: 0.00020486, Acc_Rate: 667664
04/43 - Temperature: 2.20, Magnetization: 0.06127256, Error: 0.00251489, Susceptibility: 0.00025874, Acc_Rate: 653054
05/43 - Temperature: 2.10, Magnetization: 0.06650248, Error: 0.00342411, Susceptibility: 0.00050248, Acc_Rate: 636446
06/43 - Temperature: 2.00, Magnetization: 0.07177253, Error: 0.00282175, Susceptibility: 0.00035830, Acc_Rate: 617125
07/43 - Temperature: 1.90, Magnetization: 0.07764729, Error: 0.00349707, Susceptibility: 0.00057929, Acc_Rate: 595959
08/43 - Temperature: 1.80, Magnetization: 0.07969424, Error: 0.00346145, Susceptibility: 0.00059908, Acc_Rate: 572639
09/43 - Temperature: 1.70, Magnetization: 0.10794836, Error: 0.00536923, Susceptibility: 0.00152622, Acc_Rate: 542048
10/43 - Temperature: 1.60, Magnetization: 0.12627725, Error: 0.00551466, Susceptibility: 0.00171065, Acc_Rate: 510216
11/43 - Temperature: 1.54, Magnetization: 0.23756270, Error: 0.00701833, Susceptibility: 0.00287865, Acc_Rate: 477083
12/43 - Temperature: 1.52, Magnetization: 0.19669505, Error: 0.00680513, Susceptibility: 0.00274203, Acc_Rate: 475510
13/43 - Temperature: 1.50, Magnetization: 0.26928236, Error: 0.00633434, Susceptibility: 0.00240743, Acc_Rate: 460500
14/43 - Temperature: 1.49, Magnetization: 0.28233166, Error: 0.00660430, Susceptibility: 0.00263457, Acc_Rate: 450501
15/43 - Temperature: 1.48, Magnetization: 0.31488006, Error: 0.00634694, Susceptibility: 0.00244968, Acc_Rate: 442054
16/43 - Temperature: 1.47, Magnetization: 0.27167505, Error: 0.00687240, Susceptibility: 0.00289162, Acc_Rate: 444030
17/43 - Temperature: 1.46, Magnetization: 0.28125360, Error: 0.00782904, Susceptibility: 0.00377839, Acc_Rate: 438247
18/43 - Temperature: 1.45, Magnetization: 0.31897099, Error: 0.00630256, Susceptibility: 0.00246552, Acc_Rate: 430177
19/43 - Temperature: 1.44, Magnetization: 0.33826857, Error: 0.00640851, Susceptibility: 0.00256681, Acc_Rate: 421624
20/43 - Temperature: 1.43, Magnetization: 0.33619013, Error: 0.00484151, Susceptibility: 0.00147526, Acc_Rate: 420482
21/43 - Temperature: 1.42, Magnetization: 0.40093416, Error: 0.00546004, Susceptibility: 0.00188949, Acc_Rate: 402824
22/43 - Temperature: 1.41, Magnetization: 0.40223867, Error: 0.00544269, Susceptibility: 0.00189082, Acc_Rate: 398312
23/43 - Temperature: 1.40, Magnetization: 0.41361242, Error: 0.00584992, Susceptibility: 0.00219996, Acc_Rate: 393440
24/43 - Temperature: 1.39, Magnetization: 0.42797076, Error: 0.00391937, Susceptibility: 0.00099463, Acc_Rate: 387163
25/43 - Temperature: 1.38, Magnetization: 0.45036921, Error: 0.00567672, Susceptibility: 0.00210164, Acc_Rate: 375626
26/43 - Temperature: 1.37, Magnetization: 0.45859763, Error: 0.00481411, Susceptibility: 0.00152249, Acc_Rate: 371869
27/43 - Temperature: 1.36, Magnetization: 0.44461117, Error: 0.00410840, Susceptibility: 0.00111699, Acc_Rate: 370419
28/43 - Temperature: 1.34, Magnetization: 0.46889421, Error: 0.00420211, Susceptibility: 0.00118597, Acc_Rate: 359287
29/43 - Temperature: 1.32, Magnetization: 0.51931686, Error: 0.00361887, Susceptibility: 0.00089292, Acc_Rate: 339340
30/43 - Temperature: 1.30, Magnetization: 0.54039866, Error: 0.00271472, Susceptibility: 0.00051021, Acc_Rate: 328362
31/43 - Temperature: 1.25, Magnetization: 0.59450172, Error: 0.00322554, Susceptibility: 0.00074909, Acc_Rate: 301172
32/43 - Temperature: 1.20, Magnetization: 0.61688533, Error: 0.00239843, Susceptibility: 0.00043144, Acc_Rate: 282441
33/43 - Temperature: 1.10, Magnetization: 0.70057186, Error: 0.00213217, Susceptibility: 0.00037196, Acc_Rate: 237369
34/43 - Temperature: 1.00, Magnetization: 0.75602324, Error: 0.00144626, Susceptibility: 0.00018825, Acc_Rate: 202588
35/43 - Temperature: 0.90, Magnetization: 0.79911328, Error: 0.00110576, Susceptibility: 0.00012227, Acc_Rate: 173813
36/43 - Temperature: 0.80, Magnetization: 0.83435035, Error: 0.00086618, Susceptibility: 0.00008440, Acc_Rate: 148781
37/43 - Temperature: 0.70, Magnetization: 0.86724151, Error: 0.00071909, Susceptibility: 0.00006648, Acc_Rate: 126742
38/43 - Temperature: 0.60, Magnetization: 0.89348401, Error: 0.00049948, Susceptibility: 0.00003742, Acc_Rate: 107594
39/43 - Temperature: 0.50, Magnetization: 0.91464745, Error: 0.00044023, Susceptibility: 0.00003488, Acc_Rate: 90697
40/43 - Temperature: 0.40, Magnetization: 0.93373802, Error: 0.00033215, Susceptibility: 0.00002482, Acc_Rate: 74396
41/43 - Temperature: 0.30, Magnetization: 0.95149632, Error: 0.00027094, Susceptibility: 0.00002202, Acc_Rate: 58787
42/43 - Temperature: 0.20, Magnetization: 0.96750309, Error: 0.00013980, Susceptibility: 0.00000879, Acc_Rate: 42481
43/43 - Temperature: 0.10, Magnetization: 0.98301393, Error: 0.00006677, Susceptibility: 0.00000401, Acc_Rate: 25197

In [14]:

# Set plot size and create the figure
fig, ax = plt.subplots(figsize=(9, 5))

# Ensure error_bars has the same number of elements as T_values
error_bars = error_bars[:len(T_values)]

# Plot results with error bars and lines
plt.errorbar(T_values, avg_magnetizations, yerr=error_bars, fmt='o-', label='Magnetization', 
             ecolor='red', markersize=1, elinewidth=1, linewidth=1)
plt.xscale('linear')  # Set the x-axis scale to linear
plt.xlabel('Temperature (T)')
plt.ylabel('Average Magnetization per spin')

# Add legend
legend = plt.legend()

power = int(np.log10(total_steps))
# Create the lattice_info string with dynamic steps value
lattice_info = f'Lattice Type: Cubic\nLattice Size: ${N}^3$\nSteps: $10^{power}$'

plt.text(0.78, 0.8, lattice_info, transform=ax.transAxes, fontsize=10, ha='left', va='center')

# Set legend font size
legend.get_texts()[0].set_fontsize('10')

plt.title('Classical 3D Heisenberg Model - Monte Carlo Simulation using Metropolis Algorithm')

# Customize the grid
ax.grid(color='lightgray', linestyle='--', linewidth=0.5)

# Save the figure to a file (e.g., a PNG image)
file_name_image = f"Magnetization_vs_Temperature_Simulation_results_{N}X{N}X{N}_Steps{total_steps}.png" 
plt.savefig(file_name_image, dpi=300, bbox_inches='tight')

plt.show()

In [15]:

# Plot magnetic susceptibility vs temperature
plt.figure(figsize=(9, 5))
plt.plot(T_values, susceptibility_values, 'o-', label='Susceptibility', color='blue', markersize=3, linewidth=0.7)
plt.xscale('linear')  # Set the x-axis scale to linear
plt.xlabel('Temperature (T)')
plt.ylabel('Magnetic Susceptibility')
plt.legend()

lattice_info = f'Lattice Type: Cubic\nLattice Size: ${N}^3$\nSteps: $10^{power}$'

plt.text(0.72, 0.8, lattice_info, transform=ax.transAxes, fontsize=10, ha='left', va='center')

# Set legend font size
legend.get_texts()[0].set_fontsize('10')


plt.title('Classical 3D Heisenberg Model - Monte Carlo Simulation using Metropolis Algorithm')
plt.grid(True)

# Save the figure to a file (e.g., a PNG image)
file_name_image = f"Susceptibility_vs_Temperature_Simulation_results_{N}X{N}X{N}_Steps{total_steps}.png" 
plt.savefig(file_name_image, dpi=300, bbox_inches='tight')

plt.show()

In [ ]: