cdcons.c

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

         ScaledPosition[Z]*ScaledPosition[Z];

      /*  Skip if inside cavity  */
      if (Measure <= 1.0)
         continue;


      /*  Calculate distance to the ellipse center  */
      Radius = 
         sqrt 
         (
         Position[X]*Position[X] + 
         Position[Y]*Position[Y] +
         Position[Z]*Position[Z]
         );

      /*  
      Calculate the penetration into the wall
         - NOTE: Penetration defined as 
              The distance from particle to the ellipsoid wall
              along the direction running to the ellipsoid center.
           NOT
              The distance from particle to closest point on
              ellipsoid wall.
      */
      Penetration = Radius * (1.0 - 1.0/sqrt(Measure));
      

      /*  
      Calculate normal to ellipsoid wall at intersection 
      of wall and a line joining particle and ellipsoid center
      */
      MagGradientSqr = 0.0;
      LOOP (idir, NDIR)
         {
         Gradient[idir]  = InvAxis [idir] * ScaledPosition[idir];
         MagGradientSqr += Gradient[idir] * Gradient[idir];
         }
      
		/*
		Penetration determines force magnitude
      Gradient    determines force direction
		*/
      Force = -a->CavitySpring * Penetration / sqrt(MagGradientSqr);
      TotalForce[ipart][X] += Force * Gradient[X];
      TotalForce[ipart][Y] += Force * Gradient[Y];
      TotalForce[ipart][Z] += Force * Gradient[Z];


#ifdef SPHERE
      /*  Skip if inside cavity  */
      if (Radius<= a->CavityRadius)
         continue;

      /*  Calculate penetration into the cavity wall  */
      Penetration = Radius - a->CavityRadius;

      /*  Calculate force normalized by Radius  */
      Force = -a->CavitySpring * Penetration / Radius;
      TotalForce[ipart][X] += Force * Position[X];
      TotalForce[ipart][Y] += Force * Position[Y];
      TotalForce[ipart][Z] += Force * Position[Z];
#endif

      
      }  /*  LOOP(ipart,a->np)  */
   }  


/*  APPLY EXTERNAL FORCE  */
void ApplyExternalForce (Particle_t *a)
   {
   double (*Coord     )[NDIR];
   double (*TotalForce)[NDIR];
   double *ExternalForce;
   double *ForceOrigin;
   double Displacement;
   double PotentialEnergy;
   int    ipart;
   int    idir;

   /*  Return if no external forces  */
   if (a->ForcePtrList == NULL)
      return;

   /*  Set pointer to total force and force origin  */
   Coord      = (double (*)[NDIR]) a->cur;
   TotalForce = (double (*)[NDIR]) a->f;

   /*  Loop over all particles subject external forces  */
   PotentialEnergy = 0.0;
   LOOP (ipart, a->np)
   if (a->ForcePtrList[ipart] != NULL)
      {

      /*  Store pointer to External force and origin   */
      ExternalForce = a->ForcePtrList[ipart]->Vector;
      ForceOrigin   = a->ForcePtrList[ipart]->Origin;

      /*  Loop over all components   */
      LOOP (idir, NDIR)
         {
         TotalForce[ipart][idir] += ExternalForce[idir];
         Displacement             = Coord[ipart][idir] - ForceOrigin[idir];
         PotentialEnergy         -= ExternalForce[idir] * Displacement;
         }

      } /*   if (a->ForcePtrList[ipart] != NULL)  */

   /*  Sum contributions to external energy   */
   a->epot += PotentialEnergy;
   }

