The problem
I'm trying to make a metropolis simulation of the 2D Ising model.
Basically, it's the following, for each monte-carlo step:
Visit each lattice site,
Compute energy required to flip the spin:
dEIf $exp(-\frac{dE} {kT}) > p$ where
pis a random number in [0,1], flip the spin, else leave it.
Heat capacity I retrieve from the Fluctuation dissipation theorem:
$ C_v = \frac{\sigma_E^2}{k_B T^2}$
My code reproduces some of the expected behaviour: i.e. there's a critical transition at vaguely the right temperature, but:
Energy has a smoothed-out transition. Magnetisation is a sharp step down from 1 to 0, while energy continues to rise, and doesn't show any change of gradient.
The heat capacity shows a peak, but instead of rising when I increase the lattice size, the peak gets smaller.
The heat capacity is noisy: the peak is often doubled, and fitting it to a Lorentzian is a nightmare. Worse yet, running the simulation multiple times gets me the same data (as evidenced by the small errorbar.
My project is a c++ program that simulates a single lattice, and a python script that analyses the output. As of now the script pipes the program stdout and assigns it to internal variables. It's inefficient, but that's for debugging purposes.
The code
The full code can be found here: enter link description here. The code is a c++ program that runs a simulation on a single lattice and a python script that analyses the data.
The code snippets of importance (i.e. that don't contain the irrelevant housekeeping), are
simulation.cpp
#include <iostream>
#include <math.h>
#include "include/simulation.h"
using namespace std;
/*
Advances the simulation a given number of steps, and updates/prints the statistics
into the given file pointer.
Defaults to stdout.
The number of time_steps is explcitly unsigned, so that linters/IDEs remind
the end user of the file that extra care needs to be taken, as well as to allow
advancing the simulation a larger number of times.
*/
void simulation::advance(unsigned int time_steps, FILE *output) {
unsigned int area = spin_lattice_.get_size() * spin_lattice_.get_size();
for (unsigned int i = 0; i < time_steps; i++) {
/*
To simulate an adiabatic Ising lattice, uncomment this
*/
// double temperature_delta = total_energy_/area - mean_energy_;
// if (abs(temperature_delta) < 1/area){
// cerr<<temperature_delta<<"! Reached equilibrium "<<endl;
// }
// temperature_ += temperature_delta;
if (time_ % print_interval_ == 0) {
total_magnetisation_ = spin_lattice_.total_magnetisation();
mean_magnetisation_ = total_magnetisation_ / area;
total_energy_ = compute_energy(spin_lattice_);
mean_energy_ = total_energy_ / area;
print_status(output);
}
advance();
}
}
/*
Advances the simulation a single step.
DOES NOT KEEP TRACK OF STATISTICS. Hence private.
*/
void simulation::advance() {
// #pragma omp parallel for collapse(2)
for (unsigned int row = 0; row < spin_lattice_.get_size(); row++) {
for (unsigned int col = 0; col < spin_lattice_.get_size(); col++) {
double dE = compute_dE(row, col);
double p = r_.random_uniform();
// float rnd = rand() / (RAND_MAX + 1.);
if (dE <0 || exp(-dE / temperature_) > p) {
spin_lattice_.flip(row, col);
}
}
}
time_++;
}
/*
Computes the total energy associated with spins in the spin_lattice_.
I originally used this function to test the code that tracked energy as the lattice
itself was modified, but that code turned out to be only marginally faster, and
not thread-safe. This is due to a race condition: when one thread uses a neighborhood
of a point, while another thread was computing the energy of one such point in
the neighborhood of (row, col).
*/
double simulation::compute_energy(lattice &other) {
double energy_sum = 0;
unsigned int max = other.get_size();
#pragma omp parallel for reduction(+ : energy_sum)
for (unsigned int i = 0; i < max; i++) {
for (unsigned int j = 0; j < max; j++) {
energy_sum += other.compute_point_energy(i, j);
}
}
return energy_sum/2;
}
void simulation::set_to_chequerboard(int step){
if (time_ !=0){
return;
}else{
for (unsigned int i=0; i< spin_lattice_.get_size(); ++i){
for (unsigned int j=0; j<spin_lattice_.get_size(); ++j){
if ((i/step)%2-(j/step)%2==0){
spin_lattice_.flip(i, j);
}
}
}
}
}
## investigator.py
#! /usr/local/bin/python3
from os import mkdir
from os.path import exists
from subprocess import run, PIPE
from numpy import loadtxt, linspace, log, sqrt, append, exp, \
array, shape, inf, std, mean, argmax
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
class Simulation:
def __init__(self, magnetic_field=0.0, exchange_energy=1.0,
lattice_size=64, temperature=1.0, duration=50,
grain_size=0):
self.magnetic_field = magnetic_field
self.exchange_energy = exchange_energy
self.lattice_size = lattice_size
self.temperature = temperature
self.grain_size = grain_size
self.duration = duration
file_path = 'data/' + str(self) + '.csv'
if use_disk:
self.data = self.load_from_disk(file_path)
else:
self.data = self.run()
def times(self):
return self.data[:, 0]
def mean_magnetizations(self):
return self.data[:, 1]
def mean_energies(self):
return self.data[:, 2]
def load_from_disk(self, file_path):
if not exists(path=file_path):
self.run(file_path)
data = loadtxt(file_path)
data_duration = data[-1, 0] + self.duration/200 + 10
if self.duration > data_duration:
print('duration mismatch, ' + str(self.duration) + '!=' + str(
data_duration), end='. re-', flush=True)
self.run(file_path)
data = loadtxt(file_path)
return data
def run(self, file_path=None):
if file_path is None:
file_path = 'data/' + str(self) + '.csv'
# print('running simulation... ', self.duration, end=' ', flush=True)
command = ['./main', '-d', str(self.duration), '-t',
str(self.temperature), '-n', str(self.lattice_size),
'-j', str(self.exchange_energy), '-H',
str(self.magnetic_field), '-c', str(self.grain_size)]
if self.duration > 500:
command.append('-p')
command.append(str(int(self.duration/200)))
if use_disk:
command.append('-f')
command.append(file_path)
# print('Done! ')
run(command)
else:
r = run(command, stdout=PIPE)
raw_data = loadtxt(r.stdout.decode().split(sep='\n'),
delimiter='\t')
# print('Done! ')
return raw_data
def __str__(self):
return '(H=' + str(self.magnetic_field) + ')(J=' + str(
self.exchange_energy) + ')(T=' + str(
self.temperature) + ')(N=' + str(
self.lattice_size) + ')' + str(self.grain_size)
# ------------------------------------------------------------------------
def save_plot(title):
plt.legend(loc='best')
plt.suptitle(title)
file_name = title.replace(' ', '_')
if not exists('figures/'):
mkdir('figures/')
plt.savefig('figures/%s.png'%file_name)
plt.show()
# -------------------------------------------------------------------------
def smart_duration(temperature, multiplier=1.):
return int(((10**5)*base_duration*multiplier)/(
(temperature - t_c)**2 + breadth))
def theoretical_capacity(x, a, b, c, d):
return c + (b*x)/((x - a)**2 + d)
def investigate_heat_capacity(lattice_sizes=None, temps=None, **kwargs):
if temps is None:
temps = linspace(1.5, 3, 10)
if lattice_sizes is None:
lattice_sizes = [16,32, 34, 36]
fig, axes = plt.subplots(nrows=len(lattice_sizes), sharex='col')
crit_temps = []
for l, ax in zip(lattice_sizes, axes):
print('===',l)
crit_temps.append(fit_and_plot_capacity(ax, l, temps, **kwargs))
fig.set_size_inches(10.5, 10.5)
plt.xlabel('temperature / arb. u.')
save_plot('Heat capacity')
return crit_temps
def fit_and_plot_capacity(ax, l, temps, **kwargs):
"""
Plot the heat capacity of simulations at given temperature and lattice size.
Afterwards fit a lorentzian and plot.
Parameters:
-----------
ax : pyplot.axis
what to plot to.
l : int
lattice size
temps: numpy.array
Array of temperatures where to evaluate heat capacity.
Returns:
-------
popt[0]: float
most likely critical temperature.
"""
global use_disk
use_disk = False
simulations = [Simulation(lattice_size=l, temperature=t, duration=2)
for t in temps]
# sigmas = [stdev(s.mean_energies[:]) for s in simulations]
sigmas = array([multi_run(s, **kwargs) for s in simulations]).T[3]
meta_sigmas = array([multi_run(s, **kwargs) for s in simulations]).T[4]
# print(meta_sigmas)
Cs = [sigma**2/(temp**2)*10**3 for temp, sigma in zip(temps, sigmas)]
C_errs = Cs[:]*meta_sigmas[:]
try:
# popt, pcov = curve_fit(theoretical_capacity, temps*10, Cs,
# sigma=meta_sigmas, bounds=(
# [min(temps) - .2, 0, 0, 0],
# [max(temps) + .2, inf, inf, .7]))
pass
except RuntimeError:
print('I\'m too dumb to fit')
# popt = [temps[argmax(Cs)], (max(Cs) - min(Cs))/4, min(Cs),
# (max(temps) - min(temps))/4]
# ax.axvline(x=popt[0], ls='--', color='g')
# ax.plot(temps, theoretical_capacity(temps, *popt), 'g-',
# label='N = ' + str(l) + 'd = ' + str(popt[3]) + ' fit')
ax.errorbar(temps, Cs, fmt='b.', yerr=C_errs, label='N = ' + str(l))
ax.set_ylabel(r'C $\cdot 10^3$/ arb. u.')
ax.axvline(x=t_c, ls='-.', color='k')
ax.set_xticks(
append(linspace(min(temps), t_c, 6), linspace(t_c, max(temps), 6)))
ax.legend(loc='best')
# return popt[0]
def multi_run(sim, re_runs:int=1, take_last:int=300):
"""
Re run the Ising model simulation multiple times, and gather statistics.
Parameters:
----------
sim: simulation
re_runs: int
number of times to repeat simulation
take_last: int
How many of the final points to take statistics over.
Returns:
list:
"""
global use_disk
use_disk = False
sim.duration = smart_duration(sim.temperature)
print(sim.duration)
magnetizations = []
sigma_magnetizations = []
energies = []
sigma_energies = []
for i in range(re_runs):
sim.data = sim.run()
# Make each run take 100% longer, so that we
# can see if a system is still settling
last_magnetizations = sim.mean_magnetizations()[-take_last:]
magnetizations.append(abs(mean(last_magnetizations)))
sigma_magnetizations.append(std(last_magnetizations))
last_energies = sim.mean_energies()[-take_last:]
energies.append(mean(last_energies))
sigma_energies.append(std(last_energies))
return [mean(magnetizations), mean(sigma_magnetizations),
mean(energies), mean(sigma_energies),
std(sigma_energies)/mean(sigma_energies)]
# -------------------------------------------------------------------------
def finite_size_scale(N, t_inf, a, v):
return t_inf + a*(N**(-1/v))
def investigate_finite_size_scaling(critical_temperatures, lattice_sizes,
**kwargs):
if critical_temperatures is None:
critical_temperatures = investigate_heat_capacity(lattice_sizes,
**kwargs)
args, cov = curve_fit(finite_size_scale, lattice_sizes,
critical_temperatures)
plt.plot(lattice_sizes, critical_temperatures, 'b+', label='data')
plt.plot(lattice_sizes, finite_size_scale(lattice_sizes, *args), 'r-',
label='fit')
plt.ylabel('critical temperature / arb. u.')
plt.xlabel('Lattice size')
save_plot('Finite size scaling')
return args[0], sqrt(cov[0, 0])
# -------------------------------------------------------------------------
t_c = 2/log(1 + sqrt(2))
use_disk = False
breadth = 1
base_duration = 10
sizes = [16, 32]
# investigate_time_evolution()
# investigate_temperature_dependence(lattice_sizes=sizes, re_runs=4)
critical_temps = investigate_heat_capacity(lattice_sizes=sizes,
take_last=10000)
# temp_inf = investigate_finite_size_scaling(critical_temps, sizes)
# print(temp_inf)
# print((t_c - temp_inf[0])/ temp_inf[1], ' Standard errors away')
Here's the heat capacity,
and 
