# 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:

1. Visit each lattice site,

2. Compute energy required to flip the spin: dE

3. If $exp(-\frac{dE} {kT}) > p$ where p is 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);
}
}
}

/*
Advances the simulation a single step.

DOES NOT KEEP TRACK OF STATISTICS. Hence private.
*/
// #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:
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]

if not exists(path=file_path):
self.run(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)
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)
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: 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
meta_sigmas = array([multi_run(s, **kwargs) for s in simulations]).T
# 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, ls='--', color='g')
# ax.plot(temps, theoretical_capacity(temps, *popt), 'g-',
# label='N = ' + str(l) + 'd = ' + str(popt) + ' 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

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, sqrt(cov[0, 0])

# -------------------------------------------------------------------------

t_c = 2/log(1 + sqrt(2))
use_disk = False
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)/ temp_inf, ' Standard errors away')


Here's the heat capacity, • Can you include the $c_v$ vs $T$ plots? And maybe a couple of energy time series for different lattice sizes and same $T$? – lr1985 Apr 9 '18 at 21:14

It looks like you are computing the specific heat by using the fluctuations of the energy per site rather than the total energy of the system. If you multiply the y axes of your plots by $N^2$ you should recover the expected behaviour.
By the way, you should not take into account "self energies" (that is, $e(i,i)$) in the compute_energy method.