import matplotlib.pyplot as plt
import matplotlib.animation as ani
import numpy as np
from math import *
from scipy.optimize import curve_fit
from numpy.random import default_rng
import copy
[docs]def print_progress(step, total):
"""
Prints a progress bar.
Args:
step (int): progress counter
total (int): counter at completion
"""
message = "simulation progress ["
total_bar_length = 60
percentage = int(step / total * 100)
bar_fill = int(step / total * total_bar_length)
for i in range(total_bar_length):
if i < bar_fill:
message += "|"
else:
message += " "
message += "] "+str(percentage)+" %"
if step < total:
print(message, end="\r")
else:
print(message)
[docs]def show_grid(grid):
"""
Draws the system as a :math:`N \\times N` pixel image.
Args:
grid (array): the system as an integer array
"""
plt.clf()
ax = plt.axes()
ax.set_aspect('equal')
plt.pcolormesh(grid, cmap='Greys', vmin=-1, vmax=1)
plt.show()
[docs]def save_grid_image(grid, index):
"""
Saves an image of the system as a :math:`N \\times N` pixel image.
The image is saved in a file named "ising_X.png" where
X is the given index.
Args:
grid (array): the system as an integer array
index (int): index for naming the output file
"""
plt.clf()
ax = plt.axes()
ax.set_aspect('equal')
plt.pcolormesh(grid, cmap='Greys')
plt.savefig('ising_'+str(i)+'.png', bbox_inches='tight', vmin=-1, vmax=1)
plt.close()
[docs]def draw(frame, history):
"""
Draws the system for animation.
Args:
frame (int): index of the frame to draw
history (list): list of systems at different times
"""
plt.clf()
ax = plt.axes()
ax.set_aspect('equal')
plt.pcolormesh(history[frame], cmap='Greys', vmin=-1, vmax=1)
[docs]def animate(history):
"""
Animate the simulation.
Args:
history (list): list of systems at different times
"""
nframes = len(history)
print("animating "+str(nframes)+" frames")
fig = plt.figure()
motion = ani.FuncAnimation(fig, draw, nframes, fargs=( history, ) )
plt.show()
#print("calling ffmpeg")
#filename = "L2_ising_T.mp4"
#motion.save(filename,fps=30,codec='mpeg4')
#print("animation saved as "+filename)
#plt.clf()
[docs]def generate_grid(lattice_size, randomly = False):
"""
Creates the system as :math:`N \\times N` array.
Each value in the array represents a spin with two possible values, up and down.
The value up is represented by the value 1
and the value down by -1.
Args:
lattice_size (int): system size :math:`N`
randomly (bool): If True, spins will be randomly oriented. Otherwise all spins point up.
"""
if randomly:
grid = 2*random.integers(low=0,high=2,size=(lattice_size, lattice_size))-1
else:
grid = np.ones( [ lattice_size, lattice_size] )
return grid
[docs]def calculate_flip_energy(grid, i,j):
"""
Calculates how much the total energy of the system would change,
if the spin :math:`(i,j)` was flipped.
The system is not changed in any way. The function only calculates
how much the energy would change, if a spin was flipped.
The total energy of the 2D Ising model is
.. math::
E = - \\sum_{a,b} J_{a,b} s_{a} s_{b},
where :math:`a = (i,j)` and :math:`b = (k,l)` are spin
coordinates, :math:`s_a, s_b` are the spin values at these
coordinates (either 1 or -1), and :math:`J_{a,b}` is the
coupling parameter between these two spins.
* If :math:`J > 0`, the spins interact ferromagnetically, i.e.,
they prefer to point in the same direction.
* If :math:`J < 0`, the spins interact antiferromagnetically
and prefer to point in opposite directions.
* If :math:`J = 0`, there is no interaction between the spins.
In our model, :math:`J = 1` for all spins which are directly
next to each other (vertical and horizontal neighbors),
and :math:`J = 0` for all other spins (in some suitable units).
This means that if only one spin, :math:`s_(i,j)` changes its sign,
there are exactly four terms in the energy sum that change:
the terms with :math:`s_{i,j}` and its neighbors up, down, right and left,
i.e., :math:`s_{i+1,j}, s_{i-1,j}, s_{i,j+1}, s_{i,j-1}`.
Furthermore they change by changing their sign.
For instance, if :math:`s_{i,j,\\text{before}} = 1` and :math:`s_{i+1,j} = -1`,
initially the energy of this pair is
:math:`-J_{(i,j),(i+1,j)} s_{i,j,\\text{before}} s_{i+1,j} = 1`.
After the flip, we would have :math:`s_{i,j,\\text{after}} = -1`,
and the energy of the pair would become -1.
The *change* in energy would therefore be
.. math ::
\\Delta E_1 = 2 J_{(i,j),(i+1,j)} s_{i,j,\\text{before}} s_{i+1,j} = -2.
The energies for the three other affected pairs are calculated similarly,
so substituting :math:`J = 1` we get the total energy change
.. math::
\\Delta E = 2s_{i,j,\\text{before}} ( s_{i+1,j} + s_{i-1,j} + s_{i,j+1} + s_{i,j-1}).
Args:
grid (array): the system as an integer array
i (int): row of the spin to flip
j (int): column of the spin to flip
Returns:
int: energy change :math:`\\Delta E`
"""
energy = 0.0
lattice_size = len(grid[0])
energy += 2*(grid[i,j]*grid[(i-1)%lattice_size,j])
energy += 2*(grid[i,j]*grid[(i+1)%lattice_size,j])
energy += 2*(grid[i,j]*grid[i,(j-1)%lattice_size])
energy += 2*(grid[i,j]*grid[i,(j+1)%lattice_size])
return energy
[docs]def flip(grid, i,j):
"""
Flips the spin :math:`s_{i,j}` in the system.
Since spin values up and down are represented by values of +1 and -1,
respectively, the function simply changes the sign of the given
element in the array.
Args:
grid (array): the system as an integer array
i (int): row of the spin to flip
j (int): column of the spin to flip
"""
grid[i,j] = -grid[i,j]
[docs]def metropolis_step(grid, temperature):
"""
Tries to flip one randomly chosen spin using the Metropolis algorithm.
The flip is either accepted or discarded according to a Boltzmann-like
probability:
* If the spin flip leads to *lower* total energy, it is always accepted.
* If the spin flip leads to *higher* total energy, the new configuration
is accepted with probability :math:`P = e^{- \\frac{ \Delta E }{k_B T} }`,
where :math:`\\Delta E` is the energy change, :math:`k_B` is the Boltzmann
constant and :math:`T` is the temperature.
If the spin flip is accepted, the change is applied to the given system.
The function uses a unitless temperature scale with :math:`k_B = 1`.
.. note ::
This function is incomplete!
Args:
grid (array): the system as an integer array
temperature (float): temperature
"""
# pick a random spin
lattice_size = len(grid[0])
coords = random.integers(low=0, high=lattice_size, size=2)
i = coords[0]
j = coords[1]
# if the Metropolis test succeeds, the spin is flipped
randomizer = random.random()
# todo
[docs]def calculate_magnetic_moment(grid):
"""
Calculates the total magnetization of the system.
A magnetic moment is associated with electron spins,
:math:`\\vec{\\mu} = -g\\frac{e}{2 m_e} \\vec{s}`.
In the ising model, spins only point either up or down
and get values of +1 or -1, and so
also the magnetic moment points only up or down.
Furthermore, by choosing a unitless scale :math:`g\\frac{e}{2 m_e} = 1`,
the magnetic moment becomes the same as the spin (only reversed), and
we can calculate the total magnetization as a sum over
all spins in the system,
.. math::
M = -\\frac{1}{N^2} \\sum_{i,j} s_{i,j}.
Here :math:`N \\times N` is the number of spins in the system.
Thus :math:`M` is scaled so that it only gets values between
-1 and 1. If :math:`|M| = 1`, the magnetization is complete and
all spins point in the same direction. If :math:`M = 0`, there
are equally many spins pointing up and down, and there is no
macroscopic magnetic moment.
Args:
grid (array): the system as an integer array
Returns:
float: magnetization absolute value :math:`|M|`
"""
magnetic_moment = 0
lattice_size = len(grid[0])
for i in range(lattice_size):
for j in range(lattice_size):
magnetic_moment += grid[i,j]
return np.abs(magnetic_moment)/(lattice_size**2)
[docs]def calculate_Binder_cumulant(magnetizations):
"""
Calculates the Binder cumulant.
The transition from an ordered, ferromagnetic state to a disordered,
paramagnetic state should be discontinuous at the thermodynamic limit,
but in a finite system, the transition is gradual. This makes it difficult
to pinpoint the critical temperature :math:`T_C` from a finite simulation.
This can be done very precisely using the cumulant method.
The quantity
.. math ::
U_4 = 1 - \\frac{ \\langle M^4 \\rangle }{ 3 \\langle M^2 \\rangle^2 },
known as the Binder cumulant, is also dependent on system size.
At temperatures below :math:`T_C`, it approaches 0 as the system size grows.
At temperatures above :math:`T_C`, it approaches 2/3 for the Ising model.
At :math:`T_C`, it has a non-trivial system dependent value, which is approximately
independent of system size.
Therefore one can pinpoint the transition temperature by calculating the
cumulant as a function of temperature at different system sizes
and plotting the results. The cumulants will cross approximately at :math:`T_C`.
Args:
magnetizations (array): sequence of magnetization values
Returns:
float: the cumulant :math:`U_4`
"""
m4 = 0.0
m2 = 0.0
n = 0
for m in magnetizations:
m4 += m**4
m2 += m**2
n += 1
m4 /= n
m2 /= n
return 1 - m4/(3*m2**2)
[docs]def write_data_to_file(lattice_size, temps, mags):
"""
Writes a data series in a file.
The data is written in file "ising_data_N.txt" where
N is the lattice size.
The data is written so that each column holds
a temperature and the calculated magnetization at that temperature.
Args:
lattice_size (int): system size :math:`N`
temps (list): a list of temperatures
mags (list): a list of magnetizations
"""
writelines = ""
for t, m in zip(temps, mags):
writelines += str(t)+" "+str(m)+" \n"
f = file('ising_data_'+str(lattice_size)+".txt",'w')
f.write(writelines)
f.close()
[docs]def calculate_statistics(samples):
"""
Calculates and returns the sample mean, variance,
standard deviation and error of mean for the given set of samples.
For a sequence of :math:`N` values for a random variable :math:`X`,
the mean is estimated by the average
.. math ::
\\mu_X = \\langle X \\rangle = \\frac{1}{N} \\sum_i X_i
and the variance by the sample variance
.. math ::
\\sigma^2_X = \\frac{1}{N-1} \\sum_i (X_i - \\mu_X)^2.
Standard deviation is the square root of variance
.. math ::
\\sigma_X = \\sqrt{ \\sigma^2_X }
and the standard error of mean is
.. math ::
\\Delta X = \\frac{ \\sigma_X }{ \\sqrt{ N } }.
Args:
samples (list): sequence of random results, :math:`X_i`
Returns:
float, float, float, float: mean, variance, standard deviation, error of mean
"""
expect = 0.0
variance = 0.0
expect = 0.0
variance = 0.0
n_repeats = len(samples)
for i in samples:
expect += i
expect /= float(n_repeats)
for i in samples:
variance += (i-expect)*(i-expect)
variance /= float(n_repeats-1)
std = np.array( sqrt(variance) )
error = std / sqrt(n_repeats)
return expect, variance, std, error
[docs]def run_cumulant_series(n_steps = 10000, thermalization_steps = 1000, lattice_sizes = [5, 10, 15],
min_temperature = 1.0, max_temperature = 3.0, t_steps = 11):
"""
Runs a series of simulations at different temperatures and lattice sizes.
The function calculates Binder cumulants for different systems and plots these.
If the statistics are good enough, these cumulants should intersect approximately
at the critical point :math:`T_C`.
Args:
n_steps (int): total simulation cycles (each cycle is :math:`N^2` spin flip attempts)
thermalization_steps (int): thermalization cycles in the beginning of the simulation
lattice_sizes (list): list of lattice sizes :math:`N`
min_temperature (float): the inital temperature
max_temperature (float): the final temperature (assumed larger than min_temperature)
t_steps (int): the number of different temperatures to simulate
"""
cumulants = []
for lattice_size in lattice_sizes:
print("simulating system N = "+str(lattice_size))
grid = generate_grid(lattice_size)
ts, ms, es, cs = run_temperature_series(grid,
min_temperature,
max_temperature,
t_steps,
n_steps,
thermalization_steps)
cumulants.append(cs)
for i in range( len(sizes) ):
plt.plot(ts,cumulants[i],'-',label="N = "+str(sizes[i]))
plt.legend()
plt.show()
[docs]def run_temperature_series(grid, min_temperature, max_temperature,
temp_steps, run_steps, thermalization_steps):
"""
Runs a series of simulations at different temperatures.
The function calculates the average magnetization and its error estimate (as
the standard error of mean) at each temperature and returns these as arrays.
Args:
grid (array): the system as an integer array
min_temperature (float): the inital temperature
max_temperature (float): the final temperature (assumed larger than min_temperature)
temp_steps (int): the number of different temperatures to simulate
run_steps (int): total simulation cycles (each cycle is :math:`N^2` spin flip attempts)
thermalization_steps (int): thermalization cycles in the beginning of the simulation
Returns:
array, array, array, array: temperatures, magnetizations, error estimates, cumulants
"""
# lists for storing measurements
m_averages = []
m_errors = []
cumulants = []
# list of temperatures to simulate
temperatures = np.linspace(min_temperature, max_temperature, temp_steps)
# Run a separate simulation at each temperature.
# Use the final state of the previous simulation
# as a starting point for the next one.
for T in temperatures:
# run simulation
print("starting simulation at T = "+str(T))
magnetization = run_metropolis_algorithm(grid, T, run_steps, thermalization_steps)
# calculate average magnatization and its error
avr, var, std, err = calculate_statistics(magnetization)
m_averages.append(avr)
m_errors.append(err)
# calculate cumulant
c = calculate_Binder_cumulant(magnetization)
cumulants.append(c)
return temperatures, np.array(m_averages), np.array(m_errors), np.array(cumulants)
[docs]def run_metropolis_algorithm(grid, temperature, run_steps, thermalization_steps, recording_interval = 10, animated=False):
"""
Runs the Metropolis algorithm at a contant temperature.
The algorithm repeatedly chooses random spins in the system
and tries to flip them using :meth:`metropolis_step`.
The changes are applied directly on the given grid.
After a thermalization period, the algorithm also calculates
the magnetization of the system.
By default, the recording is not done at every step because too
closely measured values will be correlated.
Args:
grid (array): the system as an integer array
temperature (float): the temperature
run_steps (int): total simulation cycles (each cycle is :math:`N^2` spin flip attempts)
thermalization_steps (int): thermalization cycles in the beginning of the simulation
recording_interval (int): number of cycles between data recording
animated (bool): if True, shows an animation of the simulation in the end
Returns:
array: sequence of magnetizations
"""
magmoms = []
history = []
n_spins = len(grid)**2 # total number of spins
# run the simulation for the required time
for i in range(run_steps):
print_progress(i,run_steps)
# for each cycle, try to flip spins as many times
# as there are spins in the system
for j in range(n_spins):
metropolis_step(grid, temperature)
# only record data after thermalization is complete
if i >= thermalization_steps and i%recording_interval == 0:
moment = calculate_magnetic_moment(grid)
magmoms.append( abs(moment) )
if animated:
history.append(copy.deepcopy(grid))
print_progress(run_steps, run_steps)
if animated:
animate(history)
return np.array(magmoms)
[docs]def main():
"""
Main program.
"""
lattice_size = 50
n_steps = 1000
thermalization_steps = 500
T = 2.20
grid = generate_grid(lattice_size, randomly=False)
mag = run_metropolis_algorithm(grid, T, n_steps, thermalization_steps, recording_interval=10, animated=True)
# calculate statistics
m_avr, m_var, m_std, m_err = calculate_statistics(mag)
# plot magnetization and the calculated average +- confidence interval
plt.plot(mag)
plt.plot(np.linspace(0, len(mag), len(mag)), m_avr*np.ones(len(mag)), 'b:')
plt.fill_between(np.linspace(0, len(mag), len(mag)), (m_avr-2*m_err)*np.ones(len(mag)), (m_avr+2*m_err)*np.ones(len(mag)), color = 'b', alpha=0.1 )
plt.xlabel("simulation step")
plt.ylabel("magnetization")
plt.show()
if __name__ == "__main__":
random = default_rng()
main()
else:
random = default_rng()