cdmd.c

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

      GearPredict (a->BoxMotion[Z]);
      a->bcur[X] = a->BoxMotion[X][0];
      a->bcur[Y] = a->BoxMotion[Y][0];
      a->bcur[Z] = a->BoxMotion[Z][0];
      }


   /*  CALCULATE FORCE  */
   ( (PotForce_f *) s->forcepnt ) (a, s); /*  calcforcmp (a, s);  */


   /*  APPLY EXTERNAL FORCE  */
   ApplyExternalForce (a);

   /*  Apply external springs  */
   ApplyExternalSpring (a);

   /*  Apply constraint  */
   ApplyConstraint (a);

   /*  Apply damping  */
   ApplyDamping (a);

	/*  Calculate thermal component of the stress  */
	if (a->stress_thermal)  CalcThermalStress(a,s);

   /*  Calculate Kinetic Energy  */
   CalcKineticEnergy (a,s);

   /*  Calculate dynamics due to box motion  */
   if (a->BoxMotionAlgorithm == BMA_ANDERSEN)
      CalcBoxDynamics   (a,s);

   /*  CORRECTION ALGORITHM  */
   m  = a->mass;
   step = 0.5 * s->dtime * s->dtime;
   LOOP (i, a->np)
   if (!IS_FIX(i))
   LOOP (j, NDIR )
      {
      ax = step * f[i][j] / m[i];
      cx = ax - c2[i][j];
      c0[i][j] += A0*cx;
      c1[i][j] += A1*cx;
      c2[i][j]  = ax;
      c3[i][j] += A3*cx;
      c4[i][j] += A4*cx;
      c5[i][j] += A5*cx;
#if 0
      /* Prevent underflow  */
      if (MINDOUBLE>fabs(c3[i][j]))  c3[i][j] = 0.0;
      if (MINDOUBLE>fabs(c4[i][j]))  c4[i][j] = 0.0;
      if (MINDOUBLE>fabs(c5[i][j]))  c5[i][j] = 0.0;
#endif
      }

   /*  Correct box motion  */
   if (a->BoxMotionAlgorithm == BMA_ANDERSEN)
      {
      GearCorrect
         (a->BoxMotion[X], a->BoxForce[X], a->BoxMass[X], s->dtime);
      GearCorrect
         (a->BoxMotion[Y], a->BoxForce[Y], a->BoxMass[Y], s->dtime);
      GearCorrect
         (a->BoxMotion[Z], a->BoxForce[Z], a->BoxMass[Z], s->dtime);
      a->bcur[X] = a->BoxMotion[X][0];
      a->bcur[Y] = a->BoxMotion[Y][0];
      a->bcur[Z] = a->BoxMotion[Z][0];
      }

   }


/*
Modify to used Clamp Sets (3 Nov 1997)
*/
void ClampVelocity (Particle_t *a, Simulation_t *s)
   {
   int    iclamp;
   double BeforeKE;
   double AfterKE;

   /*  Calculate kinetic energy before clamp  */
   CalcKineticEnergy(a,s);
   BeforeKE = a->ekin;

   /*  Apply velocity clamp  */
   LOOP (iclamp, NUM_CLAMP_TAG)
      {
      if (a->UseClamp[iclamp])
         {
         ClampVelocityTag (a, s, iclamp);
         }
      }

   /*  Calculate kinetic energy after clamp  */
   CalcKineticEnergy(a,s);
   AfterKE = a->ekin;

   /*  Calculate energy lost to heat bath  */
   a->ebath += BeforeKE - AfterKE;
   }


/*  Calcualte total kinetic energy  */
#if 0
/*  Use CalcKineticEnergy() instead  */
double CalcKE (Particle_t *a, Simulation_t *s)
   {
   int    ipart;
   double SumVelocitySqr;
   double (*Velocity)[NDIR];

   /*  Set pointer to velocity array  */
   Velocity = (double (*)[NDIR]) a->v;

   SumVelocitySqr = 0.0;
   LOOP (ipart, a->np)
      {
      SumVelocitySqr +=
         a->mass[ipart] *
         (
         Velocity[ipart][X]*Velocity[ipart][X] +
         Velocity[ipart][Y]*Velocity[ipart][Y] +
         Velocity[ipart][Z]*Velocity[ipart][Z]
         );
      }
   return 0.5 * SumVelocitySqr/(s->dtime*s->dtime);
   }
#endif


