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:
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 book1 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:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 | 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.
- Computational Physics A book by Konstantinos Anagnostopoulos ↩
