ancient ising model modernized a bit

ising_updated.py

# original: http://www.physics.rutgers.edu/~haule/681/src_MC/python_codes/ising.py
# modified for python 3 and cleaning up some (but not all) style sins
# and updating the python version with functionality from the C++ version
# and adding a numba JIT-compiled version for speed

"""
Monte Carlo simulation of the 2D Ising model
"""

import numpy as np
import matplotlib.pyplot as plt
import numba  # just to have a compiled version, not necessary at all

Nitt = 1000000  # total number of Monte Carlo steps
N = 10          # linear dimension of the lattice, lattice-size= N x N
warm = 1000     # Number of warmup steps
measure=100     # How often to take a measurement


def CEnergy(latt):
    """Energy of a 2D Ising lattice at particular configuration"""
    Ene = 0
    for i in range(len(latt)):
        for j in range(len(latt)):
            S = latt[i,j]
            WF = latt[(i+1)%N, j] + latt[i,(j+1)%N] + latt[(i-1)%N,j] + latt[i,(j-1)%N]
            Ene += -WF*S # Each neighbor gives energy 1.0
    return Ene/2. # Each par counted twice

# also compile a version for use in the jitted simulator
CEnergy_jitted = numba.njit(CEnergy)


def RandomL(N):
    """Radom lattice, corresponding to infinite temerature"""
    latt = np.zeros((N,N), dtype=int)
    for i in range(N):
        for j in range(N):
            latt[i,j] = np.sign(2*np.random.rand()-1)
    return latt

# update the python implementation with heat capacity and whatnot
def SamplePython(Nitt, latt, PW, T):
    """Monte Carlo sampling for the Ising model in Pythons"""
    Ene = CEnergy(latt)  # Starting energy
    Mn=sum(latt)         # Starting magnetization
    
    Naver=0       # Measurements
    Eaver=0.0
    Maver=0.0
    E2aver = M2aver = 0.0
    
    N2 = N*N
    for itt in range(Nitt):
        t = int(np.random.rand()*N2)
        (i,j) = (t % N, t//N)
        S = latt[i,j]
        WF = latt[(i+1)%N, j] + latt[i,(j+1)%N] + latt[(i-1)%N,j] + latt[i,(j-1)%N]
        P = PW[4+S*WF]
        if P>np.random.rand(): # flip the spin
            latt[i,j] = -S
            Ene += 2*S*WF
            Mn -= 2*S
            
        if itt>warm and itt%measure==0:
            Naver += 1
            Eaver += Ene
            Maver += Mn
            E2aver += Ene ** 2
            M2aver += Mn ** 2

    aM = Maver/Naver
    aE = Eaver/Naver
    cv = (E2aver/Naver - aE**2) / T**2
    chi = (M2aver/Naver - aM**2) / T

    return aM, aE, cv, chi


# let's just skip this for now: update the python implementation instead
# def SampleCPP(Nitt, latt, PW, T):
#     """The same Monte Carlo sampling in C++"""
#     Ene = float(CEnergy(latt))  # Starting energy
#     Mn = float(sum(latt))       # Starting magnetization
# 
#     # Measurements
#     aver = np.zeros(5,dtype=float) # contains: [Naver, Eaver, Maver]
#     
#     code="""
#     using namespace std;
#     int N2 = N*N;
#     for (int itt=0; itt<Nitt; itt++){
#         int t = static_cast<int>(drand48()*N2);
#         int i = t % N;
#         int j = t / N;
#         int S = latt(i,j);
#         int WF = latt((i+1)%N, j) + latt(i,(j+1)%N) + latt((i-1+N)%N,j) + latt(i,(j-1+N)%N);
#         double P = PW(4+S*WF);
#         if (P > drand48()){ // flip the spin
#             latt(i,j) = -S;
#             Ene += 2*S*WF;
#             Mn -= 2*S;
#         }
#         if (itt>warm && itt%measure==0){
#             aver(0) += 1;
#             aver(1) += Ene;
#             aver(2) += Mn;
#             aver(3) += Ene*Ene;
#             aver(4) += Mn*Mn;
#         }
#     }
#     """
#     weave.inline(code, ['Nitt','latt','N','PW','Ene','Mn','warm', 'measure', 'aver'],
#                  type_converters=weave.converters.blitz, compiler = 'gcc')
#     aE = aver[1]/aver[0]
#     aM = aver[2]/aver[0]
#     cv = (aver[3]/aver[0]-(aver[1]/aver[0])**2)/T**2
#     chi = (aver[4]/aver[0]-(aver[2]/aver[0])**2)/T
#     return (aM, aE, cv, chi)

@numba.njit
def SamplePython_jitted(Nitt, latt, PW, T):
    """Monte Carlo sampling for the Ising model in Pythons"""
    Ene = CEnergy_jitted(latt)  # Starting energy
    Mn=latt.sum()         # Starting magnetization
    
    Naver=0       # Measurements
    Eaver=0.0
    Maver=0.0
    E2aver = M2aver = 0.0
    
    N2 = N*N
    for itt in range(Nitt):
        t = int(np.random.rand()*N2)
        (i,j) = (t % N, t//N)
        S = latt[i,j]
        WF = latt[(i+1)%N, j] + latt[i,(j+1)%N] + latt[(i-1)%N,j] + latt[i,(j-1)%N]
        P = PW[4+S*WF]
        if P>np.random.rand(): # flip the spin
            latt[i,j] = -S
            Ene += 2*S*WF
            Mn -= 2*S
            
        if itt>warm and itt%measure==0:
            Naver += 1
            Eaver += Ene
            Maver += Mn
            E2aver += Ene ** 2
            M2aver += Mn ** 2

    aM = Maver/Naver
    aE = Eaver/Naver
    cv = (E2aver/Naver - aE**2) / T**2
    chi = (M2aver/Naver - aM**2) / T

    return aM, aE, cv, chi


if __name__ == '__main__':
    
    latt = RandomL(N)
    PW = np.zeros(9, dtype=float)

    wT = np.linspace(4,0.5,100)
    wMag=[]
    wEne=[]
    wCv=[]
    wChi=[]
    for T in wT:
        
        # Precomputed exponents
        PW[4+4] = np.exp(-4.*2/T)
        PW[4+2] = np.exp(-2.*2/T)
        PW[4+0] = np.exp(0.*2/T)
        PW[4-2] = np.exp( 2.*2/T)
        PW[4-4] = np.exp( 4.*2/T)
    
        #(aM, aE, cv, chi) = SamplePython(Nitt, latt, PW, T)  # <-- native python version
        #(aM, aE, cv, chi) = SampleCPP(Nitt, latt, PW, T)  # <-- removed C++ version using archaic weave module
        (aM, aE, cv, chi) = SamplePython_jitted(Nitt, latt, PW, T)  # <-- jit-compiled python version (compiled by numba)
        wMag.append( aM/(N*N) )
        wEne.append( aE/(N*N) )
        wCv.append( cv/(N*N) )
        wChi.append( chi/(N*N) )
        
        print(T, aM/(N*N), aE/(N*N), cv/(N*N), chi/(N*N))
        
    plt.plot(wT, wEne, label='E(T)')
    plt.plot(wT, wCv, label='cv(T)')
    plt.plot(wT, wMag, label='M(T)')
    plt.xlabel('T')
    plt.legend(loc='best')
    plt.show()
    plt.plot(wT, wChi, label='chi(T)')
    plt.xlabel('T')
    plt.legend(loc='best')
    plt.show()