cdsubs.c

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

   ASSERT(TagValue!=-1)
   ASSERT(TagValue>=0&&TagValue<NUM_CLAMP_TAG)

   /*  Set condition for this tag group  */
   if (OpenTagFound)
      {
      a->UseClamp [TagValue] = UseClamp;
      a->ClampTemp[TagValue] = InputTemperature;
      a->ClampStep[TagValue] = InputCstep;
      }

   /*  Place selected atom in tag group  */
   NumClamp = 0;
   LOOP (ipart, a->np)
      {

      /*  This particle is selected, set to itag  */
      if (!UseSelect || IS_SELECT(ipart))
         {
			APPLY_CLAMP_TAG(ipart,TagValue)
         }

      /*  This particle is in current clamp set, count it  */
      if (IS_CLAMP_TAG(ipart,TagValue))
         {
         NumClamp++;
         }
      }


   }


/*  Read COR file  */
void read_cor (char *instr)
   {
   long step;
   FILE *fin;
   char *tokstr = strhed (&instr);

   /*  OPEN FILE  */
   fin = fopen (tokstr, "rb");
   if (fin==NULL)
      {
      printf ("ERROR :  Cannot open input file %s.\n", tokstr);
      CleanAfterError();
      }

   /*  READ JUST FIRST STEP  */
   if (*instr==0)
      ReadCorFile (fin, a);
   else
      {
      step = intstrf(&instr);
      do
         {
         ReadCorFile (fin, a);
         }
         while (!feof(fin) && a->step!=step);
      }
   fclose (fin);

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }


                              /* READ CMD  */

void read_cmd (char *instr)
   {
   int nstep = intstrf (&instr);

	/*     Check for needed initializations  */
	if (!a->IsInitializedBox)
		{
		printf
			("ERROR:  Cannot peform CMD command without first setting BOX.\n");
		CleanAfterError();
		}

	if (!a->IsInitializedCoord)
		{
      printf ("ERROR:  Cannot peform CMD command without");
      printf (" first setting particle coordinates.\n");
      CleanAfterError();
      }

	if (!a->IsInitializedMass)
		{
		printf ("ERROR:  Cannot peform CMD command without");
		printf (" first setting all atoms' masses.\n");
		CleanAfterError();
		}

   if (!a->IsInitializedPotential)
      {
      printf ("WARNING:  Potential was not set");
      printf (" before performing molecular dynamics\n");
		IncrementNumberOfWarnings();
		}

	/*  Perform molecular dynamics  */
	runcmd (a, s, nstep);

	/*  Write out ending step  */
	if (s->PrintInfo)
		printf ("***   Current step is %i\n", a->step);

   }

                              /*  READ DTIME  */

void read_damp  (char *instr)
   {
   int ipart;
   BOOLEAN IsDampZero;
   char  *HeadStr    = NULL;
   double Coord[NDIR];

   HeadStr = strhed (&instr);

   if (!strcmpi(HeadStr,"ON"))
      a->UseDamp = TRUE;
   else if (!strcmpi(HeadStr, "OFF"))
      a->UseDamp = FALSE;
   else
      {
      printf ("WARNING:  Unrecognized option to DAMP command.\n");
      IncrementNumberOfWarnings();
      return;
      }

   /*  If damping is turned off, then return now  */
   if (a->UseDamp==FALSE)
      return;

   /*
   If damping requested,
   but no current damping or previous damping values,
   set damping off
   */
   HeadStr = strhed (&instr);
   if (a->Damp==NULL && HeadStr[0]=='\0')
      {
      a->UseDamp = FALSE;
      return;
      }

   /*
   If damping selected, and value specified, but no storage allocated,
   allocate now
   */
   if (a->Damp==NULL)
      {
      ALLOCATE(a->Damp, double, a->NumPartAlloc)
      }


   /*  Loop through particles  */
	CheckForNoneSelected();
   IsDampZero = FALSE;
   LOOP (ipart, a->np)
      {
      if (IS_SELECT(ipart) )
         {

         /*  Set values for X, Y, Z in formula  */
         Coord[X] = 1e8*a->cur[NDIR*ipart+X];
         Coord[Y] = 1e8*a->cur[NDIR*ipart+Y];
         Coord[Z] = 1e8*a->cur[NDIR*ipart+Z];
         EvaluateFormula
            (HeadStr, Coord, 1e+8, &a->Damp[ipart], 1.0, 1);
         }

      /*  Check if all damping coefficients are zero  */
      IsDampZero = ( IsDampZero || (a->Damp[ipart]==0.0) );
      }

   /*  All damping coeffiecients are zero  */
   if (IsDampZero)
      {

      /*  Set damp flag to false  */
      a->UseDamp = FALSE;

      /*  Free damp array  */
      FREE (a->Damp)
      }

   }
                              /*  READ DTIME  */

