CodeGuru Home VC++ / MFC / C++ .NET / C# Visual Basic VB Forums Developer.com

# Thread: Ising Model C++ Metropolis Algorithm

1. Junior Member
Join Date
Apr 2014
Posts
24

## Re: Ising Model C++ Metropolis Algorithm

Originally Posted by 2kaud
If you are using Windows 7 why don't you use the free Microsoft Visual Studio Express 2013?
You want Visual Studio Express 2013 for Windows Desktop.

As I suggested in my post #17, before you continue with 'ad hoc coding' of the program, I recommend that you take a step back and design the program before you try to code it.
Code:
`for (i < LatSize) {M+= matrix[i][j];}`
This is not valid c++ code. You also can't initialise a variable when it is defined within a class definiton (unless it is static of an integral type). If you want M to be the sum of the matrix values then you could code it as below. From your original code, B should also defined as 1?

Based upon what I understand from your postings, this is a possible program so far
Code:
```#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
using namespace std;

const int LatSize = 5;

class cSpin
{
public:
cSpin() : M(0)
{
for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
M += matrix[i][j] = (rand() % 2) * 2 - 1;
}

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[y][x];
}

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
}
}

int SpinSum()
{
return M;
}

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

int main()
{
const int J = 1;
const int B = 1;

cSpin	spin;

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

spin.display();

int M = spin.SpinSum();

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

int neighbour = 0;    // nearest neighbour count
int flip_counter = 0;
double sweep_counter = 0.0;

for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
{
neighbour += (spin(i, j) == spin(i + 1, j)) ? 1 : -1;
neighbour += (spin(i, j) == spin(i, j + 1)) ? 1 : -1;
}

int E = -J * neighbour - B * M;

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

int delta_E = -J * delta_neighbour - B * delta_M;

if (delta_E <= 0)
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
}
else
{
//if flipped, M and E is updated. If not flipped, M and E is unchanged.
}

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;
}```
I hope this helps. But you really need to do some design work for the program before you continue coding. Also I would suggest that you become more familar with the c++ language. What book(s) are you using from which to learn c++?
Hmm now when I run the program, the neighbours of (x,y) give the wrong values, I think something is wrong with the boundary conditions. Also, for a 3x3 matrix, the value of M (Sum of all spins) are wrong..

To be honest, this program is nearly complete! Once the sum of all spins, 'M' can be properly calculated and displayed and the flipping decision is done, this program is finished.

At the moment, I'm overloaded with coursework as well but I still want to finish this program and understand it.
Is there any section of C++ you would recommend to understand this program? Don't get me wrong, I find programming to be extremely useful and would definitely study it in greater detail but at the moment I'm quite tied up...

2. ## Re: Ising Model C++ Metropolis Algorithm

the neighbours of (x,y) give the wrong values
Upon looking at the code I put together from your previous postings, there appears to be inconsitencies between which element of the matix is considered to be x and which to be y. I think the code below is consistent now but I haven't tested it.
Code:
```#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
using namespace std;

const int LatSize = 5;

class cSpin
{
public:
cSpin() : M(0)
{
for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
M += matrix[j][i] = (rand() % 2) * 2 - 1;
}

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 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[j][i] << " ";  // display the current element out of the array

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

int SpinSum()
{
return M;
}

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

int main()
{
const int J = 1;
const int B = 1;

cSpin	spin;

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

spin.display();

int M = spin.SpinSum();

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

int neighbour = 0;    // nearest neighbour count
int flip_counter = 0;
double sweep_counter = 0.0;

for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
{
neighbour += (spin(j, i) == spin(j + 1, i)) ? 1 : -1;
neighbour += (spin(j, i) == spin(j, i + 1)) ? 1 : -1;
}

int E = -J * neighbour - B * M;

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

int delta_E = -J * delta_neighbour - B * delta_M;

if (delta_E <= 0)
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
}
else
{
//if flipped, M and E is updated. If not flipped, M and E is unchanged.
}

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;
}```
the value of M (Sum of all spins) are wrong..
Then you need to debug the program to find out why the value is wrong. Where does the code deviate from that expected from the program design?

