-
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;
}
}
-
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.
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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
-
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;
-
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)?
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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
-
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)?
-
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.
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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.
-
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)])
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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;
}
}
-
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;
}
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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..
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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.
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
unscientific
Doesn't seem to work, as it doesn't show the value of delta_neighbour..
Please use code tags! Your code is practically unreadable without them.
Did you step through your code with the debugger? What values are you getting for x and y?
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
D_Drmmr
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.
I agree - that's why I said in my post #10 that the easiest way would be to do text-replace. :) I only offered the #define as an option to try to get the code posted by the OP to compile as I suspect that the part of the code using spin() was obtained elsewhere and I not sure of the expertise of the OP based upon the postings.
-
Re: Ising Model C++ Metropolis Algorithm
Code:
int x = int (rand()*LatSize);
int y = int (rand()*LatSize);
What's with the '*' ? It should be % (modulo) as per my post #2!
Code:
int x = int (rand()%LatSize);
int y = int (rand()%LatSize);
Also, you don't need #include <stdlib.h> as you have #include <cstdlib>.
As D_Drmmr said in his post #15, you need to use the debugger to find out why things don't work as expected.
Before coding a program, it should first be designed (stating in English or other natural language how the program is going to work
Code:
define input requirements
define output requirements
define required algorithms
do {
design program
code program
test program
debug program
} while program_not_correct
Before you continue writing 'ad hoc' code I suggest that you first produce a program design and then code the program from this design.
What os/compiler are you using and what's your level of expertise with c++?
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
D_Drmmr
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.
I think it's more meaningful to use spin(x ,y) instead of matrix (x,y). This shouldn't screw up the code, right?
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
2kaud
Code:
int x = int (rand()*LatSize);
int y = int (rand()*LatSize);
What's with the '*' ? It should be % (modulo) as per my post #2!
Code:
int x = int (rand()%LatSize);
int y = int (rand()%LatSize);
Also, you don't need #include <stdlib.h> as you have #include <cstdlib>.
As D_Drmmr said in his post #15, you need to use the debugger to find out why things don't work as expected.
Before coding a program, it should first be designed (stating in English or other natural language how the program is going to work
Code:
define input requirements
define output requirements
define required algorithms
do {
design program
code program
test program
debug program
} while program_not_correct
Before you continue writing 'ad hoc' code I suggest that you first produce a program design and then code the program from this design.
What os/compiler are you using and what's your level of expertise with c++?
I'm using an online compiler: http://www.compileonline.com/compile_cpp_online.php
Code:
//Generate NxN lattice with -1 or +1 at each site randomly.
#include <stdlib.h>
#include <stdio.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
// Displaying the matrix on the screen
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<<matrix[i][j]; // display the current element out of the array
}
cout<<endl; // when the inner loop is done, go to a new line
}
//End of Displaying Matrix
//Boundary conditions
int x = int (rand() %LatSize);
int y = int (rand() %LatSize);
if(x<0) x += LatSize;
if(x>=LatSize) x -= LatSize;
if(y<0) y += LatSize;
if(y>=LatSize) y -= LatSize;
//End of Boundary Conditions
//Extracting value of a site at location (x,y)
printf ("The site chosen is at (x,y) = %x, %y", x, y);
cout << spin(x, y);
cout<<endl;
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;
}
-
Re: Ising Model C++ Metropolis Algorithm
Code:
if(x<0) x += LatSize;
if(x>=LatSize) x -= LatSize;
if(y<0) y += LatSize;
if(y>=LatSize) y -= LatSize;
these statements are not needed as x and y will always satisfy the condition >= 0 and < LatSize due to how these values are generated.
As you are using c++ and cout, why use printf instead of cout? The formating used is incorrect anyhow. Use
Code:
cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
Code:
int delta_neighbour = spin(x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
Note that this statement would cause the bounds of matrix to be exceeded if x or y is either 0 or LatSize -1. If x or y is either of these values, what element of matrix do you want to be used?
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
2kaud
Code:
if(x<0) x += LatSize;
if(x>=LatSize) x -= LatSize;
if(y<0) y += LatSize;
if(y>=LatSize) y -= LatSize;
these statements are not needed as x and y will always satisfy the condition >= 0 and < LatSize due to how these values are generated.
As you are using c++ and cout, why use printf instead of cout? The formating used is incorrect anyhow. Use
Code:
cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
Code:
int delta_neighbour = spin(x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
Note that this statement would cause the bounds of matrix to be exceeded if x or y is either 0 or LatSize -1. If x or y is either of these values, what element of matrix do you want to be used?
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!
Code:
//Generate NxN lattice with -1 or +1 at each site randomly.
#include <stdlib.h>
#include <stdio.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
// Displaying the matrix on the screen
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<<matrix[i][j]; // display the current element out of the array
}
cout<<endl; // when the inner loop is done, go to a new line
}
//End of Displaying Matrix
//Boundary conditions
int x = int (rand() %LatSize);
int y = int (rand() %LatSize);
if(x<0) {x += LatSize;}
if(x>=LatSize) {x -= LatSize;}
if(y<0) {y += LatSize;}
if(y>=LatSize) {y -= LatSize;}
//End of Boundary Conditions
//Extracting value of a site at location (x,y)
int delta_M = -2*matrix[x][y]; //change in magnetization energy
int delta_neighbour = matrix[x-1][y] + matrix[x+1][y]+ matrix[x][y-1] + matrix[x][y+1];
delta_neighbour = -2*matrix[x][y]* delta_neighbour; //change in neighbour product energy
printf ("The site chosen is at (x,y) = (%d, %d)", x, y);
cout<<endl;
printf ("The spin at this site is:"); cout << matrix[x][y];
cout<<endl;
printf ("The left neighbour has spin :"); cout<< matrix[x][y-1];
cout<<endl;
printf ("The right neighbour has spin:"); cout<<matrix[x][y+1];
cout<<endl;
printf ("The above neighbour has spin:"); cout<<matrix[x-1][y];
cout<<endl;
printf ("The below neighbour has spin:"); cout<<matrix[x+1][y-1];
cout<<endl;
printf ("The change in neighbour energy is:"); cout << delta_neighbour;
//End of extracting value of site at (x, y)
return 0;
}
-
Re: Ising Model C++ Metropolis Algorithm
Here's my latest code, updated with improved cout, and with printf gone! Everything works fine except the boundary conditions..
Code:
//Generate NxN lattice with -1 or +1 at each site randomly.
#include <stdlib.h>
#include <stdio.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
// Displaying the matrix on the screen
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<<matrix[i][j]; // display the current element out of the array
}
cout<<endl; // when the inner loop is done, go to a new line
}
//End of Displaying Matrix
//Boundary conditions
int x = int (rand() %LatSize);
int y = int (rand() %LatSize);
if(x<0) {x += LatSize;}
if(x>=LatSize) {x -= LatSize;}
if(y<0) {y += LatSize;}
if(y>=LatSize) {y -= LatSize;}
//End of Boundary Conditions
//Extracting value of a site at location (x,y)
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 << "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;
//End of extracting value of site at (x, y)
return 0;
}
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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.
Code:
#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;
}
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
2kaud
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.
Code:
#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;
}
It works, although this level of coding is too high for me. Let's see if I can work out some form of understanding:
1. Starting with an entire class called 'cSpin' allows us to do interesting things within the class and later the statement 'cSpin spin' makes the function 'spin' possess all the conditions, follow any iteration, tasks that are defined within that class. And those tasks may or may not be called out to perform in future.
2. sections within the class like 'void generate()' and 'void display()' are some of the tasks the computer might be doing. I'm guessing in the statement 'void generate()', generate is the 'callsign' for that set of tasks and void means you can define a function later on that uses the callsign to do those tasks.
3. Not sure what the difference between public and private is. It seems like public is where you chuck in all the things (iteration, tasks, conditions) you want the class to possess. As for private, I'm not sure what it is used for.
4. It seems the class is defined before the main program 'int main()'. I'm guessing you need to sort out the details before you execute the main program?
4. spin.generate() and spin.display() uses the callsigns 'generate' and 'display' to tell the function 'spin' to do all those tasks.
Is my understanding right?
Questions that bug me
1. What does the statement 'int operator()(int x, int y)' mean?
2. What's the difference between private and public within a class?
-
Re: Ising Model C++ Metropolis Algorithm
If you haven't met c++ classes before, you need a crash course in them.
See
http://www.learncpp.com/cpp-tutorial...class-members/
http://www.cplusplus.com/doc/tutorial/classes/
http://www.tutorialspoint.com/cplusp...es_objects.htm
What book/resource are you currently using from which to learn c++?
Very simply,
1) cSpin is the class (ie type) and the variable spin is an instance of this class.
2) generate and display are class functions (methods). They take the same format as 'normal' functions. So void generate() is a function of class cSpin that has a return type of void - ie doesn't return a value.
3) private definitions can only be accessed within the class. Public ones can be accessed outside of the class. If you try to access spin.matrix[x][y] in main() you'll get a compilation error as matrix is defined private so can't be accessed outside of the cSpin class.
4) classes (like functions, struct etc) need to be defined/declared before they are used. With more complex classes it it usual for the class to be declared in a header file which is included in the main program and for the class functions to be defined in a seperate compilation unit.
5) spin.generate() and spin.display() call the class cSpin functions generate() and display() for the class instance variable spin.
Questions
1) int operator()(int x, int y) overloads the () operator for the class cSpin so that variables of the class can be used as if they are functions. So a = spin(1, 2) calls the operator() function for the spin instance of the cSpin class and returns a value of type int.
2) see 3) above.
I hope this of help. You really need to learn about classes and the Standard Template Library if you are using c++.
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
2kaud
If you haven't met c++ classes before, you need a crash course in them.
See
http://www.learncpp.com/cpp-tutorial...class-members/
http://www.cplusplus.com/doc/tutorial/classes/
http://www.tutorialspoint.com/cplusp...es_objects.htm
What book/resource are you currently using from which to learn c++?
Very simply,
1) cSpin is the class (ie type) and the variable spin is an instance of this class.
2) generate and display are class functions (methods). They take the same format as 'normal' functions. So void generate() is a function of class cSpin that has a return type of void - ie doesn't return a value.
3) private definitions can only be accessed within the class. Public ones can be accessed outside of the class. If you try to access spin.matrix[x][y] in main() you'll get a compilation error as matrix is defined private so can't be accessed outside of the cSpin class.
4) classes (like functions, struct etc) need to be defined/declared before they are used. With more complex classes it it usual for the class to be declared in a header file which is included in the main program and for the class functions to be defined in a seperate compilation unit.
5) spin.generate() and spin.display() call the class cSpin functions generate() and display() for the class instance variable spin.
Questions
1) int operator()(int x, int y) overloads the () operator for the class cSpin so that variables of the class can be used as if they are functions. So a = spin(1, 2) calls the operator() function for the spin instance of the cSpin class and returns a value of type int.
2) see 3) above.
I hope this of help. You really need to learn about classes and the Standard Template Library if you are using c++.
Ok, here comes the interesting part. Of each spin (x,y) randomly chosen, we need to decide whether to flip it or not.
If delta_E < 0, flip. If delta_E > 0, flip with acceptance probability w = exp(-delta_E).
I think the code would be something like:
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
int flip_counter++; // Total number of sites chosen up till this point
int sweep_counter = flip_counter/(LatSize*LatSize) //Total number of sweeps till this point
return 1;
}
if (delta_E>0)
{
if flipped, M and E is updated. If not flipped, M and E is unchanged.
}
Points of concern:
1. How to ensure the same site is never picked twice in a sweep? And how do we tell when one sweep is done? For example in a 5x5 lattice, 25 different picks = 1 sweep.
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
unscientific
How to ensure the same site is never picked twice in a sweep?
you can't ( you can't choose spins randomly and indipendently *and* at the same time enforce non repetitions as this will define a dependency between spins ) why do you have such a requirement ? it's fine ( and physically sensible too ) to have repetitions as long as you realize the relationship between the simulation time and the iterations count ...
-
Re: Ising Model C++ Metropolis Algorithm
Code:
spin(x, y) *= -1; // flip spin
As the class is currently coded, you can't do this as operator() returns a value not a reference. You need to change operator() to
Code:
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];
}
From your first program
Code:
double delta_E = -J*delta_neighbour -B*delta_M;
what are J and B? What are M and E?
Code:
int flip_counter++; // Total number of sites chosen up till this point
int sweep_counter = flip_counter/(LatSize*LatSize) //Total number of sweeps till this point
No! flip_counter and sweep_counter are defined inside the block so their scope is only for the duration of the block. You can't post increment a variable that hasn't been initialised. This is basic c++ stuff. Before you go much further I would recommend that you get a good book on c++ and produce a design for your program.
Why? What os are you using?
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
2kaud
Code:
spin(x, y) *= -1; // flip spin
As the class is currently coded, you can't do this as operator() returns a value not a reference. You need to change operator() to
Code:
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];
}
By changing the operator to 'int & operator()(int x, int y)' would this fix it?
Quote:
Originally Posted by
2kaud
From your first program
Code:
double delta_E = -J*delta_neighbour -B*delta_M;
what are J and B? What are M and E?
J is a constant defined to be 1. M is a sum of all spins, so it starts from 0. I will define 'int M=0' before the public class section and as the NxN lattice is generated randomly I would insert in a line M+= matrix[i][j] here:
Code:
int M=0;
void generate()
{
for (int i = 0; i < LatSize; i++)
for (int j = 0; j < LatSize; j++)
matrix[i][j] = (rand() % 2) * 2 - 1;
for (i < LatSize) {M+= matrix[i][j];}
}
Quote:
Originally Posted by
2kaud
Code:
int flip_counter++; // Total number of sites chosen up till this point
int sweep_counter = flip_counter/(LatSize*LatSize) //Total number of sweeps till this point
No! flip_counter and sweep_counter are defined inside the block so their scope is only for the duration of the block. You can't post increment a variable that hasn't been initialised. This is basic c++ stuff. Before you go much further I would recommend that you get a good book on c++ and produce a design for your program.
Why? What os are you using?
Does 'int flip_counter' and 'int sweep_counter' go into the public class then?
At school I usually do this at the computer labs, where they use Mac OS X. But now I'm using windows 7.
-
Re: Ising Model C++ Metropolis Algorithm
If you are using Windows 7 why don't you use the free Microsoft Visual Studio Express 2013?
See http://www.visualstudio.com/download...sual-studio-vs
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++?
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
2kaud
If you are using Windows 7 why don't you use the free Microsoft Visual Studio Express 2013?
See
http://www.visualstudio.com/download...sual-studio-vs
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...
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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;
}
Quote:
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?
Quote:
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.
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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;
}
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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.
Quote:
If delta_E > 0, flip with acceptance probability w = exp(-delta_E)
?? :confused:
Quote:
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.
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
2kaud
Yes, as in maths the x-axis should go left to right and the y-axis go top to bottom.
?? :confused:
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?
-
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
}
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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."
Quote:
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;
}
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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?
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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.
-
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.
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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;
}
-
Re: Ising Model C++ Metropolis Algorithm
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?
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
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;
}
-
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.
-
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
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
2kaud
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
This would mean that the site (x, y) would be flipped, without doubt. If delta_E < 0, flip. If delta_E > 0, it means that the exponential is between 0 and 1, so flip. I need the flipping process to be such:
1. If delta_E < 0 flip.
2. Otherwise, flip with probability exp (-delta_E). If w = 0.7 = 70%, this means that if there are 10 sites with the same value of w, only 7 out of 10 times will site be flipped.
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
This would mean that the site (x, y) would be flipped, without doubt
No
If delta_E <= 0 then flip
if delta_E > 0 then flip if exp(-delta_E) > random_number_0_to_1
When evaluating logical or conditions, the left hand expression is evaluated first. If this is true then the right hand expression is not evaluated. The right hand expression is only evaluated if the left hand expression is evaluated to false. if delta_E is <= 0 then the left condition is true and the right is not evaluated. The right expression is only evaluated if the left is false.
See http://msdn.microsoft.com/en-us/library/z68fx2f1.aspx
-
Re: Ising Model C++ Metropolis Algorithm
Quote:
Originally Posted by
2kaud
No
If delta_E <= 0 then flip
if delta_E > 0 then flip if exp(-delta_E) > random_number_0_to_1
Ah I see, I missed out the > random_number_0_to_1. Ok this makes sense. If exp(-delta_E) is between 0 and 1, the chances of it being bigger than a number between 0 and 1 is also 0 and 1. This will work!