/*  Clamp velocities of particles belonging to group "itag"  */
void ClampVelocityTag (Particle_t *a, Simulation_t *s, int itag)
   {
   double (*Velocity)[NDIR];
   double CurrentTemp;
   double ClampedKE;
   double UnClampedKE;
   double NewKE;
   double SumVelocitySqr;
   double VelocityFactor;
	int    NumDegFree;
   int    ipart;

   /*  Return if no clamp implemented  */
   if (a->ClampTemp[itag]<0)
      return;

   /*  Set pointer to velocity array  */
   Velocity = (double (*)[NDIR]) a->v;

   /*  Use temperature for clamp subset of particles  */
	NumDegFree     = 0;
   SumVelocitySqr = 0.0;
   LOOP (ipart, a->np)
   if (IS_CLAMP_TAG(ipart,itag))
	if (!IS_FIX(ipart))
      {
		/*  Find number of degrees of freedom for this atom  */
		if (a->cons && a->cons[ipart]) {
			NumDegFree += (a->cons[ipart]->type==CONSTR_LINE) ? 1 : 2;
		} else {
			NumDegFree += NDIR;
		}

      SumVelocitySqr +=
         a->mass[ipart] *
         (
         Velocity[ipart][X]*Velocity[ipart][X] +
         Velocity[ipart][Y]*Velocity[ipart][Y] +
         Velocity[ipart][Z]*Velocity[ipart][Z]
         );
      }

   /*  Ignore clamp if no particles selected  */
   if (NumDegFree==0)
      return;

   ClampedKE =
      0.5 * SumVelocitySqr / (s->dtime*s->dtime);
   CurrentTemp =
      2.0 * ClampedKE / (NumDegFree * BK);


   /*  Return if current temperature 0.0 or less(?)  */
   if (CurrentTemp <= 0.0)
      return;

   /*  Calculate velocity factor  */
   ASSERT ((a->ClampTemp[itag]/CurrentTemp)>=0)

   VelocityFactor =
      pow
         (
         a->ClampTemp[itag]/CurrentTemp,
         1.0/(2*a->ClampStep[itag])
         );

   LOOP (ipart, a->np)
      {
      if (IS_CLAMP_TAG(ipart,itag))
         {
         Velocity[ipart][X] *= VelocityFactor;
         Velocity[ipart][Y] *= VelocityFactor;
         Velocity[ipart][Z] *= VelocityFactor;
         }
      }

	/* NOTE:  Global energies e->ekin, e->ebath should be set by calling
	 * function
	*/
	}


/*  REMOVE ANY VELOCITIES THAT GO UPHILL IN ENERGY  */
void QuenchVelocity (Particle_t *a, Simulation_t *s)
   {
   double dot;
   double (*f)[NDIR]  = (double (*)[NDIR]) a->f;
   double (*v)[NDIR]  = (double (*)[NDIR]) a->v;
   double *m  = a->mass;
   double eknew = 0.0;
   int i;

   LOOP (i, a->np)
      {
      dot  = f[i][X]*v[i][X] + f[i][Y]*f[i][Y] + f[i][Z]*f[i][Z];
      if (dot<0.0)
         v[i][X] = v[i][Y] = v[i][Z] = 0;
      else
         eknew += m[i] *
            (
            v[i][X]*v[i][X] +
            v[i][Y]*v[i][Y] +
            v[i][Z]*v[i][Z]
            );
      }
   eknew *= 0.5 / (s->dtime*s->dtime);
   /*  CALCULATE EBATH  */
   a->ebath += a->ekin - eknew;
   a->ekin   = eknew;
   }


