Ising Model C++ Metropolis Algorithm - Page 3
CodeGuru Home VC++ / MFC / C++ .NET / C# Visual Basic VB Forums Developer.com
Page 3 of 4 FirstFirst 1234 LastLast
Results 31 to 45 of 48

Thread: Ising Model C++ Metropolis Algorithm

  1. #31
    Join Date
    Apr 2014
    Posts
    24

    Re: Ising Model C++ Metropolis Algorithm

    Quote Originally Posted by 2kaud View Post
    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...

  2. #32
    Join Date
    Dec 2012
    Location
    England
    Posts
    2,377

    Re: Ising Model C++ Metropolis Algorithm

    the neighbours of (x,y) give the wrong values
    Upon looking at the code I put together from your previous postings, there appears to be inconsitencies between which element of the matix is considered to be x and which to be y. I think the code below is consistent now but I haven't tested it.
    Code:
    #include <iostream>
    #include <iomanip>
    #include <cstdlib>
    #include <ctime>
    using namespace std;
    
    const int LatSize = 5;
    
    class cSpin
    {
    public:
    	cSpin() : M(0)
    	{
    		for (int i = 0; i < LatSize; i++) 
    			for (int j = 0; j < LatSize; j++)
    				M += matrix[j][i] = (rand() % 2) * 2 - 1;
    	}
    
    	int& operator()(int x, int y)
    	{
    		if (x < 0)
    			x += LatSize;
    
    		if (x >= LatSize) 
    			x -= LatSize;
    
    		if (y < 0) 
    			y += LatSize;
    
    		if (y >= LatSize)
    			y -= LatSize;
    
    		return matrix[x][y];
    	}
    
    	void display()
    	{
    		for (int i = 0; i < LatSize; i++)  // loop 3 times for three lines
    		{
    			for (int j = 0; j < LatSize; j++)  // loop for the three elements on the line
    				cout << setw(2) << matrix[j][i] << " ";  // display the current element out of the array
    
    			cout << endl;  // when the inner loop is done, go to a new line
    		}
    	}
    
    	int SpinSum()
    	{
    		return M;
    	}
    
    private:
    	int matrix[LatSize][LatSize];
    	int M;
    };
    
    int main()
    {
    const int J = 1;
    const int B = 1;
    
    cSpin	spin;
    
    	srand (static_cast<unsigned> (time(0)));
    
    	spin.display();
    
    int M = spin.SpinSum();
    
    int x = int (rand() % LatSize); 
    int y = int (rand() % LatSize);
    
    int neighbour = 0;    // nearest neighbour count
    int flip_counter = 0;
    double sweep_counter = 0.0;
    
    	for (int i = 0; i < LatSize; i++)
    		for (int j = 0; j < LatSize; j++)
    		{
    			neighbour += (spin(j, i) == spin(j + 1, i)) ? 1 : -1;
    			neighbour += (spin(j, i) == spin(j, i + 1)) ? 1 : -1;
    		}
    
    int E = -J * neighbour - B * M;
    
    int delta_M = -2 * spin(x, y);	//change in magnetization energy
    
    int delta_neighbour = delta_M * (spin(x - 1, y) + spin(x + 1, y) + spin(x, y - 1) + spin(x, y + 1)); 
    
    int delta_E = -J * delta_neighbour - B * delta_M;
    
    	if (delta_E <= 0)
    	{
    		spin(x, y) *= -1;											// flip spin
    		M += delta_M;												// New Magnetization energy of flipped configuration
    		E += delta_E;												// New total energy of flipped configuration
    		flip_counter++;												// Total number of sites chosen up till this point
    		sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
    	}
    	else
    	{
    		//if flipped, M and E is updated. If not flipped, M and E is unchanged.
    	}
    
    	cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
    	cout << "The left neighbour has spin = " << spin(x, y - 1) << endl;
    	cout << "The right neighbour has spin = " << spin(x, y + 1) << endl;
    	cout << "The above neighbour has spin = " << spin(x - 1, y) << endl;
    	cout << "The below neighbour has spin = " << spin(x + 1, y) << endl;
    	cout << "The change in neighbour energy = " << delta_neighbour << endl;
    
    	return 0;
    }
    the value of M (Sum of all spins) are wrong..
    Then you need to debug the program to find out why the value is wrong. Where does the code deviate from that expected from the program design?

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

    Also try those web site references I provided in a previous post. You can't really just 'jump into' a section on c++. You need to understand the language from the beginning in a structed way.
    Last edited by 2kaud; April 12th, 2014 at 02:38 PM.
    All advice is offered in good faith only. You are ultimately responsible for effects of your programs and the integrity of the machines they run on.

  3. #33
    Join Date
    Apr 2014
    Posts
    24

    Re: Ising Model C++ Metropolis Algorithm

    Quote Originally Posted by 2kaud View Post

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

    The value of M is also right now.

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

    Code:
    #include <iostream>
    #include <iomanip>
    #include <cstdlib>
    #include <ctime>
    using namespace std;
    
    const int LatSize = 3;
    
    class cSpin
    {
    public:
        cSpin() : M(0)
    	{
    		for (int i = 0; i < LatSize; i++) 
    			for (int j = 0; j < LatSize; j++)
    				M += matrix[j][i] = (rand() % 2) * 2 - 1;
    	}
    
    	int& operator()(int x, int y)
    	{
    		if (x < 0)
    			x += LatSize;
    
    		if (x >= LatSize) 
    			x -= LatSize;
    
    		if (y < 0) 
    			y += LatSize;
    
    		if (y >= LatSize)
    			y -= LatSize;
    
    		return matrix[x][y];
    	}
    
    	void display()
    	{
    		for (int i = 0; i < LatSize; i++)  // loop 3 times for three lines
    		{
    			for (int j = 0; j < LatSize; j++)  // loop for the three elements on the line
    				cout << setw(2) << matrix[j][i] << " ";  // display the current element out of the array
    
    			cout << endl;  // when the inner loop is done, go to a new line
    		}
    	}
    
    	int SpinSum()
    	{
    		return M;
    	}
    
    private:
    	int matrix[LatSize][LatSize];
    	int M;
    };
    
    int main()
    {
    const int J = 1;
    const int B = 1;
    
    cSpin	spin;
    
    	srand (static_cast<unsigned> (time(0)));
    
    	spin.display();
    
    int M = spin.SpinSum();
    
    int x = int (rand() % LatSize); 
    int y = int (rand() % LatSize);
    
    int neighbour = 0;    // nearest neighbour count
    int flip_counter = 0;
    double sweep_counter = 0.0;
    
    	
    
    int E = -J * neighbour - B * M;
    
    int delta_M = -2 * spin(x, y);	//change in magnetization energy
    
    int delta_neighbour = delta_M * (spin(x - 1, y) + spin(x + 1, y) + spin(x, y - 1) + spin(x, y + 1)); 
    
    int delta_E = -J * delta_neighbour - B * delta_M;
    
    	if (delta_E <= 0)
    	{
    		spin(x, y) *= -1;											// flip spin
    		M += delta_M;												// New Magnetization energy of flipped configuration
    		E += delta_E;												// New total energy of flipped configuration
    		flip_counter++;												// Total number of sites chosen up till this point
    		sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
    	}
    	else
    	{
    		//if flipped, M and E is updated. If not flipped, M and E is unchanged.
    	}
    
    	cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
    	cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
    	cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
    	cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
    	cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
    	cout << "The change in neighbour energy = " << delta_neighbour << endl;
        cout << "The change in total energy dE =" << delta_E << endl;
        cout << "The number of flips = " << flip_counter <<endl;
        cout << "The magnetization Energy = " << M <<endl;
    
    	return 0;
    }

  4. #34
    Join Date
    Dec 2012
    Location
    England
    Posts
    2,377

    Re: Ising Model C++ Metropolis Algorithm

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

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

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

  5. #35
    Join Date
    Apr 2014
    Posts
    24

    Re: Ising Model C++ Metropolis Algorithm

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

    ??

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

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


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

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

  6. #36
    Join Date
    Dec 2012
    Location
    England
    Posts
    2,377

    Re: Ising Model C++ Metropolis Algorithm

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

    The easiest way to get the 100 sweeps would be to simply use a for loop
    Code:
    for (int s = 0; s < 100 * LatSize * LatSize; ++s) {
        //code for one calculation goes here
    }
    Last edited by 2kaud; April 13th, 2014 at 06:42 AM.
    All advice is offered in good faith only. You are ultimately responsible for effects of your programs and the integrity of the machines they run on.

  7. #37
    Join Date
    Apr 2014
    Posts
    24

    Re: Ising Model C++ Metropolis Algorithm

    Quote Originally Posted by 2kaud View Post
    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 View Post
    You may also need to insert #include <cmath> at the top of the code.

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

    The complete code would be this:

    Code:
    #include <math.h>
    #include <iostream>
    #include <iomanip>
    #include <cstdlib>
    #include <ctime>
    using namespace std;
    
    const int LatSize = 3;
    
    class cSpin
    {
    public:
        cSpin() : M(0)
        {
        	for (int i = 0; i < LatSize; i++) 
    			for (int j = 0; j < LatSize; j++)
    				M += matrix[j][i] = (rand() % 2) * 2 - 1;
    	}
    
    	int& operator()(int x, int y)
    	{
    		if (x < 0)
    			x += LatSize;
    
    		if (x >= LatSize) 
    			x -= LatSize;
    
    		if (y < 0) 
    			y += LatSize;
    
    		if (y >= LatSize)
    			y -= LatSize;
    
    		return matrix[x][y];
    	}
    
    	void display()
    	{
    		for (int i = 0; i < LatSize; i++)  // loop 3 times for three lines
    		{
    			for (int j = 0; j < LatSize; j++)  // loop for the three elements on the line
    				cout << setw(2) << matrix[j][i] << " ";  // display the current element out of the array
    
    			cout << endl;  // when the inner loop is done, go to a new line
    		}
    	}
    
    	int SpinSum()
    	{
    		return M;
    	}
    
    private:
    	int matrix[LatSize][LatSize];
    	int M;
    };
    
    int main()
    {
    const int J = 1;
    const int B = 1;
    
    cSpin	spin;
    
    	srand (static_cast<unsigned> (time(0)));
    
    	spin.display();
    
    int M = spin.SpinSum();
    
    int x = int (rand() % LatSize); 
    int y = int (rand() % LatSize);
    
    int neighbour = 0;    // nearest neighbour count
    int flip_counter = 0; // number of flips count
    int sweep_counter = 0; //number of sweeps count
    
    	
    
    int E = -J * neighbour - B * M;
    
    int delta_M = -2 * spin(x, y);	//change in magnetization energy
    
    int delta_neighbour = delta_M * (spin(x - 1, y) + spin(x + 1, y) + spin(x, y - 1) + spin(x, y + 1)); 
    
    int delta_E = -J * delta_neighbour - B * delta_M;
    
           for (int s = 0; s < 50 * LatSize * LatSize; ++s)
    {
    
    	if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp(-delta_E))
    	{
    		spin(x, y) *= -1;											// flip spin
    		M += delta_M;												// New Magnetization energy of flipped configuration
    		E += delta_E;												// New total energy of flipped configuration
    		flip_counter++;												// Total number of sites chosen up till this point
    		sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
    	}
    	else
    	{
    	flip_counter++;
        sweep_counter = (double)flip_counter / (LatSize * LatSize);     //if flipped, M and E is updated. If not flipped, M and E is unchanged.
    	}
    }
    	cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
    	cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
    	cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
    	cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
    	cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
    	cout << "The change in neighbour energy = " << delta_neighbour << endl;
        cout << "The change in total energy dE =" << delta_E << endl;
        cout << "The number of flips = " << flip_counter <<endl;
        cout << "The number of sweeps = " << sweep_counter <<endl;
        cout << "The magnetization Energy = " << M <<endl;
    
    	return 0;
    }

  8. #38
    Join Date
    Dec 2012
    Location
    England
    Posts
    2,377

    Re: Ising Model C++ Metropolis Algorithm

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

    What parts of the code do you want to be part of the loop and which do you want to executes only once before the loop?
    All advice is offered in good faith only. You are ultimately responsible for effects of your programs and the integrity of the machines they run on.

  9. #39
    Join Date
    Apr 2014
    Posts
    24

    Re: Ising Model C++ Metropolis Algorithm

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

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

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

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

  10. #40
    Join Date
    Dec 2012
    Location
    England
    Posts
    2,377

    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.
    All advice is offered in good faith only. You are ultimately responsible for effects of your programs and the integrity of the machines they run on.

  11. #41
    Join Date
    Apr 2014
    Posts
    24

    Re: Ising Model C++ Metropolis Algorithm

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

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

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

    Code:
    #include <math.h>
    #include <iostream>
    #include <iomanip>
    #include <cstdlib>
    #include <ctime>
    using namespace std;
    
    const int LatSize = 3;
    
    class cSpin
    {
    public:
        cSpin() : M(0)
        {
            for (int i = 0; i < LatSize; i++) 
        		for (int j = 0; j < LatSize; j++)
    				M += matrix[j][i] = (rand() % 2) * 2 - 1;
    	}
    
    	int& operator()(int x, int y)
    	{
    		if (x < 0)
    			x += LatSize;
    
    		if (x >= LatSize) 
    			x -= LatSize;
    
    		if (y < 0) 
    			y += LatSize;
    
    		if (y >= LatSize)
    			y -= LatSize;
    
    		return matrix[x][y];
    	}
    
    	void display()
    	{
    		for (int i = 0; i < LatSize; i++)  // loop 3 times for three lines
    		{
    			for (int j = 0; j < LatSize; j++)  // loop for the three elements on the line
    				cout << setw(2) << matrix[j][i] << " ";  // display the current element out of the array
    
    			cout << endl;  // when the inner loop is done, go to a new line
    		}
    	}
    
    	int SpinSum()
    	{
    		return M;
    	}
    
    private:
    	int matrix[LatSize][LatSize];
    	int M;
    };
    
    int main()
    {
    const int J = 1;
    const int B = 1;
    
    cSpin	spin;
    
    	srand (static_cast<unsigned> (time(0)));
    
    	spin.display();
    
    int M = spin.SpinSum();
    
    int x = int (rand() % LatSize);       //Randomly choose x-coordinate
    int y = int (rand() % LatSize);       //Randomly choose y-coordinate
    
    int neighbour = 0;    // nearest neighbour count
    int flip_counter = 0; // number of flips count
    int sweep_counter = 0; //number of sweeps count
    int sweep_number = 5;  // number of sweeps
    double M_avg;          // average magnetic moment
    
    	
    
    int E = -J * neighbour - B * M;   //Total Energy
    
    int delta_M = -2 * spin(x, y);	//change in Magnetic Moment
    
    int delta_neighbour = delta_M * (spin(x - 1, y) + spin(x + 1, y) + spin(x, y - 1) + spin(x, y + 1)); 
    
    int delta_E = -J * delta_neighbour - B * delta_M; //Change in total energy
    
    int M_Sum = 0;   //Cumulative Magnetic Moment starts from 0
    
           
           for (int s = 0; s < sweep_number; ++s)    //Loop for sweep_number of sweeps
    {
        
            for (int t = 0; t < LatSize * LatSize; ++t)   //Loop for 1 Sweep
            {
    
    	if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp(-delta_E))
    	{
    		spin(x, y) *= -1;											// flip spin
    		M += delta_M;												// New Magnetization energy of flipped configuration
    		E += delta_E;												// New total energy of flipped configuration
    		flip_counter++;												// Total number of sites chosen up till this point
    		sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
             
            
    	}
    	else
    	{
    	flip_counter++;   //Number of Flips
        
        sweep_counter = (double)flip_counter / (LatSize * LatSize);     //Number of Sweeps
        
         
    	}
        
        
    }
        M_Sum += M;                       //Summing up M to find M_Sum
        M_avg = (M_Sum)/sweep_counter;    //Average M
        
    
        cout << "The Magnetic Moment after "<< sweep_counter <<" sweeps = " << M << endl;
        cout << "The Cumulative Magnetic Moment after " << sweep_counter <<" sweeps = " << M_Sum << endl;
        cout << "The Average Magnetic Moment after " << sweep_counter << "sweeps = " << M_avg <<endl;
        
    }
        
    
    	//cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
    	//cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
    	//cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
    	//cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
    	//cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
    	//cout << "The change in neighbour energy = " << delta_neighbour << endl;
        //cout << "The change in total energy dE =" << delta_E << endl;
        cout << "The number of flips = " << flip_counter <<endl;
        cout << "The number of sweeps = " << sweep_counter <<endl;
        cout << "The final Magnetic Moment = " << M <<endl;
        cout << "The final average Magnetic Moment = " << M_avg <<endl;
    
    	return 0;
    }
    Last edited by unscientific; April 14th, 2014 at 08:50 AM.

  12. #42
    Join Date
    Dec 2012
    Location
    England
    Posts
    2,377

    Re: Ising Model C++ Metropolis Algorithm

    it seems to work
    Well as long as it seems then that's all right then? Have you tested it to make sure it gives the answers expected? I don't think it will as x and y are only set once at the beginning of the program and are not changed during each loop. Also you aren't calculating the initial value of neighbour. I think the code you want is more like something like this
    Code:
    #include <cmath>
    #include <iostream>
    #include <iomanip>
    #include <cstdlib>
    #include <ctime>
    using namespace std;
    
    const int LatSize = 3;
    
    class cSpin
    {
    public:
        cSpin() : M(0)
        {
            for (int i = 0; i < LatSize; i++) 
        		for (int j = 0; j < LatSize; j++)
    				M += matrix[j][i] = (rand() % 2) * 2 - 1;
    	}
    
    	int& operator()(int x, int y)
    	{
    		if (x < 0)
    			x += LatSize;
    
    		if (x >= LatSize) 
    			x -= LatSize;
    
    		if (y < 0) 
    			y += LatSize;
    
    		if (y >= LatSize)
    			y -= LatSize;
    
    		return matrix[x][y];
    	}
    
    	void display()
    	{
    		for (int i = 0; i < LatSize; i++)  // loop 3 times for three lines
    		{
    			for (int j = 0; j < LatSize; j++)  // loop for the three elements on the line
    				cout << setw(2) << matrix[j][i] << " ";  // display the current element out of the array
    
    			cout << endl;  // when the inner loop is done, go to a new line
    		}
    	}
    
    	int SpinSum()
    	{
    		return M;
    	}
    
    private:
    	int matrix[LatSize][LatSize];
    	int M;
    };
    
    int main()
    {
    const int J = 1;
    const int B = 1;
    
    cSpin	spin;
    
    	srand (static_cast<unsigned> (time(0)));
    
    	spin.display();
    
    int M = spin.SpinSum();
    
    int neighbour = 0;    // nearest neighbour count
    int flip_counter = 0; // number of flips count
    double sweep_counter = 0; //number of sweeps count
    int sweep_number = 5;  // number of sweeps
    double M_avg;          // average magnetic moment
    int M_Sum = 0;   //Cumulative Magnetic Moment starts from 0
    
    	for (int s = 0; s < sweep_number; ++s)    //Loop for sweep_number of sweeps
    	{
    		for (int t = 0; t < LatSize * LatSize; ++t)   //Loop for 1 Sweep
    		{
    		int x = int (rand() % LatSize); 
    		int y = int (rand() % LatSize);
    
    		for (int i = 0; i < LatSize; i++)
    			for (int j = 0; j < LatSize; j++)
    			{
    				neighbour += (spin(j, i) == spin(j + 1, i)) ? 1 : -1;
    				neighbour += (spin(j, i) == spin(j, i + 1)) ? 1 : -1;
    			}
    
    			int E = -J * neighbour - B * M;   //Total Energy
    
    			int delta_M = -2 * spin(x, y);	//change in Magnetic Moment
    
    			int delta_neighbour = delta_M * (spin(x - 1, y) + spin(x + 1, y) + spin(x, y - 1) + spin(x, y + 1)); 
    
    			int delta_E = -J * delta_neighbour - B * delta_M; //Change in total energy
    
    			if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp((double)-delta_E))
    			{
    				spin(x, y) *= -1;											// flip spin
    				M += delta_M;												// New Magnetization energy of flipped configuration
    				E += delta_E;												// New total energy of flipped configuration
    			}
    			flip_counter++;												// Total number of sites chosen up till this point
    			sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
    		}
    
    	    M_Sum += M;                       //Summing up M to find M_Sum
    		M_avg = M_Sum / sweep_counter;    //Average M
        
    		cout << "The Magnetic Moment after "<< sweep_counter <<" sweeps = " << M << endl;
    		cout << "The Cumulative Magnetic Moment after " << sweep_counter <<" sweeps = " << M_Sum << endl;
    		cout << "The Average Magnetic Moment after " << sweep_counter << "sweeps = " << M_avg <<endl;
    	}
        
    	//cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
    	//cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
    	//cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
    	//cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
    	//cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
    	//cout << "The change in neighbour energy = " << delta_neighbour << endl;
        //cout << "The change in total energy dE =" << delta_E << endl;
        cout << "The number of flips = " << flip_counter << endl;
        cout << "The number of sweeps = " << sweep_counter << endl;
        cout << "The final Magnetic Moment = " << M << endl;
        cout << "The final average Magnetic Moment = " << M_avg << endl;
    
    	return 0;
    }
    but you need to test/debug the code to make sure that it is giving the correct answers.

    PS Do you also need to calculate a new value of M eveytime time around the loops as well?
    Last edited by 2kaud; April 14th, 2014 at 11:14 AM.
    All advice is offered in good faith only. You are ultimately responsible for effects of your programs and the integrity of the machines they run on.

  13. #43
    Join Date
    Apr 2014
    Posts
    24

    Re: Ising Model C++ Metropolis Algorithm

    Quote Originally Posted by 2kaud View Post
    Well as long as it seems then that's all right then? Have you tested it to make sure it gives the answers expected? I don't think it will as x and y are only set once at the beginning of the program and are not changed during each loop. Also you aren't calculating the initial value of neighbour. I think the code you want is more like something like this
    Code:
    #include <cmath>
    #include <iostream>
    #include <iomanip>
    #include <cstdlib>
    #include <ctime>
    using namespace std;
    
    const int LatSize = 3;
    
    class cSpin
    {
    public:
        cSpin() : M(0)
        {
            for (int i = 0; i < LatSize; i++) 
        		for (int j = 0; j < LatSize; j++)
    				M += matrix[j][i] = (rand() % 2) * 2 - 1;
    	}
    
    	int& operator()(int x, int y)
    	{
    		if (x < 0)
    			x += LatSize;
    
    		if (x >= LatSize) 
    			x -= LatSize;
    
    		if (y < 0) 
    			y += LatSize;
    
    		if (y >= LatSize)
    			y -= LatSize;
    
    		return matrix[x][y];
    	}
    
    	void display()
    	{
    		for (int i = 0; i < LatSize; i++)  // loop 3 times for three lines
    		{
    			for (int j = 0; j < LatSize; j++)  // loop for the three elements on the line
    				cout << setw(2) << matrix[j][i] << " ";  // display the current element out of the array
    
    			cout << endl;  // when the inner loop is done, go to a new line
    		}
    	}
    
    	int SpinSum()
    	{
    		return M;
    	}
    
    private:
    	int matrix[LatSize][LatSize];
    	int M;
    };
    
    int main()
    {
    const int J = 1;
    const int B = 1;
    
    cSpin	spin;
    
    	srand (static_cast<unsigned> (time(0)));
    
    	spin.display();
    
    int M = spin.SpinSum();
    
    int neighbour = 0;    // nearest neighbour count
    int flip_counter = 0; // number of flips count
    double sweep_counter = 0; //number of sweeps count
    int sweep_number = 5;  // number of sweeps
    double M_avg;          // average magnetic moment
    int M_Sum = 0;   //Cumulative Magnetic Moment starts from 0
    
    	for (int s = 0; s < sweep_number; ++s)    //Loop for sweep_number of sweeps
    	{
    		for (int t = 0; t < LatSize * LatSize; ++t)   //Loop for 1 Sweep
    		{
    		int x = int (rand() % LatSize); 
    		int y = int (rand() % LatSize);
    
    		for (int i = 0; i < LatSize; i++)
    			for (int j = 0; j < LatSize; j++)
    			{
    				neighbour += (spin(j, i) == spin(j + 1, i)) ? 1 : -1;
    				neighbour += (spin(j, i) == spin(j, i + 1)) ? 1 : -1;
    			}
    
    			int E = -J * neighbour - B * M;   //Total Energy
    
    			int delta_M = -2 * spin(x, y);	//change in Magnetic Moment
    
    			int delta_neighbour = delta_M * (spin(x - 1, y) + spin(x + 1, y) + spin(x, y - 1) + spin(x, y + 1)); 
    
    			int delta_E = -J * delta_neighbour - B * delta_M; //Change in total energy
    
    			if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp((double)-delta_E))
    			{
    				spin(x, y) *= -1;											// flip spin
    				M += delta_M;												// New Magnetization energy of flipped configuration
    				E += delta_E;												// New total energy of flipped configuration
    			}
    			flip_counter++;												// Total number of sites chosen up till this point
    			sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
    		}
    
    	    M_Sum += M;                       //Summing up M to find M_Sum
    		M_avg = M_Sum / sweep_counter;    //Average M
        
    		cout << "The Magnetic Moment after "<< sweep_counter <<" sweeps = " << M << endl;
    		cout << "The Cumulative Magnetic Moment after " << sweep_counter <<" sweeps = " << M_Sum << endl;
    		cout << "The Average Magnetic Moment after " << sweep_counter << "sweeps = " << M_avg <<endl;
    	}
        
    	//cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
    	//cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
    	//cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
    	//cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
    	//cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
    	//cout << "The change in neighbour energy = " << delta_neighbour << endl;
        //cout << "The change in total energy dE =" << delta_E << endl;
        cout << "The number of flips = " << flip_counter << endl;
        cout << "The number of sweeps = " << sweep_counter << endl;
        cout << "The final Magnetic Moment = " << M << endl;
        cout << "The final average Magnetic Moment = " << M_avg << endl;
    
    	return 0;
    }
    but you need to test/debug the code to make sure that it is giving the correct answers.

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

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

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

    Code:
    #include <math.h>
    #include <iostream>
    #include <iomanip>
    #include <cstdlib>
    #include <ctime>
    using namespace std;
    
    const int LatSize = 30;
    
    class cSpin
    {
    public:
        cSpin() : M(0)
        {
            for (int i = 0; i < LatSize; i++) 
            	for (int j = 0; j < LatSize; j++)
    				M += matrix[j][i] = (rand() % 2) * 2 - 1;
    	}
    
    	int& operator()(int x, int y)
    	{
    		if (x < 0)
    			x += LatSize;
    
    		if (x >= LatSize) 
    			x -= LatSize;
    
    		if (y < 0) 
    			y += LatSize;
    
    		if (y >= LatSize)
    			y -= LatSize;
    
    		return matrix[x][y];
    	}
    
    	void display()
    	{
    		for (int i = 0; i < LatSize; i++)  // loop 3 times for three lines
    		{
    			for (int j = 0; j < LatSize; j++)  // loop for the three elements on the line
    				cout << setw(2) << matrix[j][i] << " ";  // display the current element out of the array
    
    			cout << endl;  // when the inner loop is done, go to a new line
    		}
    	}
    
    	int SpinSum()
    	{
    		return M;
    	}
    
    private:
    	int matrix[LatSize][LatSize];
    	int M;
    };
    
    int main()
    {
    const int J = 1;
    const int B = 1;
    
    cSpin	spin;
    
    	srand (static_cast<unsigned> (time(0)));
    
    	spin.display();
    
    int M = spin.SpinSum();
    
    
    
    int neighbour = 0;    // nearest neighbour count
    int flip_counter = 0; // number of flips count
    int sweep_counter = 0; //number of sweeps count
    int sweep_number = 520;  // number of sweeps
    double M_avg;          // average magnetic moment
    
    	
    
    
    
    
    
    
    
    
    
    int M_Sum = 0;   //Cumulative Magnetic Moment starts from 0
    
           
           for (int s = 0; s < sweep_number; ++s)    //Loop for sweep_number of sweeps
    {
        
            for (int t = 0; t < LatSize * LatSize; ++t)   //Loop for 1 Sweep
            {
    
            int x = int (rand() % LatSize);       //Randomly choose x-coordinate
            int y = int (rand() % LatSize);       //Randomly choose y-coordinate
            
           int delta_M = -2 * spin(x, y);    //change in Magnetic Moment
           int delta_neighbour = delta_M * (spin(x - 1, y) + spin(x + 1, y) + spin(x, y - 1) + spin(x, y + 1)); 
            int delta_E = -J * delta_neighbour - B * delta_M; //Change in total energy
            
            int E = -J * neighbour - B * M;   //Total Energy
             
    
    	if (delta_E <= 0 || ((double)rand() / (double)(RAND_MAX)) < (double)exp(-delta_E))
    	{
    		spin(x, y) *= -1;											// flip spin
    		M += delta_M;												// New Magnetization energy of flipped configuration
    		E += delta_E;												// New total energy of flipped configuration
    		flip_counter++;												// Total number of sites chosen up till this point
    		sweep_counter = (double)flip_counter / (LatSize * LatSize);	//Total number of sweeps till this point
             
            
    	}
    	else
    	{
    	flip_counter++;   //Number of Flips
        
        sweep_counter = (double)flip_counter / (LatSize * LatSize);     //Number of Sweeps
        
         
    	}
        
        
    }
        M_Sum += M;                       //Summing up M to find M_Sum
        M_avg = (M_Sum)/sweep_counter;    //Average M
        
    
        cout << "The Magnetic Moment after "<< sweep_counter <<" sweeps = " << M << endl;
        cout << "The Cumulative Magnetic Moment after " << sweep_counter <<" sweeps = " << M_Sum << endl;
        cout << "The Average Magnetic Moment after " << sweep_counter << "sweeps = " << M_avg <<endl;
        
    }
        
    
    	//cout << "The site chosen is at (" << x << ", " << y << ") = " << spin(x, y) << endl;
    	//cout << "The left neighbour has spin = " << spin(x-1, y) << endl;
    	//cout << "The right neighbour has spin = " << spin(x+1, y) << endl;
    	//cout << "The above neighbour has spin = " << spin(x, y-1) << endl;
    	//cout << "The below neighbour has spin = " << spin(x, y+1) << endl;
    	//cout << "The change in neighbour energy = " << delta_neighbour << endl;
        //cout << "The change in total energy dE =" << delta_E << endl;
        cout << "The number of flips = " << flip_counter <<endl;
        cout << "The number of sweeps = " << sweep_counter <<endl;
        cout << "The final Magnetic Moment = " << M <<endl;
        cout << "The final average Magnetic Moment = " << M_avg <<endl;
    
    	return 0;
    }

  14. #44
    Join Date
    Dec 2012
    Location
    England
    Posts
    2,377

    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.
    All advice is offered in good faith only. You are ultimately responsible for effects of your programs and the integrity of the machines they run on.

  15. #45
    Join Date
    Dec 2012
    Location
    England
    Posts
    2,377

    Re: Ising Model C++ Metropolis Algorithm

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

    See http://msdn.microsoft.com/en-us/libr...=vs.60%29.aspx
    Last edited by 2kaud; April 15th, 2014 at 01:57 PM.
    All advice is offered in good faith only. You are ultimately responsible for effects of your programs and the integrity of the machines they run on.

Page 3 of 4 FirstFirst 1234 LastLast

Tags for this Thread

Posting Permissions

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


Windows Mobile Development Center


Click Here to Expand Forum to Full Width

This is a CodeGuru survey question.


Featured


HTML5 Development Center