Monte Carlo Computational Physics Project

V_normed = V_norm ( x , y , N); //...the normalised potential energy. // The following piece of code will "apply the Metropolis algorithm." double a = 0.6; //...initial probability of acceptance. int counter = 0; //...count number of accepted moves. for (int j = 1; j <= mcs; j++) { int partf = Rand_partf (N); //...randomly choose a particle. double x_old = x[partf]; //...save the old position of the particle. double y_old = y[partf];
  • #1
hasan_researc
170
0
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.


Code:
// 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!
 

Related to Monte Carlo Computational Physics Project

What is Monte Carlo Computational Physics Project?

Monte Carlo Computational Physics Project is a method of using random numbers and probabilistic techniques to solve complex physics problems that are difficult to solve analytically. It involves simulating physical systems through repeated random sampling and statistical analysis.

How is Monte Carlo Computational Physics Project used in scientific research?

Monte Carlo Computational Physics Project is used in a variety of scientific research fields, including physics, chemistry, biology, and engineering. It is particularly useful for simulating complex systems with many variables and interactions, such as molecular dynamics, material properties, and quantum systems. It is also used for data analysis and optimization problems.

What are the benefits of using Monte Carlo Computational Physics Project?

One of the main benefits of using Monte Carlo Computational Physics Project is that it can solve complex problems that are difficult or impossible to solve analytically. It also allows for the simulation of physical systems that are too large or too small to be studied experimentally. It can also provide insights and predictions for systems that are not fully understood or have unknown properties.

What are the limitations of Monte Carlo Computational Physics Project?

One limitation of Monte Carlo Computational Physics Project is that it relies on random sampling and statistical analysis, which can introduce errors and uncertainties into the results. It also requires significant computational resources and can be time-consuming, especially for complex systems. Additionally, it may not be suitable for all types of physical systems and may not always provide accurate results.

What are some real-world applications of Monte Carlo Computational Physics Project?

Monte Carlo Computational Physics Project has many real-world applications, including in the fields of material science, drug discovery, climate modeling, and financial analysis. It is also used in industries such as aerospace, energy, and transportation for optimizing designs and predicting performance. Additionally, it is used in game theory and risk analysis for decision-making and strategy development.

Similar threads

  • Engineering and Comp Sci Homework Help
Replies
2
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
4
Views
2K
  • Programming and Computer Science
Replies
1
Views
772
  • Programming and Computer Science
Replies
3
Views
1K
  • Programming and Computer Science
Replies
1
Views
659
  • Engineering and Comp Sci Homework Help
Replies
1
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
6
Views
2K
  • Programming and Computer Science
Replies
5
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
2
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
2K
Back
Top