cdcor.c

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

void ReadCorFile (FILE *InputFile, Particle_t *a)
   {
   int    ititle;
   int    ipart;
   int    idir;
   int    InputNP;
   WORD2  WordCoord[NDIR];
   double Min  [NDIR];
   double Max  [NDIR];
   double Scale[NDIR];
   BOOLEAN IsReverse;
   char   InputChar;
   char   InputStr[STRING_LENGTH];
   CorFileHeader_t FileHeader;


   /*
   *********************
   Read Text Format Data
   *********************
   */

   /*  Determine if input file is a cor file  */
   GetNextStringFromFile (InputStr, STRING_LENGTH, InputFile);

   /*  Test for end of file  */
   if (feof(InputFile))
      {
      return;
      }

   if (strncmp(InputStr, CompareStr_m, strlen(CompareStr_m))!=0)
      {
      printf ("ERROR READING COR FILE: \n");
      printf ("  Input file is not a COR file.\n");
      CleanAfterError();
      }

   /*  Read input file version  */
   if (strcmp(InputStr, VersionStr_m)!=0)
      {
      printf ("ERROR READING COR FILE: \n");
      printf ("  Input file is version: %s\n", InputStr);
      printf ("  Can only read version: %s\n", VersionStr_m);
      CleanAfterError();
      }

   /*  Skip date, run, step, title, strings  */
   GetNextStringFromFile (InputStr, STRING_LENGTH, InputFile);
   GetNextStringFromFile (InputStr, STRING_LENGTH, InputFile);
   GetNextStringFromFile (InputStr, STRING_LENGTH, InputFile);
   GetNextStringFromFile (InputStr, STRING_LENGTH, InputFile);


   /*  Read title lines  */
   LOOP (ititle, 8)
      {
      GetNextStringFromFile (InputStr, STRING_LENGTH, InputFile);
      strcpy (a->title[ititle], InputStr);
      }


   /*  Next character should be DOS End-of-file character  */
   FREAD (InputChar, InputFile, 1)
   if (InputChar != '\032')
      {
      printf ("ERROR:  COR file is corrupted.\n");
      CleanAfterError();
      }


   /*
   ************************
   Read Machine Format Data
   ************************
   */


   FREAD (FileHeader.EndianWord,    InputFile, 1)
   FREAD (FileHeader.Run,           InputFile, 1)
   FREAD (FileHeader.Step,          InputFile, 1)
   FREAD (FileHeader.Time,          InputFile, 1)
   FREAD (FileHeader.BoundaryType,  InputFile, 1)
   FREAD (FileHeader.Surf[X],       InputFile, 1)
   FREAD (FileHeader.Surf[Y],       InputFile, 1)
   FREAD (FileHeader.Surf[Z],       InputFile, 1)
   FREAD (FileHeader.Box[X],        InputFile, 1)
   FREAD (FileHeader.Box[Y],        InputFile, 1)
   FREAD (FileHeader.Box[Z],        InputFile, 1)
   FREAD (FileHeader.Etot,          InputFile, 1)
   FREAD (FileHeader.Epot,          InputFile, 1)
   FREAD (FileHeader.Ekin,          InputFile, 1)
   FREAD (FileHeader.Ebath,         InputFile, 1)
   FREAD (FileHeader.Min[X],        InputFile, 1)
   FREAD (FileHeader.Min[Y],        InputFile, 1)
   FREAD (FileHeader.Min[Z],        InputFile, 1)
   FREAD (FileHeader.Max[X],        InputFile, 1)
   FREAD (FileHeader.Max[Y],        InputFile, 1)
   FREAD (FileHeader.Max[Z],        InputFile, 1)
   FREAD (FileHeader.Np,            InputFile, 1)

   /*
   Set flag for reversing byte order
   (necessary if file written on CPU of opposite Endian Type
   */
   IsReverse = (FileHeader.EndianWord != EndianWord_m);

   if (IsReverse)
      ReverseHeader (&FileHeader);

   /*  Copy header to Particle_t data  */
   a->run           =  FileHeader.Run;
   a->step          =  FileHeader.Step;
   a->time          =  FileHeader.Time;
   a->surf[X]       =  FileHeader.Surf[X];
   a->surf[Y]       =  FileHeader.Surf[Y];
   a->surf[Z]       =  FileHeader.Surf[Z];
   a->bcur[X]       =  FileHeader.Box[X];
   a->bcur[Y]       =  FileHeader.Box[Y];
   a->bcur[Z]       =  FileHeader.Box[Z];
   a->etot          =  FileHeader.Etot;
   a->epot          =  FileHeader.Epot;
   a->ekin          =  FileHeader.Ekin;
   a->ebath         =  FileHeader.Ebath;
   Min[X]           =  FileHeader.Min[X];
   Min[Y]           =  FileHeader.Min[Y];
   Min[Z]           =  FileHeader.Min[Z];
   Max[X]           =  FileHeader.Max[X];
   Max[Y]           =  FileHeader.Max[Y];
   Max[Z]           =  FileHeader.Max[Z];
   InputNP          =  (int) FileHeader.Np;

   /*  Initialize variables  */
	a->IsInitializedCoord = TRUE;
	a->IsInitializedBox   = TRUE;


   /*
   ************************
   Allocate space as needed
   ************************
   */

   /*  Allocate space for types and particles  */
   /*  Reallocate particles if new number different from previous number  */
   if (InputNP != a->np)
      ReallocateParticle (a, InputNP);
   a->np = InputNP;


   /*
   ************************
   Remove tag, set, select
   ************************
   */

	if (!a->selkeep)
      {

      /*  Remove tag  */
      if (a->tag!=NULL)
         {
         FREE (a->tag);
         a->tag = NULL;
         }

      /*  Clear select, sets  */
		LOOP (ipart, InputNP)
			CLEAR_ALL_FLAGS(ipart);


		a->nsel = 0;

      }




   /*
   **********************************
   Read particle type and coordinates
   **********************************
   */

   /*
   Read particle types (these are individual bytes ->
   don't need  to reverse byte order due to endian mis-match)
   */
   FREAD (a->type[0], InputFile, a->np)

   /*  Find scale for particles  */
   LOOP (idir, NDIR)
      {
      if (Max[idir]==Min[idir])
         Scale[idir] = 0.0;
      else
         Scale[idir] = (Max[idir] - Min[idir]) / MAX_WORD2;
      }


   /*  Read particle coordinates  */
   LOOP (ipart, a->np)
      {

      /*  Read scaled word coordinates  */
      FREAD   (WordCoord[X], InputFile, NDIR)

      /*  Reverse byte order if necessary  */
      if (IsReverse)
         {
         REVERSE (WordCoord[X])
         REVERSE (WordCoord[Y])
         REVERSE (WordCoord[Z])
         }

      /*  Scale to real coordinates  */
      LOOP (idir, NDIR)
         a->cur[NDIR*ipart+idir] = 
				Scale[idir]*(WordCoord[idir]+0.5) + Min[idir];
      }
   }
