Ising Model C++ Metropolis Algorithm

[2kaud]

I'm implying boundary conditions on this part:

Code:

int delta_neighbour = spin(x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);

by using:

Code:

if(x<0) x += LatSize;
if(x>=LatSize) x -= LatSize;
if(y<0) y += LatSize;
if(y>=LatSize) y -= LatSize;

So if you choose a site that's along the sides, it ensures continuity. For example, if Site (0,0) top most left is chosen, the right neighbour would be (0,1), the left neighbour would be (0,4), the top neighbour would be (4,0), the bottom neighbour would be (1,0).

However, when I run the code, it doesn't work, and sometimes it shows absurd values like the left neighbour having value of 25!

No. There is no implied boundary conditions in your code. That is why you are getting absurd values. The boundary conditions need to be part of spin() - so spin() needs to be function instead of just a define.

I've re-worked your code to base it upon a class.

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
using namespace std;

const int LatSize = 5;

class cSpin
{
public:
	int operator()(int x, int y)
	{
		if (x < 0)
			x += LatSize;

		if (x >= LatSize) 
			x -= LatSize;

		if (y < 0) 
			y += LatSize;

		if (y >= LatSize)
			y -= LatSize;

		return matrix[x][y];
	}

	void generate()
	{
		for (int i = 0; i < LatSize; i++) 
			for (int j = 0; j < LatSize; j++)
				matrix[i][j] = (rand() % 2) * 2 - 1;
	}


	void display()
	{
	    for (int i = 0; i < LatSize; i++)  // loop 3 times for three lines
	    {
		for (int j = 0; j < LatSize; j++)  // loop for the three elements on the line
			cout << setw(2) << matrix[i][j] << " ";  // display the current element out of the array

		cout << endl;  // when the inner loop is done, go to a new line
	    }
        }

private:
	int matrix[LatSize][LatSize];
};

int main()
{
cSpin	spin;

	srand (static_cast<unsigned> (time(0)));

	spin.generate();
	spin.display();

int x = int (rand() % LatSize); 
int y = int (rand() % LatSize);

int delta_M = -2 * spin(x, y); //change in magnetization energy

int delta_neighbour = delta_M * (spin(x - 1, y) + spin(x + 1, y) + spin(x, y - 1) + spin(x, y + 1)); 

	cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
	cout << "The left neighbour has spin = " << spin(x, y - 1) << endl;
	cout << "The right neighbour has spin = " << spin(x, y + 1) << endl;
	cout << "The above neighbour has spin = " << spin(x - 1, y) << endl;
	cout << "The below neighbour has spin = " << spin(x + 1, y) << endl;
	cout << "The change in neighbour energy = " << delta_neighbour << endl;

	return 0;
}