The Ising Model

Introduction

This is an intermediate post between the one on the Monte Carlo methods and one presenting a Monte Carlo C++ program I intend to write. My goal is to briefly expose the theory here – most of it with links – and provide a very easy JavaScript example for the Metropolis algorithm applied on the 2D Ising model.

The model

The Ising model is one of the simplest models that have a non trivial behavior and it’s very important because of the universality. For this post and the next one, I’ll consider a special case, the 2D Ising model on a square lattice. I even drop the position dependency of the magnetic field/coupling and the directional dependency of the interaction strength, so the Hamiltonian will be:

$H=-J\sum\limits_{(i,j)}s_i s_j - B\sum\limits_i s_i$

J>0 for ferromagnetic interaction and J<0 for antiferromagnetic interaction, the sum over (i, j) is summing all adjacent pairs.
There isn’t much to say about it besides what’s already in the Wikipedia pages, it’s a quite simple model. The book¹ I pointed out in last post also contains quite a bit of information on it so I won’t bother to present more theory here.

Metropolis Monte Carlo

Please take a look into my last post for the theory, here it is applied on a 2D square-lattice Ising model with periodic boundary conditions:

It’s JavaScript code I quickly wrote just to illustrate the Metropolis algorithm. Here is the code:

var monteCarlo = (function() {
 var canvas = document.getElementById("monteCarloCanvas");
 var ctx = canvas.getContext("2d");

 var isingSpins = {
    Size: 64,
    Temperature: 2.26918531421/0.95, // somewhere near the critical temperature
    spins: [], // the spin matrix
    displaySize: canvas.width / 64, // the size of the square for a spin

    // just to make the code more clear, a function to access the spin in the lattice
    Spin: function(row, col) { return this.spins[this.Size*row + col]; },

    // adjusts the 'index', taking care of periodic boundary conditions
    Index: function(index) { return index < 0 ? index + this.Size : index % this.Size; },

    // this will also take care of the periodic boundary conditions, use this one instead of 'Spin' for neighbors
    Neighbor: function(row, col) { return this.Spin(this.Index(row), this.Index(col)); },

    NeighborContribution: function(row, col) { return this.Neighbor(row - 1, col) + this.Neighbor(row + 1, col) + this.Neighbor(row, col - 1) + this.Neighbor(row, col + 1); },

    ExpMinusBetaE: function(E) { return Math.exp(-1. / this.Temperature * E); },

    EnergyDifForFlip: function(row, col) { return 2 * this.Spin(row, col) * this.NeighborContribution(row, col); },

    // initialize the matrix with random spins, this corresponds to infinite temperature
    Init: function() {
        var nr = this.Size * this.Size;
        for (var i = 0; i < nr; ++i) this.spins.push(Math.random() < 0.5 ? -1 : 1);
    },

    // the sweep over the lattice
    Sweep: function() {
        var nr = this.Size * this.Size; // number of spins in the lattice
        
        // for all spins on the lattice
        for (var i = 0; i < nr; ++i) {
            // pick a spin at random
            var row = Math.floor (Math.random() * this.Size);
            var col = Math.floor (Math.random() * this.Size);
            
            // what is the change in energy if the spin flips? 
            var energyDif = this.EnergyDifForFlip(row, col);
            if (energyDif < 0)
                // accept the spin flip
                this.spins[this.Size*row+col] *= -1;
            else
            {
                 if (Math.random() < this.ExpMinusBetaE(energyDif))
                     // accept the spin flip
                     this.spins[this.Size*row+col] *= -1;
                 // else reject it, that is, do nothing
            }
         }
    },

    Display: function(ctx) {
        for (var i = 0; i < this.Size; ++i)
            for (var j = 0; j < this.Size; ++j)
            {
               if (this.Spin(i, j) < 0) ctx.fillStyle = "#FF0000";
               else ctx.fillStyle = "#0000FF";
               ctx.fillRect(i * this.displaySize, j * this.displaySize, this.displaySize, this.displaySize);
            }
    }
 };

 isingSpins.Init();
   
 return function () { 
    // call a monte carlo sweep
    isingSpins.Sweep();
    // then display
    isingSpins.Display(ctx);
 }
})();

setInterval(monteCarlo, 100);

var monteCarlo = (function() {

var canvas = document.getElementById("monteCarloCanvas");

var ctx = canvas.getContext("2d");

var isingSpins = {

Size: 64,

Temperature: 2.26918531421/0.95, // somewhere near the critical temperature

spins: [], // the spin matrix

displaySize: canvas.width / 64, // the size of the square for a spin

// just to make the code more clear, a function to access the spin in the lattice

Spin: function(row, col) { return this.spins[this.Size*row + col]; },

// adjusts the 'index', taking care of periodic boundary conditions

Index: function(index) { return index < 0 ? index + this.Size : index % this.Size; },

// this will also take care of the periodic boundary conditions, use this one instead of 'Spin' for neighbors

Neighbor: function(row, col) { return this.Spin(this.Index(row), this.Index(col)); },

NeighborContribution: function(row, col) { return this.Neighbor(row - 1, col) + this.Neighbor(row + 1, col) + this.Neighbor(row, col - 1) + this.Neighbor(row, col + 1); },

ExpMinusBetaE: function(E) { return Math.exp(-1. / this.Temperature * E); },

EnergyDifForFlip: function(row, col) { return 2 * this.Spin(row, col) * this.NeighborContribution(row, col); },

// initialize the matrix with random spins, this corresponds to infinite temperature

Init: function() {

var nr = this.Size * this.Size;

for (var i = 0; i < nr; ++i) this.spins.push(Math.random() < 0.5 ? -1 : 1);

// the sweep over the lattice

Sweep: function() {

var nr = this.Size * this.Size; // number of spins in the lattice

// for all spins on the lattice

for (var i = 0; i < nr; ++i) {

// pick a spin at random

var row = Math.floor (Math.random() * this.Size);

var col = Math.floor (Math.random() * this.Size);

// what is the change in energy if the spin flips?

var energyDif = this.EnergyDifForFlip(row, col);

if (energyDif < 0)

// accept the spin flip

this.spins[this.Size*row+col] *= -1;

else

{

if (Math.random() < this.ExpMinusBetaE(energyDif))

// accept the spin flip

this.spins[this.Size*row+col] *= -1;

// else reject it, that is, do nothing

}

Display: function(ctx) {

for (var i = 0; i < this.Size; ++i)

for (var j = 0; j < this.Size; ++j)

{

if (this.Spin(i, j) < 0) ctx.fillStyle = "#FF0000";

else ctx.fillStyle = "#0000FF";

ctx.fillRect(i * this.displaySize, j * this.displaySize, this.displaySize, this.displaySize);

}

};

isingSpins.Init();

return function () {

// call a monte carlo sweep

isingSpins.Sweep();

// then display

isingSpins.Display(ctx);

}

})();

setInterval(monteCarlo, 100);

The code is pretty straightforward, with the help of the comments and the last post it shouldn’t be hard to understand. I tried to make it easy to understand instead of optimizing it, for example one could pre-calculate the exponentials for energy difference for a spin flip at a specific temperature, but probably the code would not be so clear so I did not bother.

Conclusion

This is just a quick post before having a C++ program implemented. It’s Sunday, so I didn’t want to spend too much time on it. Hopefully it can be useful even as it is.

If you find any mistakes, please point them out.
CodeProject

Computational Physics A book by Konstantinos Anagnostopoulos ↩

Introduction

The model

Metropolis Monte Carlo

Conclusion

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply

Introduction

The model

Metropolis Monte Carlo

Conclusion

Share this:

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply