cor0.c

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

/*

SOURCE CODE FOR READING AND WRITING RCV FILES
Uses stepdata structure, modelled after Pascal version of PLOTL

READRCV : read a single step from a rcv file.
WRITERCV: write a single step to a rcv file.

 6 Jun 1991
28 Oct 1991
 6 Nov 1991 - add code to leave types alone if selkeep flag on
15 Nov 1991 - set trans to 0 if boundrep is on
17 Nov 1991 - streamline file read
21 Nov 1991 - error with trans, was setting to zero when boundrep data
              was present rather than vice versa.
15 May 1992 - error occured in WRITERCV with line maxtbl += 1.5;.  Computer
              would hang when compiled with CLIPRCV.  Changed to maxtbl += 1;.

?? Dec 1995 - Added COR routines
23 Feb 1996 - Write elements of FileHeader struct individually
              instead of all together.  This eliminates porting
              problems on machines on RS6000 which pad the
              structure with bytes between variables, becuase
              writing the entire structure with one fwrite()
              command writes the padding bytes as well as the
              variables, which will confused programs run on
              other machines which try to read the resulting file.

*/

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

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

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



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

/*  Define default float type  */
#define FLOAT float

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

/*  Default string length  */
#define STRING_LENGTH 256

/*  Types of boundaries  */
#define BD_ORTHONORMAL  1     /*  Standard box boundary         */
#define BD_MONOCLINIC   2     /*  Parallel-piped boundary       */
#define BC_ARC          3     /*  Curved "race track" boundary  */

/*  Number of characters in file type  */
#define LEN_FILE_TYPE 3

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

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

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

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

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

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

/*  Macro to reallocate array  */
#define REALLOCATE(ARRAY,TYPE,NUM) \
   if ((ARRAY)==NULL) \
      ARRAY = (TYPE *) calloc ((NUM), sizeof(TYPE)); \
   else \
      ARRAY = realloc ((ARRAY), (NUM)*sizeof(TYPE) );

/*  Macro to free array  */
#define FREE(ARRAY) \
   if ((ARRAY)!=NULL) \
      { \
      free (ARRAY); \
      (ARRAY) = NULL; \
      }

/*  Macro for loops  */
#define LOOP(INDEX,LIMIT) \
   for (INDEX=0; INDEX<(LIMIT); INDEX++)



#define WRITEMSG \
	printf ("In File %s at Line %i: ", __FILE__, __LINE__); 

