cdmd.c

一个很好的分子动力学程序

/*
************************************************************************
Compile Switches
************************************************************************
*/
#define CALC_CORRELATION_LENGTH
#undef  CALC_CORRELATION_LENGTH

/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <values.h>
#include <time.h>
#include "cdhouse.h"
#include "cdsubs.h"
#include "cdalloc.h"
#include "cdcons.h"
#include "particle.h"
#include "cdstep.h"
#include "cdvect.h"

#ifdef __TURBOC__
#include <alloc.h>
#endif


/*
************************************************************************
Defines
************************************************************************
*/

/*  Gear coefficients  */
#define  A0    3.0/20.0
#define  A1  251.0/360.0
#define  A3   11.0/18.0
#define  A4    1.0/6.0
#define  A5    1.0/60.0

#define GETKET_TEMP 0
#define GETKET_KE   1

/*
************************************************************************
External Variables
************************************************************************
*/
extern int isig;


/*
************************************************************************
Module-Wide Variables
************************************************************************
*/
static long   NumberOfMDSteps_m = 0;
static double ElapsedTimeMD_m = 0;

/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
static void IntegrateMotion       (Particle_t *, Simulation_t *);
static void CalcBoxDynamics       (Particle_t *, Simulation_t *);
static void ClampVelocity         (Particle_t *, Simulation_t *);
static void ClampVelocityTag      (Particle_t *, Simulation_t *, int);
static void ClampPressure         (Particle_t *);
static void QuenchVelocity        (Particle_t *, Simulation_t *);
static void QuenchAllVelocities   (Particle_t *);
static void QuenchBoxVelocities   (Particle_t *);
static void CheckForSimpleErrors  (Particle_t *, Simulation_t *);
static void GearPredict           (double [NDERIV] );
static void GearCorrect           (double [NDERIV], double, double, double);
static double GetKET       (Particle_t *, Simulation_t *, BOOLEAN, int);
static void CalcThermalStress     (Particle_t *, Simulation_t *);



/*
************************************************************************
Functions
************************************************************************
*/

                  /*  perform computer molecular dynamics  */