Is there any section of C++ you would recommend to understand this program?
I would suggest you get a good book on c++. Try
http://www.amazon.co.uk/Ivor-Hortons...+visual+studio Note that there is a newer version coming 28 May.
http://www.amazon.co.uk/C-Primer-Sta...2B+programming
http://www.amazon.co.uk/How-Program-...c%2B%2B+deitel

Also try those web site references I provided in a previous post. You can't really just 'jump into' a section on c++. You need to understand the language from the beginning in a structed way.
Last edited by 2kaud; April 12th, 2014 at 02:38 PM.

3. Junior Member
Join Date
Apr 2014
Posts
24

## Re: Ising Model C++ Metropolis Algorithm

Originally Posted by 2kaud

Then you need to debug the program to find out why the value is wrong. Where does the code deviate from that expected from the program design?
Ok, I have debugged and found the reason why the neighbour values are wrong. Because now the matrix is read differently - x increases towards right, y increases towards bottom, where top left is 0. Previously it was x increasing downwards, y increasing rightwards.

The value of M is also right now.

Now, I need a boolean expression to flip it with acceptance exp(-delta_E) for delta_E > 0. Then, I would need to set the program running continuously for 100 sweeps. That would mean 100*(LatSize*LatSize) random picks.

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

const int LatSize = 3;

