Metropolis Segundos Vizinhos

Variação temporal em Python para o Google Colab

!mkdir grafs
!rm grafs/*

import numpy as np
import matplotlib.pyplot as plt
import random
from PIL import Image
import glob
import IPython.display

#Declaração de ctes
l = 100
acoplamento1 = 1
acoplamento2 = 1
campo = 0
temperatura = 0.1

beta = 1/temperatura
tf = 50*l*l
t=0

Eneg = []
Mag = []
Temp = []

Spins = np.ones(l*l)

#Randomização Inicial
for i in range(l*l):
  if random.random() <= 0.25:
    Spins[i] = -Spins[i]

X,Y = np.meshgrid(np.arange(0,l),np.arange(0,l))

#Definição da matriz de primeiros vizinhos
Vizinhos = np.zeros((l*l,4),dtype=int)

contador1 = 1
contador2 = 1
for i in range(l*l):
  if contador1 == 1:
    Vizinhos[i][0] = int(i + l - 1)
    Vizinhos[i][2] = int(i + 1)
  elif contador1 == l:
    Vizinhos[i][0] = int(i - 1)
    Vizinhos[i][2] = int(i - l + 1)
  else:
    Vizinhos[i][0] = int(i - 1)
    Vizinhos[i][2] = int(i + 1)

  if contador2 >= l*l - l:
    Vizinhos[i][3] = int(i + l - l*l)
  else:
    Vizinhos[i][3] = int(i + l)

  Vizinhos[i][1] = int(i - l)

  if contador1 == l:
    contador1 = 1
  else:
    contador1+=1

  contador2 +=1

#Definição da matriz de segundos vizinhos
Vizinhos2 = np.zeros((l*l,4),dtype=int)

for i in range(l*l):
  Vizinhos2[i][0] = Vizinhos[Vizinhos[i][0]][1]
  Vizinhos2[i][1] = Vizinhos[Vizinhos[i][1]][2]
  Vizinhos2[i][2] = Vizinhos[Vizinhos[i][2]][3]
  Vizinhos2[i][3] = Vizinhos[Vizinhos[i][3]][0]


#Definição energia inicial
energia = 0
for i in range(l*l):
  energia += -acoplamento1 * Spins[i] * (Spins[Vizinhos[i][0]] + Spins[Vizinhos[i][1]] + Spins[Vizinhos[i][2]] + Spins[Vizinhos[i][3]]) - acoplamento2 * Spins[i] * (Spins[Vizinhos2[i][0]] + Spins[Vizinhos2[i][1]] + Spins[Vizinhos2[i][2]] + Spins[Vizinhos2[i][3]]) - (campo * Spins[i])

#Loop Temporal
while t <= tf:
  #Sorteio para inversão de Spin
  inverter = random.randint(0,(l*l)-1)

  #Calculo da energia
  deltae = 2 * acoplamento1 * Spins[inverter] * (Spins[Vizinhos[inverter][0]] + Spins[Vizinhos[inverter][1]] + Spins[Vizinhos[inverter][2]] + Spins[Vizinhos[inverter][3]]) + 2 * acoplamento2 * Spins[inverter] * (Spins[Vizinhos2[inverter][0]] + Spins[Vizinhos2[inverter][1]] + Spins[Vizinhos2[inverter][2]] + Spins[Vizinhos2[inverter][3]]) + 2 * (campo * Spins[inverter])

  #Calculo da nova energia

  if deltae <= 0:
    energia += deltae
    Spins[inverter] = -Spins[inverter]
  else:
    if random.random() < np.exp(-beta*deltae):
      energia += deltae
      Spins[inverter] = -Spins[inverter]

  Mag.append(np.sum(Spins)/(l*l))
  Temp.append(1/beta)
  Eneg.append(energia)

  # Criação de figura
  if t % (l*l) == 0:
    fig = plt.figure(figsize=(5,5))
    plt.pcolormesh(X,Y,Spins.reshape([l,l]),cmap='inferno',vmin=-1,vmax=1)
    plt.title('tempo:%.2f'% (t/(l*l)))
    fig.savefig('./grafs/%c.PNG'% int(t/(l*l)+1))
    plt.close()

  t += 1

Mag = np.array(Mag)
Temp = np.array(Temp)
Eneg = np.array(Eneg)

csi = beta * l*l * (np.mean(Mag**2) - np.mean(Mag)**2)
calor = (beta**2) * (np.mean(Eneg**2) - np.mean(Eneg)**2) / (l*l)

print('Susceptibilidade: ', csi)
print('Calor específico: ', calor)

plt.plot(np.arange(0,len(Eneg))/(l*l),Eneg)
plt.xlabel('Tempo')
plt.ylabel('Energia')
plt.show()

plt.plot(np.arange(0,len(Mag))/(l*l),Mag)
# plt.plot((0,tf/(l*l)),(0,campo),label='Campo Induzido')
plt.xlabel('Tempo')
plt.ylabel('Magnetização')
plt.ylim(0,1.1)
# plt.legend()
plt.show()

# Criação do GIF
frames = (Image.open(f) for f in sorted(glob.glob("./grafs/*.PNG")))
frame_one = next(frames)
frame_one.save("animado.gif", format="GIF", append_images=frames, save_all=True, duration=150, loop=0)

IPython.display.Image(open('animado.gif','rb').read())

Variação de temperatura em Python

import numpy as np
import matplotlib.pyplot as plt
import random

#Declaração de ctes
l = 100
acoplamento1 = 1
acoplamento2 = 0
campo = 0

tf = 50*l*l
t=0

Mag = []
Temp = []


X,Y = np.meshgrid(np.arange(0,l),np.arange(0,l))

#Definição da matriz de primeiros vizinhos
Vizinhos = np.zeros((l*l,4),dtype=int)

contador1 = 1
contador2 = 1
for i in range(l*l):
  if contador1 == 1:
    Vizinhos[i][0] = int(i + l - 1)
    Vizinhos[i][2] = int(i + 1)
  elif contador1 == l:
    Vizinhos[i][0] = int(i - 1)
    Vizinhos[i][2] = int(i - l + 1)
  else:
    Vizinhos[i][0] = int(i - 1)
    Vizinhos[i][2] = int(i + 1)

  if contador2 >= l*l - l:
    Vizinhos[i][3] = int(i + l - l*l)
  else:
    Vizinhos[i][3] = int(i + l)

  Vizinhos[i][1] = int(i - l)

  if contador1 == l:
    contador1 = 1
  else:
    contador1+=1

  contador2 +=1

#Definição da matriz de segundos vizinhos
Vizinhos2 = np.zeros((l*l,4),dtype=int)

for i in range(l*l):
  Vizinhos2[i][0] = Vizinhos[Vizinhos[i][0]][1]
  Vizinhos2[i][1] = Vizinhos[Vizinhos[i][1]][2]
  Vizinhos2[i][2] = Vizinhos[Vizinhos[i][2]][3]
  Vizinhos2[i][3] = Vizinhos[Vizinhos[i][3]][0]





#Loop Temperaturas
for temperatura in np.arange(0.001,5,0.1):

  beta = 1/temperatura
  t=0
  Mag_temp = []

  #Condições iniciais
  Spins = np.ones(l*l)
    #Randomização Inicial
  for i in range(l*l):
    if random.random() <= 0.25:
      Spins[i] = -Spins[i]


  #Loop Temporal
  while t <= tf:
    #Sorteio para inversão de Spin
    inverter = random.randint(0,(l*l)-1)

    #Calculo da energia
    deltae = 2 * acoplamento1 * Spins[inverter] * (Spins[Vizinhos[inverter][0]] + Spins[Vizinhos[inverter][1]] + Spins[Vizinhos[inverter][2]] + Spins[Vizinhos[inverter][3]]) + 2 * acoplamento2 * Spins[inverter] * (Spins[Vizinhos2[inverter][0]] + Spins[Vizinhos2[inverter][1]] + Spins[Vizinhos2[inverter][2]] + Spins[Vizinhos2[inverter][3]]) + 2 * (campo * Spins[inverter])

    #Calculo da nova energia
    if deltae <= 0:
      Spins[inverter] = -Spins[inverter]
    else:
      if random.random() < np.exp(-beta*deltae):
        Spins[inverter] = -Spins[inverter]

    Mag_temp.append(np.sum(Spins)/(l*l))

    t += 1
  Mag.append(np.mean(Mag_temp[-100:]))
  Temp.append(temperatura)

plt.plot(Temp,Mag)
plt.xlabel('Temperatura')
plt.ylabel('Magnetização')
plt.ylim(-0.1,1.1)
plt.show()