4 Replies - 4373 Views - Last Post: 09 April 2014 - 03:18 PM Rate Topic: -----

#1 unscientific  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 4
  • Joined: 09-April 14

Ising Model C++ Metropolis Algorithm

Posted 09 April 2014 - 02:02 PM

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

1. Generate random NxN lattice, with each site either +1 or -1 value.

2. Select a site at random

3. 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.

4. Do like 100 sweeps.

I'm having trouble with steps 2 and 3, would appreciate any help!

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

        }



}



Is This A Good Question/Topic? 0
  • +

Replies To: Ising Model C++ Metropolis Algorithm

#2 JackOfAllTrades  Icon User is offline

  • Saucy!
  • member icon

Reputation: 6246
  • View blog
  • Posts: 24,014
  • Joined: 23-August 08

Re: Ising Model C++ Metropolis Algorithm

Posted 09 April 2014 - 02:18 PM

"Having trouble" is too non-specific.

I do note you're using drand48(), but using srand(). You should use srand48() instead.

And this code is a WTF:

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;



Are you just randomly cutting and pasting int delta_neighbour?

Here's a list of warnings and errors to fix:
g++ -Wall -pedantic -ometro metro.cpp
metro.cpp:19:19: warning: variable length arrays are a C99 feature
      [-Wvla-extension]
  int matrix[ROWS][COLS];
                  ^
metro.cpp:19:13: warning: variable length arrays are a C99 feature
      [-Wvla-extension]
  int matrix[ROWS][COLS];
            ^
metro.cpp:41:14: error: use of undeclared identifier 'N'
if(x<0) x += N;
             ^
metro.cpp:42:15: error: use of undeclared identifier 'N'
if(x>=L) x -= N;
              ^
metro.cpp:43:14: error: use of undeclared identifier 'N'
if(y<0) y += N;
             ^
metro.cpp:44:15: error: use of undeclared identifier 'N'
if(y>=L) y -= N;
              ^
