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 Ising Model C++ Metropolis Algorithm

I'm writing a code in C++ for a 2D Ising model. Here's what the code should do:

Generate random NxN lattice, with each site either +1 or -1 value.
Select a site at random
If site when flipped (+1 to -1 or -1 to +1) is a state of lower energy, flip state ie. if dE < 0, flip state. If flipped state is of higher energy, flip with acceptance rate w = e^{-b(dE)}. Where dE is the change in energy if state is flipped. 4.Do this for all NxN sites, without repetition. This is considered one sweep.
Do like 100 sweeps.
I'm having trouble with steps 2 and 3, would appreciate any help!

Code:
#include <cstdlib>
#include <ctime>
using namespace std;
#include <iostream>
int main() //random generation of spin configuration
{
int L;              //Total number of spins L = NxN
double B=1;         //magnetic field
double M;           //Total Magnetization = Sum Si
double E;           //Total Energy
int T = 1.0;
int nsweeps = 100;      //number of sweeps
int de;             //change in energy when flipped
double Boltzmann;       //Boltzmann factor
int x,y;            //randomly chosen lattice site
int i,j,a,c;            //counters
int ROWS = 5;
int COLS = 5;
int matrix[ROWS][COLS];
srand ( static_cast<unsigned> ( time ( 0 ) ) );
for ( int i = 0; i < ROWS; i++ )
{
for ( int j = 0; j < COLS; j++ )
{
matrix[i][j] = rand () % 2 *2-1;
}
}

// showing the matrix on the screen
for(int i=0;i<ROWS;i++)  // loop 3 times for three lines
{
for(int j=0;j<COLS;j++)  // loop for the three elements on the line
{
cout<<matrix[i][j];  // display the current element out of the array
}
cout<<endl;  // when the inner loop is done, go to a new line
}
return 0;  // return 0 to the OS.

//boundary conditions and range
if(x<0) x += N;
if(x>=L) x -= N;
if(y<0) y += N;
if(y>=L) y -= N;

//counting total energy of configuration
{  int neighbour = 0;    // nearest neighbour count

for(int i=0; i<L; i++)
for(int j=0; j<L; j++)
{  if(spin(i,j)==spin(i+1, j))     // count from each spin to the right and above
neighbour++;
else
neighbour--;
if(spin(i, j)==spin(i, j+1))
neighbour++;
else
neighbour--;
}

E = -J*neighbour - B*M;

//flipping spin
int x = int(drand48()*L);   //retrieves spin from randomly choosen site
int y = int(drand48()*L);

int delta_M = -2*(x, y);    //calculate change in Magnetization M
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
int delta_neighbour = -2*(x,y)* int delta_neighbour;

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

//flip or not
if (delta_E<=0)
{  (x, y) *= -1;     // flip spin and update values
M += delta_M;
E += delta_E;

}

}  Reply With Quote

2. Re: Ising Model C++ Metropolis Algorithm

I'm not sure I fully understand what you are trying to accomplish. However, looking at the code I can offer some observations which may be of help.

With reference to your NxN lattice. As the number of columns will always equal the number of rows, you don't need to specifiy seperately the number of rows and columns. I would suggest something like
Code:
const int LatSize = 5;
then L can be defined (after LatSize is defined) as
Code:
const int L = LatSize * LatSize;
as according to your comment it is supposed to be NxN.

You then don't need the variables ROWS and COLS and matrix would be defined as
Code:
int matrix[LatSize][LatSize];
and wherever in the code you use ROWS or COLUMNS these would be replaced with LatSize.

I take it that the return 0 after the code to show the matix is just for testing purposes?

You then do some testing with x and y. But x and y haven't been initialised. You use 'N' but this hasn't been defined. Assuming that 'N' is supposed to be the size of the lattice, then these would be LatSize.

If x and y are supposed to be the chosen lattice site, then they need to be set to a value in the range >=0 and < LatSize. As you have already seeded the random number generator and used it, you could use
Code:
x = rand() % LatSize;
y = rand() % LatSize;
However, looking at the following code, it seems to me that this has been taken from somewhere else as it doesn't match the definitions used at the start - ie it uses the function spin but spin isn't defined. If you want to use the code past your return statement, you will need to know what L, N, spin etc are supposed to be. I would suggest you go back to where you got this code from and see how its variables are defined and how they are supposed to be initialised.
Last edited by 2kaud; April 9th, 2014 at 05:17 PM.  Reply With Quote

3. Junior Member Join Date
Apr 2014
Posts
24

Re: Ising Model C++ Metropolis Algorithm Originally Posted by 2kaud I'm not sure I fully understand what you are trying to accomplish. However, looking at the code I can offer some observations which may be of help.
Thanks alot for your reply. I realize now my code is too messy, so I shall work on individual parts first. Let's start with:

1. Generate NxN lattice with -1 or +1 at each site randomly.
2. How to extract the value at a site at location (x, y).

//Generate NxN lattice with -1 or +1 at each site randomly.

const int LatSize = 5;
const int L = LatSize * LatSize;
int matrix[LatSize][LatSize];

int matrix[LatSize][LatSize];
srand ( static_cast<unsigned> ( time ( 0 ) ) );
for ( int i = 0; i < ROWS; i++ )
{
for ( int j = 0; j < COLS; j++ )
{
matrix[i][j] = rand () % 2 *2-1;
}
}

//End of Generating matrix  Reply With Quote

4. Re: Ising Model C++ Metropolis Algorithm

Please use code tags when you post code as you did in post #1.

As ROWS and COLS are now not defined, you could use
Code:
const int LatSize = 5;
int matrix[LatSize][LatSize];

srand(static_cast<unsigned> (time(0)));
for ( int i = 0; i < LatSize; i++)
for ( int j = 0; j < LatSize; j++)
matrix[i][j] = (rand () % 2) * 2 - 1;
Then to select a site of the lattice at random you could use
Code:
const int x = rand() % LatSize;
const int y = rand() % LatSize;  Reply With Quote

5. Junior Member Join Date
Apr 2014
Posts
24

Re: Ising Model C++ Metropolis Algorithm

Now for part 2. How to extract the value of a site at location (x, y).

//Extracting value of a site at location (x,y)

int & spin (int x, int y);
int x = int(srand48()*L);
int y = int(srand48()*L);

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

int delta_neighbour = spin(spinx-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
delta_neighbour = -2*spin(x,y)* delta_neighbour; //change in neighbour product energy

//End of extracting value of site at (x, y)

Are there any problems with this section of the code? For example Does "int H = spin(x, y);" give H the value of the site at location (x, y)?  Reply With Quote

6. Junior Member Join Date
Apr 2014
Posts
24

Re: Ising Model C++ Metropolis Algorithm Originally Posted by 2kaud Please use code tags when you post code as you did in post #1.

As ROWS and COLS are now not defined, you could use
Code:
const int LatSize = 5;
int matrix[LatSize][LatSize];

srand(static_cast<unsigned> (time(0)));
for ( int i = 0; i < LatSize; i++)
for ( int j = 0; j < LatSize; j++)
matrix[i][j] = (rand () % 2) * 2 - 1;
Then to select a site of the lattice at random you could use
Code:
const int x = rand() % LatSize;
const int y = rand() % LatSize;
oops. Here's the updated portion of code:
Code:
//Generate NxN lattice with -1 or +1 at each site randomly.

const int LatSize = 5;
const int L = LatSize * LatSize;
int matrix[LatSize][LatSize];

int matrix[LatSize][LatSize];
srand ( static_cast<unsigned> ( time ( 0 ) ) );
for ( int i = 0; i < LatSize; i++ )
{
for ( int j = 0; j < LatSize; j++ )
{
matrix[i][j] = rand () % 2 *2-1;
}
}

//End of Generating matrix  Reply With Quote

7. Junior Member Join Date
Apr 2014
Posts
24

Re: Ising Model C++ Metropolis Algorithm

For part 2 (Updated!)

[code]

//Extracting value of a site at location (x,y)

int & spin (int x, int y); //extracts value of site at location (x, y) when 'spin' is used in later parts of the code
int x = int(srand48()*L);
int y = int(srand48()*L);

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

int delta_neighbour = spin(spinx-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
delta_neighbour = -2*spin(x,y)* delta_neighbour; //change in neighbour product energy

//End of extracting value of site at (x, y)

Are there any problems with this section of the code? For example Does "int H = spin(x, y);" give H the value of the site at location (x, y)?  Reply With Quote

8. Re: Ising Model C++ Metropolis Algorithm

As I explained in my previous post #2. What is the function spin? This is being used but isn't defined.  Reply With Quote

9. Junior Member Join Date
Apr 2014
Posts
24

Re: Ising Model C++ Metropolis Algorithm Originally Posted by 2kaud As I explained in my previous post #2. What is the function spin? This is being used but isn't defined.
How do I make spin(x, y) return the value of the site at location (x, y)?

For example if I type:
Code:
int H = 2*spin(x, y);
it means that H is twice the value at site (x, y), either -2 or +2.  Reply With Quote

10. Re: Ising Model C++ Metropolis Algorithm

Possibly the easiest way is to replace spin(...) by matrix[...].

Another way would be to use a define at the top of the code
Code:
#define spin(x,y) (matrix[(x)][(y)])
Last edited by 2kaud; April 10th, 2014 at 10:38 AM.  Reply With Quote

11. Junior Member Join Date
Apr 2014
Posts
24

Re: Ising Model C++ Metropolis Algorithm Originally Posted by 2kaud Possibly the easiest way is to replace spin(...) by matrix[...].

Another way would be to use a define at the top of the code
Code:
#define spin(x,y) (matrix[(x), (y)])
Should this go in the second part or first part of the code?
Because in first section of the code to generate NxN lattice I already defined matrix to be [LatSize][LatSize]:

Code:
const int LatSize = 5;
const int L = LatSize * LatSize;
int matrix[LatSize][LatSize];

int matrix[LatSize][LatSize];
srand ( static_cast<unsigned> ( time ( 0 ) ) );
for ( int i = 0; i < LatSize; i++ )
{
for ( int j = 0; j < LatSize; j++ )
{
matrix[i][j] = rand () % 2 *2-1;
}
}  Reply With Quote

12. Re: Ising Model C++ Metropolis Algorithm

Like this
Code:
#include <iostream>
#include <cstdlib>
#include <ctime>
using namespace std;

#define spin(a, b) (matrix[(a)][(b)])

int main()
{
const int LatSize = 5;
const int L = LatSize * LatSize;

int matrix[LatSize][LatSize];

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

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

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

int H = 2 * spin(x, y);

cout << "x: " << x << "  y: " << y << "  H: " << H << endl;
return 0;
}  Reply With Quote

13. Junior Member Join Date
Apr 2014
Posts
24

Re: Ising Model C++ Metropolis Algorithm Originally Posted by 2kaud Like this
]
Ok, so here's the full code, part 1 and part 2:

//Generate NxN lattice with -1 or +1 at each site randomly.

#include <stdlib.h>
#include <iostream>
#include <cstdlib>
#include <ctime>
using namespace std;

#define spin(a, b) (matrix[(a)][(b)])

int main()
{

const int LatSize = 5;
const int L = LatSize * LatSize;
int matrix[LatSize][LatSize];

srand ( static_cast<unsigned> ( time ( 0 ) ) );
for ( int i = 0; i < LatSize; i++ )
{
for ( int j = 0; j < LatSize; j++ )
{
matrix[i][j] = rand () % 2 *2-1;
}
}

//End of Generating matrix

//Extracting value of a site at location (x,y)

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

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

int delta_neighbour = spin(x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
delta_neighbour = -2*spin(x,y)* delta_neighbour; //change in neighbour product energy

cout << delta_neighbour;

//End of extracting value of site at (x, y)

return 0;
}

Doesn't seem to work, as it doesn't show the value of delta_neighbour..  Reply With Quote

14. Re: Ising Model C++ Metropolis Algorithm Originally Posted by 2kaud Another way would be to use a define at the top of the code
Code:
#define spin(x,y) (matrix[(x)][(y)])
That's horrible. Why would you advice someone to use a macro for this? A simple text-replace can solve this without obfuscating the code.

@OP Using two different names for the same variable will make your code harder to understand. When you read your code, every time you see 'spin' you need to think "right, that actually just means matrix". If you just use a single name, you'll save yourself a lot of unnecessary thinking, freeing up your brain capacity to understand the code you're writing/debugging/reading/whatever. If spin better describes the values in the matrix, then go ahead and rename your matrix variable to spin (it's a simple search-and-replace). But don't use different names for the same thing unless you want to make your code harder to understand.  Reply With Quote

15. Re: Ising Model C++ Metropolis Algorithm Originally Posted by unscientific Doesn't seem to work, as it doesn't show the value of delta_neighbour..
Did you step through your code with the debugger? What values are you getting for x and y?  Reply With Quote

algorithm, c++  Posting Permissions

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