cdsubs.c

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

   /*  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);
      }

	/*  Allocate velocities if not already present  */
	if (a->v==NULL)
		ALLOCATE (a->v, double, NDIR*a->NumPartAlloc)

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

	
   /*  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);
      vx = dblstrf (&tokstr);
      vy = dblstrf (&tokstr);
      vz = dblstrf (&tokstr);
		CLEAR_ALL_FLAGS(ipart)
      APPLY_SELECT(ipart)
		a->type[ipart] = Type;
      Coord[ipart][X] = x  * ANGS;
      Coord[ipart][Y] = y  * ANGS;
      Coord[ipart][Z] = z  * ANGS;
      Vel  [ipart][X] = vx * s->dtime;
      Vel  [ipart][Y] = vy * s->dtime;
      Vel  [ipart][Z] = vz * s->dtime;
      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);

	/*  Warn user if input particles needed wrapping  */
   if (NumWrap>0)
      PrintWrapWarning (NumWrap);

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



                           /*  READ POTENTIAL  */

/*  potential { eam | quartic | tersoff }          */

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

   /*  SET POTENTIAL TYPES   */
	if (IsFirstToken("SET",instr))
      {

      /*  Record that potential was initialized  */
      a->IsInitializedPotential = TRUE;

		/*  Remove leading "SET" token  */
		strhed (&instr);

      /*  Get Next Token  */
      tokstr = strhed (&instr);

      /*  QUARTIC option  */
      if (!strcmpi(tokstr,"QUARTIC"))
         {
         strcpy (s->pottype, tokstr);
         s->energylistpnt = (void *) qu_energy_list;
         s->forcepnt      = (void *) qu_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) POT_Free;
         return;
         }

      /*  TERSOFF option   */
      if (!strcmpi(tokstr,"TERSOFF"))
         {
         strcpy (s->pottype, tokstr);
         s->energylistpnt = (void *) csi_energy_list;
         s->forcepnt      = (void *) csi_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) POT_Free;
         csi_read_potential ("INIT", "");
         return;
         }

      /*  Stillinger Weber  option  */
      if (!strcmpi(tokstr,"STILL"))
         {
         strcpy (s->pottype, tokstr);
         s->energylistpnt = (void *) stil_energy_list;
         s->forcepnt      = (void *) stil_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) POT_Free;
			STILL_Init (instr);
         return;
         }

      /*  EAM option  */
      if (!strcmpi(tokstr,"EAM"))
         {
         strcpy (s->pottype, tokstr);
         s->energylistpnt = (void *) em_energy_list;
         s->forcepnt      = (void *) em_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) POT_Free;
         em_read_potential ("INIT", instr);
         return;
         }

         /*  PAIR option  */
      if (!strcmpi(tokstr,"PAIR"))
         {
         strcpy (s->pottype, tokstr);
         s->energylistpnt = (void *) pr_energy_list;
         s->forcepnt      = (void *) pr_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) PAIR_Free;
			PAIR_Init(instr);
         return;
         }

      /*  Unrecognized option   */
      printf ("ERROR:  Unrecognized Potential Type (%s)\n", tokstr);
      CleanAfterError();

      }


      /*  Pass commands to appropriate potential routine  */

	/*  This routine wants a single string  */
	if (!strcmpi(s->pottype,"PAIR"))
		{
		PAIR_Parse(instr);
		}
	
	else if (!strcmpi(s->pottype,"STILL"))
		{
		STILL_Parse(instr);
		}

	/*  These routines want front token split off  */
	else
		{

		tokstr = strhed (&instr);

      if      (!strcmpi(s->pottype,"QUARTIC"))
         qu_read_potential (tokstr, instr);

      else if (!strcmpi(s->pottype,"EAM"))
         em_read_potential (tokstr, instr);

      else if (!strcmpi(s->pottype,"TERSOFF"))
         csi_read_potential (tokstr, instr);
		}

   }



/*  read PRESSURE  */

