iofunc.c

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

#include <stdio.h>
#include <stdlib.h>
#include "particle.h"
#include "cdalloc.h"
#include "cdhouse.h"
#include "direct.h"
#include "strsub.h"
#include "rcvio.h"
#include "iofunc.h"


/*
************************************************************************
Input/Output Functions
************************************************************************
*/
/*  Open file  */
BOOLEAN OpenFile (char Name[], FILE **File)
	{
	char TempStr[255];

	*File = fopen (Name, "rb");
	if (*File==NULL)
		{
		strncpy(TempStr, Name, 255);
		sprintf 
			(Name, "Cannot open file: %s", TempStr);
		return(FALSE);
		}
	else return(TRUE);
	}

/*  Read next step  */
void ReadNextStep(FILE *InputFile,
						char Name[],
						Particle_t *p,
						long *Step,
						double CenterGeom[NDIR],
						char TypeList[MAX_TYPE_NAMES][255])
	{
	double (*Coord)[NDIR];
	int    ipart;
	int    idir;
	char   temp[255];
	int    i;

	if (InputFile==NULL)
		return;

	/* Initialize TypeList */
	for (i=0;i<MAX_TYPE_NAMES; ++i)
		{
		sprintf(temp, "%d", i);
		strcpy(TypeList[i], temp);
		}
	
	if (IsCorFile(Name) == TRUE)
 	   ReadCorFile (InputFile, p);
	else if (IsRCVFile(Name) == TRUE) 
		readrcv(InputFile, p);
	else ReadXYZFile(InputFile, p, TypeList);

	*Step = p->step;

	/*  Calculate Center of Geometry  */
	CenterGeom[X] = CenterGeom[Y] = CenterGeom[Z] = 0;
	Coord = (double(*)[NDIR]) p->cur;
	for (ipart=0;ipart < (p->np); ipart++)
		{
		CenterGeom[X] += Coord[ipart][X];
		CenterGeom[Y] += Coord[ipart][Y];
		CenterGeom[Z] += Coord[ipart][Z];
		}
	if (p->np > 0)
		{
		CenterGeom[X] /= p->np;
		CenterGeom[Y] /= p->np;
		CenterGeom[Z] /= p->np;
		}
	return;
	}

/* Decide file type */
BOOLEAN IsCorFile(char Name[])
{
	int length;

	length = strlen(Name);

	if (Name[length-3] == 'c')
		return(TRUE);
 	return (FALSE);
}

BOOLEAN IsRCVFile(char Name[])
{
	int length;

	length = strlen(Name);

	if (Name[length-3] == 'r')
		return(TRUE);
 	return (FALSE);
}

void ReadXYZFile(FILE *InputFile, Particle_t *p, char TypeList[MAX_TYPE_NAMES][255])
{
	int i, np;
	char type[255];
	char title[255];
	char xPos[255], yPos[255], zPos[255];
	int j, ntypes=1;
	int filepos;

	if (InputFile == NULL)
		{
		printf("\nError:  Input XYZ file is NULL.");
		return;
		}

	/* Test to see if EOF, if not rewind to previous location */
	filepos = ftell(InputFile);

	fgetc(InputFile);
	fgetc(InputFile);
	if (fgetc(InputFile) == EOF)
		return;
	else
		{
		fseek(InputFile, filepos, SEEK_SET);
		}

	fscanf(InputFile, "%d", &np);
	fscanf(InputFile, "%s", title);

	ReallocateParticle(p, np);

	p->np = np;
	p->run = -1;
	p->step = -1;
	p->surf[X] = p->surf[Y] = p->surf[Z] = 2;
	p->bcur[X] = p->bcur[Y] = p->bcur[Z] = 0;

	for (i=0; i<np; ++i)
		{
		fscanf(InputFile, "%s", type);
		fscanf(InputFile, "%s%s%s", xPos, yPos, zPos);
		for (j=0; j<10; ++j)
			{
			if (strcmp(type, TypeList[j]) == 0)
				{
			   p->type[i] = j;
				break; /* break out of loop */
				}
	      else if (j == ntypes-1)
				{
				strcpy(TypeList[ntypes-1], type);
			   p->type[i] = ntypes-1;
			   ntypes += 1;
				break; /* break out of loop */
				}
			}	
		/* Assume reading in units of Angstroms */
		/* Convert to centimeters */
		p->cur[NDIR*i + X] = 1e-8*atof(xPos);
		p->cur[NDIR*i + Y] = 1e-8*atof(yPos);
		p->cur[NDIR*i + Z] = 1e-8*atof(zPos);
		}

	return;
}
/*
************************************************************************
Miscellaneous Functions
************************************************************************
*/
void InitializeParticle(Particle_t **p)
{
	*p = palloc(10000,0);
	return;
}

void RotateMatrix (double Rot[NDIR][NDIR], double Del[NDIR][NDIR])
{
	int    idir;
	int    jdir;
	double TempInput[NDIR][NDIR];

	LOOP (idir, NDIR)
	LOOP (jdir, NDIR)
		TempInput[idir][jdir] = Rot[idir][jdir];
	LOOP (idir, NDIR)
	LOOP (jdir, NDIR)
		{
		Rot[idir][jdir] = 
			Del[idir][X]*TempInput[X][jdir] +
			Del[idir][Y]*TempInput[Y][jdir] +
			Del[idir][Z]*TempInput[Z][jdir];
		}
}

void GetSurfaceOrient(double Orient[NDIR][NDIR],
                      double Rotation[NDIR][NDIR],
							 double RefOrient[NDIR][NDIR])
{
	int i, j;

	for (i=0; i<3; ++i)
		for (j=0; j<3; ++j)
			Orient[i][j] = Rotation[i][0]*RefOrient[0][j] +
								Rotation[i][1]*RefOrient[1][j] +
								Rotation[i][2]*RefOrient[2][j];

	ScaleDirectionToMiller(Orient[0]);
	ScaleDirectionToMiller(Orient[1]);
	ScaleDirectionToMiller(Orient[2]);

	return;
}