cor0.c

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

   /*  Reverse variables if necessary  */
   if (IsReverse)
      ReverseHeader (&FileHeader);

   /*  Copy header to stepdata  */
   a->run           =  FileHeader.Run;
   a->step          =  FileHeader.Step;
   a->time          =  FileHeader.Time;
   a->boundrep[X]   = !FileHeader.Surf[X];
   a->boundrep[Y]   = !FileHeader.Surf[Y];
   a->boundrep[Z]   = !FileHeader.Surf[Z];
   a->box[X]        =  FileHeader.Box[X];
   a->box[Y]        =  FileHeader.Box[Y];
   a->box[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;

   /*  Allocate space for types and particles  */
   if (a->selkeep)
      {
      REALLOCATE (a->sel, BYTE, InputNP)
      REALLOCATE (a->tag, BYTE, InputNP)
      }
   else
      {
      FREE (a->sel)
      FREE (a->tag)
      }

   /*  Allocate Space  */
   REALLOCATE (a->t,    BYTE,  InputNP)
   REALLOCATE (a->c[X], float, InputNP)
   REALLOCATE (a->c[Y], float, InputNP)
   REALLOCATE (a->c[Z], float, InputNP)
   if (a->c[Z]==NULL)
      {
      printf (
         "ERROR (readrcv):  No room to allocate %i new (%i total) particles.\n",
         InputNP-a->np, InputNP);
      return EOF;
      }

   a->np = InputNP;


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

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


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

      /*  Read scaled word coordinates  */
      fread (WordCoord, sizeof(WORD2), 3, InputFile);
      if (IsReverse)
         {
         REVERSE (WordCoord[X]);
         REVERSE (WordCoord[Y]);
         REVERSE (WordCoord[Z]);
         }

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

   /*  No errors - Return FALSE  */
   return 0;


   }
/*   ... ReadCorFile  */



/*
Open a coordinate file, use it analogously to fopen()

FileTypeStr can be either

"rt", "wt", "at"  - Read, Write or Append RCV (text)   file
"rb", "wb", "ab"  - Read, Write or Append COR (binary) file
"r",  "w" , "a"   - Read, Write or Append either file
                    (determine from file name extension)
*/
CFILE *OpenCoordinateFile (char *FileName, char *CFileTypeStr)
   {
   FILE  *FilePtr;
   CFILE *CFilePtr;
   int    CFileType;
   char   TempCFileTypeStr[3];

   /*  Check for illegal CFileTypeStr  */
   if (strlen(CFileTypeStr) < 1 || strlen(CFileTypeStr) > 2 )
      {
      printf 
			("INTERNAL ERROR:  (File %s Line %i) Illegal CFileTypeStr string (%s)",
			__FILE__,
			__LINE__,
         CFileTypeStr);
      printf (" sent to OpenCoordinateFile.\n");
      exit   (1);
      }

   /*  Copy Input file type string into buffer  */
   strncpy (TempCFileTypeStr, CFileTypeStr, 3);

   /*  Assign Coordinate file type from FileTypeStr  */
   if      (TempCFileTypeStr[1]=='t')
      CFileType=CF_RCV;
   else if (TempCFileTypeStr[1]=='b')
      CFileType=CF_COR;
   else
      CFileType=CF_NONE;

   /*  If no input file type specified, determine from file name.type  */
   if (CFileType==CF_NONE)
      {

      /*  Determine input file type  */
      CFileType = GetCoordFileType (FileName);
      if (CFileType==CF_NONE)
         {
         printf ("ERROR:  Input file type not recognized.  Must be\n");
         printf ("either COR or RCV.\n");
         exit   (1);
         }

      /*  Reset FileTypeStr  */
      switch (CFileType)
         {
         case CF_RCV:
            TempCFileTypeStr[1] = 't';
            break;
         case CF_COR:
            TempCFileTypeStr[1] = 'b';
         }

      /*  Make sure string terminator in place  */
      TempCFileTypeStr[2] = '\0';

      }

   /*  Open file  */
   FilePtr = fopen (FileName, TempCFileTypeStr);

   /*  Return null if file doesn't open  */
   if (FilePtr==NULL)
      return NULL;

   /*  Allocate new CFILE pointer  */
   CFilePtr = (CFILE *) calloc (1, sizeof(CFILE) );
   if (CFilePtr==NULL)
      return NULL;

   /*  Assign values to customized CFILE pointer  */
   CFilePtr->Ptr  = FilePtr;
   CFilePtr->Type = CFileType;

   return CFilePtr;
   }



/*  Returns TRUE if error  */
BOOLEAN ReadCoordinateFile (CFILE *CFilePtr, stepdata *CurStep)
   {
   BOOLEAN Status;

   /*  Read COR  format file  */
   if      (CFilePtr->Type==CF_COR)
      return ReadCorFile (CFilePtr->Ptr, CurStep);

   /*  Read RCV format file  */
   else if (CFilePtr->Type==CF_RCV)
      {

      Status = readrcv     (CFilePtr->Ptr, CurStep);

      /*  Convert units returned by readrcv() from ang to cm  */
      ScaleCoord (1e-8, CurStep);
      return Status;
      }

   /*  This statement is superfluous, but it avoids compilation warnings  */
   return TRUE;
   }


