Particle Hunter

Monte Carlo Simulation of 1d Ising Model

Introduction

The Ising model is a mathematical model used to describe the behavior of a system of ferromagnetic materials. It was first proposed by physicist Wilhelm Lenz in 1920 and has since been widely used to study phase transitions and critical phenomena in statistical mechanics.

One way to simulate the 1D Ising model in Python is through the use of the Monte Carlo method. This involves generating a large number of random configurations of the system and using these to estimate the thermodynamic properties of the system.

Monte Carlo Simulation

Monte Carlo simulation is a statistical method used to model the behavior of a system by generating a large number of random configurations and using these to estimate the probability of certain outcomes. It is commonly used in many fields to analyze complex systems and make predictions about their behavior.

To illustrate, consider the problem of estimating the value of π using a Monte Carlo simulation in Python. To do this, we can generate a large number of random points within a square and count the number that fall within a quarter circle inscribed within the square. The probability that a given point will fall within the quarter circle is equal to the ratio of the area of the quarter circle to the area of the square, which is equal to π/4. By dividing the number of points that fall within the quarter circle by the total number of points, we can estimate the value of π. This simple example illustrates the basic idea behind Monte Carlo simulation, which is to use random sampling to estimate the probability of certain outcomes

Here is an example of how this could be implemented in Python:

import random

# Number of points to generate

n = 100000

# Counter for points within the quarter circle

count = 0

# Generate the points

for i in range(n):
    x = random.uniform(-1, 1)
    y = random.uniform(-1, 1)
    
    # Check if the point is within the quarter circle

    if x**2 + y**2 <= 1:
        count += 1

# Calculate the estimated value of pi

pi_estimate = 4 * count / n
print(pi_estimate)

The first line imports the random module, which provides functions for generating random numbers. The next line sets the number of points to generate to 100,000. The following line initializes a counter for the number of points that fall within the quarter circle.

The main part of the code is a for loop that iterates over a range of n values. Within the loop, the code generates two random numbers between -1 and 1 using the random.uniform() function and assigns them to the variables x and y. These represent the coordinates of a point within the square.

The code then checks whether the point is within the quarter circle by checking if x**2 + y**2 <= 1. If this condition is true, the counter is incremented. This process is repeated for n iterations, generating n random points and counting the number that fall within the quarter circle.

Finally, the code calculates the estimated value of π by dividing the count by n and multiplying by 4. The result is printed to the console using the print() function.

1D Ising Model Simulation

To begin, we will need to import the necessary libraries and define some key variables. The first thing we will need is a way to store the state of the system. We can do this using a list, where each element represents the spin of a particular particle in the system. For simplicity, we will assume that all particles have the same spin.

import random
import math
import numpy as np

def generate_random_configuration(n, p):

    spins = []
    for i in range(n):
        if random.random() < p:
            spins.append(-1)
        else:
            spins.append(1)
    return spins

Next, we will need to define the energy of each configuration. In the 1D Ising model, the energy of a configuration is given by the sum of the interactions between neighboring particles. For a given configuration, we can calculate the energy by looping through the list of spins and adding up the interactions between each pair of neighboring spins.

Once we have a way to store the state of the system and calculate the energy, we can start using the Monte Carlo method to simulate the system. The basic idea is to generate a large number of random configurations and use these to estimate the thermodynamic properties of the system. To do this, we will need to define a function that generates a random configuration of the system. This can be done using a simple loop that flips the spin of each particle with a certain probability.

def calculate_energy(spins):
    
    n = len(spins)
    energy = 0
    for i in range(n-1):
        energy += -spins[i] * spins[i+1]
    return energy

Next, we will need to define a function that performs a single Monte Carlo step. This involves generating a random configuration, calculating the energy of the configuration, and comparing it to the energy of the current configuration. If the new configuration has a lower energy, we accept it. If the new configuration has a higher energy, we accept it with a probability given by the Boltzmann factor.

def monte_carlo_step(spins, temperature):

    n = len(spins)

    new_spins = spins.copy()
    
    # Choose a random particle

    i = random.randint(0, n-1)
    
    # Calculate the change in energy due to flipping the spin

    delta_energy = 2 * spins[i] * (spins[(i-1)%n] + spins[(i+1)%n])
    
    # Flip the spin with a probability given by the Boltzmann factor

    if delta_energy < 0 or random.random() < math.exp(-delta_energy/temperature):
        new_spins[i] = -spins[i]
    return new_spins

Finally, we will need to define a function that runs the Monte Carlo simulation for a specified number of steps. This can be done by calling the Monte Carlo step function repeatedly until the desired number of steps has been reached.

def monte_carlo_simulation(n, temperature, steps):

    # Generate a random initial configuration

    spins = generate_random_configuration(n, 0.5)
    
    # Perform the Monte Carlo steps

    for i in range(steps):
        spins = monte_carlo_step(spins, temperature)
        
    return spins

Now using the implemented function, we can run the simulation with n=10 particles, temperature=1, and steps=1000 and print the final configuration.

final_configuration = monte_carlo_simulation(10, 1, 1000)
print(final_configuration)

The output will be a list of 10 spins , where each element is 1 or -1. The spins will be arranged in a way that minimizes the energy of the system. This is because the Monte Carlo method is designed to find the configuration of the system that minimizes the energy.