/*  Check various parameters to make sure theire values are OK  */
void CheckForSimpleErrors (Particle_t *a, Simulation_t *s)
   {
   int ipart;
   double MassLimit;
   BOOLEAN IsFatalError;

   IsFatalError = FALSE;

   if (s->dtime > 1e-10)
      {
      printf ("WARNING:  Time step size (%le) may be too big.\n",
         s->dtime);
      IncrementNumberOfWarnings();
      }

   if (s->dtime < 0.0)
      {
      printf ("WARNING:  Time step size (%le) is negative).\n");
      IncrementNumberOfWarnings();
      }

   if (a->np <= 0)
      {
      printf ("WARNING: ");
      printf ("Number of particles (%i) is less or equal to zero.\n",
         a->np);
      IncrementNumberOfWarnings();
      }

   /*  Check volume settings  */
   if (a->BoxMotionAlgorithm==BMA_ANDERSEN)
      {

      /*  Check volume mass  */
      if (a->BoxMass[X]==0.0 || a->BoxMass[Y]==0.0 || a->BoxMass[Z]==0.0)
         {
			ERROR_PREFIX
         printf
         ("Cannot do volume dynamics with Box mass set to zero.\n");
         IsFatalError = TRUE;
         }

		/*  Check for hydrostatic pressure  */
		if (! (a->Pressure[X]==a->Pressure[Y] && a->Pressure[X]==a->Pressure[Z]))
			{
			ERROR_PREFIX
			printf ("Must do current implementation of Andersen algorithm with\n");
			printf ("hydrostatic pressure (that is, pressures in X,Y,Z directions equal).\n");
			IsFatalError=TRUE;
			}

      /*  Check for volume mass too low  */
      MassLimit = 0.0;
      LOOP (ipart, a->np)
         if (!IS_FIX(ipart))
            MassLimit += a->mass[ipart];

      /*  Divide total mass by ten (arbitrary limit)  */
      MassLimit *= 0.1;
      if
      (
      a->BoxMass[X] < MassLimit ||
      a->BoxMass[Y] < MassLimit ||
      a->BoxMass[Z] < MassLimit
      )
         {
         printf ("WARNING:  Box mass may be too low.\n");
         printf ("Simulation may display unstable energy or other problems.\n");
         IncrementNumberOfWarnings();
         }

      /*  Check repeating boundary conditions  */
      if (a->surf[X] || a->surf[Y] || a->surf[Z])
         {
			ERROR_PREFIX
         printf ("Cannot do volume dynamics with free surfaces on.\n");
         IsFatalError = TRUE;
         }

      }

      /*  Check for PRESSURE CLAMP without repeating boundaries  */
      if (
			a->BoxMotionAlgorithm==BMA_CLAMP && 
			(a->surf[X] && a->surf[Y] && a->surf[Z])
			)
         {
         printf ("ERROR:  Cannot do pressure clamp with free surfaces on.\n");
         IsFatalError = TRUE;
         }

      /*  Fatal Error was found, end program  */
      if (IsFatalError)
         CleanAfterError();


   }


long  GetNumberOfMDSteps(void)
   {
   return NumberOfMDSteps_m;
   }

double GetElapsedMDTime (void)
   {
   return ElapsedTimeMD_m;
   }


void CalcBoxDynamics (Particle_t *a, Simulation_t *s)
   {
   int ipart;
   int idir;
   double Factor;
   double Volume;
   double (*Position)[NDIR] = (double (*)[NDIR]) a->cur;
   double (*Force   )[NDIR] = (double (*)[NDIR]) a->f;

	/*  
	Test for pressures equal (needed for this implentation of Andersen algorithm  
	*/
	if (! (a->Pressure[X]==a->Pressure[Y] && a->Pressure[X]==a->Pressure[Z]))
		{
		ERROR_PREFIX
		printf ("When using Andersen's algorithm pressure must be equal.\n");
		CleanAfterError();
		}

   /*  Skip if no moving box  */
   if (a->BoxMotionAlgorithm!=BMA_ANDERSEN)
      return;

   /*  Calculate force on particles  */
   LOOP (idir,  NDIR)
      {

      /*  Acceleration factor for correcting force  */
      Factor =
         2.0 * a->BoxMotion[idir][2] /
         (s->dtime * s->dtime * a->BoxMotion[idir][0]);

      /*  Apply factor to all particles  */
      LOOP (ipart, a->np)
      if (!IS_FIX(ipart))
         Force[ipart][idir] +=
            a->mass[ipart] * Factor * Position[ipart][idir];

      }

   /*  Calculate force on box        */
   Volume = a->BoxMotion[X][0] * a->BoxMotion[Y][0] * a->BoxMotion[Z][0];
   a->BoxForce[X] = 2 * a->EkinBox[X] + a->Stress[XX] - a->Pressure[X]*Volume;
   a->BoxForce[Y] = 2 * a->EkinBox[Y] + a->Stress[YY] - a->Pressure[Y]*Volume;
   a->BoxForce[Z] = 2 * a->EkinBox[Z] + a->Stress[ZZ] - a->Pressure[Z]*Volume;

   /*  Adjust for isotropic box motion  (average the forces)  */
   if (a->BoxMotionGeometry==BMG_ISOTROPIC)
      {
      a->BoxForce[X] =
      a->BoxForce[Y] =
      a->BoxForce[Z] =
         (
         a->BoxForce[X] +
         a->BoxForce[Y] +
         a->BoxForce[Z]
         ) / 3;
      }

   /*  Adjust box forces by box size  */
   a->BoxForce[X] /= a->BoxMotion[X][0];
   a->BoxForce[Y] /= a->BoxMotion[Y][0];
   a->BoxForce[Z] /= a->BoxMotion[Z][0];


   /*  Calculate box   potential energy  */
   a->EpotBox = (a->Pressure[X]+a->Pressure[Y]+a->Pressure[Z])*Volume/3;

   /*  Add potential to total potential energy  */
   a->epot += a->EpotBox;
   }