void runcmd (Particle_t *a, Simulation_t *s, int nstep)
   {
   int    nfree;
   int    ipart;
   long   istep;
   time_t InitialTime;
   time_t FinalTime;
   double eold;
   double eorg;
   double etmp;
   double eavg =0;
   double evar0=0;
   double evar1=0;

   /*  Initialize time  */
   time (&InitialTime);

   /*  Check for some possible simple errors  */
   CheckForSimpleErrors (a, s);
	
   /*  Allocate storage for particle time derivatives  */
	if (a->f==NULL)
		ALLOCATE(a->f,  double, NDIR*a->NumPartAlloc)
	if (a->v==NULL)
		ALLOCATE(a->v,  double, NDIR*a->NumPartAlloc)
	if (a->c2==NULL)
		ALLOCATE(a->c2, double, NDIR*a->NumPartAlloc)
	if (a->c3==NULL)
		ALLOCATE(a->c3, double, NDIR*a->NumPartAlloc)
	if (a->c4==NULL)
		ALLOCATE(a->c4, double, NDIR*a->NumPartAlloc)
	if (a->c5==NULL)
		ALLOCATE(a->c5, double, NDIR*a->NumPartAlloc)

   /*  Open eligible periodic output files  */
   OpenStepOutput ( &(s->EnergyStepOutput)     );
   OpenStepOutput ( &(s->StressStepOutput)     );
   OpenStepOutput ( &(s->BoxStepOutput)        );
   OpenStepOutput ( &(s->TrajectoryStepOutput) );


   /*  GET ENERGY FOR ORIGIN OF AVERAGE CALCULATIONS  */
   energy_all (a, s, 0);
   eorg = a->epot;

   /*  Calculate number of free atoms  */
   nfree = a->np - a->nfix;

   /*  Test for some free atoms  */
   if (nfree==0)
      {
      printf ("ERROR (runcmd):  No free particles found.\n");
      CleanAfterError();
      }

   /*  TEST FOR EXTERNAL SIGNAL SET DURING LOOP ITERATIONS  */
   for (istep=0;istep<nstep && isig!=1;istep++)
      {

      /*  Increment step counts  */
      a->step++;
      s->nstep++;

      /*  Save previous potential energy  */
      eold = a->epot;

      /*  Integrat equations of motion  */
      IntegrateMotion (a, s);

      /*  Apply Pressure clamp  */
      ClampPressure (a);

      /*  Velocity conditioning routines  */
      if (s->quench==1) {
			if (a->epot>eold)
	         QuenchAllVelocities (a);
		} else if (s->quench==2) {
				QuenchVelocity (a,s);
		} else {
	         ClampVelocity(a,s);
		}

      /*  Box quench  */
#if 0
      if (a->BoxMotionQuench && a->epot>eold)
         QuenchBoxVelocities(a);
#endif

      /*  Total energy  */
      a->etot = a->epot + a->ekin;


		/*  Wrap particles  */
		WrapParticles(a);


      /*  Print Energy Info  */
      if (NowStepOutput( &(s->EnergyStepOutput)) )
			{
         fprintf
            (
            s->EnergyStepOutput.File,
            "%6li %+le %+le %+le %+le\n",
            a->step,
            s->eunit*a->etot/a->np,
            s->eunit*a->epot/a->np,
            s->eunit*a->ekin/a->np,
            s->eunit*(a->EkinBox[X]+a->EkinBox[Y]+a->EkinBox[Z])/a->np
            );
			}


      /*  Print Stress Info  */
      if (NowStepOutput( &(s->StressStepOutput) ) )
            fprintf
               (
               s->StressStepOutput.File,
               "%6li  %+12.5le %+12.5le %+12.5le %+12.5le %+12.5le %+12.5le\n",
               a->step,
               a->Stress[XX]/a->np,
               a->Stress[YY]/a->np,
               a->Stress[ZZ]/a->np,
               a->Stress[YZ]/a->np,
               a->Stress[ZX]/a->np,
               a->Stress[XY]/a->np
               );

      /*  Print Box Info  */
      if (NowStepOutput( &(s->BoxStepOutput) ) )
            fprintf
               (
               s->BoxStepOutput.File,
               "%6li  %13.6e %13.6e %13.6e\n",
               a->step,
               a->bcur[X],
               a->bcur[Y],
               a->bcur[Z]
               );


      /*  Print trajectory  */
      if (NowStepOutput ( &(s->TrajectoryStepOutput) ) )
         {
         LOOP (ipart,a->np)
         if (IS_SELECT(ipart))
            fprintf
               (
               s->TrajectoryStepOutput.File,
               "%+13.6e %+13.6e %+13.6e\n",
               CM*a->cur[NDIR*ipart+X],
               CM*a->cur[NDIR*ipart+Y],
               CM*a->cur[NDIR*ipart+Z]
               );
         }


      /*  CORRELATION INFO  */
      etmp   = a->epot - eorg;
      evar0 += etmp*etmp;
      evar1 += etmp*(eold-eorg);
      eavg  += etmp;
      }

   /*  Close StepOutput files  */
   CloseStepOutput ( &(s->EnergyStepOutput) );
   CloseStepOutput ( &(s->StressStepOutput) );
   CloseStepOutput ( &(s->BoxStepOutput) );
   CloseStepOutput ( &(s->TrajectoryStepOutput) );

   /*  Add to number of total MD steps  */
   NumberOfMDSteps_m += istep;

   /*  Get elapsed time  */
   time (&FinalTime);
   ElapsedTimeMD_m += difftime (FinalTime, InitialTime);

#ifdef CALC_CORRELATION_LENGTH
   if (istep>0)
      {
      double corr1;

      printf ("Number of Dynamics   steps performed %li.\n", istep);
      eavg  /= istep;
      evar0  = evar0 / istep  - eavg*eavg;
      evar1  = evar1 / istep  - eavg*eavg;
      if (evar0>0)
         corr1 = evar1 / evar0;
      if (evar0<=0.0 || corr1>=1.0)
         printf ("Correlation length: infinite or unknown.\n");
      else
         printf ("Estimated correlation length: %lf\n", 1.0 / (1.0 - corr1) );
      }
#endif   /*  CALC_CORRELATION_LENGTH  */
   }


/*  
Calculate kinetic energy according to box motion type  
- calculates a->ekin, a->EkinBoxp[], a->temp
*/
void CalcKineticEnergy (Particle_t *a, Simulation_t *s)
   {
   int ipart;
   int idir;
   int NumFree;
   double Factor;
   double VelComp;
   double (*Position)[NDIR] = (double (*)[NDIR]) a->cur;
   double (*Velocity)[NDIR] = (double (*)[NDIR]) a->v;

   /*  Initialize total kinetic energy  */
   a->ekin = 0.0;

   /*  If no particles then kinetic energy is zero  */
   if (a->np==0)
      return;


   /*  Sum "simple" Kinetic energy  */
   if (a->BoxMotionAlgorithm==BMA_NONE || a->BoxMotionAlgorithm==BMA_CLAMP)
      {
      LOOP (ipart, a->np)
		if (!IS_FIX(ipart))
			{
         a->ekin +=
            a->mass[ipart] *
            (
            Velocity[ipart][X]*Velocity[ipart][X] +
            Velocity[ipart][Y]*Velocity[ipart][Y] +
            Velocity[ipart][Z]*Velocity[ipart][Z]
            );
			}
      a->ekin /= (2.0 * s->dtime * s->dtime);
      }


   /*  Sum Kinetic energy directionally and modified by box motion  */
   else
      {

      /*  Get total kinetic energy (box and particles) in each direction  */
      LOOP (idir, NDIR)
         {

         /*  Get factor for combining particle and box motion  */
         Factor = a->BoxMotion[idir][1] / a->BoxMotion[idir][0];

         /*  Get kinetic energy of box motion  */
         a->EkinBox[idir] =
            a->BoxMass[idir] *
            a->BoxMotion[idir][1] * a->BoxMotion[idir][1];

         /*  Add kinetic energy of particles  */
         LOOP (ipart, a->np)
				if (!IS_FIX(ipart))
            {
            VelComp = Velocity[ipart][idir] - Factor * Position[ipart][idir];
            a->EkinBox[idir] += a->mass[ipart] * VelComp * VelComp;
            }

         /*  Normalize kinetic energy  */
         a->EkinBox[idir] /= (2 * s->dtime * s->dtime);

         /*  Sum kinetic energy into total  */
         a->ekin += a->EkinBox[idir];

         }


      }



   /*  Determine temperature from number of free particles  */
   NumFree = a->np - a->nfix;
   a->temp = a->ekin / (1.5 * NumFree * BK);
   }




void CalcThermalStress (Particle_t *a, Simulation_t *s) {
   int ipart;
	double v[NDIR];
   double (*Velocity)[NDIR] = (double (*)[NDIR]) a->v;

	LOOP (ipart, a->np) {
		v[X] = Velocity[ipart][X]/s->dtime;
		v[Y] = Velocity[ipart][Y]/s->dtime;
		v[Z] = Velocity[ipart][Z]/s->dtime;
		a->Stress[XX] += a->mass[ipart] * v[X] * v[X];
		a->Stress[YY] += a->mass[ipart] * v[Y] * v[Y];
		a->Stress[ZZ] += a->mass[ipart] * v[Z] * v[Z];
		a->Stress[YZ] += a->mass[ipart] * v[Y] * v[Z];
		a->Stress[ZX] += a->mass[ipart] * v[Z] * v[X];
		a->Stress[XY] += a->mass[ipart] * v[X] * v[Y];
	}
	return;
}


/*  
Returns temperature for selected particles, includes effect of
moving box on particle velocities (for Andersen? algorithm), 
but returned temperature does NOT include kinetic energy of box 
17 Jan 2000 - added degree of freedom calculation
*/
double GetTemperature (Particle_t *a, Simulation_t *s, BOOLEAN UseSelect) {
	return GetKET (a, s, UseSelect, GETKET_TEMP);
}
double GetKE          (Particle_t *a, Simulation_t *s, BOOLEAN UseSelect) {
	return GetKET (a, s, UseSelect, GETKET_KE);
}

double GetKET (Particle_t *a, Simulation_t *s, BOOLEAN UseSelect, int mode)
   {
   int ipart;
   int idir;
   int NumDegFree = 0;
	int NumPart = 0;
	double KineticEnergy = 0;
   double Factor[NDIR];
   double VelComp;
   double (*Position)[NDIR] = (double (*)[NDIR]) a->cur;
   double (*Velocity)[NDIR] = (double (*)[NDIR]) a->v;

   /*  If no particles then kinetic energy is zero  */
   if (a->np==0)
      return 0;

   /*  Sum "simple" Kinetic energy  */
   if 
	(
	a->BoxMotionAlgorithm==BMA_NONE ||
	a->BoxMotionAlgorithm==BMA_CLAMP
	)
      {
      LOOP (ipart, a->np)
		if (!IS_FIX(ipart))
		if (!UseSelect || IS_SELECT(ipart))
			{
			NumPart++;
			/*  Find number of degrees of freedom for this atom  */
			if (a->cons && a->cons[ipart]) {
				NumDegFree += (a->cons[ipart]->type==CONSTR_PLANE) ? 2 : 1;
			} else {
				NumDegFree += NDIR;
			}

			KineticEnergy +=
            a->mass[ipart] *
            (
            Velocity[ipart][X]*Velocity[ipart][X] +
            Velocity[ipart][Y]*Velocity[ipart][Y] +
            Velocity[ipart][Z]*Velocity[ipart][Z]
            );
			}
      KineticEnergy /= (2.0 * s->dtime * s->dtime);
      }


   /*  
	 *  Sum Kinetic energy directionally and modified by box motion  
	 */
   else
      {

      /*  Get factor for combining particle and box motion  */
		LOOP (idir, NDIR)
			{
         Factor[idir] = a->BoxMotion[idir][1] / a->BoxMotion[idir][0];
			}

      /*  Correct particle velocities by box motion  */
      LOOP (ipart, a->np)
		if (!IS_FIX(ipart))
		if (!UseSelect || IS_SELECT(ipart))
         {

			NumPart++;
			/*  Find number of degrees of freedom for this atom  */
			if (a->cons && a->cons[ipart]) {
				NumDegFree += (a->cons[ipart]->type==CONSTR_PLANE) ? 2 : 1;
			} else {
				NumDegFree += NDIR;
			}

	      LOOP (idir, NDIR)
				{
            VelComp = 
					Velocity[ipart][idir] - Factor[idir] * Position[ipart][idir];
            KineticEnergy += a->mass[ipart] * VelComp * VelComp;
				}
			}

		KineticEnergy /= (2 * s->dtime * s->dtime);
      }


	/*  If no free particles, temperature is zero  */
	if (NumDegFree==0 || NumPart==0)
		return 0;

	/*  Return just kinetic energy __per_particle__ */
	if (GETKET_KE==mode) 
		return KineticEnergy/NumPart;
   /*  Return temperature from kinetic energy and DOF  */
	else
		return KineticEnergy / (0.5 * NumDegFree * BK);
   }




                     /*  integrate equations of motion  */

void IntegrateMotion  (Particle_t *a, Simulation_t *s)
   {
   int i;
   int j;
   double ax, cx, step,  *m;
   double (*c0)[NDIR];
   double (*c1)[NDIR];
   double (*c2)[NDIR];
   double (*c3)[NDIR];
   double (*c4)[NDIR];
   double (*c5)[NDIR];
   double (*f )[NDIR];

   /*  PREDICTION ALGORITHM  */
   c0 = (double (*)[NDIR]) a->cur;
   c1 = (double (*)[NDIR]) a->v;
   c2 = (double (*)[NDIR]) a->c2;
   c3 = (double (*)[NDIR]) a->c3;
   c4 = (double (*)[NDIR]) a->c4;
   c5 = (double (*)[NDIR]) a->c5;
   f  = (double (*)[NDIR]) a->f;



   LOOP (i, a->np)
   if (!IS_FIX(i))
   LOOP (j, NDIR )
      {
      c0[i][j] += c1[i][j] +  c2[i][j] +  c3[i][j] +  c4[i][j] +   c5[i][j];
      c1[i][j] +=           2*c2[i][j] +3*c3[i][j] +4*c4[i][j] + 5*c5[i][j];
      c2[i][j] +=                       3*c3[i][j] +6*c4[i][j] +10*c5[i][j];
      c3[i][j] +=                                   4*c4[i][j] +10*c5[i][j];
      c4[i][j] +=                                                5*c5[i][j];
      }


   /*  Predict box motion  */
   if (a->BoxMotionAlgorithm == BMA_ANDERSEN)
      {
      a->BoxMotion[X][0] = a->bcur[X];
      a->BoxMotion[Y][0] = a->bcur[Y];
      a->BoxMotion[Z][0] = a->bcur[Z];
      GearPredict (a->BoxMotion[X]);
      GearPredict (a->BoxMotion[Y]);