Monte Carlo simulation of 2D Ising model

ising2d-montecarlo.py

import numpy as np

# Set parameters
J = 1.0  # Interaction energy
h = 0.0  # External magnetic field
kT = 1.0  # Temperature times Boltzmann's constant
num_steps = 10000  # Number of steps in the simulation
N = 10  # Size of the lattice, N x N

# Initialize lattice with random spins
lattice = 2*np.random.randint(2, size=(N, N))-1

# Function to calculate the energy of the system
def energy(lattice):
    E = -J * (lattice * np.roll(lattice, 1, axis=0) + lattice * np.roll(lattice, 1, axis=1)).sum() - h * lattice.sum()
    return E

# Function to update the lattice using the Metropolis algorithm
def metropolis_step(lattice):
    for _ in range(N*N):
        x, y = np.random.randint(0, N, 2)  # Choose a random spin
        spin = lattice[x, y]
        dE = 2 * J * spin * (lattice[(x+1)%N, y] + lattice[(x-1)%N, y] + lattice[x, (y+1)%N] + lattice[x, (y-1)%N]) + 2 * h * spin
        if dE <= 0 or np.random.rand() < np.exp(-dE/kT):  # Metropolis condition
            lattice[x, y] = -spin  # Flip the spin
    return lattice

# Run the simulation
for _ in range(num_steps):
    lattice = metropolis_step(lattice)

print(lattice)