xmd Code

Brought to you by: jrifkin
[r6]: / trunk / src / cdcor.c Maximize Restore History

748 lines (601 with data), 19.0 kB

/*
xmd - molecular dynamics for metals and ceramics
By Jonathan Rifkin <jon.rifkin@uconn.edu>
Copyright 1995-2004 Jonathan Rifkin

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/


/*
************************************************************************
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.
2009-09-26
   Suppress compiler warnings by assigning fread return value.
*/



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

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "cdsubs.h"
#include "cdboun.h"

#include "cdcor.h"

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

/*
Import
   GetCurrentTimeStr()
   CleanAfterError()
   FREE()
*/

#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_return_m = 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";

size_t fread_return_m;



/*
************************************************************************
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 WriteCorFile (FILE *OutputFile, Particle_t *a, BOOLEAN UseSelect,
BOOLEAN UseTag)
   {
   int    ititle;
   int    ipart;
   int    idir;
   int    OutputNP;
   WORD2  WordCoord[NDIR];
   double Min  [NDIR];
   double Max  [NDIR];
   double Scale[NDIR];
   char   OutputStr[STRING_LENGTH];
   char  *LineFeed = "\n\r";
   CorFileHeader_t FileHeader;

   /*  Check for no selected particles  */
   if ( (UseSelect && a->nsel==0) || a->np==0)
      return;


   /*  Wrap coordinates first  */
   WrapParticles (a);

   /*
   ****************************
   Determine range of particles
   ****************************
   */
   GetMinMaxCoord (a, Min, Max, UseSelect, TRUE);


   /*  Determine number of output particles  */
   if (UseSelect)
      OutputNP = a->nsel;
   else
      OutputNP = a->np;


   /*
   **********************
   Write Text Format Data
   **********************
   */

   /*  Write file type and version  */
   FWRITESTR (VersionStr_m, OutputFile);
   FWRITESTR (LineFeed, OutputFile);

   FWRITESTR ("DATE: ", OutputFile);
   FWRITESTR (GetCurrentTimeStr(), OutputFile);
   FWRITESTR (LineFeed, OutputFile);

   /*  Write run and step  */
   sprintf   (OutputStr, "RUN: %li", a->run);
   FWRITESTR (OutputStr, OutputFile)
   FWRITESTR (LineFeed,  OutputFile)

   sprintf   (OutputStr, "STEP: %li", a->step);
   FWRITESTR (OutputStr, OutputFile)
   FWRITESTR (LineFeed,  OutputFile)

   FWRITESTR ("TITLE:", OutputFile)
   FWRITESTR (LineFeed,  OutputFile)


   /*  Movie file comments go here, at header end  */
   LOOP (ititle, 8)
      {
      FWRITESTR (a->title[ititle], OutputFile);
      FWRITESTR (LineFeed, OutputFile);
      }

   /*  Write header end (DOS End-of-file character)  */
   FWRITESTR ("\032", OutputFile);


   /*
   *************************
   Write Machine Format Data
   *************************
   */


   /*  Set values of file header  */
   FileHeader.EndianWord   = EndianWord_m;
   FileHeader.Run          = a->run;
   FileHeader.Step         = a->step;
   FileHeader.Time         = a->time;
   FileHeader.BoundaryType = a->BoundaryType;
   FileHeader.Surf[X]      = a->surf[X];
   FileHeader.Surf[Y]      = a->surf[Y];
   FileHeader.Surf[Z]      = a->surf[Z];
   FileHeader.Box[X]       = a->bcur[X];
   FileHeader.Box[Y]       = a->bcur[Y];
   FileHeader.Box[Z]       = a->bcur[Z];
   FileHeader.Etot         = a->etot ;
   FileHeader.Epot         = a->epot ;
   FileHeader.Ekin         = a->ekin ;
   FileHeader.Ebath        = a->ebath;
   FileHeader.Min[X]       = Min[X];
   FileHeader.Min[Y]       = Min[Y];
   FileHeader.Min[Z]       = Min[Z];
   FileHeader.Max[X]       = Max[X];
   FileHeader.Max[Y]       = Max[Y];
   FileHeader.Max[Z]       = Max[Z];
   FileHeader.Np           = OutputNP;

   /*  Write output header  */
   FWRITE (FileHeader.EndianWord,    OutputFile, 1)
   FWRITE (FileHeader.Run,           OutputFile, 1)
   FWRITE (FileHeader.Step,          OutputFile, 1)
   FWRITE (FileHeader.Time,          OutputFile, 1)
   FWRITE (FileHeader.BoundaryType,  OutputFile, 1)
   FWRITE (FileHeader.Surf[X],       OutputFile, 1)
   FWRITE (FileHeader.Surf[Y],       OutputFile, 1)
   FWRITE (FileHeader.Surf[Z],       OutputFile, 1)
   FWRITE (FileHeader.Box[X],        OutputFile, 1)
   FWRITE (FileHeader.Box[Y],        OutputFile, 1)
   FWRITE (FileHeader.Box[Z],        OutputFile, 1)
   FWRITE (FileHeader.Etot,          OutputFile, 1)
   FWRITE (FileHeader.Epot,          OutputFile, 1)
   FWRITE (FileHeader.Ekin,          OutputFile, 1)
   FWRITE (FileHeader.Ebath,         OutputFile, 1)
   FWRITE (FileHeader.Min[X],        OutputFile, 1)
   FWRITE (FileHeader.Min[Y],        OutputFile, 1)
   FWRITE (FileHeader.Min[Z],        OutputFile, 1)
   FWRITE (FileHeader.Max[X],        OutputFile, 1)
   FWRITE (FileHeader.Max[Y],        OutputFile, 1)
   FWRITE (FileHeader.Max[Z],        OutputFile, 1)
   FWRITE (FileHeader.Np,            OutputFile, 1)

   /*
   ***************************
   Write Types and Coordinates
   ***************************
   */

   /*  Write particle types  */
   LOOP (ipart, a->np)
      if (!UseSelect || IS_SELECT(ipart))
			if (UseTag)
	         FWRITE (a->tag[ipart], OutputFile, 1)
			else
	         FWRITE (a->type[ipart], OutputFile, 1)

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

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

      /*  Write particle if selected  */
      if (!UseSelect || IS_SELECT(ipart))
         {

         LOOP (idir, NDIR)
            {

            /*  Scale floating point coordinate to 2-byte integer  */
            WordCoord[idir] =
               Scale[idir] * (a->cur[NDIR*ipart+idir] - Min[idir]);
            }

         /*  Write 2-byte integer  */
         FWRITE (WordCoord[0], OutputFile, NDIR)
         }


   }
/*   ... WriteCorFile  */


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