cdwrite.c

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

/*  NOTES:

    WARNING   if run cmd then add particles then re-run, will not have
              storage allocated for f,c0, etc.
*/


/*
***********************************************************************
File Includes
***********************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "iomngr.h"
#include "parse.h"
#include "strsub.h"
#include "cdhouse.h"
#include "cdsubs.h"
#include "cdalloc.h"
#include "cdcor.h"
#include "cdboun.h"
#include "rcvio.h"
#include "neigh.h"


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

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


/*
***********************************************************************
Macros
***********************************************************************
*/
/*  If pointer p is null, return "0", else return value string */
#define PFMT(p,v) (p ? sprintf(buffer,fmt,(v)),buffer : "0")



/*
***********************************************************************
External Variables
***********************************************************************
*/

extern Particle_t   *a;
extern Simulation_t *s;
extern FILE         *out;
extern LIST         *inlist;



/*
***********************************************************************
Module Wide Variables
***********************************************************************
*/
BOOLEAN	UseTypes_m;
int		NumTable_m;
int		Type1_m;
int      Type2_m;
long		*RdfTable_m;
BYTE		*Type_m;
double   MinRadius_m;
double   MaxRadius_m;
double	RdfTableFactor_m;


/*
***********************************************************************
Local Function Prototypes
***********************************************************************
*/
void PrintExternalVectorList
(
FILE             *OutputFile,
ExternalVector_t **VectorList,
BYTE             *Type,
SEL_T            *Select,
int               NumList,
BOOLEAN           UseSelect
);

void WriteXMOL (FILE *, Particle_t *, BOOLEAN);
void WritePDB  (FILE *, Particle_t *, BOOLEAN);

void CalcAvgMinMax
(
double *Array,
int    NumValue,
int    NumDim,
double *Avg,
double *Min,
double *Max,
BOOLEAN UseSelect
);

void CalcRdf (WORD, WORD *, double *, WORD);

char *get_pptr (Particle_t *a, char *symbol, int index, char *fmt);
char *get_pfmt (char *symbol);

/*
***********************************************************************
Exported Functions
***********************************************************************
*/


