# Monte Carlo Computational Physics Project

1. Jun 12, 2010

### hasan_researc

Hi,

I am doing an undergraduate physics project by writing a C++ code to implement the Metropolis algorithm to a simple 2d one component plasma. In short, I have to determine the equilibrium configuration at each temperature by means of the Metropolis algorithm and then compute ensemble averages (such as heat capacity) over the microstates.

I have written the code, but when it comes to interpreting values of the heat capacity, I have no clue what mistakes in the code have lead to the supposedly erroneous values.

Code (Text):

// 2d Monte Carlo simulation of a one component plasma

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

using namespace std;

int rand_no (int, int);                                    // generates a random number in between two integer limits.
double rand_num ( );                                       // generates a decimal random number between 0 and 1.
vector <double> initial ( int , int);                      // initialises the positions of the particles.
int Rand_partf ( int );                                    // chooses a particle at random.
double V_norm ( vector <double> , vector <double> , int ); // calculates the normalised potential energy.
double Vf ( vector <double> , vector<double> , int );      // calculates the potential energy for a given configuration.
double V_charf ( );                                        // calculates the characteristic potential energy.
double pbc_pos (double);                                   // applies pbc to the positions of the particles.
double pbc_dist ( double, double, double, double, int);    // applies pbc to calculate the potential energy.

