cdstate.c

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

/*  READ AND WRITE CMD STATE	*/

/*
************************************************************************
Compiler Switches
************************************************************************
*/


/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "particle.h"
#include "cdalloc.h"
#include "cdhouse.h"
#include "cdsubs.h"

/*
************************************************************************
Defines
************************************************************************
*/
#define BOOLEAN int
#define TRUE 1
#define FALSE 0

#define ENDIAN_TEST_VALUE 1

/*  Version and release of state file	*/
#define STATE_VERSION 1
#define STATE_RELEASE 5


/*
************************************************************************
Macros
************************************************************************
*/
#define WRITE_SCALAR(FILE,SCALAR) \
	fwrite (&SCALAR, sizeof(SCALAR), 1, FILE);

#define WRITE_ARRAY(FILE,ARRAY, N) \
	fwrite (ARRAY, sizeof(ARRAY[0]), N, FILE);


#define READ_SCALAR(FILE,SCALAR) \
	fread (&SCALAR, sizeof(SCALAR), 1, FILE);

#define READ_ARRAY(FILE, ARRAY, N) \
	fread (ARRAY, sizeof(ARRAY[0]), N, FILE);



/*
************************************************************************
Global Variables
************************************************************************
*/
extern int	XMD_Version_g;
extern int	XMD_Release_g;

/*
************************************************************************
Module-Wide variables
************************************************************************
*/
static int StateVersion_m = STATE_VERSION;
static int StateRelease_m = STATE_RELEASE;
static int InputStateVersion_m;
static int InputStateRelease_m;

/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void writearray  (size_t, void  *, size_t, FILE *);
void readarray   (void **, FILE *);
void ReadParticleStructure  (Particle_t *, FILE *);
void WriteParticleStructure (Particle_t *, FILE *);
void ReadSimulationStructure	(Simulation_t *, FILE *);
void WriteSimulationStructure (Simulation_t *, FILE *);
void WriteVersion (FILE *);
void ReadVersion	(FILE *);
void WriteUserInfo(FILE *);
void ReadUserInfo (FILE *);
void WriteStringToBinaryFile	(FILE *, char *);
void ReadStringFromBinaryFile (FILE *, char *, int);
void WriteTitle	(FILE *, Particle_t *);
void ReadTitle 	(FILE *, Particle_t *);
void WriteDelimiter (FILE *);
void ReadDelimiter  (FILE *);

void WriteExternalVector (int NumVector, ExternalVector_t **Vector, FILE *fout);
void ReadExternalVector  (ExternalVector_t ***Vector, FILE *fin);


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

