Este manual está elaborado con el fin de proporcionar una guía práctica para la elaboración de un programa Monte Carlo para la simulación de un sistema Ising 2D. El programa está construido en python3. Es cierto que la complejidad del programa es mínima, pero la eficiencia de los ciclos en python es mínima. Es por ello que cabe resaltar que es recomendable emplear otro lenguaje, sea C++ o Fortran, para la implementación del paso Monte Carlo y el algoritmo de Metrópolis.
El sistema que se desea simular es un Ising 2D, cuyo Hamiltoniano consiste únicamente en el término de intercambio
$$\mathcal{H}=-J\sum_{\left\langle i,j\right\rangle }\sigma_{i}\sigma_{j}$$donde $J$ es la constante de intercambio y $\sigma_{i}$ y $\sigma_{j}$ son los espines de la posición $i$ y $j$, respectivamente.
Para calcular la magnetización se emplea
$$M=\sum_{i}\sigma_{i}$$y para el cálculo de la susceptibilidad magnética y el calor específico se sigue $$\chi = \frac{\left(\left\langle M^{2}\right\rangle -\left\langle M\right\rangle ^{2}\right)}{k_{B}T}$$
$$C_{v} = \frac{\left(\left\langle E^{2}\right\rangle -\left\langle E\right\rangle ^{2}\right)}{k_{B}T^{2}}$$
In [19]:
import numpy
from collections import defaultdict
from matplotlib import pyplot
import itertools
%matplotlib inline
In [10]:
length = 10
J = 1.0
kB = 1.0
In [6]:
sites = list()
spins = dict()
nbhs = defaultdict(list)
defaultdict(list), es decir, los values por defecto serán listas vacías. En este diccionario las keys son las parejas $(i, j)$ y los values son una lista de tuplas, donde cada pareja corresponde a un vecino.
In [7]:
for x, y in itertools.product(range(length), range(length)):
sites.append((x,y))
Debemos crear la red cuadrada de lado length y agregar a la lista sites las parejas $(i, j)$ como tuplas. Se puede observar que sites es una lista donde cada elemento es una tupla que corresponde a un sitio:
In [8]:
print(sites)
In [23]:
def random_configuration():
for spin in sites:
spins[spin] = numpy.random.choice([-1, 1])
La función random_configuration recorre todos los sitios y les asigna un valor de 1 o -1 aleatoriamente. Se empleará para darle el estado inicial al sistema a temperatura alta.
Si se ejecuta la función random_configuration, el diccionario spins tomará valores de -1 o 1 para cada key:
In [25]:
random_configuration()
print(spins)
In [28]:
def plot_spins():
pyplot.figure()
colors = {1: "red", -1: "blue"}
for site, spin in spins.items():
x, y = site
pyplot.quiver(x, y, 0, spin, pivot="middle", color=colors[spin])
pyplot.xticks(range(-1,length+1))
pyplot.yticks(range(-1,length+1))
pyplot.gca().set_aspect("equal")
pyplot.grid()
pyplot.show
La función consiste en recorrer cada pareja $(i, j)$ y en cada punto colocar una flecha de color rojo si el sitio tiene valor 1 o azul si tiene valor -1. Se ejecuta la función plot_spins para visualizar el estado del sistema:
In [29]:
plot_spins()
In [38]:
nbhs = defaultdict(list)
for site in spins:
x, y = site
if x + 1 < length:
nbhs[site].append(((x + 1) % length, y))
if x - 1 >= 0:
nbhs[site].append(((x - 1) % length, y))
if y + 1 < length:
nbhs[site].append((x, (y + 1) % length))
if y - 1 >= 0:
nbhs[site].append((x, (y - 1) % length))
Recorremos todos los sitios y agregamos en su lista de nbhs los sitios vecinos. Recordemos que el sistema tiene condiciones periódicas de frontera.
In [40]:
def energy_site(site):
energy = 0.0
for nbh in nbhs[site]:
energy += spins[site] * spins[nbh]
return -J * energy
def total_energy():
energy = 0.0
for site in sites:
energy += energy_site(site)
return 0.5 * energy
def magnetization():
mag = 0.0
for spin in spins.values():
mag += spin
return mag
Si se imprime la magnetización del sistema, debe ser consecuente con la suma de los espines en la figura:
In [42]:
plot_spins()
print("magnetization = ", magnetization())
In [43]:
def metropolis(site, T):
oldSpin = spins[site]
oldEnergy = energy_site(site)
spins[site] *= -1
newEnergy = energy_site(site)
deltaE = newEnergy - oldEnergy
if deltaE <= 0:
pass
else:
if numpy.random.uniform(0, 1) <= numpy.exp(-deltaE/(kB*T)):
pass
else:
spins[site] *= -1
Para mayor información, se puede consultar el algoritmo de Metrópolis en http://hua-zhou.github.io/teaching/st758-2014fall/top10/metropolis.pdf
In [45]:
def monte_carlo_step(T):
for i in range(len(sites)):
int_rand_site = numpy.random.randint(0, len(sites))
rand_site = sites[int_rand_site]
metropolis(rand_site, T)
Un paso Monte Carlo consiste en recorrer la cantidad de sitios que tenga la red aleatoriamente y en cada elección aplicar el algoritmo de Metrópolis.
In [46]:
amount_mcs = 10000
T_high = 5.0
T_low = 0.01
step = -0.1
Se definen los parámetros de la simulación. Entre ellos está la cantidad de pasos Monte Carlo amount_mcs, la temperatura inicial T_high y la temperatura final T_low. Recordemos que $T_{high} \geq T_{low}$. A su vez, se define el paso entre temperaturas step.
In [49]:
%%time
temps = numpy.arange(T_high, T_low, step)
energies = numpy.zeros(shape=(len(temps), amount_mcs))
magnetizations = numpy.zeros(shape=(len(temps), amount_mcs))
random_configuration()
for ind_T, T in enumerate(temps):
for i in range(amount_mcs):
monte_carlo_step(T)
energies[ind_T, i] = total_energy()
magnetizations[ind_T, i] = magnetization()
Se crea el arreglo de temperaturas con los valores establecidos anteriormente. Se aloja espacio para almacenar la energía del sistema y la magnetización. Se recorre las temperaturas y los pasos Monte Carlo. En cada paso de temperatura y en cada paso Monte Carlo se almacena la energía total y la magnetización del sistema. Los arreglos energies y magnetizations se emplearán para calcular los promedios estadísticos.
Una simulación con 10000 pasos Monte Carlo, desde $T_{high}=5.0$ hasta $T_{low}=0.01$ con una paso de $\Delta T = -0.1$ en un computador con Intel® Core™ i7-3612QM CPU @ 2.10GHz tardó 20min 45s.
Podrías observar el estado final de la muestra a $T=0.01$:
In [55]:
plot_spins_spinst_spins()
In [60]:
tau = amount_mcs // 2
energy_mean = numpy.mean(energies[:, tau:], axis=1)
magnetization_mean = abs(numpy.mean(magnetizations[:, tau:], axis=1))
Se dejará la mitad de pasos Monte Carlo para relajación. La otra mitad se ha de promediar.
In [61]:
pyplot.figure()
pyplot.plot(temps, energy_mean, label="Energy")
pyplot.legend()
pyplot.xlabel(r"$T$")
pyplot.ylabel(r"$\left<E\right>$")
pyplot.grid()
pyplot.show()
pyplot.figure()
pyplot.plot(temps, magnetization_mean, label="Magnetization")
pyplot.legend()
pyplot.xlabel(r"$T$")
pyplot.ylabel(r"$\left<M\right>$")
pyplot.grid()
pyplot.show()
In [62]:
magnetization_std = numpy.std(numpy.abs(magnetizations[:, tau:]), axis=1)
susceptibility = magnetization_std ** 2 / (kB * temps)
pyplot.figure()
pyplot.plot(temps, susceptibility, label="Susceptibility")
pyplot.legend()
pyplot.xlabel(r"$T$")
pyplot.ylabel(r"$\chi$")
pyplot.grid()
pyplot.show()
In [63]:
energy_std = numpy.std(energies[:, tau:], axis=1)
specific_heat = energy_std ** 2 / (kB * temps * temps)
pyplot.figure()
pyplot.plot(temps, specific_heat, label="Specific heat")
pyplot.legend()
pyplot.xlabel(r"$T$")
pyplot.ylabel(r"$C_v$")
pyplot.grid()
pyplot.show()
Te invitamos a continuar con tu aprendizaje desde este punto. Podrías intentar obtener curvas de relajación o graficar otro tipo de cantidades como cumulantes de Binder.