/*   ... ReadCorFile  */



/*
************************************************************************
Local Functions
************************************************************************
*/

/*  Reverse bytes of word (works for any length WORD)  */
void ReverseVars (void *VarPtr, int SizeVar, int NumVar)
   {
   BYTE *BytePtr;
   BYTE TempByte;
   int  FirstIndex;
   int  LastIndex;
   int  ivar;

   /*  Set BYTE pointer to WORD array  */
   BytePtr = (BYTE *) VarPtr;

   LOOP (ivar, NumVar)
      {

      /*  Find start and end of current word in list  */
      FirstIndex = ivar * SizeVar;
      LastIndex  = FirstIndex + SizeVar - 1;

      while (FirstIndex<LastIndex)
         {
         TempByte = BytePtr[FirstIndex];
         BytePtr[FirstIndex] = BytePtr[LastIndex];
         BytePtr[LastIndex ] = TempByte;
         FirstIndex++;
         LastIndex--;
         }

      }

   }
/* .. Reverse Vars  */


/*  Read binary file into string until LineFeed encountered  */
void GetNextStringFromFile (char *InputString, int Length, FILE *InputFile)
   {
   int iread;
   char InputChar;

   iread = 0;
   do
      {
      FREAD (InputChar, InputFile, 1);
      if (InputChar!='\r' && InputChar!='\n' && iread < Length-1)
         {
         InputString[iread] = InputChar;
/*test*/
#if 0
printf ("iread InputChar %4i %c <%3i>\n", iread, InputChar, InputChar);
#endif
/*end test*/
         iread++;
         }
      }
#if 0
      while (InputChar != '\n' && !feof(InputFile));
#endif
      while (InputChar != '\r' && !feof(InputFile));

   InputString[iread] = '\0';

   }

/*  Return time and date string  */
#if 0
static char *GetCurrentTimeStr (void)
   {
   struct tm *timeptr;
   static char TimeString[23];
   time_t CurrentTime;

   /*  Get current time  (in seconds from 1900)  */
   time (&CurrentTime);

   /*  Convert to date, time  */
   timeptr = localtime (&CurrentTime);

   /*  Form string  */
   sprintf (TimeString, "%02d/%02d/%04d   %02d:%02d:%02d",
      timeptr->tm_mon+1, timeptr->tm_mday, 1900+timeptr->tm_year,
      timeptr->tm_hour, timeptr->tm_min, timeptr->tm_sec);

   /*  Return string  */
   return TimeString;
   }
#endif


/*  Reverse bytes of variable  */
void ReverseVar (void *VarPtr, int SizeVar)
   {
   BYTE *BytePtr;
   BYTE TempByte;
   int  FirstIndex;
   int  LastIndex;

   /*  Return if no bytes to switch  */
   if (SizeVar<2)
      return;

   /*  Set BYTE pointer to WORD array  */
   BytePtr = (BYTE *) VarPtr;

   /*  Find start and end of current word in list  */
   FirstIndex = 0;
   LastIndex  = FirstIndex + SizeVar - 1;

   while (FirstIndex<LastIndex)
      {
      TempByte = BytePtr[FirstIndex];
      BytePtr[FirstIndex] = BytePtr[LastIndex];
      BytePtr[LastIndex ] = TempByte;
      FirstIndex++;
      LastIndex--;
      }


   }
/* .. Reverse Var  */

void ReverseHeader (CorFileHeader_t *FileHeader)
   {
   REVERSE (FileHeader->Run)
   REVERSE (FileHeader->Step)
   REVERSE (FileHeader->Time)
   REVERSE (FileHeader->BoundaryType)
   REVERSE (FileHeader->Surf[X])
   REVERSE (FileHeader->Surf[Y])
   REVERSE (FileHeader->Surf[Z])
   REVERSE (FileHeader->Box[X])
   REVERSE (FileHeader->Box[Y])
   REVERSE (FileHeader->Box[Z])
   REVERSE (FileHeader->Etot)
   REVERSE (FileHeader->Epot)
   REVERSE (FileHeader->Ekin)
   REVERSE (FileHeader->Ebath)
   REVERSE (FileHeader->Min[X])
   REVERSE (FileHeader->Min[Y])
   REVERSE (FileHeader->Min[Z])
   REVERSE (FileHeader->Max[X])
   REVERSE (FileHeader->Max[Y])
   REVERSE (FileHeader->Max[Z])
   REVERSE (FileHeader->Np)
   }