void read_pressure   (char *instr)
   {
   char *LeadingToken;

   while (!IsBlank (instr))
      {

      LeadingToken = strhed (&instr);

      if (!strcmpi("OFF",LeadingToken))
         a->BoxMotionAlgorithm = BMA_NONE;

      else if (!strcmpi("ANDERSEN",LeadingToken))
			{
			/*  Set pressure algorithm to Andersen  */
         a->BoxMotionAlgorithm = BMA_ANDERSEN;

         /*  Read in box masses  */
         a->BoxMass[X] = dblstrf (&instr)/AVAG;
         a->BoxMass[Y] = dblstrf (&instr)/AVAG;
         a->BoxMass[Z] = dblstrf (&instr)/AVAG;

         /*  Quit if X mass is invalid  */
         if (a->BoxMass[X]<=0)
            {
            printf ("ERROR:  Box mass is less or equal to zero.\n");
            CleanAfterError();
            }

         /*  Set Y and Z to X mass if they are invalid  */
         if (a->BoxMass[Y]<=0)
            a->BoxMass[Y] = a->BoxMass[X];
         if (a->BoxMass[Z]<=0)
            a->BoxMass[Z] = a->BoxMass[Y];
         }

		else if (!strcmpi("CLAMP",LeadingToken))
			{

			/*  Set pressure algorithm to Andersen  */
         a->BoxMotionAlgorithm = BMA_CLAMP;

			/*  Convert from MBAR to erg/cm^3  */
			a->BulkModulus = dblstrf (&instr)/MBAR;

			/*  Test Bulkmodulus  */
			if (a->BulkModulus <= 0.0)
				{
				printf
					(
					"ERROR:  Bulk Modulus (%le) equal to or less than zero\n",
					a->BulkModulus*MBAR
					);
				CleanAfterError();
				}

			/*  Read in cstep if present  */
			if (!IsBlank(instr))
				a->VolClampStep = dblstrf(&instr);
			}

		/*  
		Read X,Y,Z external pressures
		if only one pressure,  then X,Y,Z pressure equal
		if two      pressures, then Y=Z
		*/
      else if (!strcmpi("EXTERNAL",LeadingToken))
			{
         a->Pressure[X] = dblstrf (&instr)/MBAR;
			if (IsBlank(instr))
				{
				a->Pressure[Y] = a->Pressure[Z] = a->Pressure[X];
				}
			else
				{
				a->Pressure[Y] = dblstrf (&instr)/MBAR;
				if (IsBlank(instr))
					{
					a->Pressure[Z] = a->Pressure[Y];
					}
				else
					{
					a->Pressure[Z] = dblstrf (&instr)/MBAR;
					}
				}
			}

		/*  Box dimensions expands/contracts uniformly  */
      else if (!strcmpi("ISOTROPIC",LeadingToken))
         a->BoxMotionGeometry = BMG_ISOTROPIC;

		/*  Box dimensions expand/contract independently  */
      else if 
			(
			!strcmpi("ORTHONORMAL",LeadingToken) ||
			!strcmpi("ORTHORHOMBIC",LeadingToken)
			)
         a->BoxMotionGeometry = BMG_ORTHORHOMBIC;

      else
         {
         printf ("ERROR:  Unrecognized option to PRESSURE command.\n");
         CleanAfterError();
         }

      }

   }



                           /*  READ PSTRAIN  */

void read_pstrain (char *instr)
   {
   int  i;
   double dot, eps, dx,dy,dz, nx,ny,nz;
   double (*Coord)[NDIR] = (double (*)[NDIR]) a->cur;

   eps  = dblstrf (&instr);

   dx   = dblstrf (&instr);
   dy   = dblstrf (&instr);
   dz   = dblstrf (&instr);

   nx   = dblstrf (&instr);
   ny   = dblstrf (&instr);
   nz   = dblstrf (&instr);

	CheckForNoneSelected();
   LOOP (i, a->np)
   if (IS_SELECT(i) )
      {
      dot = eps * (nx*Coord[i][0] + ny*Coord[i][1] + nz*Coord[i][2]);
      Coord[i][0] += dot*dx;
      Coord[i][1] += dot*dy;
      Coord[i][2] += dot*dz;
      }
   WrapParticles (a);
   }

                           /*  READ QUENCH  */
void read_quench (char *instr)
	{
   int nstep = intstrf (&instr);


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

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

	if (!a->IsInitializedMass)
		{
		printf ("ERROR:  Cannot peform QUENCH 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();
		}

   /*  TURN QUENCH FLAG ON  */
   if (*instr==0)
      s->quench=1;   

	/*  Read non-standard quench value  */
   else
      s->quench= intstrf (&instr);

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

	/*  TURN QUENCH FLAG OFF */
   s->quench=0;   

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

	}

                           /*  READ RCV FILE  */

/*  ADD FILE TO STACK in[]  */
void read_rcv (char *instr)  {
   long step;
   FILE *fin;
   char *tokstr = strhed (&instr);

   /*  OPEN FILE  */
   fin = fopen (tokstr, "rt");
   if (fin==NULL) {
      printf ("ERROR (read_rcv):  Cannot open input file %s.\n", tokstr);
      return;
   }
   /*  READ JUST FIRST STEP  */
   if (*instr==0)
      readrcv (fin, a);
   else {
      step = intstrf(&instr);
      do {
         readrcv (fin, a);
      } while (!feof(fin) && a->step!=step);
   }
   fclose (fin);

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


                           /*  READ FILE  */

/*  ADD FILE TO STACK in[]  */
void read_read (char *InputStr)
   {
   FILE *InputFile;
   char *LeadingTokenStr;
   char  BufferStr[NSTR];
   BOOLEAN QuitFlag;

   /*  Get leading token  */
   LeadingTokenStr = strhed (&InputStr);


   /*  Open file  */
   InputFile = fopen (LeadingTokenStr, "rt");
   if (InputFile==NULL)
      {
      printf ("ERROR:  Cannot open file %s.\n", LeadingTokenStr);
      CleanAfterError();
      }

   /*  Add file to input list  */
   m_add_list (&inlist, InputFile, "f");


   /*  Read file until end  */
   QuitFlag = FALSE;
   while (!QuitFlag && m_gets_list(BufferStr,NSTR,inlist)!=NULL)
      {
      read_command (BufferStr, &QuitFlag);
      }

   /*  Remove file from input list (also closes file)   */
   m_del_list(&inlist);

   }

                              /*  READ REFSTEP  */
void read_refstep (char *instr)
   {
   double *cp;
   double (*c)[3], (*r)[3];
   double b[3];
   int   i,j;
   unsigned char *ti;
   char *tokstr = NULL;

   tokstr = strhed (&instr);

   /*  TEST DATA INTEGRITY  */
   if (a->ref == NULL)
      {
      if      (!strcmpi(tokstr,"SWAP"))
         {
         printf ("ERROR (read_refstep):  Cannot swap particles because\n");
         printf ("   there are no reference particles.\n");
         CleanAfterError();
         }
      else if (!strcmpi(tokstr,"COPY"))
         {
         printf ("ERROR (read_refstep):  Cannot copy particles because\n");
         printf ("   there are no reference particles.\n");
         CleanAfterError();
         }
      }


   /*  Swap current and reference step  */
   if      (!strcmpi(tokstr,"SWAP"))
      {

      /*  Save current particle info  */
      cp   = a->cur;
      b[0] = a->bcur[0];
      b[1] = a->bcur[1];
      b[2] = a->bcur[2];
      ti   = a->type;

      /*  Copy reference info into current step  */
      a->cur     = a->ref;
      a->bcur[0] = a->bref[0];
      a->bcur[1] = a->bref[1];
      a->bcur[2] = a->bref[2];
      a->type    = a->rtype;

      /*  Copy saved current info into reference step  */
      a->ref     = cp;
      a->bref[0] = b[0];
      a->bref[1] = b[1];
      a->bref[2] = b[2];
      a->rtype   = ti;
      }

   /*  COPY REFERENCE PARTICLES TO CURRENT STEP  */
   else if (!strcmpi(tokstr,"COPY"))
      {
      /*  Copy box  */
      for (i=0;i<3;i++)
         a->bcur[i] = a->bref[i];

      /*  Copy types  */
      for (j=0;j<a->np;j++)
         a->type[j] = a->rtype[j];

      /*  Copy coordinates  */
      c = (double (*)[3]) a->cur;
      r = (double (*)[3]) a->ref;
      for (i=0;i<a->np;i++)
         for (j=0;j<3;j++)
            c[i][j] = r[i][j];
      }

   /*  CLEAR REFERENCE STEP  */
   else if (!strcmpi(tokstr,"CLEAR"))
      {
      FREE (a->ref)
      FREE (a->rtype)
      }

   /*  COPY CURRENT PARTICLES TO REFERENCE STEP  */
   else if (*instr==0)
      {

      /*  Allocate reference particle array  */
      if (a->ref == NULL)
			{
			ALLOCATE (a->ref, double, NDIR*a->NumPartAlloc)
			}

      /*  Allocate reference type array  */
      if (a->rtype == NULL)
         ALLOCATE (a->rtype, BYTE, a->NumPartAlloc)

      /*  Copy Particles  */
      c = (double (*)[3]) a->cur;
      r = (double (*)[3]) a->ref;
      for (i=0;i<3;i++)
         {
         a->bref[i] = a->bcur[i];
         for (j=0;j<a->np;j++)
            r[j][i] = c[j][i];
         }

      /*  Copy reference types  */
      for (j=0;j<a->np;j++)