Ising Model in C++

[unscientific]

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;

        }



}