Hamiltoniano:
$$ H = - \sum_{\langle i, j \rangle} J s_i s_j - \sum_{i}(H(t) + h_i)s_i $$El campo aleatorio se saca de la distribución:
$$ P(h) = \frac{1}{\sqrt{2 \pi} R} e^{-h^2/2R^2} $$
In [1]:
using PyPlot
In [2]:
type MicroEstado
σ::Array{Int,2}
h::Array{Float64,2}
#Vamos a suponer que todas las configuraciones son cuadradas
L::Int
end
import Base.show
show(io::IO, m::MicroEstado) = print(io, m.σ)
import PyPlot.imshow
imshow(m::MicroEstado) = imshow(m.σ, interpolation="none")d
Out[2]:
In [3]:
function edo_inicial(L::Int)
σ = -ones(Int, (L,L))
h = Array(Float64, (L,L))
for i in 1:L^2
# Temporal, hay que cambiar la distribución
h[i] = randn()
end
MicroEstado(σ, h, L)
end
Out[3]:
In [4]:
function energia_espin(m::MicroEstado, i::Int, j::Int)
σ = m.σ ; L = m.L
#Tal vez deberíamos de cambiar el signo de la energía
σ[i,j]*( σ[mod1(i-1, L),j] + σ[mod1(i+1, L),j] + σ[i,mod1(j-1, L)] + σ[i,mod1(j+1, L)] ) +
m.h[i,j]
end
Out[4]:
In [5]:
function voltea_espin!(m::MicroEstado, i::Int, j::Int)
if m.σ[i,j] == -1
m.σ[i,j] *= -1
end
end
Out[5]:
In [25]:
function max_energia_abajo(m::MicroEstado,min_energia)
energias = Array(Float64,(m.L,m.L))
for i = 1:m.L, j = 1:m.L
if m.σ[i,j] == -1
energias[i,j] = energia_espin(m,i,j)
else
energias[i,j] = min_energia
end
end
#Si cambiamos signo de la energía sería findmin
f = findmax(energias) # Da el máximo y el ínidice lineal del máximo
i = mod1(f[2], m.L)
j = int(ceil(f[2]/m.L))
-f[1], i, j
end
Out[25]:
In [7]:
function espines_vecinos_abajo(m::MicroEstado,i::Int,j::Int)
vecinos = [ (mod1(i-1,m.L),j), (mod1(i+1,m.L),j), (i,mod1(j-1,m.L)), (i,mod1(j+1,m.L)) ]
vecinos_abajo = (Int,Int)[]
for k in vecinos
if m.σ[k...] == -1
push!(vecinos_abajo, k)
end
end
vecinos_abajo
end
Out[7]:
In [26]:
function avalancha(m::MicroEstado, i::Int, j::Int, H::Float64)
voltea_espin!(m,i,j)
#puede haber un problema, si se escoje de entrada un espín que ya está volteado
espines_volteados = [(i,j)]
candidatos = espines_vecinos_abajo(m,i,j)
while length(candidatos) > 0
nuevos_candidatos = (Int,Int)[]
for k in candidatos
if energia_espin(m,k...) < H
voltea_espin!(m,k...)
push!(espines_volteados,k)
nuevos_candidatos = vcat(nuevos_candidatos, espines_vecinos_abajo(m,k...))
end
end
candidatos = nuevos_candidatos
end
m, espines_volteados
end
Out[26]:
In [9]:
T = edo_inicial(6)
energia_minima = minimum(T.h)-5
max_energia_abajo(T,energia_minima)
Out[9]:
In [18]:
T.h
Out[18]:
In [19]:
while sum(T.σ) < T.L^2
# figure(figsize=(4,4))
# imshow(T)
H,i,j = max_energia_abajo(T,energia_minima)
@show H-4
T, volteados = avalancha(T,i,j,H)
@show volteados
println(T)
println()
end
In [27]:
magnetizacion(m::MicroEstado) = sum(m.σ)
function magnetizacion_aumenta_H(m::MicroEstado)
energia_minima = minimum(m.h)-5
N = m.L^2
mag = [magnetizacion(m)]
hs = [0.]
while mag[end] < N
H,i,j = max_energia_abajo(m,energia_minima)
push!(hs, H)
m, volteados = avalancha(m,i,j,H)
ΔM = 2*length(volteados)
#Ojo: esto sólo sirve porque volteamos espines hacia arriba
push!(mag, mag[end] + ΔM)
end
mag, hs
end
Out[27]:
In [30]:
L = 50
N = L^2
m = edo_inicial(L)
@time mag, hs = magnetizacion_disminuye_H(m);
In [31]:
figure(figsize=(5,3))
xlabel(L"M/N")
ylabel(L"H/J") #J=1
plot(mag/N,hs);
In [82]:
function microEstados_aumenta_H(m::MicroEstado, num_pasos::Int)
energia_minima = minimum(m.h)-5
edos = Array{Int,2}[copy(m.σ)] #Si no se hace ésto, nada más se copia el apuntador
sizehint(edos, num_pasos)
for i in 1:num_pasos-1
H,i,j = max_energia_abajo(m,energia_minima)
avalancha(m,i,j,H)
push!(edos, copy(m.σ))
end
edos
end
Out[82]:
In [85]:
L = 10
m = edo_inicial(L)
num_pasos = 50
@time edos = microEstados_aumenta_H(m, num_pasos);
In [87]:
for estado in edos
figure(figsize=(4,4))
imshow(estado, interpolation="none")
end
In [ ]: