cdsubs.c

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

         a->rtype[j] = a->type[j];
      }

   /*  Unrecognized Option  */
   else
      {
      printf ("WARNING:  Unrecognized option to REFSTEP command.\n");
      IncrementNumberOfWarnings();
      }

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

                              /*  REMOVE ATOMS  */

#pragma argsused
void read_remove (char *instr)
   {
   int ipart;
   int isave;

   /*
   Copy saved particles (the unselected particles)
   to cover up "hole" that is removed
   */
   isave = 0;
	CheckForNoneSelected();
   LOOP (ipart, a->np)
      {
      if (!IS_SELECT(ipart) )
         {
         if (isave != ipart)
            {
            CopyParticle (a, isave, ipart);
            }
         isave++;
         }
      }

   /*  Set new number of particles  */
   a->np = isave;

   /*  Set number of selected particles to zero  */
   a->nsel = 0;

   /*  Realloate particle info  */
   ReallocateParticle (a, isave);

   /*  Mark neighbor list as invalid  */
   a->InvalidNeighborList = TRUE;

   /*  Print message  */
   printf ("***   NUMBER REMAINING %i\n", isave);

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


                              /*  READ ROTATE  */
void read_rotate (char *InputString)
   {
   double Axis   [NDIR];
   double Center [NDIR];
   double Cosine;
   double Sine;
   double Rotation[NDIR][NDIR];
   double Asym    [NDIR][NDIR];
   double Magnitude;
   double Angle;
   int    idir;
   int    jdir;
	BOOLEAN UseSelect;


   /*  Read paramters  */
   Axis  [X] = dblstrf (&InputString);
   Axis  [Y] = dblstrf (&InputString);
   Axis  [Z] = dblstrf (&InputString);
   Center[X] = 1e-8 * dblstrf (&InputString);
   Center[Y] = 1e-8 * dblstrf (&InputString);
   Center[Z] = 1e-8 * dblstrf (&InputString);
   Angle     = M_PI * dblstrf (&InputString)/180.0;

   /*  No rotation is angle is zero  */
   if (Angle==0.0)
      return;

   /*  Normalize axis  */
   Magnitude = sqrt (Axis[X]*Axis[X] + Axis[Y]*Axis[Y] + Axis[Z]*Axis[Z]);
   if (Magnitude==0.0)
      {
      printf ("ERROR (read_rotate):  Null axis.\n");
      CleanAfterError();
      }
   Axis[X] /= Magnitude;
   Axis[Y] /= Magnitude;
   Axis[Z] /= Magnitude;

   /*  SET UP ASYMMETRIC MATRIX (REPRESENTS CROSS PRODUCT)  */
   Asym[X][X] = Asym[Y][Y] = Asym[Z][Z] = 0;
   Asym[X][Y] = -Axis[Z];
   Asym[Y][Z] = -Axis[X];
   Asym[Z][X] = -Axis[Y];
   Asym[Y][X] = -Asym[X][Y];
   Asym[Z][Y] = -Asym[Y][Z];
   Asym[X][Z] = -Asym[Z][X];

   /*  SET UP MATRIX  */
   Sine   = sin(Angle);
   Cosine = cos(Angle);
   for (idir=0; idir<NDIR; idir++)
   for (jdir=0; jdir<NDIR; jdir++)
      if (idir==jdir)
         {
         Rotation[idir][jdir] =
            Axis[idir]*Axis[jdir] +
            Cosine * (1 - Axis[idir]*Axis[jdir]) +
            Sine   * Asym[idir][jdir];
         }
      else
         Rotation[idir][jdir] =
            Axis[idir]*Axis[jdir] * (1 - Cosine) +
            Sine * Asym[idir][jdir];

   /*  TRANSFORM COORDINATES  */
	UseSelect = TRUE;
	CheckForNoneSelected();
   RotateCoordinates (Rotation, Center, UseSelect, a->cur, a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->v , a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->c2, a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->c3, a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->c4, a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->c5, a->np);
   RotateVectorList
      (Rotation, Center, UseSelect, a->ForcePtrList,  a->np);
   RotateVectorList
      (Rotation, Center, UseSelect, a->SpringPtrList, a->np);
   }


                           /*  READ RUN  */

void read_run (char *instr)  {
   a->run = lngstrf (&instr);
}


                            /*  READ SCALE  */

void read_scale    (char *instr)
   {
   double scale[3];
   int   idir, ipart;
   double (*Coord)[3] = (double (*)[3]) a->cur;


   /*  PARSE SCALE  */
   LOOP (idir, NDIR)
       scale[idir] = dblstrf (&instr);

   /*  Test for zero scale  */
   if (scale[X]==0.0)
      {
      printf ("ERROR:  Cannot scale by zero.\n");
      CleanAfterError();
      }

   if (scale[Y]==0.0)
       scale[Y] = scale[X];
   if (scale[Z]==0.0)
       scale[Z] = scale[Y];

   /*  Scale Box  */
   LOOP (idir, NDIR)
      {
      a->bcur[idir] *= scale[idir];
      }

   /*  Scale cavity contraint center and axis  */
   LOOP (idir, NDIR)
      {
      a->CavityCenter[idir] *= scale[idir];
      a->CavityAxis  [idir] *= scale[idir];
      }


   /*  Test for particles present  */
   if (!a->IsInitializedCoord)
      {
      printf ("WARNING:  Cannot scale particle without initializing them.\n");
      IncrementNumberOfWarnings();
      return;
      }

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

   }





void read_screw    (char *instr)
   {
   int    idir;
   int    ipart;
   double (*Coord)[NDIR] = (double (*)[NDIR]) a->cur;
   double BurgersVector[NDIR];
   double Origin    [NDIR];
   double Axis   [NDIR];
   double Reference [NDIR];
   double Perpendicular[NDIR];
   double Position  [NDIR];
   double Angle;
   double Xpos;
   double Ypos;
   double ScrewFactor;

   /*  Parse Screw Burgers Vector  */
   LOOP (idir, NDIR)
      BurgersVector[idir] = 1e-8*dblstrf (&instr);
   LOOP (idir, NDIR)
      Origin[idir] = 1e-8*dblstrf (&instr);
   LOOP (idir, NDIR)
      Reference[idir] = 1e-8*dblstrf (&instr);


   /*
   If magnitude of burgers vector is zero,
   then no need to move atoms
   */
   if (GetMagnitude(BurgersVector)==0)
      {
      return;
      }

   /*  Set screw axis as normalized Burgers Vector  */
   LOOP (idir, NDIR)
      {
      Axis[idir] = BurgersVector[idir];
      }
   Normalize(Axis);


   /*  If no angle reference vector specified, create one  */
   if (GetMagnitude(Reference)==0.0)
      {
      Reference[X] =  Axis[Y];
      Reference[Y] = -Axis[Z];
      Reference[Z] =  Axis[X];
      }

   /*  Orthogonalize reference vector to axis  */
   Orthogonalize (Reference, Axis);
   Normalize(Reference);
   if (GetMagnitude(Reference)==0.0)
      {
      printf ("INPUT ERROR: Reference vector is parallel to Burgers vector.\n");
      CleanAfterError();
      }

   /*  Find vector perpendicular to Axis[] and Reference[]  */
   CalcCross (Perpendicular, Axis, Reference);


   /*
   For each particle,
      1)  find position vector relative to Origin[]
      2)  orthgonalize position vector relative to Axis[]
      3)  obtain angle between orthoganilized position vector and Reference[]
      4)  displace atom by factor proportional to burgers vector and
          sin() of angle between Reference[] and orthogonalized position vector
   */
	CheckForNoneSelected();
   LOOP (ipart, a->np)
   if (IS_SELECT(ipart))
      {

      /*  Find atom position relative to origin  */
      Position[X] = Coord[ipart][X] - Origin[X];
      Position[Y] = Coord[ipart][Y] - Origin[Y];
      Position[Z] = Coord[ipart][Z] - Origin[Z];

      /*  Orthoginalize position vector to Axis[]  */
      Orthogonalize (Position, Axis);

      /*
      Find X and Y componenent of position given that Axis[] is Z direction
      */
      Xpos = GetDot (Position, Reference);
      Ypos = GetDot (Position, Perpendicular);
      Angle = atan2 (Ypos, Xpos);
      ScrewFactor = Angle/(2.0*M_PI);

      LOOP (idir, NDIR)
         {
         Coord[ipart][idir] += ScrewFactor * BurgersVector[idir];
         }

      }

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


                           /*  READ SEED  */
void read_seed (char *instr)
   {
   long l = lngstrf (&instr);
   l &= 0xffff;
   if (l!=0)
      {
      s->seed = (unsigned) l;
      initrand (s->seed);
      }
   }


                           /*  READ SIZE */

void read_size   (char *instr)
   {
   s->size  = ANGS*dblstrf (&instr);
   }

                           /*  READ SSAVE FLAG  */

void read_ssave  (char *instr)
   {
   ReadStepOutput ( &(s->StressStepOutput), instr);
   }

                           /*  READ STATE */

void read_state   (char *instr)
   {
   readstate (&a, s, instr);

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

void read_step   (char *instr)
   {
   char *tokstr = strhed(&instr);

   if (*tokstr=='+')
      {
      tokstr++;
      a->step += intstrf (&tokstr);
      }
   else
      a->step  = intstrf  (&tokstr);
   }


                            /*  READ STRESS */

void read_stress (char *instr) {
	char *tokstr = strhed (&instr);
	if (!strcmpi("THERMAL", tokstr)) {
		tokstr = strhed (&instr);
		if (!strcmpi("ON", tokstr)) {
			a->stress_thermal = TRUE;
		} else if (!strcmpi("OFF", tokstr)) {
			a->stress_thermal = FALSE;
		} else {
			printf ("ERORR:  Unrecognized option to STRESS THERMAL command.\n");
      	CleanAfterError();
		}
	} else {
		printf ("ERORR:  Unrecognized option to STRESS command.\n");
     	CleanAfterError();
	}
}


                           /*  READ SURFACE  */

void read_surface  (char *instr)
   {
   char *tokstr = NULL;
   BOOLEAN OnFlag;
   int  NumWrap;

   /*  GET FIRST "OFF" or "ON"  */
   tokstr = strhed (&instr);
   if      (!strcmpi(tokstr,"ON"))
      OnFlag = TRUE;
   else if (!strcmpi(tokstr,"OFF"))
      OnFlag = FALSE;
   else
      {
      printf ("        ***  ERRROR:  UNKNOWN OPTION");
      return;
      }

   /*  PARSE REMAINING OPTIONS  */
   tokstr = strhed (&instr);
   while (strcmpi(tokstr,""))
      {
      if      (!strcmpi(tokstr,"ON"))
         OnFlag = TRUE;
      else if (!strcmpi(tokstr,"OFF"))
         OnFlag = FALSE;
      else if (!strcmpi(tokstr,"X"))
         a->surf[0] = OnFlag;
      else if (!strcmpi(tokstr,"Y"))
         a->surf[1] = OnFlag;
      else if (!strcmpi(tokstr,"Z"))
         a->surf[2] = OnFlag;
      else
         printf ("        ***  ERRROR:  UNKNOWN OPTION \"%s\"\n",tokstr);
      tokstr = strhed (&instr);
      }

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

   if (NumWrap>0)
      PrintWrapWarning (NumWrap);

   /*  Mark neighbor list as invalid  */
   a->InvalidNeighborList = TRUE;

	/*  Make surface as initialized  */
	a->IsInitializedSurface = TRUE;

   }



                           /*  READ TSAVE FLAG  */

void read_tsave  (char *instr)
   {
   ReadStepOutput ( &(s->TrajectoryStepOutput), instr);
   }

                           /*  READ TYPE  */

void read_type     (char *instr)
   {
   int i;
   int curtype;

   curtype = intstrf (&instr)-1;
   LOOP (i, a->np)
      if (IS_SELECT(i))
			a->type[i] = curtype;
   }


void read_typelist (char *instr)
   {
   char *tokstr = NULL;
   char bufstr[NSTR];
   int i, nread, ntype, sflag;
   int *typelist = NULL;

   tokstr = strhed (&instr);
   if (!strcmpi(tokstr,"SEL"))
      {
      sflag = TRUE;
      ntype = intstrf (&instr);
      }
   else
      {
      sflag = FALSE;
      ntype = intstrf (&tokstr);
      }

   if (ntype==0)
      printf ("ERROR:  Number of types is 0\n");

      /*  READ TYPES  */
   else
      {
      /*  ALLOCATE ARRAY  */
      ALLOCATE (typelist, int, ntype)

      nread = 0;
      while (   nread<ntype
                && NULL!=m_gets_list(bufstr,NSTR,inlist))
         {
         tokstr = bufstr;
         while (*tokstr)
            typelist[nread++] = intstrf(&tokstr)-1;
         }
      if (nread<ntype)
         {
         printf ("Number of particles expected is %d.\n", ntype);
         printf ("Number of particles read     is %d.\n", nread);
         CleanAfterError();
         }

		/*  Check for no atoms selected  */
		if (sflag)
			CheckForNoneSelected();

      /*  ASSIGN TYPES  */
      nread = 0;
      LOOP (i, a->np)
         {
         if (!sflag || IS_SELECT(i) )
            a->type[i] = typelist[nread++];
         }
      /*  FREE ALLOCATION  */
      FREE (typelist)
      }
   }


/*  Read name string to associate with each type (used to WRITE XMOL command)  */
void read_typename (char *instr)
   {
   int   Type;
   char *NameStr;

   /*  Raad type name pairs on line  */
   while (*instr != '\0')
      {

      /*  Read type and name string   */
      Type    = intstrf (&instr)-1;
      NameStr = strhed  (&instr);

      /*  Insure that type number in range  */
      if (Type>MAX_TYPE_NAMES)
         {
         printf ("ERROR:  Maximum type for TYPENAME command is %i\n",
            MAX_TYPE_NAMES);
         CleanAfterError();
         }

      /*  Store type name string  */
      ALLOCATE (a->TypeNames[Type], char, strlen(NameStr)+1)
      strcpy   (a->TypeNames[Type], NameStr);

      }
   }


void read_verbose  (char *instr)
	{
	char *tokstr = NULL;

	tokstr = strhed (&instr);

	/*  VERBOSE FILE { ON | OFF }  */
	if (!strcmpi("FILE", tokstr))
		{
		tokstr = strhed (&instr);
		if (!strcmpi("ON",tokstr))
			{
			s->UseVerboseOutputFile = TRUE;
			}
		else if (!strcmpi("OFF",tokstr))
			{
			s->UseVerboseOutputFile = FALSE;
			}
		else
			{
   		printf ("WARNING:  Unknown option to VERBOSE FILE.\n");
  			IncrementNumberOfWarnings();
			}
		}
	/*  VERBOSE