metro.cpp:51:11: error: use of undeclared identifier 'spin'
    {  if(spin(i,j)==spin(i+1, j))     // count from each spin to the ri...
          ^
metro.cpp:51:22: error: use of undeclared identifier 'spin'
    {  if(spin(i,j)==spin(i+1, j))     // count from each spin to the ri...
                     ^
metro.cpp:55:15: error: use of undeclared identifier 'spin'
           if(spin(i, j)==spin(i, j+1))
              ^
metro.cpp:55:27: error: use of undeclared identifier 'spin'
           if(spin(i, j)==spin(i, j+1))
                          ^
metro.cpp:61:10: error: use of undeclared identifier 'J'
    E = -J*neighbour - B*M;
         ^
metro.cpp:67:19: warning: expression result unused [-Wunused-value]
int delta_M = -2*(x, y);    //calculate change in Magnetization M
                  ^
metro.cpp:68:25: warning: expression result unused [-Wunused-value]
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
                       ~^~
metro.cpp:68:34: error: use of undeclared identifier 'spin'
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
                                 ^
metro.cpp:68:48: error: use of undeclared identifier 'spin'
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
                                               ^
metro.cpp:68:63: error: use of undeclared identifier 'spin'
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
                                                              ^
metro.cpp:69:5: error: redefinition of 'delta_neighbour'
int delta_neighbour = -2*(x,y)* int delta_neighbour;
    ^
metro.cpp:68:5: note: previous definition is here
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
    ^
metro.cpp:69:27: warning: expression result unused [-Wunused-value]
int delta_neighbour = -2*(x,y)* int delta_neighbour;
                          ^
metro.cpp:69:37: error: expected
      '(' for function-style cast or type construction
int delta_neighbour = -2*(x,y)* int delta_neighbour;
                                ~~~ ^
metro.cpp:71:19: error: use of undeclared identifier 'J'
double delta_E = -J*delta_neighbour -B*delta_M;
                  ^
metro.cpp:76:9: warning: expression result unused [-Wunused-value]
    {  (x, y) *= -1;     // flip spin and update values
        ^
metro.cpp:86:1: error: expected '}'
^
metro.cpp:6:1: note: to match this '{'
{
^
6 warnings and 16 errors generated.

Was This Post Helpful? 0
  • +
  • -

#3 unscientific  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 4
  • Joined: 09-April 14

Re: Ising Model C++ Metropolis Algorithm

Posted 09 April 2014 - 02:22 PM

View PostJackOfAllTrades, on 09 April 2014 - 02:18 PM, said:

"Having trouble" is too non-specific.

I do note you're using drand48(), but using srand(). You should use srand48() instead.

And this code is a WTF:

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;



Are you just randomly cutting and pasting int delta_neighbour?

Here's a list of warnings and errors to fix:
g++ -Wall -pedantic -ometro metro.cpp
metro.cpp:19:19: warning: variable length arrays are a C99 feature
      [-Wvla-extension]
  int matrix[ROWS][COLS];
                  ^
metro.cpp:19:13: warning: variable length arrays are a C99 feature
      [-Wvla-extension]
  int matrix[ROWS][COLS];
            ^
metro.cpp:41:14: error: use of undeclared identifier 'N'
if(x<0) x += N;
             ^
metro.cpp:42:15: error: use of undeclared identifier 'N'
if(x>=L) x -= N;
              ^
metro.cpp:43:14: error: use of undeclared identifier 'N'
if(y<0) y += N;
             ^
metro.cpp:44:15: error: use of undeclared identifier 'N'
if(y>=L) y -= N;
              ^
metro.cpp:51:11: error: use of undeclared identifier 'spin'
    {  if(spin(i,j)==spin(i+1, j))     // count from each spin to the ri...
          ^
metro.cpp:51:22: error: use of undeclared identifier 'spin'
    {  if(spin(i,j)==spin(i+1, j))     // count from each spin to the ri...
                     ^
metro.cpp:55:15: error: use of undeclared identifier 'spin'
           if(spin(i, j)==spin(i, j+1))
              ^
metro.cpp:55:27: error: use of undeclared identifier 'spin'
           if(spin(i, j)==spin(i, j+1))
                          ^
metro.cpp:61:10: error: use of undeclared identifier 'J'
    E = -J*neighbour - B*M;
         ^
metro.cpp:67:19: warning: expression result unused [-Wunused-value]
int delta_M = -2*(x, y);    //calculate change in Magnetization M
                  ^
metro.cpp:68:25: warning: expression result unused [-Wunused-value]
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
                       ~^~
metro.cpp:68:34: error: use of undeclared identifier 'spin'
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
                                 ^
metro.cpp:68:48: error: use of undeclared identifier 'spin'
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
                                               ^
metro.cpp:68:63: error: use of undeclared identifier 'spin'
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
                                                              ^
metro.cpp:69:5: error: redefinition of 'delta_neighbour'
int delta_neighbour = -2*(x,y)* int delta_neighbour;
    ^
metro.cpp:68:5: note: previous definition is here
int delta_neighbour = (x-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
    ^
metro.cpp:69:27: warning: expression result unused [-Wunused-value]
int delta_neighbour = -2*(x,y)* int delta_neighbour;
                          ^
metro.cpp:69:37: error: expected
      '(' for function-style cast or type construction
int delta_neighbour = -2*(x,y)* int delta_neighbour;
                                ~~~ ^
metro.cpp:71:19: error: use of undeclared identifier 'J'
double delta_E = -J*delta_neighbour -B*delta_M;
                  ^
metro.cpp:76:9: warning: expression result unused [-Wunused-value]
    {  (x, y) *= -1;     // flip spin and update values
        ^
metro.cpp:86:1: error: expected '}'
^
metro.cpp:6:1: note: to match this '{'
{
^
6 warnings and 16 errors generated.


Thanks alot for replying! I forgot to define it at the beginning - I'll change that. delta_neighbour measures the energy the site shares with its neighbours - the product of it and its neighbours.

I'm having trouble with steps 1, 2 and 3. For step 1, I managed to create and display a lattice, but I can't seem to extract the value of a site at location (x, y). Steps 2 and 3, how do I use a boolean expression of some sort to flip according to acceptance probability?
Was This Post Helpful? 0
  • +
  • -

#4 unscientific  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 4
  • Joined: 09-April 14

Re: Ising Model C++ Metropolis Algorithm

Posted 09 April 2014 - 02:30 PM

I have updated the code, defining N and delta_neighbour and the section of the code- how do i make the expression "spin(x, y)" extract the value of that site at location (x, y)?

#include <cstdlib>
#include <ctime>
using namespace std;
#include <iostream>
int main() //random generation of spin configuration
{
int L;				//Total number of spins L = NxN
int N = 30			//A square lattice of length 30
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(srand48()*L);	//retrieves spin from randomly choosen site
int y = int(srand48()*L);

int delta_M = -2*spin(x, y);	//calculate change in Magnetization M
int delta_neighbour = (spinx-1, y) + spin(x+1, y)+ spin(x, y-1) + spin(x, y+1);
int delta_neighbour = -2*spin(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;
       
        }

	

}


Was This Post Helpful? 0
  • +
  • -

#5 unscientific  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 4
  • Joined: 09-April 14

Re: Ising Model C++ Metropolis Algorithm

Posted 09 April 2014 - 03:18 PM

Ok, I've decided that this code is too messy. I shall start off by working on these 2 parts first:


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

Part 1:
//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



Part 2:

//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)?
Was This Post Helpful? 0
  • +
  • -

Page 1 of 1