/*  WRITE STATE  */
void writestate (Particle_t *a, Simulation_t *s, char *fname)
	{
	int i;
	int TestWord = ENDIAN_TEST_VALUE;
	BOOLEAN TempBoolean;

	/*  DECLARE AND SET ARRAY SIZES	*/
	size_t np = a->np;
	size_t ng = a->ng;

	/*  OPEN OUTPUT FILE  */
	FILE *fout = fopen (fname, "wb");
	if (fout==NULL)
		return;


	/*  Write user readable info	*/

	/*  Write user readable version information	*/
	WriteUserInfo (fout);

	/* Write title */
	WriteTitle (fout, a);

	/*  Write ascii text delimiter  */
	WriteDelimiter (fout);


	/*  Write machine info	*/

	/*
	Write test word
	(to be used by computer reading file to check "endian" type)
	*/
	WRITE_SCALAR (fout, TestWord);

	/*  Write version information  */
	WriteVersion (fout);

	/*  WRITE ARRAY SIZES  */
	fwrite (&np, sizeof(np), 1, fout);
	fwrite (&ng, sizeof(ng), 1, fout);

	/*  STORE STRUCTURES  */
	WriteParticleStructure	 (a, fout);
	WriteSimulationStructure (s, fout);

	/*  WRITE PARTICLE ARRAYS	*/
	writearray	(3 * np, (void  *) a->cur, 	sizeof(a->cur	[0]), fout);
	writearray	(3 * np, (void  *) a->ref, 	sizeof(a->ref	[0]), fout);
	writearray	(3 * np, (void  *) a->ngh, 	sizeof(a->ngh	[0]), fout);
	writearray	(3 * np, (void  *) a->v,		sizeof(a->v 	[0]), fout);
	writearray	(3 * np, (void  *) a->f,		sizeof(a->f 	[0]), fout);
	writearray	(3 * np, (void  *) a->c2 , 	sizeof(a->c2	[0]), fout);
	writearray	(3 * np, (void  *) a->c3 , 	sizeof(a->c3	[0]), fout);
	writearray	(3 * np, (void  *) a->c4 , 	sizeof(a->c4	[0]), fout);
	writearray	(3 * np, (void  *) a->c5 , 	sizeof(a->c5	[0]), fout);
	writearray	(	  np, (void  *) a->type,	sizeof(a->type [0]), fout);
	writearray	(	  np, (void  *) a->rtype,	sizeof(a->rtype[0]), fout);
	writearray	(	  np, (void  *) a->mass,	sizeof(a->mass [0]), fout);
	writearray	(	  np, (void  *) a->eatom,	sizeof(a->eatom[0]), fout);
	writearray	(	  np, (void  *) a->tag, 	sizeof(a->tag	[0]), fout);
	writearray	(	  np, (void  *) a->Selection, sizeof(a->Selection[0]), fout);


	/*  Write ExternalVector_t structures	*/
/*test*/
#if 0
	WriteExternalVector ( np, a->ForcePtrList,	fout);
	WriteExternalVector ( np, a->SpringPtrList, fout);
#endif

	/*  Write Temp CLAMP info  */
	LOOP (i, NUM_CLAMP_TAG)
		WRITE_SCALAR(fout,a->UseClamp[i])
	LOOP (i, NUM_CLAMP_TAG)
		WRITE_SCALAR(fout,a->ClampTemp[i])
	LOOP (i, NUM_CLAMP_TAG)
		WRITE_SCALAR(fout,a->ClampStep[i])

	/*  
	Write Volume CLAMP info 
	( a->UseVolClamp replaced by a->BoxMotionAlgorithm==BMA_CLAMP )
	 */
	TempBoolean = (a->BoxMotionAlgorithm==BMA_CLAMP);
	WRITE_SCALAR(fout, TempBoolean);
	WRITE_SCALAR(fout, a->VolClampStep)
	WRITE_SCALAR(fout, a->BulkModulus)

#if 0
	writearray
		(NUM_CLAMP_TAG, (void  *) &(a->UseClamp[0]), sizeof(a->UseClamp [0]), fout);
	writearray
		(NUM_CLAMP_TAG, (void  *) &(a->ClampTemp[0]),sizeof(a->ClampTemp[0]), fout);
	writearray
		(NUM_CLAMP_TAG, (void  *) &a->ClampStep,sizeof(a->ClampStep[0]), fout);
#endif


	/*  Write pressure info  */
	WRITE_SCALAR(fout, a->BoxMotionGeometry)
	WRITE_SCALAR(fout, a->BoxMotionAlgorithm)
	WRITE_SCALAR(fout, a->EpotBox)
	WRITE_SCALAR(fout, a->BulkModulus)
	WRITE_SCALAR(fout, a->VolClampStep)
		
	WRITE_ARRAY (fout, a->Pressure,  NDIR)
	WRITE_ARRAY (fout, a->EkinBox,   NDIR)
	WRITE_ARRAY (fout, a->BoxMotion, NDIR*NDERIV)
	WRITE_ARRAY (fout, a->BoxForce,  NDIR)
	WRITE_ARRAY (fout, a->BoxMass,   NDIR)

	/*  Close output file  */
	fclose(fout);
}



int  
readstate
(
Particle_t **InputParticleStructure,
Simulation_t *s,
char *fname
)
	{
	size_t np;
	size_t ng;
	FILE *fin;
	int  TestWord;
	int  i;
	BOOLEAN IsInitializedPotential;
	BOOLEAN TempBoolean;
	Particle_t *a = *InputParticleStructure;

	/*  Save potential status  */
	IsInitializedPotential = a->IsInitializedPotential;

	/*  Open input file	*/
	fin = fopen (fname, "rb");
	if (fin==NULL)
		{
		ERROR_PREFIX
		printf ("Cannot open file <%s>\n", fname);
		CleanAfterError();
		}

	/*  Read user readable info  */

	/*  Read user readable version information  */
	ReadUserInfo (fin);

	/* Read title */
	ReadTitle (fin, a);

	/*  Read ascii text delimiter  */
	ReadDelimiter (fin);


	/*  Machine readable info	*/

	/*  Check "endian" type  */
	READ_SCALAR (fin, TestWord);

	if (TestWord != ENDIAN_TEST_VALUE)
		{
		ERROR_PREFIX
		printf ("State file was written by incompatible computer.\n");
		printf ("  (Wrong endian type.\n");
		CleanAfterError();
		}

	/*  Read Version Information	*/
	ReadVersion (fin);

	/*  Free arrays inside particle structure  */
	pfree (a);

	/*  Free particle structure itself	*/
	FREE (a)

	/*  READ ARRAY SIZES  */
	fread (&np, sizeof(np), 1, fin);
	fread (&ng, sizeof(ng), 1, fin);

	/*
	Allocate new particle structure
	(and set global value of limits of number of particles and neighbors)
	*/
	/*  Set initial neighbor array to small value  */
	a = palloc (np, 4);

	/*  Read particle record  */
	ReadParticleStructure (a, fin);

	/* Read simulation record */
	ReadSimulationStructure (s, fin);

	/*  Original data for stat file version 1 release 0  */

	/*  READ  PARTICLE ARRAYS	*/
	readarray  ( (void  **) &a->cur, 	fin);
	readarray  ( (void  **) &a->ref, 	fin);
	readarray  ( (void  **) &a->ngh, 	fin);
	readarray  ( (void  **) &a->v,		fin);
	readarray  ( (void  **) &a->f,		fin);
	readarray  ( (void  **) &a->c2,		fin);
	readarray  ( (void  **) &a->c3,		fin);
	readarray  ( (void  **) &a->c4,		fin);
	readarray  ( (void  **) &a->c5,		fin);
	readarray  ( (void  **) &a->type,	fin);
	readarray  ( (void  **) &a->rtype,	fin);
	readarray  ( (void  **) &a->mass,	fin);
	readarray  ( (void  **) &a->eatom,	fin);

	/*  Determine if masses are initialized  */
	a->IsInitializedMass = (a->mass != NULL);

	/*  Attempt to read data previous to version.release 1.2  */
	if (InputStateVersion_m==1 && InputStateRelease_m<2)
		{
		int i;
		BYTE *sel = NULL;
		BYTE *set = NULL;
		BYTE *fix = NULL;
		
		/*  Read temp and perm arrays  */
		readarray  ( (void  **) &sel	    , fin);
		readarray  ( (void  **) &a->tag	 , fin);
		readarray  ( (void  **) &set	    , fin);
		readarray  ( (void  **) &fix	    , fin);

		/*  Combine sel, set and fix into *Selection  */
		/*  Read set data into Selection[]  */
		if (set!=NULL)
			{
			LOOP (i, np)
				a->Selection[i] = set[i];
			}

		/*  Merge sel data into Selection[]  */
		if (sel!=NULL)
			{
			LOOP (i, np)
				{
				REMOVE_SELECT(i)
				if (sel[i])
					APPLY_SELECT(i)
				}
			}

		/*  Merge fix data into Selection[]  */
		if (fix!=NULL)
			{
			LOOP (i, np)
				{
				REMOVE_FIX(i)
				if (fix[i])
					APPLY_FIX(i)
				}
			}


		/*  Free temporary storage  */
		FREE(sel)
		FREE(set)
		FREE(fix)
		}

	else
		{
		readarray  ( (void  **) &a->tag	 , fin);
		readarray  ( (void  **) &a->Selection	 , fin);
		}


	/*  Following code added for state file version 1 release 1  */
	if (InputStateRelease_m > 0)
		{
		/*  Read ExternalVector_t structures  */
/*test*/
#if 0
		ReadExternalVector ( &a->ForcePtrList,  fin);
		ReadExternalVector ( &a->SpringPtrList, fin);
#endif
		}

	/*  read CLAMP info (state version 1.4)  */
	if (InputStateRelease_m > 3)
		{
		
		LOOP (i, NUM_CLAMP_TAG)
			READ_SCALAR(fin, a->UseClamp[i])
		LOOP (i, NUM_CLAMP_TAG)
			READ_SCALAR(fin, a->ClampTemp[i])
		LOOP (i, NUM_CLAMP_TAG)
			READ_SCALAR(fin, a->ClampStep[i])

#if 0
	READ_SCALAR(fin, a->UseVolClamp)
#endif
		READ_SCALAR(fin, TempBoolean);
		READ_SCALAR(fin, a->VolClampStep)
		READ_SCALAR(fin, a->BulkModulus)
		if (TempBoolean)
			{
			a->BoxMotionAlgorithm = BMA_CLAMP;
			}

#if 0
		readarray ( (void **) &(a->UseClamp[0]) , fin);
		readarray ( (void **) &(a->ClampTemp), fin);
		readarray ( (void **) &(a->ClampStep), fin);
#endif
		}


	/*  read PRESSURE info (state version 1.5)  */
	if (InputStateRelease_m > 4)
		{

		READ_SCALAR(fin, a->BoxMotionGeometry)
		READ_SCALAR(fin, a->BoxMotionAlgorithm)
		READ_SCALAR(fin, a->EpotBox)
		READ_SCALAR(fin, a->BulkModulus)
		READ_SCALAR(fin, a->VolClampStep)
			
		READ_ARRAY (fin, a->Pressure,  NDIR)
		READ_ARRAY (fin, a->EkinBox,   NDIR)
		READ_ARRAY (fin, a->BoxMotion, NDIR*NDERIV)
		READ_ARRAY (fin, a->BoxForce,  NDIR)
		READ_ARRAY (fin, a->BoxMass,   NDIR)

		}


	/*  Close input file  */
	fclose(fin);

	/*  Pass address of new pointer back to calling routine	*/
	*InputParticleStructure = a;

	/*  Restore potential status  */
	a->IsInitializedPotential = IsInitializedPotential;


	return 0;
	}


/*
************************************************************************
Local routines
************************************************************************
*/


/*  Write array and its length in bytes  */
void writearray
(
size_t NumElements,
void	 *ArrayPointer,
size_t ElementSize,
FILE	 *OutputFile
)
	{
	size_t NumChars;

	/*  If array holds data write it out  */
	if (NumElements>0 && ArrayPointer!=NULL && OutputFile!=NULL)
		{
		NumChars = NumElements * ElementSize;
		fwrite (&NumChars,	 sizeof(NumChars), sizeof(char), OutputFile);
		fwrite (ArrayPointer, ElementSize,		 NumElements,	OutputFile);
		}

	/*  .. else write a zero  */
	else
		{
		NumChars = 0;
		fwrite (&NumChars, sizeof(NumChars), sizeof(char), OutputFile);
		}
	}



/*  Read array and its length in bytes  */
void readarray (void **ArrayPointer, FILE *InputFile)
	{
	size_t NumChars;

	/*  Read array size	*/
	fread (&NumChars, sizeof(NumChars), 1, InputFile);
#ifdef DEBUG
	printf ("File %s Line %i : NumChars = %i\n",
		__FILE__, __LINE__, NumChars);
#endif

	/*  If not null array allocate memory and read array	*/