cor0.c

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

   for (i=0;i<NTITLE;i++)
      if (s->title[i]==NULL)
         fprintf (fout, "* title *\n");
      else
         fprintf (fout, "%s\n", s->title[i]);

   fprintf (fout, "%14e  %15e  %10.4f  %10.4f %13e %13e\n",
            s->mass, s->vibp, s->dtime, s->cutoff, s->strain, s->stress);
   fprintf (fout, "%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %2i %2i %2i\n",
            boxsc[0], boxsc[1], boxsc[2],
            tr[0], tr[1], tr[2], bx[0], bx[1], bx[2]);



   /*  SEND CODED PARTICLE COORDINTES  */
   nchar = 0;
   for (i=0;i<s->np;i++) {
      if (!selfl || s->sel==NULL || s->sel[i]) {
         for (j=0;j<3;j++) {
            code = res * (s->c[j][i] + tr[j]) / boxsc[j];
            line[nchar++] = code % base + off;
            line[nchar++] = code / base + off;
         }
         if (nchar==78) {
            line[nchar++] = '.';
            for (k=0;k<nchar;k++)
               putc (line[k], fout);
            putc ('\n', fout);
            nchar = 0;
          }
          else if (nchar>78) {
             fprintf (stderr,
                "INTERNAL ERROR (File %s Line %i): nchar (%i) > 78.\n", 
						__FILE__, __LINE__, nchar);
             exit(1);
          }


      }
   }

   /*  SEND REMAINING CHARACTERS  */
   if (nchar>0) {
      for (k=0;k<nchar;k++)
         putc (line[k], fout);
      putc ('\n', fout);
   }


   /*  RDF TABLE  */
   nbin  = 100;
   maxtbl = s->rdftable[0];
   for (i=1;i<nbin;i++)
      if (maxtbl<s->rdftable[i])
         maxtbl = s->rdftable[i];
   /* adjust maxtbl to avoid overflow  */
   maxtbl += 1;
   if (maxtbl==0)
      maxtbl = 1;
   fprintf (fout, "%li\n", maxtbl);
   nchar=0;
   for (i=0;i<nbin;i++) {
      code = res * s->rdftable[i] / maxtbl + 0.5;
      line[nchar++] = code % base + off;
      line[nchar++] = code / base + off;
      if (nchar==78) {
         line[nchar++] = 0;
         for (k=0;k<nchar;k++)
            putc (line[k], fout);
         putc ('\n', fout);
         nchar = 0;
      }
      else if (nchar>78) {
         fprintf (stderr,
            "INTERNAL ERROR 2 (File %s Line %i): nchar (%i) > 78.\n", 
					__FILE__, __LINE__, nchar);
         exit(1);
      }
   }

   /*  SEND REMAINING CHARACTERS  */
   if (nchar>0)  {
      for (k=0;k<nchar;k++)
         putc (line[k], fout);
      putc ('\n', fout);
   }

}
/*  END WRITERCV  */




/*  Free Step Data  */
void FreeStep (stepdata *cur)
   {
   /*  Free old Variables  */
   if (cur->c[X]!=NULL)
      free (cur->c[X]);
   if (cur->c[Y]!=NULL)
      free (cur->c[Y]);
   if (cur->c[Z]!=NULL)
      free (cur->c[Z]);
   if (cur->t    !=NULL)
      free (cur->t    );
   if (cur->cd[X]!=NULL)
      free (cur->cd[X]);
   if (cur->cd[Y]!=NULL)
      free (cur->cd[Y]);
   if (cur->cd[Z]!=NULL)
      free (cur->cd[Z]);
   if (cur->radius!=NULL)
      free (cur->radius);
   if (cur->color!=NULL)
      free (cur->color );
   if (cur->symbol!=NULL)
      free (cur->symbol);
   if (cur->fill!=NULL)
      free (cur->fill  );
   if (cur->dispfl!=NULL)
      free (cur->dispfl);
   }




/*  Allocate Step Data  (not coordinates and related arrays)  */
stepdata *make_step ()
   {
   return (stepdata *) calloc ((unsigned) 1, sizeof(stepdata) );
   }




void readtype (FILE *fin, stepdata *cur)
   {
   int i,n,nalloc;
   int InputType;
   fscanf (fin, "TYPELIST%i", &n);
   if (n>cur->np)
      {
      nalloc  = n;
      cur->np = n;
      }
   else
      {
      nalloc = cur->np;
      n      = cur->np;
      }

   /*  RE ALLOCATE STORAGE  */
   if (cur->t!=NULL)
      cur->t = (unsigned char *)
                  realloc (cur->t, nalloc*sizeof(unsigned char));

   /*  ALLOCATE NEW STORAGE  */
   else
      cur->t = (unsigned char *) calloc (nalloc, sizeof(unsigned char));

   /*  CHECK ALLOCATION  */
   if (cur->t==NULL)
      {
      fprintf (stderr, "ERROR in readtype:  Insufficient memory for %i types.\n",
         nalloc);
      exit(1);
      }

   /*  READ TYPES  */
   i = 0;
   while (EOF!=fscanf(fin,"%i", &InputType)  &&  i<n)
      {
      cur->t[i] = InputType-1;
      i++;
      }

   /*  Test for premature end of file  */
   if (i<n)
      {
      fprintf (stderr, "ERROR in readtype:  premature end of file.\n");
      exit(1);
      }
   }




