ising2d-reach-equilibrium.py

import matplotlib.pyplot as plt

# Additional lists to store properties at each step
magnetizations = []
energies = []
heat_capacities = []

# Run the simulation
for i in range(num_steps):
    lattice = metropolis_step(lattice)
    
    # Calculate and store magnetization
    magnetization = np.sum(lattice)
    magnetizations.append(magnetization)
    
    # Calculate and store energy
    E = energy(lattice)
    energies.append(E)
    
    # Calculate and store heat capacity if not the first step
    if i > 0:
        heat_capacity = (energies[i] - energies[i-1]) / kT
        heat_capacities.append(heat_capacity)
    else:
        heat_capacities.append(0)
        
# Plot magnetization
plt.figure()
plt.plot(magnetizations)
plt.title('Magnetization')

# Plot energy
plt.figure()
plt.plot(energies)
plt.title('Energy')

# Plot heat capacity
plt.figure()
plt.plot(heat_capacities)
plt.title('Heat Capacity')

plt.show()