class cSpin
{
public:
cSpin() : M(0)
{
for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
M += matrix[j][i] = (rand() % 2) * 2 - 1;
}

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 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[j][i] << " ";  // display the current element out of the array

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

int SpinSum()
{
return M;
}

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

int main()
{
const int J = 1;
const int B = 1;

cSpin	spin;

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

spin.display();

int M = spin.SpinSum();

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

int neighbour = 0;    // nearest neighbour count
int flip_counter = 0;
double sweep_counter = 0.0;

int E = -J * neighbour - B * M;

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

int delta_E = -J * delta_neighbour - B * delta_M;

if (delta_E <= 0)
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
}
else
{
//if flipped, M and E is updated. If not flipped, M and E is unchanged.
}

cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
cout << "The change in neighbour energy = " << delta_neighbour << endl;
cout << "The change in total energy dE =" << delta_E << endl;
cout << "The number of flips = " << flip_counter <<endl;
cout << "The magnetization Energy = " << M <<endl;

return 0;
}```

4. ## Re: Ising Model C++ Metropolis Algorithm

Because now the matrix is read differently - x increases towards right, y increases towards bottom, where top left is 0. Previously it was x increasing downwards, y increasing rightwards.
Yes, as in maths the x-axis should go left to right and the y-axis go top to bottom.

If delta_E > 0, flip with acceptance probability w = exp(-delta_E)
??

That would mean 100*(LatSize*LatSize) random picks.
This comes back to what is meant by a sweep? During one sweep can the same element be chosen more than once or not? If yes, then some cells could get chosen multiple times and some cells not chosen at all. If no then during one sweep only previously unpicked cells can be chosen. This means that for each sweep you need to keep a record of which cells have already been picked in that sweep.

5. Junior Member
Join Date
Apr 2014
Posts
24

## Re: Ising Model C++ Metropolis Algorithm

Originally Posted by 2kaud
Yes, as in maths the x-axis should go left to right and the y-axis go top to bottom.

??

This comes back to what is meant by a sweep? During one sweep can the same element be chosen more than once or not? If yes, then some cells could get chosen multiple times and some cells not chosen at all. If no then during one sweep only previously unpicked cells can be chosen. This means that for each sweep you need to keep a record of which cells have already been picked in that sweep.
If delta_E > 0, it would mean w = exp(-delta_E) would be between 0 and 1. Let's say w = 0.7, so the system has 70% chance of flipping the spin at location (x,y) that was randomly chosen. If delta_E < 0, it would mean w > 1, which means we flip.

Yes, the same site can be chosen more than once. We pick a site at random, and decide whether to flip it or not, depending on value of w. Once I have done this LatSize*LatSize times, we call this one sweep.

I did some research and found this; maybe we do something like this:

Code:
```if (delta_E<=0 || drand48()<exp(-delta_E/T))
{  spin(x, y) *= -1;
M += delta_M;
E += delta_E;
return 1;
}```
Do you happen to know what this code mean? Is this an example of a boolean 'if' function? what does the || in the brackets mean?
Last edited by unscientific; April 12th, 2014 at 06:06 PM.

6. ## Re: Ising Model C++ Metropolis Algorithm

What is the variable T? || in c++ conditional expression means logical OR. ie if either or both of the conditions are met then the condition is evaluated to TRUE. drand48() returns a random number in the range 0.0 to 1.0 and is used with some Unix-like systems but isn't available as standard with Windows. The quick way to get the same result in Windows is to use rand() / double(RAND_MAX). So change
Code:
```        if (delta_E <= 0)
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
}
else
{
//if flipped, M and E is updated. If not flipped, M and E is unchanged.
}```
to
Code:
```        if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < exp((double)-delta_E))
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
}```
You may also need to insert #include <cmath> at the top of the code.

The easiest way to get the 100 sweeps would be to simply use a for loop
Code:
```for (int s = 0; s < 100 * LatSize * LatSize; ++s) {
//code for one calculation goes here
}```
Last edited by 2kaud; April 13th, 2014 at 06:42 AM.

7. Junior Member
Join Date
Apr 2014
Posts
24

## Re: Ising Model C++ Metropolis Algorithm

Originally Posted by 2kaud
What is the variable T? || in c++ conditional expression means logical OR. ie if either or both of the conditions are met then the condition is evaluated to TRUE. drand48() returns a random number in the range 0.0 to 1.0 and is used with some Unix-like systems but isn't available as standard with Windows. The quick way to get the same result in Windows is to use rand() / double(RAND_MAX). So change

Code:
```        if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < exp((double)-delta_E))
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
}```
T is set equal to 1. Does this 'if' function satisfy the boolean flip condition? "If delta_E < 0, flip. If delta_E > 0, the system flips with a w = exp(-delta_E) chance."

Originally Posted by 2kaud
You may also need to insert #include <cmath> at the top of the code.

The easiest way to get the 100 sweeps would be to simply use a for loop
Code:
```for (int s = 0; s < 100 * LatSize * LatSize; ++s) {
//code for one calculation goes here
}```
So this part of the code goes before the 'boolean delta_E' section?

The complete code would be this:

Code:
```#include <math.h>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
using namespace std;

const int LatSize = 3;

class cSpin
{
public:
cSpin() : M(0)
{
for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
M += matrix[j][i] = (rand() % 2) * 2 - 1;
}

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 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[j][i] << " ";  // display the current element out of the array

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

int SpinSum()
{
return M;
}

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

int main()
{
const int J = 1;
const int B = 1;

cSpin	spin;

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

spin.display();

int M = spin.SpinSum();

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

int neighbour = 0;    // nearest neighbour count
int flip_counter = 0; // number of flips count
int sweep_counter = 0; //number of sweeps count

int E = -J * neighbour - B * M;

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

int delta_E = -J * delta_neighbour - B * delta_M;

for (int s = 0; s < 50 * LatSize * LatSize; ++s)
{

if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp(-delta_E))
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
}
else
{
flip_counter++;
sweep_counter = (double)flip_counter / (LatSize * LatSize);     //if flipped, M and E is updated. If not flipped, M and E is unchanged.
}
}
cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
cout << "The change in neighbour energy = " << delta_neighbour << endl;
cout << "The change in total energy dE =" << delta_E << endl;
cout << "The number of flips = " << flip_counter <<endl;
cout << "The number of sweeps = " << sweep_counter <<endl;
cout << "The magnetization Energy = " << M <<endl;

return 0;
}```

8. ## Re: Ising Model C++ Metropolis Algorithm

So this part of the code goes before the 'boolean delta_E' section?
Not really. If you want x and y to selected differently for each loop then the for needs to be before the int x. Do you want neighbour etc to be reset each time around the loop or values covering all the loop iterations?

What parts of the code do you want to be part of the loop and which do you want to executes only once before the loop?

9. Junior Member
Join Date
Apr 2014
Posts
24

## Re: Ising Model C++ Metropolis Algorithm

Originally Posted by 2kaud
Not really. If you want x and y to selected differently for each loop then the for needs to be before the int x. Do you want neighbour etc to be reset each time around the loop or values covering all the loop iterations?

What parts of the code do you want to be part of the loop and which do you want to executes only once before the loop?
Within one loop, these things happen:

1. Pick site (x, y) at random
2. calculate delta_E, decide whether to flip or not ---> If delta_E > 0 : Is it possible to flip only with w = exp(-delta_E) probability? Say if w=0.6, the computer has a 60% chance of flipping it.

When this loop is done LatSize*LatSize times, that is one sweep. We do 100 sweeps. The same site may or may not be chosen again.
Last edited by unscientific; April 13th, 2014 at 06:35 PM.

10. ## Re: Ising Model C++ Metropolis Algorithm

Fine. So just do the quick change to the code so that it matches what you want to do and then you'll have it.

11. Junior Member
Join Date
Apr 2014
Posts
24

## Re: Ising Model C++ Metropolis Algorithm

Originally Posted by 2kaud
Fine. So just do the quick change to the code so that it matches what you want to do and then you'll have it.
Here's the complete code, it seems to work. I have edited the code such that after each sweep, it calculates the magnetic moment, cumulative magnetic moment (sum of magnetic moments till date) and average magnetic moment (cumulative divided by sweep number).

I have seen people store returned values into a .txt file. Is it possible to export (M, M_Sum, M_Avg, sweep_counter) after each sweep into a .txt file so I can plot a graph later?

About the flipping process - I'm not too sure about the flipping process when delta_E > 0. Does the system choose to flip with only a w% chance? (w = exp(-delta_E)) is the probability between 0 and 1.

Code:
```#include <math.h>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
using namespace std;

const int LatSize = 3;

class cSpin
{
public:
cSpin() : M(0)
{
for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
M += matrix[j][i] = (rand() % 2) * 2 - 1;
}

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 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[j][i] << " ";  // display the current element out of the array

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

int SpinSum()
{
return M;
}

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

int main()
{
const int J = 1;
const int B = 1;

cSpin	spin;

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

spin.display();

int M = spin.SpinSum();

int x = int (rand() % LatSize);       //Randomly choose x-coordinate
int y = int (rand() % LatSize);       //Randomly choose y-coordinate

int neighbour = 0;    // nearest neighbour count
int flip_counter = 0; // number of flips count
int sweep_counter = 0; //number of sweeps count
int sweep_number = 5;  // number of sweeps
double M_avg;          // average magnetic moment

int E = -J * neighbour - B * M;   //Total Energy

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

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

int delta_E = -J * delta_neighbour - B * delta_M; //Change in total energy

int M_Sum = 0;   //Cumulative Magnetic Moment starts from 0

for (int s = 0; s < sweep_number; ++s)    //Loop for sweep_number of sweeps
{

for (int t = 0; t < LatSize * LatSize; ++t)   //Loop for 1 Sweep
{

if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp(-delta_E))
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point

}
else
{
flip_counter++;   //Number of Flips

sweep_counter = (double)flip_counter / (LatSize * LatSize);     //Number of Sweeps

}

}
M_Sum += M;                       //Summing up M to find M_Sum
M_avg = (M_Sum)/sweep_counter;    //Average M

cout << "The Magnetic Moment after "<< sweep_counter <<" sweeps = " << M << endl;
cout << "The Cumulative Magnetic Moment after " << sweep_counter <<" sweeps = " << M_Sum << endl;
cout << "The Average Magnetic Moment after " << sweep_counter << "sweeps = " << M_avg <<endl;

}

//cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
//cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
//cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
//cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
//cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
//cout << "The change in neighbour energy = " << delta_neighbour << endl;
//cout << "The change in total energy dE =" << delta_E << endl;
cout << "The number of flips = " << flip_counter <<endl;
cout << "The number of sweeps = " << sweep_counter <<endl;
cout << "The final Magnetic Moment = " << M <<endl;
cout << "The final average Magnetic Moment = " << M_avg <<endl;

return 0;
}```
Last edited by unscientific; April 14th, 2014 at 08:50 AM.

12. ## Re: Ising Model C++ Metropolis Algorithm

it seems to work
Well as long as it seems then that's all right then? Have you tested it to make sure it gives the answers expected? I don't think it will as x and y are only set once at the beginning of the program and are not changed during each loop. Also you aren't calculating the initial value of neighbour. I think the code you want is more like something like this
Code:
```#include <cmath>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
using namespace std;

const int LatSize = 3;

class cSpin
{
public:
cSpin() : M(0)
{
for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
M += matrix[j][i] = (rand() % 2) * 2 - 1;
}

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 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[j][i] << " ";  // display the current element out of the array

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

int SpinSum()
{
return M;
}

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

int main()
{
const int J = 1;
const int B = 1;

cSpin	spin;

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

spin.display();

int M = spin.SpinSum();

int neighbour = 0;    // nearest neighbour count
int flip_counter = 0; // number of flips count
double sweep_counter = 0; //number of sweeps count
int sweep_number = 5;  // number of sweeps
double M_avg;          // average magnetic moment
int M_Sum = 0;   //Cumulative Magnetic Moment starts from 0

for (int s = 0; s < sweep_number; ++s)    //Loop for sweep_number of sweeps
{
for (int t = 0; t < LatSize * LatSize; ++t)   //Loop for 1 Sweep
{
int x = int (rand() % LatSize);
int y = int (rand() % LatSize);

for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
{
neighbour += (spin(j, i) == spin(j + 1, i)) ? 1 : -1;
neighbour += (spin(j, i) == spin(j, i + 1)) ? 1 : -1;
}

int E = -J * neighbour - B * M;   //Total Energy

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

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

int delta_E = -J * delta_neighbour - B * delta_M; //Change in total energy

if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp((double)-delta_E))
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
}
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
}

M_Sum += M;                       //Summing up M to find M_Sum
M_avg = M_Sum / sweep_counter;    //Average M

cout << "The Magnetic Moment after "<< sweep_counter <<" sweeps = " << M << endl;
cout << "The Cumulative Magnetic Moment after " << sweep_counter <<" sweeps = " << M_Sum << endl;
cout << "The Average Magnetic Moment after " << sweep_counter << "sweeps = " << M_avg <<endl;
}

//cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
//cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
//cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
//cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
//cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
//cout << "The change in neighbour energy = " << delta_neighbour << endl;
//cout << "The change in total energy dE =" << delta_E << endl;
cout << "The number of flips = " << flip_counter << endl;
cout << "The number of sweeps = " << sweep_counter << endl;
cout << "The final Magnetic Moment = " << M << endl;
cout << "The final average Magnetic Moment = " << M_avg << endl;

