# ————————————————————————————————————————————————————————————————————————————————————————————————
# Import Libraries
# ————————————————————————————————————————————————————————————————————————————————————————————————
import numpy as np
from scipy.optimize import curve_fit
from scipy.ndimage import convolve, generate_binary_structure
import matplotlib as mpl
import matplotlib.pyplot as plt
from ipywidgets import interactive, fixed
import matplotlib.animation
from matplotlib.animation import FuncAnimation
from IPython.core.display import HTML
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
plt.rcParams.update({
'mathtext.fontset':'cm',
'font.family':'serif',
'font.sans-serif':'Times New Roman',
'figure.figsize':[10,8],
'figure.titlesize':22,
"figure.dpi":90,
'savefig.dpi':90,
'axes.titlesize':20,
'axes.labelsize':18,
'axes.titley': 1.0,
'axes.titlepad': 5.0,
'axes.edgecolor':'black',
'axes.grid': True,
'grid.alpha': .5,
'xtick.labelsize':14,'ytick.labelsize':14,
'xtick.major.size':6,'ytick.major.size':6,
'xtick.major.width':1.25, 'ytick.major.width':1.25,
'xtick.direction':'inout','ytick.direction':'inout',
'xtick.top':False, 'ytick.right':False,
'legend.title_fontsize':14, 'legend.fontsize':14,
'legend.borderaxespad': 1, 'legend.borderpad': 0.5,
'legend.framealpha': 1,
'legend.handleheight': 0.5, 'legend.handlelength': 2.0, 'legend.handletextpad':0.5,
'legend.labelspacing': 0.25,
'legend.fancybox':False,
'legend.edgecolor': '0',
'legend.frameon': True,
'legend.markerscale': 1.25,
'animation.embed_limit':2**128,
'animation.html': 'jshtml'})
def totalEnergy(spins, J):
neighbors = generate_binary_structure(2,1)
neighbors[1][1] = False
E = - J * np.sum(spins * convolve(spins, neighbors, mode = 'constant'))
return E
def metropolis(init_spins, steps, J, beta, E):
σ = init_spins.copy()
E_record = np.zeros(steps - 1)
σ_record = np.zeros((steps-1,) + σ.shape)
Nx, Ny = σ.shape
for t in range(steps-1):
i = np.random.randint(0, Nx)
j = np.random.randint(0, Ny)
original_energy = 0
proposed_energy = 0
if i > 0:
original_energy += -σ[i,j] * σ[i - 1, j]
proposed_energy += σ[i,j] * σ[i - 1, j]
if i < (Nx-1):
original_energy += -σ[i,j] * σ[i + 1, j]
proposed_energy += σ[i,j] * σ[i + 1, j]
if j > 0:
original_energy += -σ[i,j] * σ[i, j - 1]
proposed_energy += σ[i,j] * σ[i, j - 1]
if j < (Ny-1):
original_energy += -σ[i,j] * σ[i, j + 1]
proposed_energy += σ[i,j] * σ[i, j + 1]
δE = J * (proposed_energy - original_energy)
φ = np.exp(-beta*δE)
θ = np.random.random()
if (δE > 0) and (θ < φ):
σ[i,j] *= -1
E += J*δE
elif δE <= 0:
σ[i,j] *= -1
E += J*δE
E_record[t] = E
σ_record[t] = σ.copy()
fig, [ax1,ax2] = plt.subplots(1, 2, figsize = (10,5))
cmap = plt.cm.get_cmap('tab20b', 2)
config1 = ax1.imshow(σ_record[0], cmap = cmap)
config2 = ax2.imshow(σ_record[-1], cmap = cmap)
ax2.set_aspect('equal')
ax1.set_aspect('equal')
ax1.axis(False)
ax2.axis(False)
fig.tight_layout()
plt.show()
return σ_record, E_record
up = 0.5
σ0 = np.random.random((100,100))
σ0[σ0 >= up] = 1
σ0[σ0 < up] = -1
β = 1
J = 10
E0 = totalEnergy(σ0, J)
N = 1_000_000
spin_record, net_energy = metropolis(
init_spins = σ0,
steps = N,
J = J,
beta = β,
E = E0
)
%config InlineBackend.figure_format = 'retina'
fig, ax = plt.subplots(figsize = (6,6))
cmap = plt.cm.get_cmap('tab20b', 2)
config = ax.imshow(spin_record[-1], cmap = cmap)
# fig.colorbar(config, ticks=[-1,1], fraction = 0.046, pad = 0.04, label = r'$\sigma$')
ax.set_aspect('equal')
ax.axis(False)
def animate(i):
config.set_array(spin_record[i])
# if i%1000 == 0:
# ax.set_title(r'Algorithm Step $ = {}$'.format(i))
fig.tight_layout()
plt.close()
FPS = 60
N_FRAMES = int(60*5)
anim = FuncAnimation(fig, animate, interval = int(1e3/FPS), frames = np.arange(0, N, int(N/N_FRAMES)), repeat = False)
anim
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
writergif = mpl.animation.PillowWriter(fps = FPS/2)
anim.save('ising.gif', writer = writergif)