1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

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



Similar Discussions: Monte Carlo Computational Physics Project
  1. Monte-Carlo Method (Replies: 2)

Loading...