cdcorlt.c

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

/*
************************************************************************
History
************************************************************************
*/
/*
8 Jan 1997
   Added test in ReadCorFile to return if end-of-file encountered
when reading first string.  This was a problem encountered
when trying to read a COR file with more than one step.
*/



/*
************************************************************************
Compile Switches
************************************************************************
*/

/*
************************************************************************
File Includes
************************************************************************
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "cdcor.h"

/*
Import
   Structure Particle_t
*/
#include "particle.h"

/*
Import
   CleanAfterError()
*/

#include "cdhouse.h"

/*
Import
  ReallocateParticle()
*/
#include "cdalloc.h"



/*
************************************************************************
Defines
************************************************************************
*/

/*  Boolean values  */
#define BOOLEAN int
#define TRUE    1
#define FALSE   0

/*
Define types for standardizing file reads between machines
  - Following definitions valid for PC, Linux and AIX
*/
#define INT2    short
#define INT4    int
#define WORD1   unsigned char
#define WORD2   unsigned short
#define WORD4   unsigned long
#define REAL4   float
#define REAL8   double


/*  Maximum value of word  */
#define MAX_WORD2 0xffff

/*  Default string length  */
#define STRING_LENGTH 256


/*
************************************************************************
Macros
************************************************************************
*/


/*  Macro for reversing byte order  */
#define REVERSE(VAR) \
   ReverseVar (&(VAR), sizeof(VAR));

/*  Macro for block write  */
#define FWRITE(VAR,FILE,N) \
   fwrite (&(VAR), sizeof(VAR), (N), (FILE));

/*  Macro for block read  */
#define FREAD(VAR,FILE,N) \
   fread (&(VAR), sizeof(VAR), (N), (FILE));

/*  Macro for string write  */
#define FWRITESTR(STRING,FILE) \
   fwrite (STRING, strlen(STRING), (1), (FILE));



/*
************************************************************************
Type Definitions
************************************************************************
*/
typedef struct
   {
   INT2   EndianWord;
   INT4   Run;
   INT4   Step;
   REAL8  Time;
   WORD2  BoundaryType;
   WORD1  Surf[NDIR];
   REAL8  Box [NDIR];
   REAL8  Etot;
   REAL8  Epot;
   REAL8  Ekin;
   REAL8  Ebath;
   REAL8  Min [NDIR];
   REAL8  Max [NDIR];
   INT4   Np;
   }
   CorFileHeader_t;


/*
************************************************************************
Module-wide Variables
************************************************************************
*/
INT2   EndianWord_m         = 0x0001;
INT2   SwitchedEndianWord_m = 0x0100;

char *VersionStr_m = "COR FILE: Version 1.0.0";
char *CompareStr_m = "COR FILE";



/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void GetNextStringFromFile (char *InputString, int Length, FILE *InputFile);
void ReverseVar (void *WordPtr, int SizeVar);
void ReverseHeader (CorFileHeader_t *FileHeader);


/*
************************************************************************
Exported Functions
************************************************************************
*/
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)
WRITEVAR(FileHeader.Run,%i)
WRITEVAR(FileHeader.Step,%i)
WRITEVAR(FileHeader.Time,%le)
WRITEVAR(FileHeader.Box[X],%le)
WRITEVAR(FileHeader.Box[Y],%le)
WRITEVAR(FileHeader.Box[Z],%le)

   /*
   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->coordflag     = TRUE;
   a->boxflag       = 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);

   /*  Reallocate tag info if keeping it  */
   if (a->selkeep)
      {
      if (a->tag==NULL)
         a->tag  = (BYTE *) calloc   ( (BYTE) InputNP, sizeof(BYTE) );
      }

   /* .. else remove tag info, clear select  */
   else
      {

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

      /*  Clear select  */
		if (a->Selection != NULL)
         LOOP (ipart, InputNP)
				REMOVE_SELECT(ipart);

      }

   a->np = InputNP;



   /*
   **********************************
   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)
      {
WRITEVAR(Min[idir],%le)
WRITEVAR(Max[idir],%le)
      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] + 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;
         iread++;
         }
      }
      while (InputChar != '\r' && !feof(InputFile));

   InputString[iread] = '\0';

   }



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