TP Ising

ising2018.py

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

#%%

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


#%%

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('Estado final', S)

#%%

nterm = 50000

Temperaturas = np.linspace(2,3.0,num=20)
betas = 1/Temperaturas

Magmed = []
Magstd = []
Enermed = []
Enerstd = []

for t in betas:
    Mag = np.zeros(1001)
    Ener = np.zeros(1001)
    S = 2*(np.random.rand(L,L)>0.5) -1
        
    for i in range(nterm):
        S, dE, dM = metropolis(S,t)
    
    Mag[0] = calcMagnet(S)
    Ener[0] = calcEnergia(S)
        
    for i in range(1000):
        S, dE, dM = metropolis(S,t)
        M[i] = M[i-1] + dM
        E[i] = E[i-1] + dE
    Magmed.append(np.average(Mag))
    Enermed.append(np.average(Ener))
    Magstd.append(np.std(Mag))
    Enerstd.append(np.std(Ener))

       
plt.plot(Temperaturas, Magmed,'r.')
    
#np.savetxt('Ising_11_6', np.c_[betas, Magmed, Magstd,Enermed,Enerstd], header='Betas, Magmed, Magenermed, Enerstd', delimiter='')