MD程序

//
// lj_md.C
//
// Molecular dynamics simulation code for a Lennard-Jones system.
//
// Written by: Yanting Wang                      October 22, 2009
//

#include <iostream>
#include <fstream>
using namespace std;

#include <stdlib.h>

#include "vector.h"


//
// Define reduced units
//

double m = 1.0;         // particle mass
double epsilon = 1.0;   // LJ energy coefficient
double sigma = 1.0;     // LJ distance coefficient
double kB = 1.0;        // Boltzmann constant


//
// Define global variables
//

Vector *r;               // particle positions
Vector *v;               // particle velocities
Vector *a;               // particle accelerations
Vector *pa;              // previous accelerations

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

double Ek;               // total kinetic energy
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
double dt;               // integration time interval
bool scale;              // whether scale temperature
double T = 0.01;         // temperature


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


//
// Update particle positions with the Velocity Verlet algorithm
//

void position()
{
  for( int i=0; i<np; ++i )
  {
    r[i] += v[i] * dt + 0.5 * a[i] * dt * dt;   // Velocity Verlet integration

    r[i] = pbc( r[i] );    // put back into the box if out of boundary
  }
}


//
// Calculate forces and potentials according to the current positions
//

void force()
{
  Ep = 0.0;

  for( int i=0; i<np; ++i )
  {
    pa[i] = a[i];
    a[i].clear();        // set all three components to be 0
  }

  //
  // Calculate pair forces and update the system potential
  //

  for( int i=0; i<np-1; ++i )
  {
    for( int j=i+1; j<np; ++j )
    {
      Vector dr = pbc( r[i] - r[j] );

      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;
	
	// Below is actually f/r to save computational time
	double f = 24.0 * epsilon * i2 * i6 * ( 2.0 * i6 - 1.0 );
	
	a[i] += f * dr / m;
	a[j] -= f * dr / m;
	
	Ep += 4.0 * epsilon * i6 * ( i6 - 1.0 ) - Ecut;
      }
    }
  }
}


//
// Update velocities with the velocity Verlet algorithm
//

void velocity()
{
  Ek = 0.0;

  for( int i=0; i<np; ++i )
  {
    v[i] += 0.5 * dt * ( a[i] + pa[i] );

    Ek += v[i] * v[i];
  }

  Ek *= 0.5;
}


//
// Scale temperature with the isokinetics (Evans) thermostat
//

void scale_T()
{
  double Tc = 2.0 / 3.0 * Ek / np / kB;   // current temperature
  double fs = sqrt( T / Tc );             // scaling factor

  for( int i=0; i<np; ++i ) v[i] *= fs;   // scale velocity of each particle

  Ek *= fs * fs;                          // update kinetic energy
}


//
// System and variable initialization
//

void init()
{
  //
  // Read initial configuration
  //

  ifstream fc( "md_init.xyz" );

  if( !fc.is_open() )       // failed to open the file
  {
    cerr << "File md_init.xyz can not be opened for reading." << endl;
    exit( -4 );
  }

  fc >> np >> L;

  r = new Vector[np];
  v = new Vector[np];
  a = new Vector[np];
  pa = new Vector[np];

  string pname;                // particle name;

  for( int i=0; i<np; ++i )
  {
    fc >> pname >> r[i] >> v[i] >> a[i] >> pa[i];
  }

  fc.close();


  //
  // Read simulation parameters
  //

  ifstream fp( "md_para.dat" );

  if( !fp.is_open() )        // failed to open the file
  {
    cerr << "File md_para.dat can not be opened for reading." << endl;
    exit( -4 );
  }

  fp >> step >> sample_int >> dt >> rc >> scale;

  if( scale ) fp >> T;

  fp.close();

  if( step <= 0 || sample_int <= 0 || dt <= 0.0 || rc <= 0.0 || T <= 0.0 )
  {
    cerr << "Error: input parameter is less than 0." << endl;
    cerr << "Input parameters: " << step << " " << sample_int << " " 
	 << dt << " " << rc << " " << T << endl;

    exit( -3 );
  }

  //
  // 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( "md_out.dat" );
  od.close();

  ofstream oc( "md_out.xyz" );
  oc.close();
}


//
// Sample and dump thermodynamic data and configurations
//
// Input: cstep -- current step
//

void sample( int cstep )
{
  //
  // Output thermodynamic data
  //

  ofstream od( "md_out.dat", ios::app );

  od << cstep * dt             // current simulation time
     << " " << (Ek + Ep) / np  // average total energy per particle
     << " " << Ek / np         // average kinetic energy per particle
     << " " << Ep / np         // average potential energy per particle
     << endl;

  od.close();

  //
  // Output an instantaneous configuration
  //

  ofstream oc( "md_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
       << " " << v[i]          // velocities
       << " " << a[i]          // accelerations
       << " " << pa[i]         // previous accelerations
       << endl;
  }

  oc.close();
}


//
// Termination process
//

void end()
{
  //
  // Release memory allocated for the arrays
  //

  delete []r;
  delete []v;
  delete []a;
  delete []pa;
}


int main()
{
  init();

  for( int i=0; i<step; ++i )
  {
    position();
    force();
    velocity();

    if( scale ) scale_T();

    if( i % sample_int == 0 ) sample( i+1 );
  }

  end();

  return 0;
}