-
April 9th, 2014, 04:08 PM
#1
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;
}
}
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
-
Forum Rules
|
Click Here to Expand Forum to Full Width
|