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]:
[<matplotlib.lines.Line2D at 0x7fa39dcd47f0>]

In [11]:
doctest.testmod()


Out[11]:
TestResults(failed=0, attempted=6)