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.
Does anyone have any clue at all? Please reply.
// 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;
}Please help!
