Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Monte Carlo Computational Physics Project

  1. Jun 12, 2010 #1
    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 (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;
    }

           




     
    Please help!
     
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted