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!

And how the hell are we supposed to tell what the errors are?

Are you assuming we know what the correct values should be?

Re-word your question.

Firstly, thanks for posting neat, organized, well commented code as that is definitely appreciated around here. It's helpful for you to write out your questions instead of leaving them in the comments, though.

In looking over your program, one error that sticks out (there may be more that I missed) is on line 311. Two problems with it: firstly it is an assignment rather than a comparison operator and the second is that it is difficult to compare any floating point value for equality since your representation of 1.0 in memory could be 1.00001 or 0.99998 which would not be equal to 1.0. Instead, establish the absolute value of (p - 1.0) is less than some preset epsilon (say 0.00002).

See if that change affects your problem on line 216 at all.

This article has been dead for over six months. Start a new discussion instead.