/*
**    READ WRITE INSTRUCTIONS
*/
void read_write   (char *instr)
	{
   int i;
   unsigned NumPartUsed;
   char *tokstr;
   char *OutputFileName;
   FILE *OutputFile;
	BOOLEAN UseAvg;
	BOOLEAN UseMin;
	BOOLEAN UseMax;
   BOOLEAN UseSelect;
   BOOLEAN UseOutputFile;
   BOOLEAN Continue;
   int     cnt;
	double  Avg[NDIR];
	double  Min[NDIR];
	double  Max[NDIR];
	double  Value;
	char    *FormatPtr = NULL;
	BOOLEAN UseFormatPtr = FALSE;
	BOOLEAN UseWrap = FALSE;
	BOOLEAN UseTag  = FALSE;
	char    buffer[NSTR];
	char    *bptr = NULL;
	char    *vptr = NULL;

   /*  Initialize flags  */
   UseSelect       = FALSE;
   UseOutputFile   = FALSE;
	UseAvg          = FALSE;
	UseMin          = FALSE;
	UseMax          = FALSE;

   /*  Initialize number of particle to use  */
   NumPartUsed = a->np;


   /*
   Check for FILE or SEL tokens
   Strip and parse leading token until its not SEL or FILE
   */

   do
      {
      /*  Initialize continue loop flag to false  */
      Continue = FALSE;

      /*  Get next token from  instr[]  */
      tokstr = strhed (&instr);

      /*  Test SELECT option  */
      if (!strcmpi (tokstr, "SEL"))
         {
         UseSelect = TRUE;
         Continue = TRUE;
         NumPartUsed = a->nsel;
			CheckForNoneSelected();
         }

		else if (!strcmpi(tokstr,"AVG"))
			{
			UseAvg   = TRUE;
			Continue = TRUE;
			}

		else if (!strcmpi(tokstr,"MIN"))
			{
			UseMin   = TRUE;
			Continue = TRUE;
			}

		else if (!strcmpi(tokstr,"MAX"))
			{
			UseMax   = TRUE;
			Continue = TRUE;
			}

      /*  Test DISK file option  */
      else if (!strcmpi (tokstr, "FILE"))
         {
         /*  Set disk file flag  */
         UseOutputFile = TRUE;
         /*  Get pointer to file name  */
         OutputFileName = strhed (&instr);
         /*  Try for another flag  */
         Continue = TRUE;
         }

		/*  Read optional format string  */
		else if (!strcmpi(tokstr,"FORMAT"))
			{
			FormatPtr = GetQuotedString (&instr, '"');
			ReplaceSpecialChar (FormatPtr);
			UseFormatPtr = TRUE;
			Continue = TRUE;
			}

      } while (Continue);

   /*
   Open output file
      some commands such as WRITE RCV, TYPELIST, ILIST don't rely on this
      this output file variable OutputFile, but instead rely on a local
      fout variable
   */

   /*  If output file specified, open disk file ... */
   if (UseOutputFile)
      {
      /*  Test for file append  */
      if (OutputFileName[0]=='+')
         OutputFile = fopen (OutputFileName+1, "at");
      else
         OutputFile = fopen (OutputFileName,   "wt");

      /*  Test for no file  */
      if (OutputFile==NULL)
         {
         printf ("ERROR:  Cannot open output file %s.\n", OutputFileName);
         }
   }
   /*  ... else use stardard out  */
   else
      OutputFile = stdout;




	/*  
	================================
	Examine tokstr for write command
	================================
	*/
   if (!strcmpi("BOX", tokstr))
      fprintf (OutputFile, "BOX %14.8f %14.8f %14.8f\n",
               1e8*a->bcur[0], 1e8*a->bcur[1], 1e8*a->bcur[2]);

   else if (!strcmpi("NEIGH", tokstr))
      {
      fprintf (OutputFile, "NEIGH %i\n", a->TotalNumNeighbors);
      }

   else if (!strcmpi("NP", tokstr))
      {
      fprintf (OutputFile, "NP %i\n", a->np);
      }

   else if (!strcmpi("PARTICLE", tokstr) || !strcmpi("POSITION", tokstr))
      {

		/*  Wrap particle positions  */
		WrapParticles(a);

		/*  Print position statistics  */
		if (UseAvg || UseMin || UseMax)
			{
			CalcAvgMinMax (a->cur, a->np, NDIR, Avg, Min, Max, UseSelect);

			if (UseAvg)
				fprintf
					(
					OutputFile,
					"AVG_POSITION %12.4le %12.4le %12.4le\n",
					Avg[X]*1e8,
					Avg[Y]*1e8,
					Avg[Z]*1e8
					);

			if (UseMin)
				fprintf
					(
					OutputFile,
					"MIN_POSITION %12.4le %12.4le %12.4le\n",
					Min[X]*1e8,
					Min[Y]*1e8,
					Min[Z]*1e8
					);

			if (UseMax)
				fprintf
					(
					OutputFile,
					"MAX_POSITION %12.4le %12.4le %12.4le\n",
					Max[X]*1e8,
					Max[Y]*1e8,
					Max[Z]*1e8
					);
			}


		/*  Print individual velocities  */
		else
			{
	      fprintf (OutputFile, "PARTICLE %8i\n", NumPartUsed);
			if (!UseFormatPtr)
				FormatPtr = "%6i %10.4f %10.4f %10.4f\n";
			LOOP (i, a->np)
				if (!UseSelect || IS_SELECT(i))
	         fprintf (OutputFile, FormatPtr, 
	                  a->type[i]+1,
	                  1e8*a->cur[X+3*i],
	                  1e8*a->cur[Y+3*i],
	                  1e8*a->cur[Z+3*i]);
			}
      }

	else if (!strcmpi("MASS", tokstr)) {
      fprintf (OutputFile, "MASS %8i\n", NumPartUsed);
		LOOP (i, a->np)
			if (!UseSelect || IS_SELECT(i))
				fprintf (OutputFile, "%6i %12.4e\n", a->type[i]+1, a->mass[i]*AVAG);
	}

   else if (!strcmpi("VELOCITY", tokstr))
      {
		if (a->v==NULL)
			{
			printf
				("WARNING:  Cannot print velocities, they are not intialized!\n");
         IncrementNumberOfWarnings();
			return;
			}

      /*  Check for s->dtime zero  */
      if (s->dtime==0.0)
         {
         printf ("ERROR:  Time step is zero, cannot calculate velocity.\n");
         CleanAfterError();
         }

		/*  Print velocity statistics  */
		if (UseAvg || UseMin || UseMax)
			{
			CalcAvgMinMax (a->v, a->np, NDIR, Avg, Min, Max, UseSelect);

			if (UseAvg)
				fprintf
					(
					OutputFile,
					"AVG_VELOCITY %12.4le %12.4le %12.4le\n",
					Avg[X]/s->dtime,
					Avg[Y]/s->dtime,
					Avg[Z]/s->dtime
					);

			if (UseMin)
				fprintf
					(
					OutputFile,
					"MIN_VELOCITY %12.4le %12.4le %12.4le\n",
					Min[X]/s->dtime,
					Min[Y]/s->dtime,
					Min[Z]/s->dtime
					);

			if (UseMax)
				fprintf
					(
					OutputFile,
					"MAX_VELOCITY %12.4le %12.4le %12.4le\n",
					Max[X]/s->dtime,
					Max[Y]/s->dtime,
					Max[Z]/s->dtime
					);
			}


		/*  Print individual velocities  */
		else
			{

	      fprintf (OutputFile, "VELOCITY %8i\n", NumPartUsed);
			LOOP (i, a->np)
				if (!UseSelect || IS_SELECT(i))
	            fprintf (OutputFile, "%4i %12.4le %12.4le %12.4le\n",
	               a->type[i]+1,
	               a->v[X+3*i]/s->dtime,
	               a->v[Y+3*i]/s->dtime,
	               a->v[Z+3*i]/s->dtime);
	      }
		
		}

   else if (!strcmpi("POSVEL", tokstr))
      {
		if (a->v==NULL)
			{
			printf
				("WARNING:  Cannot print velocities, they are not intialized!\n");
         IncrementNumberOfWarnings();
			return;
			}

      /*  Check for s->dtime zero  */
      if (s->dtime==0.0)
         {
         printf ("ERROR:  Time step is zero, cannot calculate velocity.\n");
         CleanAfterError();
         }

		/*  Print statistics  */
		if (UseAvg || UseMin || UseMax)
			{

			/*  Print position statistics  */
			CalcAvgMinMax (a->cur, a->np, NDIR, Avg, Min, Max, UseSelect);

			if (UseAvg)
				fprintf
					(
					OutputFile,
					"AVG_POSVEL %12.4le %12.4le %12.4le",
					Avg[X]*1e+8,
					Avg[Y]*1e+8,
					Avg[Z]*1e+8
					);

			if (UseMin)
				fprintf
					(
					OutputFile,
					"MIN_POSVEL %12.4le %12.4le %12.4le",
					Min[X]*1e+8,
					Min[Y]*1e+8,
					Min[Z]*1e+8
					);

			if (UseMax)
				fprintf
					(
					OutputFile,
					"MAX_POSVEL %12.4le %12.4le %12.4le",
					Max[X]*1e+8,
					Max[Y]*1e+8,
					Max[Z]*1e+8
					);

			/*   Velocity statistics  */
			CalcAvgMinMax (a->v, a->np, NDIR, Avg, Min, Max, UseSelect);

			if (UseAvg)
				fprintf
					(
					OutputFile,
					" %12.4le %12.4le %12.4le\n",
					Avg[X]/s->dtime,
					Avg[Y]/s->dtime,
					Avg[Z]/s->dtime
					);

			if (UseMin)
				fprintf
					(
					OutputFile,
					" %12.4le %12.4le %12.4le\n",
					Min[X]/s->dtime,
					Min[Y]/s->dtime,
					Min[Z]/s->dtime
					);

			if (UseMax)
				fprintf
					(
					OutputFile,
					" %12.4le %12.4le %12.4le\n",
					Max[X]/s->dtime,
					Max[Y]/s->dtime,
					Max[Z]/s->dtime
					);
			}


		/*  Print individual positions and velocities  */
		else
			{

	      fprintf (OutputFile, "POSVEL %8i\n", NumPartUsed);
			LOOP (i, a->np)
				if (!UseSelect || IS_SELECT(i))
	            fprintf (OutputFile, 
					"%4i %14.7le %14.7le %14.7le %14.7le %14.7le %14.7le\n",
	               a->type[i]+1,
						a->cur[X+3*i]*1e8,
						a->cur[Y+3*i]*1e8,
						a->cur[Z+3*i]*1e8,
	               a->v  [X+3*i]/s->dtime,
	               a->v  [Y+3*i]/s->dtime,
	               a->v  [Z+3*i]/s->dtime);
	      }
		
		}

   else if (!strcmpi("DISP", tokstr))
      {
		if (a->disp==NULL)
			{
			printf ("WARNING:  No displacements.  ");
			printf ("Did you run the DISP command first?\n");
         IncrementNumberOfWarnings();
			}
      fprintf (OutputFile, "DISP READ %8i\n", NumPartUsed);
		LOOP (i, a->np)
			{
			if (!UseSelect || IS_SELECT(i))
				{
				if (a->disp==NULL)
					fprintf (OutputFile, "0.0 0.0 0.0\n");
				else
					fprintf