Ising Baby

ising2018_2.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 = 32
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)


#%%
#LONG DE CORRELACIÓN 
    
def correlation(S):
    i = int(L/2) 
    j = int(L/2)
    #vector con el valor del spin elegido
    c1 = S[i,j]*np.ones(L) 
    #vector con el valor de los spines de la columna (media columna)
    c2 = [S[i,(j+r)%L] for r in range(int(L/2))] 
    #vector con la multiplicación el spin con los demás de la columna
    c3 = [ c1[j]*c2[j] for j in range(int(L/2))]
    return c1, c2 ,c3

def meancorr(S,T,npre,measures,stop):
    beta=1/T
    #pretermalizo la red
    for i in range(npre):
        S=metropolis(S,beta)[0]
    #quiero c1, c2 y c3 para cada medición
    #los voy a guardar en los vectores C mayuscula
    C1=[]
    C2=[]
    C3=[]
    #tomo las mediciones
    for i in range(measures):    
        wait=0
        while wait<stop:
            S=metropolis(S,beta)[0]
            wait+=1
        a,b,c=correlation(S)    
        C1.append(a); C2.append(b); C3.append(c)
    #armo el vector correlación 
    cr=[]
    for r in range(int(L/2)):
        suma1=0
        suma2=0
        suma3=0
        for i in range(measures):
            c1=C1[i]; suma1=suma1+c1[r]
            c2=C2[i]; suma2=suma1+c2[r]
            c3=C3[i]; suma3=suma3+c3[r]
        promS1=suma1/measures
        promS2=suma2/measures
        prom_prod=suma3/measures
        cr.append(prom_prod-(promS1*promS2))
    return cr

#%% GRAFICO VECTOR CR PARA DISTINTAS TEMP
    
t0=datetime.now()
stop=1000
npre = 40000
measures = 600
L=32
J=1.00

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

cr_T=[]

for t in T:
    S = 2*(np.random.rand(L,L)>0.5) -1
    cr_T.append(meancorr(S,t,npre,measures,stop))
    print(datetime.now()-t0)


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

#%% FITEO PARA LONG CORR Y GRAFICO

from scipy.optimize import curve_fit

t0=datetime.now()

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

leng = []

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

print(datetime.now()-t0)  

plt.figure(9)  
plt.plot(T,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)