/*  Apply external springs  */
void ApplyExternalSpring (Particle_t *a)
   {
   double (*Coord)[NDIR];
   double (*Force)[NDIR];
   double *SpringConstant;
   double PotentialEnergy;
   double SpringMag;
   double Dot;
   int ipart;
   int idir;

   /*  Return if no external springs  */
   if (a->SpringPtrList == NULL)
      return;

   /*  Get pointer to coordinates  */
   Coord = (double (*)[NDIR]) a->cur;
   Force = (double (*)[NDIR]) a->f;


   /*  Initialize potential energy   */
   PotentialEnergy = 0.0;

   /*  Loop over all particles  */
   LOOP (ipart, a->np)

   /*  If spring constants assigned to current particle   */
   if (a->SpringPtrList[ipart] != NULL)
      {

      /*  Calculate force components  */
      LOOP (idir, NDIR)
         {

         /*  Store spring constant pointer */
         SpringConstant = a->SpringPtrList[ipart]->Vector;

         /*  Calculate magnitude squared of spring constant  */
         SpringMag =
            (
            SQR(SpringConstant[X]) +
            SQR(SpringConstant[Y]) +
            SQR(SpringConstant[Z])
            );

         /*  Test for zero spring  */
         if (SpringMag==0)
            continue;

         /*  Get magnitude  */
         SpringMag    = sqrt(SpringMag);

         /*  Calculate Dot Product   */
         Dot = 0.0;
         LOOP (idir, NDIR)
            {
            Dot +=
               SpringConstant[idir] *
               (Coord[ipart][idir] - a->SpringPtrList[ipart]->Origin[idir]);
            }

         /*  Scale dot product by magnitude of spring vector  */
         Dot /= SpringMag;

         /*  Calculate force and Potential */
         PotentialEnergy += SQR(Dot);
         LOOP (idir, NDIR)
            Force[ipart][idir] -= SpringConstant[idir] * Dot;

         }


      }

      /*  Add current potential energy to total  */
      a->epot += 0.5 * SpringMag * PotentialEnergy;
   }



void ApplyDamping (Particle_t *a)
	{
	int ipart;
	int idir;
	double Coeff;
	double (*Force   )[NDIR] = (double (*)[NDIR]) a->f;
	double (*Velocity)[NDIR] = (double (*)[NDIR]) a->v;


	/*  Check if damping is one  */
	ASSERT(!(a->UseDamp && a->Damp==NULL))
	if (a->UseDamp==FALSE)
		return;

	/*  Apply damping  */
	LOOP (ipart, a->np)
		{
		Coeff = a->Damp[ipart] / a->dtime;
		LOOP (idir,  NDIR)
			{
			Force[ipart][idir] -= Coeff * Velocity[ipart][idir];
			}
		}

	}



/*
************************************************************************
Local functions
************************************************************************
*/
void ReadCavity (char *InputStr)
   {
   char *TokenStr;
   
   TokenStr = strhed (&InputStr);
   
   /*  Read sphere option  */
   if (!strcmpi(TokenStr, "SPHERE"))
      ReadCavityEllipsoid (InputStr);

   /*  Read ellipsoid option  */
   else if (!strcmpi(TokenStr, "ELLIPSOID"))
      ReadCavityEllipsoid (InputStr);

   /*  Unrecognized option  */
   else
      {
      printf ("ERROR:  Unknown option (%s) to CAVITY.\n",
         TokenStr);
      CleanAfterError();
      }
   }


void ReadCavityEllipsoid (char *InputStr)
   {
   char *TokenStr;

   TokenStr = strhed (&InputStr);

   if (!strcmpi(TokenStr,"OFF"))
      {
      a->UseCavity = FALSE;
      return;
      }

   a->UseCavity = TRUE;

   a->CavitySpring    = dblstrf (&TokenStr);

   a->CavityCenter[X] = 1e-8*dblstrf (&InputStr);
   a->CavityCenter[Y] = 1e-8*dblstrf (&InputStr);
   a->CavityCenter[Z] = 1e-8*dblstrf (&InputStr);

   a->CavityAxis[X]   = 1e-8*dblstrf (&InputStr);
   a->CavityAxis[Y]   = 1e-8*dblstrf (&InputStr);
   a->CavityAxis[Z]   = 1e-8*dblstrf (&InputStr);

   /*  Handle Sphere case  */
   if (a->CavityAxis[Y]==0.0)  a->CavityAxis[Y] = a->CavityAxis[X];
   if (a->CavityAxis[Z]==0.0)  a->CavityAxis[Z] = a->CavityAxis[Y];

      
   /* Sanity Check  */
   if 
   (
   (a->CavityAxis[X] <= 0.0) ||
   (a->CavityAxis[Y] <= 0.0) ||
   (a->CavityAxis[Z] <= 0.0) 
   )
      {
      printf ("ERROR:  Cavity axis (%le,%le,%le) cannot be less than zero.\n",
         a->CavityAxis[X], a->CavityAxis[Y], a->CavityAxis[Z]);
      CleanAfterError();
      }
   if (a->CavitySpring <= 0.0)
      {
      printf ("ERROR:  Cavity spring (%le) cannot be less than zero.\n",
         a->CavitySpring);
      CleanAfterError();
      }

   }