return 0;
}```
but you need to test/debug the code to make sure that it is giving the correct answers.

PS Do you also need to calculate a new value of M eveytime time around the loops as well?
Last edited by 2kaud; April 14th, 2014 at 11:14 AM.

13. Junior Member
Join Date
Apr 2014
Posts
24

## Re: Ising Model C++ Metropolis Algorithm

Originally Posted by 2kaud
Well as long as it seems then that's all right then? Have you tested it to make sure it gives the answers expected? I don't think it will as x and y are only set once at the beginning of the program and are not changed during each loop. Also you aren't calculating the initial value of neighbour. I think the code you want is more like something like this
Code:
```#include <cmath>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
using namespace std;

const int LatSize = 3;

class cSpin
{
public:
cSpin() : M(0)
{
for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
M += matrix[j][i] = (rand() % 2) * 2 - 1;
}

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 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[j][i] << " ";  // display the current element out of the array

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

int SpinSum()
{
return M;
}

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

int main()
{
const int J = 1;
const int B = 1;

cSpin	spin;

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

spin.display();

int M = spin.SpinSum();

int neighbour = 0;    // nearest neighbour count
int flip_counter = 0; // number of flips count
double sweep_counter = 0; //number of sweeps count
int sweep_number = 5;  // number of sweeps
double M_avg;          // average magnetic moment
int M_Sum = 0;   //Cumulative Magnetic Moment starts from 0

for (int s = 0; s < sweep_number; ++s)    //Loop for sweep_number of sweeps
{
for (int t = 0; t < LatSize * LatSize; ++t)   //Loop for 1 Sweep
{
int x = int (rand() % LatSize);
int y = int (rand() % LatSize);

for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
{
neighbour += (spin(j, i) == spin(j + 1, i)) ? 1 : -1;
neighbour += (spin(j, i) == spin(j, i + 1)) ? 1 : -1;
}

int E = -J * neighbour - B * M;   //Total Energy

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

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

int delta_E = -J * delta_neighbour - B * delta_M; //Change in total energy

if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp((double)-delta_E))
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
}
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
}

M_Sum += M;                       //Summing up M to find M_Sum
M_avg = M_Sum / sweep_counter;    //Average M

cout << "The Magnetic Moment after "<< sweep_counter <<" sweeps = " << M << endl;
cout << "The Cumulative Magnetic Moment after " << sweep_counter <<" sweeps = " << M_Sum << endl;
cout << "The Average Magnetic Moment after " << sweep_counter << "sweeps = " << M_avg <<endl;
}

//cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
//cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
//cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
//cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
//cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
//cout << "The change in neighbour energy = " << delta_neighbour << endl;
//cout << "The change in total energy dE =" << delta_E << endl;
cout << "The number of flips = " << flip_counter << endl;
cout << "The number of sweeps = " << sweep_counter << endl;
cout << "The final Magnetic Moment = " << M << endl;
cout << "The final average Magnetic Moment = " << M_avg << endl;

return 0;
}```
but you need to test/debug the code to make sure that it is giving the correct answers.

PS Do you also need to calculate a new value of M eveytime time around the loops as well?
Yeah I'm supposed to find the cumulative magnetic moment after each sweep and take that value, divided by sweep number to get average magnetic moment. It should converge to some value.

I'm more bothered by this section of the code. What does this actually mean? :

Code:
```	if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp(-delta_E))
else {}```
Here's the full updated code, with random site picks (x, y) put into each loop for 1 sweep. How do I store values of "M, M_Sum, M_Avg, sweep_counter" into a .txt file or something like that so I can plot a graph later?

Code:
```#include <math.h>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
using namespace std;

const int LatSize = 30;

