MC

//// lj_mc.C//// Monte Carlo simulation code for a Lennard-Jones system.//// Written by: Yanting Wang                      October 24, 2009//#include <iostream>#include <fstream>using namespace std;#include <stdlib.h>#include <math.h>#include <time.h>#include "vector.h"//// Define reduced units//double epsilon = 1.0;   // LJ energy constantdouble sigma = 1.0;     // LJ distance constantdouble kB = 1.0;        // Boltzmann constant//// Define global variables//Vector *r;               // particle positionsint np;                  // number of particlesdouble L;                // cubic simulation box sizedouble Ep;               // total potential energydouble Ecut;             // potential energy contribution beyond cutoff double rc;               // cutoff distance of the potentialint step;                // total simulation stepsint sample_int;          // sampling intervalint acc_count = 0;       // acceptance countdouble dx;               // maximum displacementdouble T;                // temperaturedouble beta;             // 1/( kB*T )long seed;               // seed for generating random numbers//// Periodic boundary condition (PBC)//Vector pbc( Vector r ){  //  // PBC correction  //  double hL = L / 2.0;  if( r.x > hL )  r.x -= L;  if( r.x < -hL ) r.x += L;  if( r.y > hL )  r.y -= L;  if( r.y < -hL ) r.y += L;  if( r.z > hL )  r.z -= L;  if( r.z < -hL ) r.z += L;  //  // Check if the vector is now inside the box  //  if( r.x > hL || r.x < -hL )  {    cerr << "r.x = " << r.x << " is out of simulation box." << endl;    exit( -1 );  }    if( r.y > hL || r.y < -hL )  {    cerr << "r.y = " << r.y << " is out of simulation box." << endl;    exit( -1 );  }  if( r.z > hL || r.z < -hL )  {    cerr << "r.z = " << r.z << " is out of simulation box." << endl;    exit( -1 );  }  return r;}//// Random number generator for a uniform distribution//// "Minimmal" random number generator of Park and Miller with Bays-Durham // shuffle and added safeguards. Returns a uniform random deviate between// 0.0 and 1.0 ( exclusive of the endpoint values ). Call with seed a// negative integer to initialize; thereafter, do not alter seed between// successive deviates in a sequence. RNMX should approximate the largest// floating value that is less than 1.//// From Numerical Recipies in C by Press et al. with slight modifications.//double gen_uni_rand(){  const long IA = 16807;  const long IM = 2147483647;  const double AM = 1.0 / IM;  const long IQ = 127773;  const long IR = 2836;  const long NTAB = 32;  const long NDIV = 1 + (IM-1) / NTAB;  const double EPS = 1.2e-7;  const double RNMX = 1.0 - EPS;  static long iy = 0;  static long iv[NTAB];  long k;  int j;  //  // Initialize. Be sure to prevent nSeed = 0.  //  if( seed <= 0 || !iy )  {    if( -seed < 1 )  seed = 1;    else  seed = -seed;    //    // Load the shuffle table ( after 8 warm-ups ).    //    for(j= NTAB+7;j>=0;j-- )    {      k = seed / IQ;      seed = IA * ( seed - k * IQ ) - IR * k;      if( seed < 0 ) seed += IM;      if( j < NTAB ) iv[ j ] = seed;    }    iy = iv[ 0 ];  }  //  // start here when not initializing.  //  k = seed /IQ;  //  // Compute idum=(IA*idum)%IM without overflows by Schrage's method.  //  seed = IA * (seed - k * IQ ) - IR * k;  if( seed < 0 ) seed += IM;  //  // Will be in the range 0..NTAB-1.  //  j = iy / NDIV;  //  // Output previously stored value and refill the shuffle table.  //    double temp, ran;  iy = iv[ j ];  iv[ j ] = seed;  if( (temp = AM*iy) > RNMX ) ran = RNMX;  else ran = temp;  return ran;}//// Calculate potential with the j particle at the p position//double potential( int j, Vector p ){  double V = 0.0;  for( int i=0; i<np; ++i )  {    if( i != j )    {      Vector dr = pbc( r[i] - p );      double d = dr.r();          // modulus of vector dr      if( d < rc )       // within cutoff distance      {double id = sigma / d;double i2 = id * id;double i6 = i2 * i2 * i2;V += 4.0 * epsilon * i6 * ( i6 - 1.0 ) - Ecut;      }    }  }  return V;}//// Make one trial of Monte Carlo movement//void mc_move(){  //  // Randomly pick up one particle  //  int j = int( gen_uni_rand() * np );  //  // Calculate its old potential energy  //  double Eo = potential( j, r[j] );  //  // Make a trial move  //  Vector dr( ( gen_uni_rand() - 0.5 ) * dx,      ( gen_uni_rand() - 0.5 ) * dx,      ( gen_uni_rand() - 0.5 ) * dx );  Vector rn = pbc( r[j] + dr );   double En = potential( j, rn );  //  // Determine if accept this trial move  //  if( En < Eo || gen_uni_rand() < exp( -beta * (En-Eo) ) )  {    r[j] = rn;    Ep += En - Eo;    ++ acc_count;  }}//// System and variable initialization//void init(){  //  // Read initial configuration  //  ifstream fc( "mc_init.xyz" );  if( !fc.is_open() )       // failed to open the file  {    cerr << "File mc_init.xyz can not be opened for reading." << endl;    exit( -4 );  }  fc >> np >> L;  r = new Vector[np];  string pname;                // particle name;  for( int i=0; i<np; ++i )  {    fc >> pname >> r[i];  }  fc.close();  //  // Read simulation parameters  //  ifstream fp( "mc_para.dat" );  if( !fp.is_open() )        // failed to open the file  {    cerr << "File mc_para.dat can not be opened for reading." << endl;    exit( -4 );  }  fp >> step >> sample_int >> dx >> rc >> T;  fp.close();  if( step <= 0 || sample_int <= 0 || dx <= 0.0 || rc <= 0.0 || T <= 0.0 )  {    cerr << "Error: input parameter is less than 0." << endl;    cerr << step << " " << sample_int << " " << dx << " "  << rc << " " << T << endl;    exit( -3 );  }  beta = 1.0 / kB / T;  //  // Determine if rc is valid  //  if( rc > 0.5 * L )  {    cerr << "Error: rc=" << rc  << " is larger than half of the box size L=" << L << endl;    exit( -2 );  }  //  // Calculate Ecut  //  double id = sigma / rc;  double i2 = id * id;  double i6 = i2 * i2 * i2;  Ecut = 4.0 * epsilon * i6 * ( i6 - 1.0 );  //  // Refresh output files  //  ofstream od( "mc_out.dat" );  od.close();  ofstream oc( "mc_out.xyz" );  oc.close();  //  // Initialize seed for the random generator  //  seed = -abs( time( NULL ) );  //  // Calculate initial total potential energy  //  Ep = 0.0;  for( int i=0; i<np; ++i )  Ep += potential( i, r[i] );  Ep /= 2.0;}//// Sample and dump thermodynamic data and configurations//// Input: cstep -- current step//void sample( int cstep ){  //  // Output thermodynamic data  //  ofstream od( "mc_out.dat", ios::app );  od << cstep                      // current simulation step     << " " << Ep / np             // average potential energy per particle     << " " << double( 1.0 * acc_count / cstep )   // acceptance rate     << endl;  od.close();  //  // Output an instantaneous configuration  //  ofstream oc( "mc_out.xyz", ios::app );  oc << np << endl             // number of particles     << L  << endl;            // simulation box size  for( int i=0; i<np; ++i )  {    oc << "He"                 // particle name       << " " << r[i]          // positions       << endl;  }  oc.close();}//// Termination process//void end(){  //  // Release memory allocated for the arrays  //  delete []r;}int main(){  init();  for( int i=0; i<step; ++i )  {    mc_move();    if( i % sample_int == 0 ) sample( i+1 );  }  end();  return 0;}