void read_dtime (char *instr)
   {
   double NewTimeStep;
   double OldTimeStep;
   double TimeScale;
   double CumulativeTimeScale;
   int    i;
   int    j;

   /*  No tim step entered - leave tim step unchanged  */
   if (*instr==0)
      return;

   /*  Read new time step  */
   NewTimeStep = dblstrf (&instr);

   /*  Check time step  */
   if (NewTimeStep<=0)
      {
      printf ("ERROR: DTIME (%f) <= 0\n", NewTimeStep);
      return;
      }

   /*  Save old time step  */
   OldTimeStep    = s->dtime;

   /*  Store new time step  */
   s->dtime = NewTimeStep;

   /*  RESCALE TIME DERIVATIVES  */
   if (OldTimeStep != NewTimeStep)
      {

      /*  Scaling Factor  */
      TimeScale  = NewTimeStep/OldTimeStep;

      /*  Rescale velocities  */
      CumulativeTimeScale = TimeScale;
      if (a->v != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->v[i + 3*j] *= CumulativeTimeScale;

      /*  Rescale accelerations  */
      CumulativeTimeScale *= TimeScale;
      if (a->c2 != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->c2[i + 3*j] *= CumulativeTimeScale;

      /*  Rescale 3rd time derivatives  */
      CumulativeTimeScale *= TimeScale;
      if (a->c3 != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->c3[i + 3*j] *= CumulativeTimeScale;

      /*  Rescale 4th time derivatives  */
      CumulativeTimeScale *= TimeScale;
      if (a->c4 != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->c4[i + 3*j] *= CumulativeTimeScale;

      /*  Rescale 5th time derivatives  */
      CumulativeTimeScale *= TimeScale;
      if (a->c5 != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->c5[i + 3*j] *= CumulativeTimeScale;
      }
   }

                              /*  READ DEBUG  */

void read_debug (char *instr)
   {
   char *HeadStr;

   HeadStr = strhed (&instr);

   if      (!strcmpi(HeadStr, "ON" ))
      s->debug = TRUE;
   else if (!strcmpi(HeadStr, "OFF"))
      s->debug = FALSE;

   else
      printf ("ERROR:  Unrecognized option to DEBUG command.\n");
   Debug_g = s->debug;
   }

                            /*  READ DUP  */

void read_dup      (char *instr)
   {
   int OldNumPart;
   int NewNumPart;
   int idir;
   int idup;
   int iold;
   int inew;
   int indexnew;
   int indexold;
   int ndup;
   int NumWrap;
   double dx;
   double dy;
   double dz;
   double tx;
   double ty;
   double tz;
   double (*Coord)[NDIR] = NULL;

   ndup = intstrf (&instr);
   dx = ANGS*dblstrf (&instr);
   dy = ANGS*dblstrf (&instr);
   dz = ANGS*dblstrf (&instr);
   tx = ty = tz = 0;

   /*  RETURN IF NO DUPLICATIONS SPECIFIED  */
   if (ndup==0)
      return;

   OldNumPart = a->np;
   NewNumPart = a->np + ndup * a->nsel;

   /*  REALLOCATE Particles  */
   ReallocateParticle (a, NewNumPart);

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

   /*  DUPLICATE PARTICLES  */
   inew = OldNumPart;
   LOOP(idup,ndup)
      {
      tx += dx;
      ty += dy;
      tz += dz;
      LOOP(iold,OldNumPart)
      if (IS_SELECT(iold) )
         {
         Coord[inew][X] = Coord[iold][X] + tx;
         Coord[inew][Y] = Coord[iold][Y] + ty;
         Coord[inew][Z] = Coord[iold][Z] + tz;
         a->type[inew] = a->type [iold];
         if (a->mass!=NULL) a->mass[inew] = a->mass[iold];

         /*  Copy velocity and higher derivatives  */
         LOOP (idir, NDIR)
            {
            indexnew = inew * NDIR + idir;
            indexold = iold * NDIR + idir;
            if (a->v  != NULL)  a->v [indexnew] = a->v [indexold];
            if (a->c2 != NULL)  a->c2[indexnew] = a->c2[indexold];
            if (a->c3 != NULL)  a->c3[indexnew] = a->c3[indexold];
            if (a->c4 != NULL)  a->c4[indexnew] = a->c4[indexold];
            if (a->c5 != NULL)  a->c5[indexnew] = a->c5[indexold];
            }
         APPLY_SELECT(inew)
         a->nsel++;
         inew ++;
         }
      }

   /*  Set new number of particles  */
   ASSERT(inew==NewNumPart)
   a->np = NewNumPart;

   /*  Wrap particles   */
   NumWrap = WrapParticles(a);

   if (NumWrap>0)
      PrintWrapWarning (NumWrap);

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }
                           /*  READ ENERGY UNITS  */

void read_eunit  (char *instr)
   {
   double tunit;
   char *tokstr = strhed (&instr);

   /*  NO UNITS SPECIFIED - PRINT OUT CURRENT UNIT  */
   if (!*tokstr)
      {
      printf ("   EUNIT %s %e\n", s->eulabel, s->eunit);
      return;
      }

   /* KELVIN */
   if (!strcmpi(tokstr,"K"))
      {
      strcpy (s->eulabel, "K");
      s->eunit = 1/BK;
      return;
      }

   /* ERGS   */
   if (!strcmpi(tokstr,"ERG"))
      {
      strcpy (s->eulabel, "erg");
      s->eunit = 1;
      return;
      }

   /* JOULES  */
   if (!strcmpi(tokstr,"JOULE"))
      {
      strcpy (s->eulabel, "joule");
      s->eunit = JOULE;
      return;
      }

   /* ELECTRON VOLTS  */
   if (!strcmpi(tokstr,"EV"))
      {
      strcpy (s->eulabel, "eV");
      s->eunit = EV;
      return;
      }

   /* ELECTRON VOLTS  */
   if (!strcmpi(tokstr,"HARTREE"))
      {
      strcpy (s->eulabel, "Hartree");
      s->eunit = HARTREE;
      return;
      }

   /*  USER SPECIFIED UNITS  */
   tunit = dblstrf (&instr);
   if (tunit==0)
      printf ("ERROR (read_eunit): energy unit 0 or unreadable.");
   else
      {
      strcpy (s->eulabel, tokstr);
      s->eunit   = tunit;
      }
   }

                          /*  READ EXTERNAL FORCE  */
void read_extforce (char *InputString)
   {
   double InputForce[NDIR];
   double (*Coord)[NDIR];
   int ipart;
   int idir;
   int NumActiveForces;

   if (!strcmpi(InputString,"CLEAR"))
      {
      /*  Free force constant and origin storage for individual particles  */
      if (a->ForcePtrList != NULL)
         LOOP (ipart, a->np)
            FREE (a->ForcePtrList[ipart])

      /*  Free storage pointer storage of entire list  */
      FREE (a->ForcePtrList);

      return;
      }

#if 0
   /*  Read input force components  */
   LOOP (idir, NDIR)
      InputForce[idir] = dblstrf (&InputString) / s->eunit;

   /*  Set flag if no force  */
   ZeroForce =
      (
      InputForce[X]==0.0 &&
      InputForce[Y]==0.0 &&
      InputForce[Z]==0.0
      );

   /*  If zero force and currently no list, than don't create one  */
   if (ZeroForce && a->ForcePtrList==NULL)
      return;
#endif

   /*  ALLOCATE ARRAY  */
   if (a->ForcePtrList==NULL)
      ALLOCATE (a->ForcePtrList, ExternalVector_t *, a->NumPartAlloc)

   /*
   For each selected particle, allocate and store for force
   (NOTE:  If all remaining force constants are set to zero,
   then list will be removed from memory below)
   */
	CheckForNoneSelected();
   Coord = (double (*)[NDIR]) a->cur;
   NumActiveForces = 0;
   LOOP (ipart, a->np)
      {
      if (IS_SELECT(ipart) )
         {

         EvaluateFormula
            (
            InputString,
            Coord[ipart],
            1e8,
            InputForce,
            1.0,
            NDIR
            );

         /*  Free Storage if applied force constant is zero  */
         if (InputForce[X]==0.0 && InputForce[Y]==0.0 && InputForce[Z]==0.0)
            {
            if (a->ForcePtrList[ipart] != NULL)
               FREE (a->ForcePtrList[ipart]);
            }

         else
            {
            /*  Allocate new storage if needed  */
            if (a->ForcePtrList[ipart] == NULL)
					ALLOCATE (a->ForcePtrList[ipart], ExternalVector_t, 1);

            /*  Copy origin and force constant  */
            LOOP (idir, NDIR)
               {
               a->ForcePtrList[ipart]->Origin   [idir] =
                  Coord[ipart][idir];
               a->ForcePtrList[ipart]->Vector  [idir] =
                  InputForce[idir];
               }
            }
         }

      /*  Count number of active forces  */
      if (a->ForcePtrList[ipart] != NULL)
         NumActiveForces++;
      }

   /*  If no active forces remain, free storage for list  */
   if (NumActiveForces==0)
      {
      FREE (a->ForcePtrList)
      }

   }


                          /*  READ EXTERNAL SPRING  */
void read_extspring (char *InputString)
   {
   double  (*Coord)[NDIR];
   double InputSpringConstant[NDIR];
   int NumActiveSprings;
   int ipart;
   int idir;

   /*  Remove storage for external spring  */
   if (!strcmpi(InputString,"CLEAR"))
      {

      /*  Free spring constant and origin storage for individual particles  */
      if (a->SpringPtrList != NULL)
         LOOP (ipart, a->np)
            FREE (a->SpringPtrList[ipart])

      /*  Free storage pointer storage of entire list  */
      FREE (a->SpringPtrList);

      /*  Return to calling function  */
      return;
      }

#if 0
   /*  Read input spring constant  */
   LOOP (idir, NDIR)
      InputSpringConstant[idir] = dblstrf (&InputString) / s->eunit;

   /*  Set flag if no spring  */
   ZeroSpring =
      (
      InputSpringConstant[X]==0.0 &&
      InputSpringConstant[Y]==0.0 &&
      InputSpringConstant[Z]==0.0