cdsubs.c

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

      );

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


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

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

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

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

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

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

            }

         }

      /*  Count number of active springs  */
      if (a->SpringPtrList[ipart] != NULL)
         NumActiveSprings++;
      }

   /*  If no active springs remain, free storage for list  */
   if (NumActiveSprings==0)
      {
      FREE (a->SpringPtrList)
      }

   }

                           /*  READ FIX FLAG    */
void read_fix (char *instr)
   {
   int i;
   char *tokstr = NULL;

   tokstr = strhed (&instr);

   if (!strcmpi(tokstr,"ON"))
      {

      /*  
		Make sure that fix is allocated  
		(Am pretty sure this ALLOCATION unnecessary, because
      Selection is allocated during program inialization May 6, 1998)
		*/
      if (a->Selection==NULL)
         ALLOCATE(a->Selection,SEL_T,a->NumPartAlloc)

      /*  Set selected particles as fixed  */
		CheckForNoneSelected();
      LOOP(i,a->np)
         {

         /*  If selected ..  */
         if (IS_SELECT(i) )
            {

            /*  Set fix flag  */
            APPLY_FIX(i)

            /*  Zero derivatives  */
            ZeroDerivatives (i, a);
            }
         }
      }
   else if (!strcmpi(tokstr,"OFF"))
      {
		CheckForNoneSelected();
      if (a->Selection!=NULL)
         {
         LOOP (i, a->np)
            {
            if (IS_SELECT(i) )
               REMOVE_FIX(i)
            }
         }
      }
   else
      {
      printf ("WARNING:  Unrecognized  option to FIX command.\n");
      IncrementNumberOfWarnings();
      }

   /*  Count number of fixed particles */
   a->nfix = 0;
   if (a->Selection!=NULL)
      {
      LOOP (i, a->np)
         {
         if (IS_FIX(i))
            a->nfix++;
         }
      }

	/*  Update number of degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }
                           /*  READ ECHO FLAG   */
void read_echo (char *instr)
   {
   char *tokstr = NULL;

   tokstr = strhed (&instr);

   if      (!strcmpi(tokstr,"ON" )) s->echo = TRUE;
   else if (!strcmpi(tokstr,"OFF")) s->echo = FALSE;
   }


/*  Erase file  */
void read_erase (char *instr)
   {
   int result;

   result = remove(instr);
   if (result!=0)
      {
      printf ("WARNING:  Could not erase file %s\n", instr);
      IncrementNumberOfWarnings();
      }

   }

                           /*  READ ESAVE FLAG  */

void read_esave  (char *instr)
   {
   ReadStepOutput ( &(s->EnergyStepOutput), instr);
   }

                               /*  READ ITEMP  */

/*
NOTE:
To intialize Temp for a subset of particle,
use FIX command prior to calling ITEMP.
*/
void read_itemp (char *instr)
   {
   double Temp;
   char *tokstr = NULL;
   BOOLEAN  xflag;
   BOOLEAN  yflag;
   BOOLEAN  zflag;
   BOOLEAN  OptionRead;
   BOOLEAN  UseSelect      = FALSE;
   BOOLEAN  UseOldVelocity = FALSE;
   int  ndim;


   /*  Check for SEL and OLD options  */
   do
      {
      OptionRead = FALSE;
      tokstr = strhed (&instr);
      if (!strcmpi(tokstr, "SEL"))
         {
         UseSelect = TRUE;
         OptionRead = TRUE;
			CheckForNoneSelected();
         }
      if (!strcmpi(tokstr, "OLD"))
         {
         UseOldVelocity = TRUE;
         OptionRead = TRUE;
         }
      }
      while (OptionRead);

   /*  Read input temperature  */
   Temp = dblstrf (&tokstr);

   /*  Find Dimensions to Initialize  */
   if (*instr==0)
      {
      xflag = TRUE;
      yflag = TRUE;
      zflag = TRUE;
      ndim  = 3;
      }

   else
      {
      xflag = FALSE;
      yflag = FALSE;
      zflag = FALSE;
      tokstr = strhed (&instr);
      while (*tokstr!=0)
         {
         if      (!strcmpi(tokstr, "X"))   xflag = TRUE;
         else if (!strcmpi(tokstr, "Y"))   yflag = TRUE;
         else if (!strcmpi(tokstr, "Z"))   zflag = TRUE;
         else
            {
            printf ("ERROR IN TEMPERATURE INITIALIZATION.\n");
            printf ("    Invalid dimension specified \"%s\".\n", tokstr);
            }
         tokstr = strhed (&instr);
         }
      ndim = 0;
      if (xflag) ndim++;
      if (yflag) ndim++;
      if (zflag) ndim++;
      if (ndim==0)
         {
         printf ("ERROR (read_itemp) IN TEMPERATURE INITIALIZATION.\n");
         printf ("    No valid dimension specified.\n");
         CleanAfterError();
         }
      }

   /*  Call routine to initialize velocities  */
   InitialTemp (a, s, Temp, UseSelect, UseOldVelocity, xflag,yflag,zflag);
   }


                               /*  READ LABEL  */
void read_label (char *instr)
{
   int i,ntitle;
   char bufstr[NTITLESTR];

   ntitle = intstrf(&instr);
   i = 0;
   while (   i<ntitle
          && NULL!=m_gets_list(bufstr,NTITLESTR,inlist))
   {
      if (i<8)
         strcpy (a->title[i], bufstr);
      i++;
   }
}


                           /*  READ MASS  */

void read_mass (char *instr)
   {
   int ipart;
   double Mass;

   Mass = dblstrf (&instr) / AVAG;
   a->IsInitializedMass = TRUE;

	CheckForNoneSelected();
   LOOP (ipart, a->np)
      {

      /*  Assign mass to selected particles  */
      if (IS_SELECT(ipart) )
         a->mass[ipart] = Mass;

      /*  Test if all masses initialized  */
      a->IsInitializedMass =
         (a->IsInitializedMass && a->mass[ipart]!=0.0);
      }

   }

                           /*  READ MC  */

void read_mc (char *instr) {
   int nstep;
   double temp;
   nstep = intstrf (&instr);
   temp  = dblstrf (&instr);
   runmc (a, s, nstep, temp);
}

                           /*  READ MOVE  */

void read_move     (char *instr)
   {
   int    ipart;
   double (*Coord)[NDIR] = (double (*)[NDIR]) a->cur;
   double Values[NDIR];


	CheckForNoneSelected();
   LOOP (ipart, a->np)
   if (IS_SELECT(ipart) )
      {

      /*  Map formula string into Values[] array  */
      EvaluateFormula (instr, Coord[ipart], 1e+8, Values, 1e-8, NDIR);

/*test*/
#if 0
{
char  buffer[160];
char* tokstr=buffer;

sprintf(tokstr, "X=%le Y=%le Z=%le  X Y Z", 0.1, 0.1, 0.1);
dblstrf(&tokstr);
dblstrf(&tokstr);
dblstrf(&tokstr);
Values[X] = dblstrf(&tokstr)*1e-8;
Values[Y] = dblstrf(&tokstr)*1e-8;
Values[Z] = dblstrf(&tokstr)*1e-8;
}
#endif
/*end test*/

      /*  Increment position by Values[]  */
      Coord[ipart][X] += Values[X];
      Coord[ipart][Y] += Values[Y];
      Coord[ipart][Z] += Values[Z];

      }

   /*  Wrap displaced particles throught repeating boundary conditions  */
   WrapParticles (a);
   }


                           /*  READ NRANGE */

void read_nrange   (char *instr)
   {
   double SearchCutoffRatio;

   /*  Read ratio of neighbor search cutoff to potential cutoff  */
   SearchCutoffRatio = dblstrf (&instr);


   /*  Check value of ratio  */
   if (SearchCutoffRatio < 1)
      {
      printf ("ERROR (read_nrange):  NRANGE must be 1 or greater.\n");
      CleanAfterError();
      }

   if (SearchCutoffRatio > 2)
      {
      printf ("WARNING:  Are you sure you want NRANGE greater than 2?\n");
      printf ("  This may result in a large and expensive neighbor list.\n");
      IncrementNumberOfWarnings();
      }

   /*  Save value  */
   s->uratio  = SearchCutoffRatio;
   }


                         /*  READ PARTICLE  */

void read_particle (char *instr)
   {
   char *tokstr=NULL;
   char bufstr[NSTR];
   int  LastNewPart;
   int  LastOldPart;
   int  FirstNewPart;
   int  NumRead;
   int  ipart;
   int  Type;
   int  NumWrap;
   BOOLEAN addflag;
   int    i;
   double x,y,z;
   double (*Coord)[NDIR];

   addflag = FALSE;

   /*  Parse first token  */
   tokstr = strhed (&instr);

   /*  ADD option  */
   if (!strcmpi(tokstr,"ADD"))
      {
      addflag = TRUE;
      tokstr = strhed (&instr);
      }

   /*  No particle number entered  */
   if (tokstr==NULL)
      {
      printf ("ERROR (read_particle):  No number specified for PARTICLE.");
      CleanAfterError();
      }

   /*  Number of particles to input  */
   NumRead = intstrf (&tokstr);

   /*  Invalid number  */
   if (NumRead==0)
      {
      printf ("ERROR (read_particle): ");
      printf ("Number for PARTICLE is 0 or unrecognizable.\n");
      CleanAfterError();
      }

   /*  Save old number of particles  */
   LastOldPart = a->np;

   /*  Clear selected particles  */
   LOOP (i, a->np)
		REMOVE_SELECT(i)

   /*  READ COORDINATES  */
   if (addflag)
      {
      FirstNewPart = a->np;
      }
   else
      {
      FirstNewPart = 0;
      }

	/*  Last particle after new particles are read  */
   LastNewPart  = FirstNewPart + NumRead;

   /*  Reallocate particles if new number different from previous number  */
   if (LastNewPart != LastOldPart)
      {
      ReallocateParticle (a, LastNewPart);
      }

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

	
   /*  Save only new particles  */
   ipart = FirstNewPart;
   while (   ipart < LastNewPart
          && NULL!=m_gets_list(bufstr,NSTR,inlist))
      {
      tokstr = bufstr;
      Type = intstrf (&tokstr)-1;
      x = dblstrf (&tokstr);
      y = dblstrf (&tokstr);
      z = dblstrf (&tokstr);
		CLEAR_ALL_FLAGS(ipart)
      APPLY_SELECT(ipart)
		if (a->tag != NULL)
			a->tag [ipart] = (BYTE) 0;
		a->type[ipart] = (BYTE) Type;
      Coord[ipart][X] = x * ANGS;
      Coord[ipart][Y] = y * ANGS;
      Coord[ipart][Z] = z * ANGS;
      ipart++;
      }

   if (ipart != LastNewPart)
      {
      printf ("ERROR (read_particle):  Could not find all particles.");
      printf ("Number of particles expected is %d.\n", LastNewPart);
      printf ("Number of particles read     is %d.\n",
               ipart-LastOldPart);
      CleanAfterError();
      }

   /*  Set number of particles  */
   a->np   = LastNewPart;
   a->nsel = NumRead;

   /*  Set Coordinate flag  */
   a->IsInitializedCoord = TRUE;


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


   if (NumWrap>0)
      PrintWrapWarning (NumWrap);

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


                         /*  READ POSVEL */

void read_posvel (char *instr)
   {
   char *tokstr=NULL;
   char bufstr[NSTR];
   int  LastNewPart;
   int  LastOldPart;
   int  FirstNewPart;
   int  NumRead;
   int  ipart;
   int  Type;
   int  NumWrap;
   BOOLEAN addflag;
   int    i;
   double x,y,z;
   double vx,vy,vz;
   double (*Coord)[NDIR];
   double (*Vel  )[NDIR];

   addflag = FALSE;

   /*  Parse first token  */
   tokstr = strhed (&instr);

   /*  ADD option  */
   if (!strcmpi(tokstr,"ADD"))
      {
      addflag = TRUE;
      tokstr = strhed (&instr);
      }

   /*  No particle number entered  */
   if (tokstr==NULL)
      {
      printf ("ERROR (read_posvel):  No number specified for POSVEL.");
      CleanAfterError();
      }

   /*  Number of particles to input  */
   NumRead = intstrf (&tokstr);

   /*  Invalid number  */
   if (NumRead==0)
      {
      printf ("ERROR (read_posvel): ");
      printf ("Number for POSVEL is 0 or unrecognizable.\n");
      CleanAfterError();
      }

   /*  Save old number of particles  */
   LastOldPart = a->np;

   /*  Clear previous particle selection  */
   LOOP (i, a->np)
		REMOVE_SELECT(i)