Ising Baby

ising2018.py

#%%
import numpy as np
from matplotlib import pyplot as plt

from datetime import datetime


#%%

def calcMagnet(S):
    
    M = np.sum(S)
    return M


def esitio(S,i,j):
    
    c = -J*S[i, j]*(S[i, (j+1)%L] + S[i,(j-1)%L] + S[(i-1)%L, j] + S[(i+1)%L, j])
    return c

def calcEnergia(S):
    
    E_spin = 0
    for i in range(0,L):
        for j in range(0,L):
            E_spin = E_spin + esitio(S,i,j)
    return E_spin/4.00


def newstate(S):
    
    Snew = S.copy()
    i , j = np.random.randint(L), np.random.randint(L)        
    Snew[i,j] = -S[i,j]
    
    dE = esitio(Snew,i,j) - esitio(S,i,j)
    dM = 2*Snew[i,j]
    return Snew, dE, dM

def metropolis(S, beta):
    
    Snew, deltaE, deltaM = newstate(S)
    prob = np.exp(-beta*deltaE)
    ran = np.random.rand()
    
    if deltaE <= 0: #Si disminuye la energía, cambio
        S = Snew
        dE = deltaE
        dM = deltaM
    elif ran <= prob: #Si aumenta la energía, cambio con proba prob
        S = Snew
        dE = deltaE
        dM = deltaM
    else: #Descarto el cambio
        S = S
        dE = 0
        dM = 0
    return S, dE, dM


#%%

startTime = datetime.now()   #para timear las iteraciones

L = 16
beta = 1/1.2
J = 1.00

S = 2*(np.random.rand(L,L)>0.5) -1

n = 100000 #iteraciones

M = np.zeros(n)
E = np.zeros(n)
M[0] = calcMagnet(S)
E[0] = calcEnergia(S)

plt.figure(1)
plt.imshow(S,interpolation='none') #estado inical
plt.show(block=False)

for i in range(1,n):
    S, dE, dM = metropolis(S,beta)
    M[i] = M[i-1] + dM
    E[i] = E[i-1] + dE

plt.figure(2)
plt.imshow(S,interpolation='none') #estado final
plt.show(block=False)

Msitio = [i/(L*L) for i in M]


plt.figure(1) 
plt.plot(Msitio,'.')
plt.ylabel('Magnetización por sitio')    
plt.xlabel('Iteracioines')
plt.show(1)


plt.figure(2) 
plt.plot(E,'.')
plt.ylabel('Energía')    
plt.xlabel('Iteracioines')
plt.show()      

print(datetime.now() - startTime)

#%%
##### Este es el que va 2.0

startTime = datetime.now()   #para timear las iteraciones

L= 16
J=1
nterm = 100000
measures = 10000
correlation = 200    

Temperaturas = np.linspace(0.8,3.5,num=20)
betas = 1/Temperaturas

Magmed = []
Enermed = []
 
for t in betas:

    S = 2*(np.random.rand(L,L)>0.5) -1
    
    Mag = []
#    Ener = []
    
    for i in range(nterm):
        S = metropolis(S,t)[0]
    
    Mag.append(abs(calcMagnet(S)))
#    Ener.append(abs(calcEnergia(S)))
    
    for i in range(measures):
        steps = 0
        while steps < correlation:    
            S = metropolis(S,t)[0]
            steps += 1
            if steps == correlation:
                Mag.append(abs(calcMagnet(S))/(L*L))
#                Ener.append(abs(calcEnergia(S))/(L*L))
                
    Magmed.append(np.average(Mag))
#    Enermed.append(np.average(Ener))
    
Tcrit = 2.2691
T_Tcrit = [ i/Tcrit for i in Temperaturas]
    
plt.plot(Temperaturas, Magmed,'r.')

print(datetime.now() - startTime)

#%%
#%%

#LONGITUD DE CORRELACION
#LONGITUD DE CORRELACION

def correlation(S):
    i = int(L/2)
    j = int(L/2)
    c1 = S[i,j]*np.ones(L)
    c2 = [(S[i,(j+r)%L] for r in range(int(L/2))]
    c3 = [ c1[j]*c2[j] for j in range(int(L/2))]
    return c1, c2 ,c3


def expo(x,e,y,C0):
    return y + C0*np.exp(- x / e)

def meancorr(S,beta, npre, measures, stop):
    for i in range(npre):
        S = metropolis(S,beta)[0]
    
    c1 = np.zeros(int(L/2))
    c2 = np.zeros(int(L/2))
    c3 = np.zeros(int(L/2))
    
    for i in range(measures):
        wait = 0
        while wait < stop:
            S = metropolis(S,beta)[0]
            wait += 1
            if wait == stop:
                c = correlation(S)
                for i in range(int(L/2)):
                    c1[i] += c[0][i]
                    c2[i] += c[1][i]
                    c3[i] += c[2][i]
    cr = [np.abs((c3[i]-c1[i]*c2[i]))/(measures**2) for i in range(int(L/2))]
    return cr

#%%
    
t0 = datetime.now()
stop=1000
npre = 25000
measures = 2000
L =16

temperatures = np.linspace(1.5,3,num=20)
betas = [1/i for i in temperatures]


cr = []

for b in betas:
    S = 2*(np.random.rand(L,L)>0.5) -1
    cr.append(meancorr(S,b, npre, measures, stop))

#
for i in range(len(betas)):
    plt.plot(range(int(L/2)), cr[i], label ='T %s'%temperatures[i])
    plt.xlabel('Distancia <r>',fontsize=12)
    plt.ylabel('Correlación',fontsize=12)
    plt.legend()  


leng = []

for i in range(len(betas)):
    popt, pcov = curve_fit(expo, np.arange(int(L/2)),cr[i])
    leng.append(popt[1])

print(datetime.now() -t0)  
plt.figure(9)  
plt.plot(temperatures,np.abs(leng),'r.')
plt.rc('xtick',labelsize=12)
plt.rc('ytick',labelsize=12)
plt.xlabel('Temperaturas',fontsize=12)
plt.ylabel('Longitud de correlación',fontsize=12)     
plt.grid(True)
plt.grid(color='k', linewidth=.5, linestyle=':') 
plt.show(9)