void  GearCorrect
(
double Coord[NDERIV],
double Force,
double Mass,
double Dtime
)
   {
   double AccDisp;
   double Corr;

   AccDisp = 0.5 * Dtime * Dtime * Force / Mass;
   Corr    = AccDisp - Coord[2];
   Coord[0] += A0*Corr;
   Coord[1] += A1*Corr;
   Coord[2]  = AccDisp;
   Coord[3] += A3*Corr;
   Coord[4] += A4*Corr;
   Coord[5] += A5*Corr;

   }



void GearPredict (double Coord[NDERIV])
   {
   Coord[0] += Coord[1] +   Coord[2] +   Coord[3] +   Coord[4] +    Coord[5];
   Coord[1] +=            2*Coord[2] + 3*Coord[3] + 4*Coord[4] +  5*Coord[5];
   Coord[2] +=                         3*Coord[3] + 6*Coord[4] + 10*Coord[5];
   Coord[3] +=                                      4*Coord[4] + 10*Coord[5];
   Coord[4] +=                                                    5*Coord[5];
   }



void QuenchAllVelocities (Particle_t *a)
   {
   int i;
   /*  Set particle velocities to zero  */
   LOOP (i, NDIR*a->np)
      {
      a->v [i] = 0.0;
      a->c2[i] = 0.0;
      a->c3[i] = 0.0;
      a->c4[i] = 0.0;
      a->c5[i] = 0.0;
      }

   /*  Set box motion to zero  */
   QuenchBoxVelocities(a);
   }



void QuenchBoxVelocities (Particle_t *a)
   {
   int i;
   int idir;


   /*  Set box motion to zero  */
   if (a->BoxMotionAlgorithm==BMA_ANDERSEN)
      {
      LOOP (idir, NDIR)
         for (i=1; i<NDERIV; i++)
            a->BoxMotion[idir][i] = 0.0;
      }

   }


void ClampPressure(Particle_t *a)
   {
   int    idir;
   int    ipart;
   double Volume;
   double Factor;
   double Scale [NDIR];
   double Strain[NDIR];
   double (*Coord)[NDIR] = (double (*)[NDIR]) a->cur;


   if (a->BoxMotionAlgorithm!=BMA_CLAMP)
      return;
   /*  Calculate volume (to correct stress)  */
   Volume = a->bcur[X]*a->bcur[Y]*a->bcur[Z];

   /*  Factor  to obtain strain from stress  */
   Factor = 1.0 / (3.0 * a->BulkModulus * a->VolClampStep);

   /*
   Stress is internal stress, + means expansion.
   Pressure is external pressure, + means compression
   */
   Strain[X] = Factor * (a->Stress[XX]/Volume - a->Pressure[X]);
   Strain[Y] = Factor * (a->Stress[YY]/Volume - a->Pressure[Y]);
   Strain[Z] = Factor * (a->Stress[ZZ]/Volume - a->Pressure[Z]);

	/*  Correct Strain if using ISOTROPIC geometry  */
	if (a->BoxMotionGeometry==BMG_ISOTROPIC)
		{
		Strain[X] = 
		Strain[Y] = 
		Strain[Z] = 
			(Strain[X] + Strain[Y] + Strain[Z])/3.0;
		}

   /*  Convert from Strain[] to Transformation  */
   Scale[X] = 1.0 + Strain[X];
   Scale[Y] = 1.0 + Strain[Y];
   Scale[Z] = 1.0 + Strain[Z];

   /*  Scale particles  */
   LOOP (idir, NDIR)
   if (!a->surf[idir])
      {
      LOOP (ipart, a->np)
         Coord[ipart][idir] *= Scale[idir];
      a->bcur[idir] *= Scale[idir];
      }

   }