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 massdouble epsilon = 1.0;   // LJ energy coefficientdouble sigma = 1.0;     // LJ distance coefficientdouble kB = 1.0;        // Boltzmann constant//// Define global variables//Vector *r;               // particle positionsVector *v;               // particle velocitiesVector *a;               // particle accelerationsVector *pa;              // previous accelerationsint np;                  // number of particlesdouble L;                // cubic simulation box sizedouble Ek;               // total kinetic energydouble Ep;               // total potential energydouble Ecut;             // potential energy contribution beyond cutoff double rc;               // cutoff distance of the potentialint step;                // total simulation stepsint sample_int;          // sampling intervaldouble dt;               // integration time intervalbool scale;              // whether scale temperaturedouble 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 timedouble 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;}