I'm very new to python so my code may be messy. I was trying to do this physics homework, using the Ising model, where an array of 16x16 is set with a spin value (-1 or 1) for each cell. Then, using mcmove(config,iT), the equilibrium situation is found for a given temperature, using calcMag(config). The program also computes the total magnetization for each situation.
My problem is that I have achieved the animation of the different arrays (using matplotlib.animation.ArtistAnimation), but I would like to animate the magnetization too. My ultimate goal would be to have the six subplots animated and with every frame syncronized: the array for each temperature and its magnetization value at the same frame for the three temperatures.
Here's my code (the previous definitions can be skipped, I think my problem is inside the loop when plotting the magnetization):
import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt
import matplotlib.animation as animation
def initialstate(N):
'''Generates a random spin configuration for initial conditions'''
state=2*np.random.randint(2,size=(N,N))-1
return state
def mcmove(config,beta):
'''Monte Carlo move using Metropolis algorithm'''
for i in range(N):
for j in range(N):
#select random spin for NxN system
a=np.random.randint(0,N)
b=np.random.randint(0,N)
s=config[a,b]
#calculate energy cost of this new configuration (the % is for calculation of periodic boundary conditions)
nb=config[(a+1)%N,b] + config[a,(b+1)%N]+config[(a-1)%N,b]+config[a,(b-1)%N]
cost = 2*s*nb
#flip spin or not, depending on the cost and its Boltzmann factor
## (acceptance probability is given by Boltzmann factor with beta=1/kBT)
if cost<0:
s*=-1
elif rand() < np.exp(-cost*beta):
s*=-1
config[a,b]=s
return config
def calcEnergy(config):
'''Energy of a given configuration'''
energy = 0
for i in range(len(config)):
for j in range(len(config)):
S = config[i,j]
nb=config[(i+1)%N,j]+config[i,(j+1)%N] + config[(i-1)%N,j] + config[i,(j-1)%N]
energy += -nb*S
return energy/4.
def calcMag(config):
'''Magnetization of a given configuration'''
mag = np.sum(config)
return mag
nt = 64 # number of temperature points
N = 16 # size of the lattice, NxX
eqSteps = 10 # number of MC sweeps for equilibration
n1= 1.0/(eqSteps*N*N)
Magnetization = np.zeros(eqSteps)
n=0
tem='Temperature = '
c=0
f=0
id=1
images1 = []
images2 = []
images3 = []
#images4 = []
#images5 = []
#images6 = []
fig = plt.figure()
ax1 = fig.add_subplot(2,3,1)
ax2 = fig.add_subplot(2,3,2)
ax3 = fig.add_subplot(2,3,3)
ax4 = fig.add_subplot(2,3,4)
ax5 = fig.add_subplot(2,3,5)
ax6 = fig.add_subplot(2,3,6)
list=[[0.5, ax1, images1, ax4],
[1.5, ax2, images2, ax5],
[3.0, ax3, images3, ax6]]
for T, c, p, g in list:
E1 = M1 = E2 = M2 = 0
config = initialstate(N)
iT = 1.0/T
iT2 = iT*iT
print ("Simulation at temperature T=",T)
titles = tem + str(T)
for i in range(eqSteps): # equilibrate
mcmove(config,iT) # Monte Carlo moves
im = c.imshow(config, animated=True)
c.set_title(titles)
p.append([im])
Mag = calcMag(config)
M1 = M1 + Mag
Magnetization[i] = n1*M1
g.plot(i, abs(Magnetization[i]),'.',color="#348ABD")
g.set_title('Magnetization')
ani = animation.ArtistAnimation(fig, images1, interval=50, blit=True, repeat_delay=0)
ani = animation.ArtistAnimation(fig, images2, interval=50, blit=True, repeat_delay=0)
ani = animation.ArtistAnimation(fig, images3, interval=50, blit=True, repeat_delay=0)
THANK YOU SO MUCH in advance.
About more trivial things I would also like some help in:
- Can I remove the axis from the arrays?
- Can I have the same axis for the magnetization plots for every frame?
- How do I put xlabels and ylabels on the magnetization plots?
- Can I put a label in the figure that updates with every frame? I want this for enumarating them