Need to evolve magnetisation and energy per spin with my monte carlo steps

[wsroyle]

Hey guys, I could do with some help to get my program working. In the int main section of my program I could do with some help to evolve my program with monte carlo steps so the magnetisation and energy per spin changes with step. I think I need to add something in the last for loop.
I use microsoft visual studio 2010.

// Ising Model in 2D
// Nearest Neighbours only

#include <cmath>
#include <cstdlib>
#include <iostream>
#include <fstream>

using namespace std;

inline double std_rand()
{
return rand() / (RAND_MAX + 1.0);
}

// define and intialise variables

double J = +1; // ferromagnetic coupling
int Lx, Ly; // number of spins in and y
int N; // number of spins
int **s; // the spins
double T; // temperature
double H; // External Field

//double prob[5];
double w[17][3]; //Boltzman Factors

void computeBoltzmannFactors () {

/*int i;
for (int i = 2; i < 5; i += 2){
prob[i] = exp((-2*i)/T);
}*/

for (int i = -8; i <= 8; i += 4){

w[i + 8][0] = exp( - (i*J + 2*H) / T); //w = exp(-E/kT)
w[i + 8][2] = exp( - (i*J - 2*H) / T);
}
}

int steps = 0; // initalise steps, number of steps so far

void initialise () {
s = new int*[Lx];

for (int i = 0;i < Lx; i++)
s[i] = new int[Ly];

for (int i = 0;i < Lx; i++)
for (int j = 0; j < Ly; j++)

s[i][j] = (std_rand() < 0.5) ? +1 : -1; //hot start. ? is conditional operator, condition ? result1 : result2

computeBoltzmannFactors();
steps = 0;

}

// One of N soins chosen at random, Metropolis method applied. If spin is at boundary
// periodic boundary conditions are applied

bool MetropolisStep(){ //true or false data type

//choose random spin

int i = int(Lx*std_rand());
int j = int(Ly*std_rand());

//find its neighbours using periodic boundary conditions

int iPrev = i == 0 ? Lx - 1 : i - 1;
int iNext = i == Lx - 1 ? 0 : i + 1;
int jPrev = j == 0 ? Ly - 1 : j - 1;
int jNext = j == Ly - 1 ? 0 : j + 1;

//find sum of neighbours

int sumNeighbours = s[j][iPrev] + s[j][iNext] + s[i][jPrev] + s[i][jNext];
int delta_ss = 2*s[i][j]*sumNeighbours;

//ratio of Boltzmann factors

double ratio = w*[delta_ss + 8][1 + s[i][j]];
if (std_rand() < ratio) {
s[i][j] = -s[i][j];
return true;
}
else return false;
}

//take at least N Metropolis Steps to generate next config, gives opportunity to change state

double acceptanceRatio;

void oneMonteCarloStepPerSpin () {
int accepts = 0;
for (int i = 0; i < N; i++)
if (MetropolisStep())
++accepts;
acceptanceRatio = accepts / double (N);
++steps;
}

//compute average magnetisation per spin of lattice

double magnetisationPerSpin() {
int sSum = 0;
for (int i = 0; i < Lx; i++)
for (int j = 0; j < Ly; j++){

sSum += s[i][j];

}

return sSum / double (N);
}

//energy per spin, accounting for boundary conditons

double energyPerSpin () {
int sSum = 0, ssSum = 0;
for (int i = 0; i < Lx; i++)
for (int j = 0; j < Ly; j++){

sSum += s[i][j];

int iNext = i == Lx - 1 ? 0 : i + 1;
int jNext = j == Ly - 1 ? 0 : j + 1;

ssSum += s[i][j]*(s[iNext][j] + s[i][jNext]);
}

return - (J*ssSum + H*sSum)/N;
}

//main function, user input

int main(int argc, char *argv[]){

cout << " Two Dimensional Ising Model\n"
<< "----------------------------\n"
<< " Enter Number of spins L in each direction: ";
cin >> Lx;
Ly = Lx;
N = Lx*Ly;
cout << " Enter Temperature T: ";
cin >> T;
cout << " Enter Magnetic Field H: ";
cin >> H;
cout << " Enter number of Monte Carlo Steps: ";
int MCSteps;
cin >> MCSteps;
initialise ();

for (int i = 0; i < Lx; i++)
for (int j = 0; j < Ly; j++)

cout << s[i][j] << endl;

//perform thermalization steps first, allows system to reach thermal equilibrium at specified temperature

int thermSteps = int(0.2 * MCSteps);
cout << " Performing " << thermSteps
<< " Steps to thermalise ..." << endl;

for (int s = 0; s < thermSteps; s++)
oneMonteCarloStepPerSpin();

//production steps, take data, compute averages.

cout << "Done\n Performing Production Steps ..." << endl;

double mAv = 0, m2Av = 0, eAv = 0, e2Av = 0;

//ofstream fileuni("H:\\4th Year\\Project\\Ising Data.csv");
ofstream filehome("C:\\Users\\Will\\Documents\\University\\"Ising Data.csv")

Physics\\Ising Data.csv");

for (int s = 0; s <MCSteps; s++)
{

double m = magnetisationPerSpin();
double e = energyPerSpin();

mAv += m; m2Av = m*m;
eAv += e; e2Av = e*e;

filehome << m << '\t' << e << '\n';

}

filehome.close();

mAv /= MCSteps; m2Av /= MCSteps;
eAv /= MCSteps; e2Av /= MCSteps;

cout << "\n Magnetisation and Energy per spin is written in file " << endl;
cout << "\\ising.data\\" << endl;

cout << " <m> = " << mAv << endl;
cout << " <e> = " << eAv << endl;

system("PAUSE");

return 0;

}