int main ()
{
srand (time (NULL));

double error;                                          /* error is the limiting fluctuation
which would signal that equilibrium has been reached.*/
cout << "Limiting fluctuation on equilibrium = ";
cin >> error;

ofstream stream1 ("Data sheet T against Cv");
stream1 << "T" << "\t" << "Cv" << endl;
ofstream stream2 ("Data Sheet 2");
stream2 << "T" << "\t" << "x" <<  "\t" << "y" << endl;

for (int i = 0; i < 50; i++)                          // This loop calculates ensemble averages for 0 < T =< 50 (T in K).

{
// The following piece of code will "initialise the particles on a 2d square grid."
int N = 10, mcs = 100000;                          /* N is the number of particles on the square grid.
mcs is the number of Monte Carlo Steps per particle.*/
vector <double>  x(N);
x = initial ( N , mcs );                            // x[n] is the x- component of the n'th particle.
vector <double>  y(N);
y = initial ( N , mcs);                             // y[n] is the y- component of the n'th particle.

// The following piece of code will "calculate the energy of the system."
double V_unnorm = Vf( x , y , N) ;                   // V_unnorm is the actual potential energy.
double V_char = V_charf();                           // V_char is the characteristic potential energy
double V = V_unnorm / V_char;                        // V is the normalized potential energy.
double Vbefore = V;                          /* V is stored in Vbefore because it will be compared to
the potential energy after the particle configuration is altered. */

// The following piece of code will "repeat steps 2 through 7 until the system has settled down to an equilibrium."
double Vbefore1 = 0; // Vbefore1 set to an arbitrary number; will be used later in the condition statement.
double Vafter = 0;   // Vafter set to an arbitrary number; will be used later in the condition statement.
double compare;
double T = ( i + 1 );                         // The simulation will run from T = 1K to 500K.
double invbeta = (  ( 1.38*pow(10.0,-23.0) ) * T  );
double beta = 1 / invbeta;
double gamma = V_char * beta;
double iterations = 0;                          // number of iterations of the loop.

double f;                                       // f is Boltzmann's constant.
stream2 << T << "\t" << iterations << "\t" << V << endl;

do
{
iterations = iterations + 1;

// The following piece of code will "randomly choose a particle".
int Rand_part = rand_no ( 1, N );                  // Rand_part is the randomly chosen particle.

// The following piece of code will "move the particle by a small random offset."
double dx = ( 2*rand_num() - 1 );  // dx is the offset in the x-component of the chosen particle.
double dy = ( 2*rand_num() - 1 );  // dy is the offset in the y-component of the chosen particle.

Rand_part = Rand_part - 1;                // this is done so that the vector subscript does not get out of range.
double x_old = x[Rand_part];              // x_old stores the previous x-component of the chosen particle.
double y_old = y[Rand_part];
x[Rand_part] = x[Rand_part] + dx;         // x[Rand_part] is the new x-component of the particle.
x[Rand_part] = pbc_pos ( x[Rand_part] );  // Pbc applied to x[Rand_part].
y[Rand_part] = y[Rand_part] + dy;
y[Rand_part] = pbc_pos ( y[Rand_part]);

// The following piece of code will "calculate the energy of the system."
double V_unnorm = Vf( x , y , N) ;        // V_unnorm is the actual potential energy.
double V_char = V_charf();                // V_norm is the characteristic potential energy.
double V = V_unnorm / V_char;             // V is the normalized potential energy.
double Vafter = V;                        // V is stored in Vafter because it will be compared to Vbefore.

stream2 << T << "\t" << iterations << "\t" << V << endl;

if ( T = 1)
{ stream2 <<  T << << x << y <<
// The following piece of code will "accept the new microstate unconditionally if the energy decreases" or
// "calculate the boltzmann factor f and generate a new random number r if the energy increases.
// If r =< f," the code will "accept the new microstate, otherwise" it will "retain previous microstate."
if (Vafter > Vbefore)
{
double r = rand_num();
f = exp ( -(V_unnorm * beta) );      // f is Boltzmann's factor.

if (r > f)
{
x[Rand_part] = x_old;
y[Rand_part] = y_old;                // These retain the previous microstate.
}
else
{                                         // These retain the new microstate.
}
}
else
{                                             // These retain the new microstate.
}

double Vbefore1 = Vbefore;
Vbefore = Vafter;                          // this is done in preparation for the next iteration.
compare = Vbefore1 - Vafter;
}
while ( abs(compare) >= error );

stream2 << endl << endl;

// The following piece of code will "calculate the desired physical quantities."
double Z = f;                                 // Z is the partition function.
double ener_f = V_unnorm * f;                 // ener_f is the numerator of the formula for mean energy.
double ener2_f = (V_unnorm*V_unnorm) * f;     // ener2_f is the numerator of the formula for mean (energy squared).

double step_number = 0;                          // number of iterations of the loop.
for ( int i = 0; i < 100; i++)
{
step_number = step_number + 1;

// The following piece of code will "randomly choose a particle".
int Rand_part = rand_no ( 1, N );                  // Rand_part is the randomly chosen particle.

// The following piece of code will "move the particle by a small random offset."
double dx = ( 2*rand_num() - 1 );  // dx is the offset in the x-component of the chosen particle.
double dy = ( 2*rand_num() - 1 );  // dy is the offset in the y-component of the chosen particle.

Rand_part = Rand_part - 1;                // this is done so that the vector subscript does not get out of range.
double x_old = x[Rand_part];              // x_old stores the previous x-component of the chosen particle.
double y_old = y[Rand_part];
x[Rand_part] = x[Rand_part] + dx;         // x[Rand_part] is the new x-component of the particle.
x[Rand_part] = pbc_pos ( x[Rand_part] );  // Pbc applied to x[Rand_part].
y[Rand_part] = y[Rand_part] + dy;
y[Rand_part] = pbc_pos ( y[Rand_part]);

// The following piece of code will "calculate the energy of the system."
double V_unnorm = Vf( x , y , N) ;        // V_unnorm is the actual potential energy.
double V_char = V_charf();                // V_norm is the characteristic potential energy
double V = V_unnorm / V_char;             // V is the normalized potential energy.
double Vafter = V;                        // V is stored in Vafter because it will be compared to Vbefore.

stream2 << T << "\t" << iterations << "\t" << V << endl;

// The following piece of code will "accept the new microstate unconditionally if the energy decreases" or
// "calculate the boltzmann factor f and generate a new random number r if the energy increases.
// If r =< f," the code will "accept the new microstate, otherwise" it will "retain previous microstate."
if (Vafter > Vbefore)
{
double r = rand_num();
f = exp ( -(V_unnorm * beta) );      // f is Boltzmann's factor.

if (r > f)
{
x[Rand_part] = x_old;
y[Rand_part] = y_old;                // These retain the previous microstate.
}
else
{                                         // These retain the new microstate.
}
}
else
{                                             // These retain the new microstate.
}

Z = Z + f;
ener_f = ener_f + (V_unnorm*f);
ener2_f = (V_unnorm*V_unnorm) * f;
}

double Eave = ener_f / Z;
double E2ave = ener2_f / Z;
cout << Eave << E2ave << endl;
double Cv0 =  ( beta/T ) * ( E2ave - (Eave*Eave) );               // Cv0 is the heat capacity at constant volume.
double Cv = Cv0 / ( 10.0 / (6.02*pow(10.0,23.0)) );               // Cv is the molar heat capacity at constant volume.
stream1 << T << "\t" << Cv << endl;
cout << T << "\t" << Cv << endl;
}

return 0;

}

vector <double> initial ( int N , int mcs)
{
vector <double> r(N);
for ( int n = 0; n < N; n++)                 // This loop gives the N particles their random positions.
{
double rcum = 0;                         // rcum is the cumulative sum to be calculated from the mcs iterations.
for (int j = 0; j < mcs; j++)            // Calculates r[n] mcs times per particle.
{
int factor = 10000;                  // 1/factor is the precision to which the positions will be specified.
r[n] = rand_no (0, factor);          // Gives r[n] a number between 1 and factor inclusive.
r[n] = r[n]/factor;                  // Changes r[n] to a number in between 0 and 1 inclusive.
rcum = r[n] + rcum;
}
r[n] = rcum/mcs;                         // Calculates the average r[n] from the mcs iterations.
// ALL r[n] CLOSE TO 0.50000. MEANS THERE'S A PROBLEM.
}
return r;
}

double V_charf ( )
{
double n = pow (10.0, 19.0);              // n is the number density
const double pi  = acos (-1.0);
double base = 3 / ( 4.0 * pi * n );
double expo = 1.0 / 3.0;
double a = pow (base, expo);                // a is the characteristic distance

const double q = - 1.60 * pow (10.0,-19.0);
const double epsilon_naught = 8.85 * pow(10.0,-12.0);
double Vc = ( q * q ) / ( 4.0 * pi * epsilon_naught * a );
return Vc;
}

double Vf ( vector <double> x , vector<double> y , int N)
{
double V = 0;                         /* V is the potential energy of the system.
V is initialised to 0 as the system has just been given N particles.*/
const double q = - 1.60 * pow (10.0,-19.0);
const double epsilon_naught = 8.85 * pow(10.0,-12.0);
const double pi  = acos (-1.0);
for (int i = 0; i < N; i++)             // This loop is an implementation of equation 5b (Page 3 of the Project Outline.)
{
for (int j = 0; j < N; j++)
{
double separation = pbc_dist (x[i], y[i], x[j], y[j], N);
if (separation > 0)
{
V = V + (q*q)/(4*pi*epsilon_naught*separation);
}
}
}
V = 0.5*V;
return V;

}

double pbc_pos ( double r)
{
if ( r < 0.0)
{
r = r + 1.0;
}
if (r >= 1.0)
{
r = r - 1.0;
}
else
{
}
return r;
}

double pbc_dist (double xi , double yi ,double xj, double yj, int N)
{
double separation = sqrt (  ( (xi-xj)*(xi-xj) ) + ( (yi-yj)*(yi-yj) )  );
for ( int s = -1; s < 2; s++)
{
double x_trial = xj + s;
for ( int t = -1; t < 2; t++)
{
double y_trial = yj + t;
double separation_trial = sqrt (  ( (xi-x_trial)*(xi-x_trial) ) + ( (yi-y_trial)*(yi-y_trial) )  );
if ( (separation_trial < separation) && (separation_trial != 0) )
{
separation = separation_trial;
}
else {}
}
}
return separation;
}

int rand_no (int lowest, int largest)     // lowest is the smallest random number that will be generated.
// largest is the maximum random number that will be generated.
{
int range = ( largest - lowest ) + 1;
double a = RAND_MAX;
double b = rand();
double q = b / a;
int r = 0;
if ( q < 1.0)                       // This makes sure r = 2 ( if rand()  = RAND_MAX ) is excluded as it puts the vector subscript out of range later.
{
r = ( range * q ) + 1;        // r is the random number to be sent to the main function.
}
else if ( q = 1.0)
{
r = (range * q);
}
else {}
return r;
}

double rand_num()
{
double a = RAND_MAX;
double b = rand();
double r = b/a;
return r;
}