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;
}