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 constant
double sigma = 1.0;     // LJ distance constant
double kB = 1.0;        // Boltzmann constant


//
// Define global variables
//

Vector *r;               // particle positions

int np;                  // number of particles
double L;                // cubic simulation box size

double Ep;               // total potential energy

double Ecut;             // potential energy contribution beyond cutoff 
double rc;               // cutoff distance of the potential

int step;                // total simulation steps
int sample_int;          // sampling interval
int acc_count = 0;       // acceptance count
double dx;               // maximum displacement
double T;                // temperature
double 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;
}