[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

To summerize the code:

Finally, the code calls the monte_carlo_simulation() function with n=10, temperature=1, and steps=1000 to run the simulation with 10 particles, a temperature of 1, and 1000 Monte Carlo steps. It prints the final configuration to the console using the print() function.

Ising Model Animation

To animate the 1D Ising model , we can use a library that provides tools for creating animations, such as matplotlib, mayavi, or plotly. I’m going to use the famous matplotlib to do this. First we need to import the library

import matplotlib.pyplot as plt
import matplotlib.animation as animation

Now we set up the figure and axes

fig, ax = plt.subplots()

After that we need to Initialize the data for the plot using the generate_random_configuration() function defined in the previous section. It creates a scatter plot using the scatter function, with the spins as the y-coordinates and the particle indices as the x-coordinates

n = 10
spins = generate_random_configuration(n, 0.5)
scatter = ax.scatter(range(n), spins, c=spins, cmap='coolwarm')

We then create an update() function, which updates the plot at each frame of the animation. This function performs a single Monte Carlo step using the monte_carlo_step() function and then updates the scatter plot with the new spin values using the set_offsets() and set_array() methods. I was getting AttributeError: 'list' object has no attribute 'ndim' error because the scatter.set_array() method expects a numpy array as its argument, but I was passing it a Python list. This causes the ndim attribute to be accessed on the list object, which does not exist and raises the AttributeError. To fix this, I converted the spins list to a numpy array using the np.array() function before passing it to the set_array() method:

# Function to update the plot at each frame

def update(i):
    global spins
    spins = monte_carlo_step(spins, 1)
    scatter.set_offsets(np.c_[range(n), spins])
    scatter.set_array(np.array(spins))
    return scatter

Finally, we create the animation using the FuncAnimation() function, specifying the figure, the update function, the number of frames, and the blit parameter. It then displays the animation using the show() function.

# Create the animation

anim = animation.FuncAnimation(fig, update, frames=100, blit=False)
plt.show()

PS: If you are using Google colab, it will not show you the animation. The easiest way is to ‘jshtml’ to display matplotlib animation. You need to add the following line of code before the animation code:

from matplotlib import rc
rc('animation', html='jshtml')

Then, just type your animation object. It will display itself, i.e anim in this case.

anim

Alternatively, you can use anim.to_html5_video from IPython.display.HTML using the following code:

from IPython.display import HTML
HTML(anim.to_html5_video())

The outout will be like this:

Magnetization and heat capacity

To calculate the magnetization and heat capacity of the 1D Ising model, we can define two additional functions in our code:

def calculate_magnetization(spins):
    """Calculate the magnetization of a spin configuration."""
    return np.mean(spins)

The calculate_magnetization() function takes a spin configuration (a list or array of 1s and -1s) and returns the mean value of the spins, which is a measure of the overall magnetization of the system.

and for the heat capacity:

def calculate_heat_capacity(spins, temperature):
    """Calculate the heat capacity of a spin configuration at a given temperature."""
    energy = calculate_energy(spins)
    magnetization = calculate_magnetization(spins)
    heat_capacity = (energy**2 + magnetization**2) / temperature**2
    return heat_capacity

The calculate_heat_capacity() function takes a spin configuration and a temperature, and returns the heat capacity of the system. It calculates the energy and magnetization of the system using the calculate_energy() and calculate_magnetization() functions defined earlier, and then uses these values to compute the heat capacity using the formula: C = (E^2 + M^2) / T^2 or heat_capacity = (energy2 + magnetization2) / temperature**2.

We can now use these functions to calculate the magnetization and heat capacity of the system at each step of the animation by calling them in the update() function:

def update(i):
    global spins
    spins = monte_carlo_step(spins, 1)
    magnetization = calculate_magnetization(spins)
    heat_capacity = calculate_heat_capacity(spins, temperature)
    scatter.set_offsets(np.c_[range(n), spins])
    scatter.set_array(np.array(spins))
    return scatter

We can then use these values to update a plot of the magnetization and heat capacity as the animation progresses. You can create these plots using additional axes on the same figure, or you can create separate figures for each plot. We need to write our animation as a function of temperature. This can be done as follows:

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np

def plot_ising_model(n, temperature):

    # Set up the figure and axes

    fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
    ax1.set_title('1D Ising Model')
    ax1.set_ylabel('Magnetization')
    ax2.set_ylabel('Heat Capacity')

    # Initialize the data

    spins = generate_random_configuration(n, 0.5)
    scatter = ax1.scatter(range(n), spins, c=spins, cmap='coolwarm')

    # Set up the magnetization and heat capacity plots

    magnetization_line, = ax1.plot([], [], color='blue')
    heat_capacity_line, = ax2.plot([], [], color='red')

    # Function to update the plots at each frame

    def update(i):
        global spins
        spins = monte_carlo_step(spins, temperature)
        magnetization = calculate_magnetization(spins)
        heat_capacity = calculate_heat_capacity(spins, temperature)
        scatter.set_offsets(np.c_[range(n), spins])
        scatter.set_array(np.array(spins))
        magnetization_line.set_data(range(i), magnetization)
        heat_capacity_line.set_data(range(i), heat_capacity)
        return scatter, magnetization_line, heat_capacity_line

    # Create the animation

    anim = animation.FuncAnimation(fig, update, frames=100, blit=False)
    #plt.show()

The plot_ising_model() function takes the number of spins and the temperature as arguments, and creates a figure with two axes. It then creates a scatter plot of the spin configuration, and two line plots for the magnetization and heat capacity. It then defines an update() function that updates the scatter plot and the two line plots at each frame. Finally, it creates the animation and displays it using the show() function. To use this function, we can call it with different values of the temperature:

# Set the number of particles and the temperature (in units of J/k)

n = 10
temperature = 100

# Plot the 1D Ising model

plot_ising_model(n, temperature)

Conclusion

Ising model is a simple model that can be used to study the phase transition in a system. In this post, we have studied the 1D Ising model using the Monte Carlo method. We have also seen how to create an animation of the system using matplotlib. We have also seen how to calculate the magnetization and heat capacity of the system at each step of the animation. In the next post, we will study the 2D Ising model.

#Monte Carlo #Ising Model #Python