In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import itertools as itr
import doctest
In [2]:
N = 10
M = 100
relax = 10*N*N
J = 1
s = np.ones((N,N), dtype=int)
In [3]:
def energy(s, J):
'''
Returns the total energy of s.
>>> energy(np.ones((10,10)), J=1)
-400.0
>>> energy(np.ones((10,10)), J=-1)
400.0
'''
E = 0
for i,j in itr.product(range(N),range(N)):
E += s[i,j]*(s[(i+1)%N,j] + s[(i-1),j] + s[i,(j+1)%N] + s[i,(j-1)])
return -J * E
In [4]:
def flip(s, i, j, beta, J):
'''
Flips spin i,j of matrix s given beta and J.
>>> flip(np.ones((10,10)), 1,1, 0, 1)
True
>>> flip(np.ones((10,10)), 1,1, np.inf, 1)
False
>>> flip(np.ones((10,10)), 1,1, np.inf, -1)
True
>>> flip(np.ones((10,10)), 1,1, 0, -1)
True
'''
E_old = -J * s[i,j]*(s[(i+1)%N,j] + s[(i-1),j] + s[i,(j+1)%N] + s[i,(j-1)])
s[i,j] *= -1
E_new = -J * s[i,j]*(s[(i+1)%N,j] + s[(i-1),j] + s[i,(j+1)%N] + s[i,(j-1)])
delta_E = E_new - E_old
if delta_E > 0:
if np.random.random() > np.exp(-beta*delta_E):
s[i,j] *= -1 # rejecting the move
return False
return True
In [5]:
def run(n, s, beta, J):
'''
Run the algorithm for n times.
'''
for k in range(n):
i = np.random.randint(N)
j = np.random.randint(N)
flip(s, i, j, beta, J)
In [6]:
def cv(E, beta):
return beta**2*E.std()
def jack_knife(a,beta,f):
c = np.zeros_like(a)
for i in range(1,len(a)+1):
idx = list(range(i-1))+list(range(i, len(a)))
c[i-1] = f(a[idx],beta)
c0 = f(a,beta)
return(np.sqrt(((c-c0)**2).sum()))
In [7]:
run(relax, s, 0, J)
beta = 1.5**np.arange(-12,3)
En = np.zeros_like(beta)
En_error = np.zeros_like(beta)
Mag = np.zeros_like(beta)
C = np.zeros_like(beta)
C_error = np.zeros_like(C)
for i in range(len(beta)):
run(relax, s, beta[i], J)
E = np.zeros(M)
mag = np.zeros(M)
for m in range(M):
run(relax, s, beta[i], J)
E[m] = energy(s,J)/(N*N)
mag[m] = s.mean()
En[i] = E.mean()
Mag[i] = mag.mean()
C[i] = cv(E,beta[i])
En_error[i] = E.std()/np.sqrt(M)
C_error[i] = jack_knife(E,beta[i],cv)
In [8]:
ax = plt.subplot(1,1,1)
ax.errorbar(1/beta, En, yerr=En_error)
ax.set_xscale('log')
In [9]:
ax = plt.subplot(1,1,1)
ax.errorbar(1/beta, C, yerr=C_error)
ax.set_xscale('log')
In [10]:
plt.semilogx(1/beta,Mag)
Out[10]:
In [11]:
doctest.testmod()
Out[11]: