rcvio.c

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

/*

SOURCE CODE FOR READING AND WRITING RCV FILES
Uses particle 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;.
29 Jul 1992 - modified RCVIO0 (stepdata based) to particle based for PMD
 4 Feb 1993 - modify NSTR value (from 81 defined in PARTICLE.H to 256).
	Dec 1993 - change cur to double *, array of x y z x y z .. interleaved
 9 Dec 1993 - make default type 0
 4 Dec 1996 - readrcv() - correct box calculation when using free surfaces.
*/

/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include "strsub.h"
#include "cdhouse.h"
#include "cdalloc.h"
#include "rcvio.h"

/*
************************************************************************
Defines
************************************************************************
*/
#define NSTR 256
#define ANG_TO_CM 1e8

#define X 0
#define Y 1
#define Z 2
#define NDIR 3

#define BOOLEAN int
#define FALSE 0
#define TRUE  1


/*
************************************************************************
Module 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 				 };


/*
************************************************************************
Exported Functions
************************************************************************
*/
int readrcv (FILE *fin, Particle_t *a)
	{
	int		 Resolution  = 10000;
	double	 scale[NDIR];
	double    CoordRange[NDIR];
	double	 tablefactor;
	unsigned  np;
	unsigned  i;
	unsigned  n;
	unsigned  iline;
	unsigned  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 %lf %lf %lf %i %i\n",
		&a->run, &a->step, &a->time, &a->tempc, &a->temp,
		&np, &a->ndim);
	if (teof==EOF)
		return(EOF);
	/*  AJUST NUMBER OF POINTS  */
	np--;
	LOOP (i, 8)
		{
		fgets (instr, NSTR, fin);
		/*  COPY FIRST NTITLESTR CHARACTERS TO TITLE  */
		n = 0;
		while (instr[n] && instr[n]!='\n' && n<NTITLESTR-1)
			{
			a->title[i][n] = instr[n];
			n++;
			}
		a->title[i][n] = 0;
      }
						  /*	READ MISC DATA LINE	*/
	fgets  (instr, NSTR, fin);
	tstr = instr;
	/*  Skip mass (obsolete)  */
	dblstr (&tstr);
	a->vibp	 = dblstr (&tstr);
	a->dtime  = dblstr (&tstr);
	a->cutoff = dblstr (&tstr);
	/*  Don't read in strain or stress (obsolete)  */

	/*  Read Coordinate Range	*/
	fgets  (instr, NSTR, fin);
	tstr = instr;
	LOOP (i, NDIR)
		CoordRange[i] = 1e-8*dblstr(&tstr);

	/*  READ Trans IF PRESENT	*/
	if (*tstr!=0)
		LOOP (i, NDIR)
			a->trans[i] = 1E-8*dblstr(&tstr);
	else
		LOOP (i, NDIR)
			 a->trans[i] = 0;

	/*  BOUNDARY IF PRESENT  */
	if (*tstr!='\0')
		{
		LOOP (i, NDIR)
			{
			a->surf[i] = intstr(&tstr);
			a->surf[i] = !a->surf[i];
			}
		}
	else
		LOOP (i, NDIR)
			{
			a->surf [i] = TRUE;
			a->trans[i] = 0.0;
			}

	/*  Set box size	*/
	LOOP (i, NDIR)
		{
		if (a->surf[i])
			a->bcur[i] = CoordRange[i] - a->trans[i];
		else
			a->bcur[i] = CoordRange[i];
		}


	/*  Calculate scale for converting from code to coordinates  */
	LOOP (i, NDIR)
		scale[i] = CoordRange[i] / Resolution;

	/*  Reallocate particles if new number different from previous number  */
	if (np != a->np)
		{
		ReallocateParticle (a, np);
		a->np = np;
		}

	/*  Reallocate tag info if keeping it	*/
	if (a->selkeep)
		{
		if (a->tag==NULL)
			ALLOCATE (a->tag, BYTE, a->NumPartAlloc)
		}

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

		/*  Remove tag  */
		FREE (a->tag)

		/*  Clear select	*/
		LOOP (i, a->np)
			REMOVE_SELECT(i);

		}


	/*  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))
					{
					printf ("ERROR:  Premature end of file in step %i.\n", a->step);
					exit(1);
					}
				/*  TEST FOR N TOO BIG	*/
				if (n==np)
					{
					printf ("ERROR (readrcv): Data is corrupted.\n");
					printf (
						" Number of particles in header (%i) and body (%i)\n", np, n);
					printf (
						" for run (%i) step (%i) do not agree.\n", a->run, a->step);
					exit(1);
					}
				a->cur[X + 3*n] =
					scale[X] * (table[xlo]+table2[xhi]-off2+0.5) - a->trans[X];
				a->cur[Y + 3*n] =
					scale[Y] * (table[ylo]+table2[yhi]-off2+0.5) - a->trans[Y];
				a->cur[Z + 3*n] =
					scale[Z] * (table[zlo]+table2[zhi]-off2+0.5) - a->trans[Z];
				if (!a->selkeep)
					a->type[n] = 0;
				n++;
				i++;
				}
			}
			while (i!=13 && xlo !=0 && xlo!='\n' && n!=np);
		iline++;
		}
		while (xlo !=0 && xlo!='\n' && n!=np);

	if (n<np)
		{
		printf ("ERROR (readrcv): Only %i of %i points read for step %i.\n",
								n, np, a->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)
					a->rdftable[ntable++] =
							  tablefactor * (table[xlo]+table2[xhi]-off2);
				}
			done = (done || ntable==100);
			}
			while (iline<39 && !done);
		while (getc(fin)!='\n');
		}
		while (!done);

	/*  EVERYTHING OK:  RETURN 0	*/
	return 0;

	}	/* readrcv()  */



void writercv (FILE *fout, Particle_t *a, BOOLEAN selfl)
	{
	long maxtbl;
	int nchar;
	int code;
	int ipart;
	int idir;
	int ibin;
	int ichar;
	int ititle;
	int line[79];
	int base = 100;
	int off	=	27;
	int res	= 10000;
	int nbin = 100;
	double  Component;
	double  MinCoord	[NDIR];
	double  MaxCoord	[NDIR];
	double  CoordRange[NDIR];
	double  CoordTrans[NDIR];
	BOOLEAN RepeatingBoundary[NDIR];
	double  (*Coord)[NDIR];
	double  *Box;


	/*  Set Coord	*/
	Coord = (double (*)[NDIR])  a->cur;
	Box	= a->bcur;

	/*  FIND MINIMUM AND MAXIMUM COORDINATE  */
	for (idir=X; idir<=Z; idir++)
	{
		MinCoord[idir] = MaxCoord[idir] = Coord[0][idir];
		for (ipart=1;ipart<a->np;ipart++)
		{
			if (MinCoord[idir]>Coord[ipart][idir])
				MinCoord[idir] = Coord[ipart][idir];
			if (MaxCoord[idir]<Coord[ipart][idir])
				MaxCoord[idir] = Coord[ipart][idir];
		}
	}

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

	for (idir=X; idir<=Z; idir++)
		{
		if (a->surf[idir])
			{
			CoordRange[idir] = (MaxCoord[idir]-MinCoord[idir]) * 1.01;
			CoordTrans[idir] = -MinCoord[idir];
			RepeatingBoundary[idir] = FALSE;
			}
		else
			{
			CoordRange[idir] = Box[idir];
			CoordTrans[idir] = 0;
			RepeatingBoundary[idir] = TRUE;
			}
		}

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

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

#ifdef OLD
	fprintf (fout, "%14e  %15e  %10.4f  %10.4f %13e %13e\n",
				1.0, a->vibp, a->dtime, a->cutoff, a->strain, a->stress);
#endif
	/*  Write out dummy value for mass, stress and strain (1.0, 0, 0)  */
	fprintf (fout, "1.0 %15le %10.4lf %10.4lf 0.0 0.0\n",
		a->vibp, a->dtime, a->cutoff);

	fprintf (fout, "%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %2i %2i %2i\n",
		ANG_TO_CM*CoordRange[X],
		ANG_TO_CM*CoordRange[Y],
		ANG_TO_CM*CoordRange[Z],
		ANG_TO_CM*CoordTrans[X],
		ANG_TO_CM*CoordTrans[Y],
		ANG_TO_CM*CoordTrans[Z],
		RepeatingBoundary[X],
		RepeatingBoundary[Y],
		RepeatingBoundary[Z]
		);



	/*  SEND CODED PARTICLE COORDINTES	*/
	nchar = 0;
	LOOP (ipart, a->np)
		{
		if (!selfl || IS_SELECT(ipart))
			{

			for (idir=X; idir<=Z; idir++)
				{

				/*  Isolate current position component  */
				Component = Coord[ipart][idir];

				/*  Wrap component if repeating boundary	*/
				if (RepeatingBoundary[idir])
					{
					if (Component < 0)
						Component += Box[idir];
					else if (Component >= Box[idir])
						Component -= Box[idir];
					}

				/*  Map component into integer (between 0 and 10000  */
				code = res *
					(Component + CoordTrans[idir]) / CoordRange[idir];

				/*  Code integer into ASCII character	*/
				line[nchar++] = code % base + off;
				line[nchar++] = code / base + off;
				}

			if (nchar==78)
				{
				line[nchar++] = '.';
				for (ichar=0; ichar<nchar; ichar++)
					putc (line[ichar], fout);
				putc ('\n', fout);
				nchar = 0;
				}

			ASSERT(nchar<=78)
			}
		}

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

	/*  RDF TABLE	*/
	nbin	= 100;
	maxtbl = a->rdftable[0];
	for (ibin=1; ibin<nbin ;ibin++)
		if (maxtbl<a->rdftable[ibin])
			maxtbl = a->rdftable[ibin];

	/* adjust maxtbl to avoid overflow	*/
	maxtbl += 1;
	if (maxtbl==0)
		maxtbl = 1;
	fprintf (fout, "%li\n", maxtbl);
	nchar=0;
	for (ibin=0;ibin<nbin;ibin++)
		{
		code = res * a->rdftable[ibin] / maxtbl + 0.5;
		line[nchar++] = code % base + off;
		line[nchar++] = code / base + off;
		if (nchar==78)
			{
			line[nchar++] = 0;
			for (ichar=0; ichar<nchar; ichar++)
				putc (line[ichar], fout);
			putc ('\n', fout);
			nchar = 0;
			}
		ASSERT(nchar<=78)
		}

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

	}
/*  END WRITERCV	*/


void readtype (FILE *fin, Particle_t *a)
	{
	int itype;
	int NumType;
	int InputType;
	char InputString[NSTR];
	char *TempString;

	/*  Read header  TYPELIST n	*/
	fgets (InputString, NSTR, fin);
	TempString = InputString;

	/* Strip leading token	*/
	strhed (&TempString);

	/*  Read number  */
	NumType = intstr (&TempString);

	/*  READ TYPES  */
	itype = 0;
	while
	(
	EOF != fscanf(fin,"%i", &InputType)  &&
	itype<NumType
	)
		{
		a->type[itype] = InputType;
		itype++;
		}

	/*  Check for corret number of types read  */
	if (itype<NumType)
		{
		printf ("ERROR (readtype):  premature end of file.\n");
		exit(1);
		}
	}

void readilist(FILE *fin, Particle_t *a)
	{
	unsigned i,n,is,nalloc;
	unsigned *tind = NULL;

	/*  Read number of indices from file  */
	fscanf (fin, "ILIST%i", &n);

	/*  Allocate space for indices  */
	ALLOCATE (tind, unsigned, n)

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

	/*  FIND MAXIMUM INDEX	*/
	nalloc = tind[0];
	LOOP (i, n)
		if (tind[i]>nalloc)
			nalloc = tind[i];

	/*  Check for maximum index value too large	*/
	if (nalloc > a->np)
		{
		printf ("ERROR: ");
		printf ("Value of largest index exceeds current number of particles.\n");
		CleanAfterError();
		exit(1);
		}

	/*  Clear select  */
	LOOP (i, a->np)
		REMOVE_SELECT(i);

	/*  Set select for indices that were read  */
	LOOP (i, n)
		APPLY_SELECT(tind[i]-1);

	a->nsel = n;

	/*  FREE STORAGE	*/
	FREE(tind)
	}

/*
************************************************************************
Local functions
************************************************************************
*/
int translate (int lo, int hi) {
	return (table[lo] + base*table[hi]- off2);
}