class cSpin
{
public:
cSpin() : M(0)
{
for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
M += matrix[j][i] = (rand() % 2) * 2 - 1;
}

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 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[j][i] << " ";  // display the current element out of the array

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

int SpinSum()
{
return M;
}

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

int main()
{
const int J = 1;
const int B = 1;

cSpin	spin;

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

spin.display();

int M = spin.SpinSum();

int neighbour = 0;    // nearest neighbour count
int flip_counter = 0; // number of flips count
int sweep_counter = 0; //number of sweeps count
int sweep_number = 520;  // number of sweeps
double M_avg;          // average magnetic moment

int M_Sum = 0;   //Cumulative Magnetic Moment starts from 0

for (int s = 0; s < sweep_number; ++s)    //Loop for sweep_number of sweeps
{

for (int t = 0; t < LatSize * LatSize; ++t)   //Loop for 1 Sweep
{

int x = int (rand() % LatSize);       //Randomly choose x-coordinate
int y = int (rand() % LatSize);       //Randomly choose y-coordinate

int delta_M = -2 * spin(x, y);    //change in Magnetic Moment
int delta_neighbour = delta_M * (spin(x - 1, y) + spin(x + 1, y) + spin(x, y - 1) + spin(x, y + 1));
int delta_E = -J * delta_neighbour - B * delta_M; //Change in total energy

int E = -J * neighbour - B * M;   //Total Energy

if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp(-delta_E))
{
spin(x, y) *= -1;											// flip spin
M += delta_M;												// New Magnetization energy of flipped configuration
E += delta_E;												// New total energy of flipped configuration
flip_counter++;												// Total number of sites chosen up till this point
sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point

}
else
{
flip_counter++;   //Number of Flips

sweep_counter = (double)flip_counter / (LatSize * LatSize);     //Number of Sweeps

}

}
M_Sum += M;                       //Summing up M to find M_Sum
M_avg = (M_Sum)/sweep_counter;    //Average M

cout << "The Magnetic Moment after "<< sweep_counter <<" sweeps = " << M << endl;
cout << "The Cumulative Magnetic Moment after " << sweep_counter <<" sweeps = " << M_Sum << endl;
cout << "The Average Magnetic Moment after " << sweep_counter << "sweeps = " << M_avg <<endl;

}

//cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
//cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
//cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
//cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
//cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
//cout << "The change in neighbour energy = " << delta_neighbour << endl;
//cout << "The change in total energy dE =" << delta_E << endl;
cout << "The number of flips = " << flip_counter <<endl;
cout << "The number of sweeps = " << sweep_counter <<endl;
cout << "The final Magnetic Moment = " << M <<endl;
cout << "The final average Magnetic Moment = " << M_avg <<endl;

return 0;
}```

14. ## Re: Ising Model C++ Metropolis Algorithm

To store data in a file, you first need to open the file for writing and then output data to the file. Data can be written to the file using stream insertion (just like << for cout).
See http://www.cplusplus.com/reference/ostream/ostream/

For c++ stream i/o you will need to have #include <fstream> at the beginning of the program. To open a file for output use
Code:
```ofstream outs("spins.txt");

if (!outs.is_open()) {
cout << "Problem opening output file\n";
return 1;
}```
where spins.txt is the name of the output file to be created. This code would typically be placed near the beginning of the main() function.

To write to the file when required use
Code:
`outs << "data to write to the file\n";`
you can use the same method of writing to the file as with cout.

15. ## Re: Ising Model C++ Metropolis Algorithm

Code:
`if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp((double)-delta_E))`
the if condition is a logical or. If either or both of the seperate conditons (left and right of the ||) is true then the if condition is true.
The first condition is easy. For the second, (double) means cast (convert) the result of what follows to be of type double. As rand() and RAND_MAX are both integers, they are first cast to double so that the result is of type double because in c/c++ an integer divided by an integer gives an integer (3/4 = 0). So the expression to the left of < gives a number of type double in the range 0 to 1 (for the probability required). The expression to the right of the < is what you stated was needed. The (double) casts are needed again to make sure that the types are as required.

See http://msdn.microsoft.com/en-us/libr...=vs.60%29.aspx
Last edited by 2kaud; April 15th, 2014 at 01:57 PM.

#### Posting Permissions

• You may not post new threads
• You may not post replies
• You may not post attachments
• You may not edit your posts
•