TP Ising

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] + S[i-1, 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()

L = 32
beta = 0.1
J = 1.00

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

n = 1000000 #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)
plt.show(1)

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)
plt.show(2)

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


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


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

print(datetime.now() - startTime)

#print('Estado final', S)

#%%
#%%
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)

#%%
start = datetime.now()

Temperaturas = np.linspace(0.8,3,num=10)
betas = 1/Temperaturas
n=30
term = 50000
corr=200
measures = 1000

longitudes = []


for t in betas:
    
    for i in range(n):
        S = 2*(np.random.rand(L,L)>0.5) -1
        
        for i in range(term):
            metropolis(S,t)
            
        longmean = []
        for i in range(measures):
            t = 0
            while t <corr:    
                metropolis(S,t)
                t +=1
                
            s = int(np.random.random()*L*L)
            i = int(s/L)
            j = s%L
            mean1 = (S[i,j] * S[i,(j+1)%L])/2.
            mean2 = S[i,j] * S[i,(j+1)%L]
            longmean.append(mean1-mean2)
    longitudes.append(np.average(longmean))            

plt.figure(7)
plt.plot(Temperaturas, longitudes,'g.')        
plt.plot(Temperaturas, longitudes,'y-')
plt.grid(True)
plt.xlabel('Temperaturas')
plt.ylabel('Longitud de Correlación')
plt.show(7)

print(datetime.now()-start)