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


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

Advances the simulation a single step.

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);

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){
    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,
        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)
            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):
        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)
            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:
        if use_disk:
            # print('Done! ')
            r = run(command, stdout=PIPE)
            raw_data = loadtxt(r.stdout.decode().split(sep='\n'),
            # 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):
    file_name = title.replace(' ', '_')
    if not exists('figures/'):

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

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):
        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.

    ax : pyplot.axis
    what to plot to.

    l : int
    lattice size

    temps: numpy.array
    Array of temperatures where to evaluate heat capacity.

    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[:]
        # 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]))
    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')
        append(linspace(min(temps), t_c, 6), linspace(t_c, max(temps), 6)))
    # 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.

    sim: simulation
    re_runs: int
    number of times to repeat simulation
    take_last: int
    How many of the final points to take statistics over.


    global use_disk
    use_disk = False
    sim.duration = smart_duration(sim.temperature)
    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:]

        last_energies = sim.mean_energies()[-take_last:]

    return [mean(magnetizations), mean(sigma_magnetizations),
            mean(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,
    if critical_temperatures is None:
        critical_temperatures = investigate_heat_capacity(lattice_sizes,
    args, cov = curve_fit(finite_size_scale, lattice_sizes,
    plt.plot(lattice_sizes, critical_temperatures, 'b+', label='data')
    plt.plot(lattice_sizes, finite_size_scale(lattice_sizes, *args), 'r-',
    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,

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

heat cap

and mean magnetisation

  • $\begingroup$ Can you include the $c_v$ vs $T$ plots? And maybe a couple of energy time series for different lattice sizes and same $T$? $\endgroup$
    – lr1985
    Apr 9, 2018 at 21:14

2 Answers 2


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.


I have done this simulation before. I think you can solve your problem by first increasing the number averaging (Monte-Carlo method) and then increasing the size of the lattice. Since when you are doing in small sizes, the model feels the effect of size. You can take a look at the following paper: http://iopscience.iop.org/article/10.1088/1742-5468/2009/07/P07030/pdf

Which talks about the effects of size.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.