void WriteCoordinateFile
(
CFILE *CFilePtr,
stepdata *CurStep,
BOOLEAN SelFlag
)
   {

   if      (CFilePtr->Type==CF_COR)
      WriteCorFile (CFilePtr->Ptr, CurStep, SelFlag);

   else if (CFilePtr->Type==CF_RCV)
      writercv     (CFilePtr->Ptr, CurStep, SelFlag);

   }



void CloseCoordinateFile (CFILE *CFilePtr)
   {
   fclose (CFilePtr->Ptr);
   free   (CFilePtr);
   }



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

/*  Reverse bytes of word (works for any length WORD)  */
void ReverseWords (WORD *WordPtr, int NumWords)
   {
   BYTE *BytePtr;
   BYTE TempByte;
   int  FirstIndex;
   int  LastIndex;
   int  iword;

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

   LOOP (iword, NumWords)
      {

      /*  Find start and end of current word in list  */
      FirstIndex = iword * sizeof(WORD);
      LastIndex  = FirstIndex + sizeof(WORD) - 1;

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

      }

   }
/* .. Reverse Word  */


/*  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, 1, 1, InputFile);
      if (InputChar!='\r' && InputChar!='\n' && iread < Length-1)
         {
         InputString[iread] = InputChar;
         iread++;
         }
      }
      while (InputChar != '\r' && iread < Length-1);

   InputString[iread] = '\0';

   }

/*  Return time and date string  */
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/19%02d   %02d:%02d:%02d",
      timeptr->tm_mon+1, timeptr->tm_mday, timeptr->tm_year,
      timeptr->tm_hour, timeptr->tm_min, timeptr->tm_sec);

   /*  Return string  */
   return TimeString;
   }


/*  Return code for current file type (either CF_COR or CF_RCV)  */
int GetCoordFileType (char *FileName)
   {
   if (!strcmp(GetFileType(FileName),"COR"))
      return CF_COR;
   if (!strcmp(GetFileType(FileName),"RCV"))
      return CF_RCV;
   return CF_NONE;
   }


/*  Get file type from file name  */
char *GetFileType (char *FileName)
   {
   static char FileTypeBuffer[LEN_FILE_TYPE+1];
   int ichar;

   /*  Find period working from end of string  */
	ichar = strlen(FileName);
	while (ichar>=0 && FileName[ichar]!='.')
		ichar--;

	/*  No . suffix found  */
	if (ichar==-1)
		{
		FileTypeBuffer[0] = '\0';
		return FileTypeBuffer;
		}

   /*  Copy file type to buffer  */
   strncpy (FileTypeBuffer, FileName+ichar+1, LEN_FILE_TYPE+1);

   /*  Convert to uppercase  */
   LOOP (ichar, LEN_FILE_TYPE)
      {
      if (FileTypeBuffer[ichar]>='a' && FileTypeBuffer[ichar]<='z')
         FileTypeBuffer[ichar] += 'A' - 'a';
      }

   /*  Return pointer to buffer  */
   return FileTypeBuffer;
   }


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


void CalcMinMaxCoord
(
stepdata *a,
BOOLEAN UseSelect,
BOOLEAN UseBoundary
)
   {
   int FirstPart;
   int idir;
   int ipart;
   float TempCoord;


   /*  Find first selected particle  */
   FirstPart = 0;
   if (UseSelect)
      while (!a->sel[FirstPart] && FirstPart<a->np)
         FirstPart++;

   /*  Test for no selected partices, then use all  */
   if (FirstPart==a->np)
      {
      FirstPart = 0;
      UseSelect = FALSE;
      }

   LOOP (idir, NDIR)
      {

      /*  Check boundary conditions  */
      if (UseBoundary && a->boundrep[idir])
         {
         a->MinCoord[idir] = 0.0;
         a->MaxCoord[idir] = a->box[idir];
         }

      /*  use particles positions  */
      else
         {
         /*  Initialize min and maximum particle  */
         a->MinCoord[idir] = a->MaxCoord[idir] = a->c[idir][FirstPart];

         /*  Search for max amd min coordinate in each direction  */
         for (ipart=FirstPart+1; ipart<a->np; ipart++)
            if (!UseSelect || a->sel[ipart])
               {
               TempCoord = a->c[idir][ipart];
               if (TempCoord < a->MinCoord[idir])
                  a->MinCoord[idir] = TempCoord;
               else if (TempCoord > a->MaxCoord[idir])
                  a->MaxCoord[idir] = TempCoord;
               }
         }
      }
   }


void ScaleCoord (double Scale, stepdata *Step)
   {
   int idir;
   int ipart;

   LOOP (idir, NDIR)
   LOOP (ipart, Step->np)
      if (Step->c[idir] != NULL)
         {
         Step->c[idir][ipart] *= Scale;
         }

   LOOP (idir, NDIR)
      {
      Step->box[idir] *= Scale;
      }

   }