#define WRITEVAR(VAL,FORMAT) \
	printf ("File %s Line %i: ", __FILE__, __LINE__); \
   printf ("%s = ", #VAL); \
   printf (#FORMAT, VAL); \
   printf ("\n");

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

/* TRANSLATE: CONVERTS TWO Code byte VALUES INTO integer FROM 0 TO 9999 */
/* THIS TABLE TRANSLATES BETWEEN charACTERS SENT FROM PURDUE
 CDC 205 AND THE CRAY AT LLL TO THE  IBM PC.  THE SEND PROGRAM
 TRANSLATES THE FOLLOWING charACTERS
 CHARACTER      REASON
 ---------      -------
 27 -> 17       CRAY
 28 -> 14
 30 -> 15
 31 -> 25       CYBER LINE FEED CONTROL CHARACTER
 92 -> 16
(94 -> 20)      ASCII TO EBCDIC CONVERSION  )
(124-> 21)      ASCII TO EBCDIC CONVERSION  )
 126-> 20       VAX MAIL UTILITY
*/

int off  =   27;
int base =  100; /* NOTE:  Base is fixed in values of table2  */
int off2 = 2727; /* off2 = off * (base + 1)  */
int table[128] =
   {  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,
          11,  12,  13,  28,  30,  92,  27,  18,  19, 126,
          21,  22,  23,  24,  31,  26,  27,  28,  29,  30,
          31,  32,  33,  34,  35,  36,  37,  38,  39,  40,
          41,  42,  43,  44,  45,  46,  47,  48,  49,  50,
          51,  52,  53,  54,  55,  56,  57,  58,  59,  60,
          61,  62,  63,  64,  65,  66,  67,  68,  69,  70,
          71,  72,  73,  74,  75,  76,  77,  78,  79,  80,
          81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
          91,  92,  93,  94,  95,  96,  97,  98,  99, 100,
         101, 102, 103, 104, 105, 106, 107, 108, 109, 110,
         111, 112, 113, 114, 115, 116, 117, 118, 119, 120,
         121, 122, 123, 124, 125, 126, 127                };

int table2[128] =
 {  000,   100,   200,   300,   400,   500,   600,   700,   800,   900,  1000,
          1100,  1200,  1300,  2800,  3000,  9200,  2700,  1800,  1900, 12600,
          2100,  2200,  2300,  2400,  3100,  2600,  2700,  2800,  2900,  3000,
          3100,  3200,  3300,  3400,  3500,  3600,  3700,  3800,  3900,  4000,
          4100,  4200,  4300,  4400,  4500,  4600,  4700,  4800,  4900,  5000,
          5100,  5200,  5300,  5400,  5500,  5600,  5700,  5800,  5900,  6000,
          6100,  6200,  6300,  6400,  6500,  6600,  6700,  6800,  6900,  7000,
          7100,  7200,  7300,  7400,  7500,  7600,  7700,  7800,  7900,  8000,
          8100,  8200,  8300,  8400,  8500,  8600,  8700,  8800,  8900,  9000,
          9100,  9200,  9300,  9400,  9500,  9600,  9700,  9800,  9900, 10000,
         10100, 10200, 10300, 10400, 10500, 10600, 10700, 10800, 10900, 11000,
         11100, 11200, 11300, 11400, 11500, 11600, 11700, 11800, 11900, 12000,
         12100, 12200, 12300, 12400, 12500, 12600, 12700              };

INT2  EndianWord_m = 0x0001;
INT2  SwitchedEndianWord_m = 0x0100;
INT4  BoundaryType_m = BD_ORTHONORMAL;

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);
static char *GetCurrentTimeStr (void);
int GetCoordFileType (char *FileName);
char *GetFileType (char *FileName);
void ReverseVar (void *WordPtr, int SizeVar);
void ReverseHeader (CorFileHeader_t *);

void ScaleCoord (double, stepdata *);


/*
************************************************************************
Exported Routines
************************************************************************
*/


int readrcv (FILE *fin, stepdata *cur)
   {
   int res  = 10000;
   float     scale[3], tablefactor;
   unsigned  np,i, n,iline, ntable;
   long      maxtable;
   int       teof,xlo,xhi,ylo,yhi,zlo,zhi;
   int done;
   char instr[NSTR], *tstr;

   /* READ FILE HEADER                       */
   xlo = getc (fin); /*  GOBBLE UP FIRST CHARACTER  */
   if (xlo==EOF)
      return EOF;
   fgets (instr, NSTR, fin);
   teof = sscanf (instr, "%li %li %f %f %f %i %i\n",
      &cur->run, &cur->step, &cur->time, &cur->tempc, &cur->temp,
      &np, &cur->ndim);
   if (teof==EOF)
      return EOF;
   /*  AJUST NUMBER OF POINTS  */
   np--;
   for (i=0;i<NTITLE;i++)
      {
      fgets (cur->title[i], NSTR, fin);
      strtok(cur->title[i], "\n");
      }

   /*  READ MISC DATA LINE  */
   fgets  (instr, NSTR, fin);
   tstr = instr;
   cur->mass   = dblstr (&tstr);
   cur->vibp   = dblstr (&tstr);
   cur->dtime  = dblstr (&tstr);
   cur->cutoff = dblstr (&tstr);

   /*  READ Strain AND Stress  */
   if (*tstr!=0)
      {
      cur->strain = dblstr (&tstr);
      cur->stress = dblstr (&tstr);
      }
   else
      cur->stress = cur->strain = 0.0;

   /*  READ BOX  */
   fgets  (instr, NSTR, fin);
   tstr = instr;
   for (i=0;i<3;i++)
      cur->box[i] = dblstr(&tstr);

   /*  READ Trans IF PRESENT  */
   if (*tstr!=0)
      for (i=0;i<3;i++)
         cur->trans[i] = dblstr(&tstr);
   else
      for (i=0;i<3;i++)
          cur->trans[i] = 0;

   /*  BOUNDARY IF PRESENT  */
   if (*tstr!=0)
      {
      for (i=0;i<3;i++)
         cur->boundrep[i] = intstr(&tstr);
      }
   else
      for (i=0;i<3;i++)
         {
         cur->boundrep[i] = 0;
         cur->trans[i] = 0.0;
         }

   for (i=0;i<3;i++)
      scale[i] = cur->box[i] / res;

   /*  Allocate Select Info  */
   if (cur->selkeep)
      {
      REALLOCATE (cur->sel, BYTE, np)
      REALLOCATE (cur->tag, BYTE, np)
      }
   else
      {
      FREE (cur->sel)
      FREE (cur->tag)
      }

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

   /*  Read COORDINATES IN CodeD FORM   */
   n = 0;
   iline = 0;
   do
      /* UNTIL END-OF-FILE  OR  UNTIL ASCII 0 Code FOUND */
      {
      i = 0;
      fgets (instr, NSTR, fin);
      tstr = instr;
      do
         /* UNTIL 13 POSITIONS Read IN   OR  END-OF-FILE  OR  0 FOUND*/
         {
         xlo = *(tstr++);
         if (xlo!='\n' && xlo!=0)
            {
            xhi = *(tstr++);
            ylo = *(tstr++);
            yhi = *(tstr++);
            zlo = *(tstr++);
            zhi = *(tstr++);
            if (feof(fin))
               {
               fprintf (stderr, "ERROR (readrcv):  Premature end of file in step %i.\n", cur->step);
               exit(1);
               }
            /*  TEST FOR N TOO BIG  */
            if (n==np)
               {
               fprintf (stderr, "ERROR (readrcv): Data is corrupted.\n");
               fprintf (stderr, "                 Number of particles in header (%i) and body (%i)\n", np, n);
               fprintf (stderr, "                 for run (%i) step (%i) do not agree.\n", cur->run, cur->step);
               exit(1);
               }
            cur->c[0][n] = scale[0] * (table[xlo]+table2[xhi]-off2) - cur->trans[0];
            cur->c[1][n] = scale[1] * (table[ylo]+table2[yhi]-off2) - cur->trans[1];
            cur->c[2][n] = scale[2] * (table[zlo]+table2[zhi]-off2) - cur->trans[2];
            if (!cur->selkeep)
               cur->t[n] = 1;
            n++;
            i++;
            }
         } while (i!=13 && xlo !=0 && xlo!='\n' && n!=np);
      iline++;
      } while (xlo !=0 && xlo!='\n' && n!=np);

   if (n<np)
      {
      fprintf (stderr, "ERROR (readrcv): Only %i of %i points read for step %i.\n",
                        n, np, cur->step);
      exit(1);
      }


   /*  Read RDF table  */
   fscanf (fin, "%li\n", &maxtable);
   tablefactor = maxtable/10000.0;
   ntable = 0;
   do
      {
      /* UNTIL END-OF-FILE  OR  UNTIL ASCII 0 Code FOUND */
      iline = 0;
      do
       /*     UNTIL 13 POSITIONS Read IN
                OR  END-OF-FILE  OR  0 FOUND    */
         {
         xlo = getc(fin);
         done = (xlo==13 || xlo==0 || xlo==EOF);
         if (!done)
            {
            xhi = getc(fin);
            iline++;
            if (xhi!=0 && xhi!=13 && xhi!=EOF)
               cur->rdftable[ntable++] =
                       tablefactor * (table[xlo]+table2[xhi]-off2);
            }
         done = (done || ntable==100);
         } while (iline<39 && !done);

      /*   Goto end of line  */
      while (getc(fin)!='\n');

      } while (!done);


   /*  Find min and maximum coordinate  */
   CalcMinMaxCoord (cur, FALSE, FALSE);


   /*  EVERYTHING OK:  RETURN not eof  */
   return 0;
   }
/*  END READRCV  */




void writercv (FILE *fout, stepdata *s, int selfl)
   {
   long maxtbl;
   int i,j,k, nchar, code;
   int line[79];
   int base = 100;
   int off  =  27;
   int res  = 10000;
   int nbin = 100;
   float boxsc[3];
   float tr[3];
   int   bx[3];

   /*
   Find min and max coordinate
    - use only selected particles if selfl
    - use repeating boundary conditions
   */
   CalcMinMaxCoord (s, selfl, TRUE);

   /*  ADJUST BOX IF FREE SURFACE      */
   /*  ADJUST PARTICLE TRANSLATION     */
   /*  SET BOUNDARY FLAG (FOR OUTPUT)  */

   for (i=0;i<3;i++) {
      if (!s->boundrep[i])
         {
         boxsc[i] = ( s->MaxCoord[i] - s->MinCoord[i] ) * 1.01;
         tr   [i] = -s->MinCoord[i];
         bx   [i] = 0;
         }
      else {
         boxsc[i] = s->box[i];
         tr   [i] = 0;
         bx   [i] = 1;
      }
   }

   /*  PLACE REPEATING PARTICLES INSIDE StepData  */
   for (i=0;i<3;i++)
      if (s->boundrep[i])
         for (j=0;j<s->np;j++)
            if (s->c[i][j] > s->box[i])
               s->c[i][j] -= s->box[i];
            else if (s->c[i][j]<0)
               s->c[i][j] += s->box[i];


   /*  SEND HEADER DATA  */
   if (selfl)
      fprintf (fout, "  %li %li %f %f %f %i 3\n", s->run, s->step, s->time,
                     s->tempc, s->temp, s->nsel+1);
   else
      fprintf (fout, "  %li %li %f %f %f %i 3\n", s->run, s->step, s->time,
                     s->tempc, s->temp, s->np+1);