void readilist(FILE *fin, stepdata *cur)
   {
   unsigned i,n,is,nalloc,*tind;
   fscanf (fin, "ILIST%i", &n);
   /*  ALLOCATE TEMPORARY INDEX STORAGE  */
   tind = (unsigned *)  calloc (n, sizeof(unsigned));
   if (tind==NULL)
      {
      fprintf (stderr,
      "ERROR in readilist:  Insufficient memory for allocating temporary indices.\n");
      exit(1);
      }

   /*  READ ILIST  */
   i = 0;
   while (EOF!=fscanf(fin,"%i", &is)  &&  i<n)
      {
      tind[i] = is;
      i++;
      }
   if (i<n)
      {
      fprintf (stderr, "ERROR in readilist:  premature end of file.\n");
      exit(1);
      }


   /*  FIND MAXIMUM INDEX  */
   nalloc = tind[0];
   for (i=1;i<n;i++)
      if (tind[i]>nalloc)
         nalloc = tind[i];


   /*  ADJUST STORAGE FOR SELECT FLAGS  */
   if (cur->np>nalloc)
      nalloc = cur->np;
   if (cur->sel==NULL)
      cur->sel = (unsigned char *) calloc (nalloc, sizeof(unsigned));
   else
      cur->sel = (unsigned char *) realloc(cur->sel, nalloc*sizeof(unsigned));
   if (cur->sel==NULL) {
      fprintf (stderr,
      "ERROR in readilist:  Insufficient memory for allocating select flags.\n");
      exit(1);
   }


   /*  SET SELECT FLAGS */
   for (i=0;i<n;i++)
      cur->sel[tind[i]-1] = 1;
   cur->nsel = n;


   /*  FREE STORAGE  */
   free(tind);
   }





/*
************************************************************************
Cor Routines
************************************************************************
*/

void WriteCorFile (FILE *OutputFile, stepdata *a, BOOLEAN UseSelect)
   {
   CorFileHeader_t FileHeader;
   int    ititle;
   int    ipart;
   int    idir;
   FLOAT  Scale[NDIR];
   char   OutputStr[STRING_LENGTH];
   char  *LineFeed = "\n\r";
   INT4   OutputNP;
   WORD2  WordCoord[NDIR];


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

   /*
   Calc min and max coordinates
    - heed the Select flag
    - do not use boundary conditions
   */
   CalcMinMaxCoord (a, UseSelect, FALSE);


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


   /*
   **********************************
   Write Text portions of File Header
   **********************************
   */

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

   /*  Set values of file header  */
   FileHeader.EndianWord   = EndianWord_m;
   FileHeader.Run          = a->run;
   FileHeader.Step         = a->step;
   FileHeader.Time         = a->time;
   FileHeader.BoundaryType = BoundaryType_m;
   FileHeader.Surf[X]      = !a->boundrep[X];
   FileHeader.Surf[Y]      = !a->boundrep[Y];
   FileHeader.Surf[Z]      = !a->boundrep[Z];
   FileHeader.Box[X]       = a->box[X];
   FileHeader.Box[Y]       = a->box[Y];
   FileHeader.Box[Z]       = a->box[Z];
   FileHeader.Etot         = a->etot ;
   FileHeader.Epot         = a->epot ;
   FileHeader.Ekin         = a->ekin ;
   FileHeader.Ebath        = a->ebath;
   FileHeader.Min[X]       = a->MinCoord[X];
   FileHeader.Min[Y]       = a->MinCoord[Y];
   FileHeader.Min[Z]       = a->MinCoord[Z];
   FileHeader.Max[X]       = a->MaxCoord[X];
   FileHeader.Max[Y]       = a->MaxCoord[Y];
   FileHeader.Max[Z]       = a->MaxCoord[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)

#if 0
   /*  Write file header  */
   fwrite   (&FileHeader, sizeof(FileHeader), 1, OutputFile);
#endif


   /*  Write particle types  */
   LOOP (ipart, a->np)
      if (!UseSelect || a->sel[ipart])
         fwrite (&a->t[ipart], sizeof(WORD1), 1, OutputFile);

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

   /*  Write particle coordinates  */
   LOOP (ipart, a->np)
      if (!UseSelect || a->sel[ipart])
         {
         LOOP (idir, NDIR)
            WordCoord[idir] =
               Scale[idir] * (a->c[idir][ipart] - a->MinCoord[idir]);
         fwrite (WordCoord, sizeof(WORD2), NDIR, OutputFile);
         }


   }
/*   ... WriteCorFile  */


int ReadCorFile (FILE *InputFile, stepdata *a)
   {
   int    ititle;
   int    ipart;
   int    idir;
   int    InputNP;
   double Scale[NDIR];
   char   InputChar;
   char   InputStr[STRING_LENGTH];
   CorFileHeader_t FileHeader;
   BOOLEAN IsReverse;
   WORD2   WordCoord[NDIR];
   double  Min[NDIR];
   double  Max[NDIR];

   /*  Determine if input file is a cor file  */
   GetNextStringFromFile (InputStr, STRING_LENGTH, InputFile);

   /*  Test for end-of-file  */
   if (feof(InputFile))
      return EOF;

   if (strncmp(InputStr, CompareStr_m, strlen(CompareStr_m))!=0)
      {
      printf ("ERROR READING COR FILE: \n");
      printf ("  Input file is not a COR file.\n");
      return EOF;
      }

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

   /*  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, sizeof(WORD1), 1, InputFile);
   if (InputChar != '\032')
      {
      printf ("ERROR:  COR file is corrupted.\n");
      return EOF;
      }

   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)

#if 0
   /*  Read file header  */
   fread (&FileHeader, sizeof(FileHeader), 1, InputFile);
#endif

   /*  Determine if data from different endian machine  */
   IsReverse = (FileHeader.EndianWord != EndianWord_m);