cdcor.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 "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 (&(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 WriteCorFile (FILE *OutputFile, Particle_t *a, BOOLEAN UseSelect,
BOOLEAN UseTag)
   {
   int    ititle;
   int    ipart;
   int    idir;
   int    FirstPart;
   int    OutputNP;
   double TempCoord;
   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
   ****************************
   */

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

   /*  Test for internal error  */
   if (FirstPart==a->np)
      {
      printf ("INTERNAL ERROR:  Number of selected particles is in error.\n");
      CleanAfterError();
      }

   /*  Initialize min and maximum particle  */
   Min[X] = Max[X] = a->cur[NDIR*FirstPart+X];
   Min[Y] = Max[Y] = a->cur[NDIR*FirstPart+Y];
   Min[Z] = Max[Z] = a->cur[NDIR*FirstPart+Z];

   /*  Search for max amd min coordinate in each direction  */
   for (ipart=FirstPart+1; ipart<a->np; ipart++)
      if (!UseSelect || IS_SELECT(ipart))
         {
         LOOP (idir, NDIR)
            {
            TempCoord = a->cur[NDIR*ipart+idir];
            if (TempCoord < Min[idir])
               Min[idir] = TempCoord;
            else if (TempCoord > Max[idir])
               Max[idir] = TempCoord;
            }
         }

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