# Ising Model C++ Metropolis Algorithm

Page 1 of 1

## 4 Replies - 4593 Views - Last Post: 09 April 2014 - 03:18 PMRate Topic: //<![CDATA[ rating = new ipb.rating( 'topic_rate_', { url: 'http://www.dreamincode.net/forums/index.php?app=forums&module=ajax&section=topics&do=rateTopic&t=344377&amp;s=83f26c1705888a4c9db4808c169386e7&md5check=' + ipb.vars['secure_hash'], cur_rating: 0, rated: 0, allow_rate: 0, multi_rate: 1, show_rate_text: true } ); //]]>

### #1 unscientific

• New D.I.C Head

Reputation: 0
• 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

• Saucy!

Reputation: 6246
• 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.
```

### #3 unscientific

• New D.I.C Head

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

## Re: Ising Model C++ Metropolis Algorithm

Posted 09 April 2014 - 02:22 PM

JackOfAllTrades, 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?

### #4 unscientific

• New D.I.C Head

Reputation: 0
• 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;

}

}

```

### #5 unscientific

• New D.I.C Head

Reputation: 0
• 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)?