xmd Code

/*
xmd - molecular dynamics for metals and ceramics
By Jonathan Rifkin <jon.rifkin@uconn.edu>
Copyright 1995-2004 Jonathan Rifkin

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/


/*  NOTES:

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


/*
***********************************************************************
About
***********************************************************************
*/
/*
2006-03-01  Added NOHEAD option for Bin Li.
2006-09-11  Added CTABLE command.
*/
/*
***********************************************************************
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;
	BOOLEAN WriteHead = TRUE;
	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;
			}

		else if (!strcmpi(tokstr,"NOHEAD"))
			{
			WriteHead = FALSE;
			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");

      /*  Quit if cannot open output file.  */
      if (OutputFile==NULL)
         {
         printf ("ERROR:  Cannot open output file %s.\n", OutputFileName);
		 CleanAfterError();
         }
   }
   /*  ... 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
			{
			if (WriteHead) 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("CTABLE", tokstr) )
      {

      if (IsBlank(instr)) {
			tokstr = "CTABLE";
		} else {
	      tokstr = strhed (&instr);
		}

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

		if (!UseFormatPtr)
			FormatPtr = "%6i %10.4f %10.4f %10.4f\n";
		LOOP (i, a->np)
			if (!UseSelect || IS_SELECT(i)) {
				fprintf (OutputFile, "%s %ld %e %d ", tokstr, a->step, a->time, i+1);
				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)) {
      if (WriteHead) 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
			{

	      if (WriteHead) 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
			{

	      if (WriteHead) 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();
			}
      if (WriteHead) 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
						(OutputFile,
						"%10.4f %10.4f %10.4f\n",
						1e8*a->disp[X+3*i],
						1e8*a->disp[Y+3*i],
						1e8*a->disp[Z+3*i]
					);
				}
			}
      }


   else if (!strcmpi("EPOT", tokstr) || !strcmpi("ENERGY", tokstr) )
      {
		double Epot = 0.0;
		double NumPart = 0;

		/*  calculate latest energy  */
      energy_all (a, s, UseSelect);

		/*  test for mass and velocity present  */
		if (a->eatom!=NULL && a->np!=0)
			{
			LOOP (i, a->np)
				{
				if (!UseSelect || IS_SELECT(i))
					{
					Epot += a->eatom[i];
					NumPart++;
					}
				}
			}

		/*  prepare energy  */
		if (NumPart==0)
			Epot = 0.0;
		else
   		Epot = s->eunit*Epot/NumPart;

  		fprintf (OutputFile, "EPOT %17.9e\n", Epot);
      }

   else if (!strcmpi("EATOM", tokstr))
      {
		double ekin;
		double (*Vel)[NDIR] = (double (*)[NDIR]) a->v;

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

		/*  Get latest energy  */
      energy_all (a, s, UseSelect);

		/*  Print results  */
      if (WriteHead) fprintf (OutputFile, "EATOM %i\n", NumPartUsed);

		/*  Print verbose comments  */
		if (WriteHead && s->UseVerboseOutputFile)
			{
			fprintf (OutputFile, 
   	    "# type          x          y          z epot (%6s) ekin (%6s)\n",
				s->eulabel, s->eulabel);
			fprintf (OutputFile, "#----- ---------- ----------");
			fprintf (OutputFile, " ---------- ------------- -------------\n");
			}

		/*  Print energies  */
		LOOP(i, a->np)
		if (!UseSelect || IS_SELECT(i))
         {
         fprintf
				(
				OutputFile, "%6i %10.4lf %10.4lf %10.4lf %+13.6e",
            a->type[i]+1,
            1e8*a->cur[X+3*i],
            1e8*a->cur[Y+3*i],
            1e8*a->cur[Z+3*i],
				a->eatom[i]*s->eunit
				);

			/*  Calculate kinetic energy if appropriate  */
			if (s->dtime==0 || Vel==NULL || IS_FIX(i))
				fprintf (OutputFile, "\n");
			else
				{
				ekin = 0.5 * a->mass[i] *
					(SQR(Vel[i][X]) + SQR(Vel[i][Y]) + SQR(Vel[i][Z]))
					/ SQR(s->dtime);
				fprintf (OutputFile, " %+13.6e\n", ekin*s->eunit);
				}
         }
      }

   else if (!strcmpi("ETOT", tokstr))
      {
      fprintf (OutputFile, "ETOT %14.6e\n", (a->epot+a->ekin)*s->eunit/a->np);
      }

	/*  Temperature of particles which are both SELECTED and FREE  */
   else if (!strcmpi("TEMP", tokstr) )
		{
 	   fprintf 
			( 
			OutputFile, 
			"TEMP %14.6e\n", 
			GetTemperature(a,s,UseSelect)
			);
		}

	else if (!strcmpi("EKIN", tokstr) )
      {
		
      /*  KE of selected particle per atom  */
		if (UseSelect) {
			Value = s->eunit*GetKE(a, s, UseSelect);
		/*  Kinetic Energy of PARTICLES and BOX divided by ALL Part */
		} else {
			CalcKineticEnergy(a,s);
			Value = s->eunit*a->ekin/a->np;
		}

  	   fprintf (OutputFile, "EKIN %14.6e\n", Value);
      }

   else if (!strcmpi("DTIME", tokstr))
      {
      fprintf (OutputFile, "DTIME %14.6e\n", s->dtime);
      }

   else if (!strcmpi("TIME", tokstr))
      {
      fprintf (OutputFile, "TIME %14.6e\n", a->time);
      }

   else if (!strcmpi("CSTEP", tokstr))
      {
      fprintf (OutputFile, "CSTEP %14.6e\n", s->cstep);
      }

   else if (!strcmpi("COR", tokstr))
      {
      FILE *fout;

      tokstr = strhed (&instr);

		/*  Check for TAG option */
		if (!strcmpi("TAG",tokstr)) {
			tokstr = strhed (&instr);
			UseTag = TRUE;
			}

      if (*tokstr=='+')
         fout = fopen (tokstr+1, "ab");
      else
         fout = fopen (tokstr, "wb");
      if (fout==NULL)
         {
         printf ("ERROR:  Cannot open file %s\n", tokstr);
         return;
         }
      WriteCorFile (fout, a, UseSelect, UseTag);
      fclose (fout);
      }


   else if (!strcmpi("STRESS", tokstr))
      {
      BOOLEAN  DeallocForce;
		BOOLEAN  OldStressState;

      /*  Determine if force should be deallocated  */
      DeallocForce = (a->f == NULL);

		/*  Save previous stress output state  */
		OldStressState = s->StressStepOutput.Save;
		s->StressStepOutput.Save = TRUE;

		/*  Set stress select flag  */
		s->UseSelectStress = UseSelect;


      /*  Allocate space if needed  */
		if (a->f==NULL)
			{
			ALLOCATE (a->f, double, NDIR*a->NumPartAlloc)
			}

      /*  Calculate force  */
      ( (PotForce_f *) s->forcepnt ) (a, s);

      /*  Print Stress */
		fprintf
			(
			OutputFile,
			"STRESS %10.3le %10.3le %10.3le %10.3le %10.3le %10.3le\n",
			a->Stress[XX]/a->np,
			a->Stress[YY]/a->np,
			a->Stress[ZZ]/a->np,
			a->Stress[YZ]/a->np,
			a->Stress[ZX]/a->np,
			a->Stress[XY]/a->np
			);

      /*  Deallocate force if needed  */
      if (DeallocForce)
         FREE(a->f)

		/*  Reset old stress state  */
		s->StressStepOutput.Save = OldStressState;

		/*  Reset stress select flag  */
		s->UseSelectStress = FALSE;
      }


   else if (!strcmpi("FORCE", tokstr))
      {
      BOOLEAN  DeallocForce;

      /*  Determine if force should be deallocated  */
      DeallocForce = (a->f == NULL);

      /*  Allocate space if needed  */
		if (a->f==NULL)
			{
			ALLOCATE (a->f, double, NDIR*a->NumPartAlloc)
			}

      /*  Calculate force  */
      ( (PotForce_f *) s->forcepnt ) (a, s);

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

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

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

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


		/*  Print individual velocities  */
		else
			{
	      /*  Print force  */
	      if (WriteHead) fprintf (OutputFile, "FORCE %i\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->f[X + 3*i],
	                    a->f[Y + 3*i],
	                    a->f[Z + 3*i]);
			}

      /*  Deallocate force if needed  */
      if (DeallocForce)
         FREE(a->f)
      }

   else if (!strcmpi("EXTFORCE", tokstr))
      {
      /*  Test for existance of external force  */
      if (a->ForcePtrList== NULL)
         {
         printf ("WARNING: There is no applied external force.\n");
         IncrementNumberOfWarnings();
         }
      else
         {
			fprintf (OutputFile, "#  Force units are dynes\n");
         fprintf (OutputFile, "EXTFORCE %i\n", NumPartUsed);
         PrintExternalVectorList
            (
            OutputFile,
            a->ForcePtrList,
            a->type,
            a->Selection,
            a->np,
            UseSelect
            );
         }
      }


   else if (!strcmpi("EXTSPRING", tokstr))
      {
      /*  Test for existance of external spring  */
      if (a->SpringPtrList== NULL)
         {
         printf ("WARNING: There is no applied external spring.\n");
         IncrementNumberOfWarnings();
         }
      else
         {
			fprintf (OutputFile, "#  Spring constant units are erg/cm^2\n");
         fprintf (OutputFile, "EXTSPRING %i\n", NumPartUsed);
         PrintExternalVectorList
            (
            OutputFile,
            a->SpringPtrList,
            a->type,
            a->Selection,
            a->np,
            UseSelect
            );
         }
      }

   /*  Write data in PDB (Protein Data Bank) format  */
   else if (!strcmpi("PDB", tokstr))
      {
      FILE *fout;

      tokstr = strhed (&instr);
      if (*tokstr=='+')
         fout = fopen (tokstr+1, "at");
      else
         fout = fopen (tokstr, "wt");
      if (fout==NULL)
         {
         printf ("ERROR:  Cannot open file %s\n", tokstr);
         return;
         }

      /*  Write XMOL format  */
      WritePDB (fout, a, UseSelect);
      fclose (fout);
      }


   else if (!strcmpi("STATE", tokstr))
		{
		/*  Get first token for state file name  */
		tokstr = strhed (&instr);
		
		/*  Write state file  */
      writestate (a, s, tokstr);

		/*  
		Update Neighbor List so that subsequent trajectory
		will agree with trajectory that follows from this 
		state file
		*/
		UpdateNeighborList(a);
		}

   else if (!strcmpi("STEP", tokstr))
      fprintf (OutputFile, "STEP %li\n", a->step);

   else if (!strcmpi("RCV", tokstr))
      {
      FILE *fout;

      tokstr = strhed (&instr);
      if (*tokstr=='+')
         fout = fopen (tokstr+1, "at");
      else
         fout = fopen (tokstr, "wt");
      if (fout==NULL)
         {
         printf ("ERROR:  Cannot open file %s\n", tokstr);
         return;
         }
      writercv (fout, a, UseSelect);
      fclose (fout);
      }

   else if (!strcmpi("RDF", tokstr))
      {
		double   fac;
		double   r;
      BOOLEAN rep[NDIR];
      long   CumulativeCount;

      NumTable_m  = intstrf (&instr);
      MinRadius_m = 1e-8 * dblstrf (&instr);
      MaxRadius_m = 1e-8 * dblstrf (&instr);

      if (IsBlank(instr))
         UseTypes_m = FALSE;
      else
         {
         UseTypes_m = TRUE;
         Type1_m = intstrf (&instr)-1;
         Type2_m = intstrf (&instr)-1;
         }

      /*  Test for type in range  */
      if (UseTypes_m && (Type1_m < 0  ||  Type2_m < 0))
         {
         printf ("ERROR:  One or more types are out of range (<0).\n");
         CleanAfterError();
         }

      /*  ALLOCATE TABLE  */
      if (NumTable_m<=0)
         return;
		ALLOCATE (RdfTable_m, long, NumTable_m)
      rep[X] = !a->surf[X];
      rep[Y] = !a->surf[Y];
      rep[Z] = !a->surf[Z];

#if 0
/*
Old code for use with module bondsub.c
*/
      /*  Call RDF without types if type==0  */
      if (UseTypes_m)
         CalcTypeRDF
            (
            a->cur, a->type, a->np, rmin, rmax,
            type1, type2, a->bcur, rep, table, ntable
            );
      else
         CalcRDF
            (
            a->cur, a->np, rmin, rmax, a->bcur, rep, table, ntable
            );
#endif

/*
New code (dec 12, 1997) for use with CalcNeighborPerAtom()
*/

		/*
		Store values for recall via CallBack Function
		*/
		Type_m      = a->type;
		RdfTableFactor_m = NumTable_m / (MaxRadius_m-MinRadius_m);

		/*  Call neighbor routine (which calls callback function)  */
		if (s->nsearch==NEIGH_SEARCH_SORT)
		CalcNeighborPerAtom_Sort
			(
			(double (*)[NDIR]) a->cur,
			a->Selection,
			a->np,
			UseSelect,
			MinRadius_m,
			MaxRadius_m,
			a->bcur,
			rep,
			TRUE,
			CalcRdf
			);
		else if (s->nsearch==NEIGH_SEARCH_CELL) 
		CalcNeighborPerAtom_Cell
			(
			(double (*)[NDIR]) a->cur,
			a->Selection,
			a->np,
			UseSelect,
			MinRadius_m,
			MaxRadius_m,
			a->bcur,
			rep,
			TRUE,
			CalcRdf
			);
		else if (s->nsearch==NEIGH_SEARCH_CELL2) 
		CalcNeighborPerAtom_Cell2
			(
			(double (*)[NDIR]) a->cur,
			a->Selection,
			a->np,
			UseSelect,
			MinRadius_m,
			MaxRadius_m,
			a->bcur,
			rep,
			TRUE,
			CalcRdf
			);
		else {
			printf ("INTERNAL ERROR (read_write):  Invalid value for s->nsearch (%d)\n", s->nsearch);
			CleanAfterError();
		}
			

      if (WriteHead) fprintf
			(
			OutputFile,
			"RDF  %i  %lf  %lf\n",
			NumTable_m,
			1e8*MinRadius_m,
			1e8*MaxRadius_m
			);
		fac = (MaxRadius_m - MinRadius_m) / NumTable_m;
      CumulativeCount = 0;
		LOOP (i, NumTable_m)
         {
         r = 1e8 * (MinRadius_m + (i+0.5)*fac);
         CumulativeCount += RdfTable_m[i];
         fprintf
				(
				OutputFile,
				"%8.2lf %8li %8li\n",
				r,
				RdfTable_m[i],
				CumulativeCount
				);
         }

      /*  FREE STORAGE  */
		FREE (RdfTable_m)
      }

   else if (!strcmpi("RUN", tokstr))
      fprintf (OutputFile, "RUN %li\n", a->run);

   else if (!strcmpi("ILIST", tokstr))
      {
      FILE *fout;
      BOOLEAN IsFile;

      tokstr = strhed (&instr);

      /*  Is output to a file?  */
      IsFile = (tokstr[0]!='\0');


      /*  Open file  */
      if (IsFile)
         {

         /*  Open file  */
         if (*tokstr=='+')
            fout = fopen (tokstr+1, "at");
         else
            fout = fopen (tokstr, "wt");

         /*  Error opening file  */
         if (fout==NULL)
            {
            printf ("ERROR:  Cannot open file %s\n", tokstr);
            return;
            }

         }

      else
         fout = stdout;


      if (WriteHead) fprintf (fout,"ILIST %i\n", NumPartUsed);
      cnt = 0;
		LOOP (i, a->np)
         {
         /*  If select flag, print only selected particles  */
			if (!UseSelect || IS_SELECT(i))
            {
            fprintf (fout,"%4i", i+1);
            cnt++;

            /*  If end of line, print new line  */
            if (cnt==16)
               {
               cnt=0;
               fprintf (fout,"\n");
               }

            /*  .. else print blank for separator  */
            else
               fprintf (fout, " ");
            }
         }
      if (cnt)
         fprintf (fout,"\n");

      /*  Close file  */
      if (IsFile)
         fclose(fout);
      }

   else if (!strcmpi("TYPELIST", tokstr))
      {
      FILE *fout;
      BOOLEAN IsFile;

      tokstr = strhed (&instr);

      /*  Is output to a file?  */
      IsFile = (tokstr[0]!='\0');

      /*  Open file  */
      if (IsFile)
         {

         /*  Open file  */
         if (*tokstr=='+')
            fout = fopen (tokstr+1, "at");
         else
            fout = fopen (tokstr, "wt");

         /*  Error opening file  */
         if (fout==NULL)
            {
            printf ("ERROR:  Cannot open file %s\n", tokstr);
            return;
            }
         }

      else
         fout = stdout;


      if (WriteHead) fprintf (fout, "TYPELIST %i\n", NumPartUsed);
      cnt = 0;
		LOOP (i, a->np)
         {
			if (!UseSelect || IS_SELECT(i))
            {
            fprintf (fout," %4i", a->type[i]+1);
            cnt++;
            }
         if (cnt==16)
            {
            cnt=0;
            fprintf (fout,"\n");
            }
         }
      if (cnt)
         fprintf (fout,"\n");

      /*  Close file  */
      if (IsFile)
        fclose(fout);
      }


   /*  Write data for XMOL program  */
   else if (!strcmpi("XMOL", tokstr))
      {
      FILE *fout;

      tokstr = strhed (&instr);
      if (*tokstr=='+')
         fout = fopen (tokstr+1, "at");
      else
         fout = fopen (tokstr, "wt");
      if (fout==NULL)
         {
         printf ("ERROR:  Cannot open file %s\n", tokstr);
         return;
         }

      /*  Write XMOL format  */
      WriteXMOL (fout, a, UseSelect);
      fclose (fout);
      }

   /*  Write VAR format data  */
   else if (!strcmpi("VAR", tokstr)) {
	   if (WriteHead) fprintf (OutputFile, "VAR %s %8i\n", instr, NumPartUsed);
		LOOP (i, a->np)
			if (!UseSelect || IS_SELECT(i)) {
				/*  Run through list of variables  */
				strncpy (buffer, instr, NSTR);
				bptr = &buffer[0];
				vptr = strhed(&bptr);
				while (*vptr!=0) {
					fprintf (OutputFile, " %s", get_pptr(a,vptr,i,get_pfmt(vptr)) );
				vptr = strhed (&bptr);
				}
			fprintf (OutputFile, "\n");
			}
		}

	/*  Write random number seed  */
	else if (!strcmpi("SEED",tokstr)) {
		fprintf (OutputFile, "SEED %u\n", s->seed);
	}


   /*  TREAT AS NUMBER  */
   else
      {
      char *hedstr = tokstr;
      double x = dblstrf (&tokstr);
      fprintf (OutputFile, "%s %g\n", hedstr, x);
      }

   /*  If file opened then close it  */
   if (UseOutputFile)
      fclose (OutputFile);
	}



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

   /*  Empty list - print nothing  */
   if (VectorList==NULL)
      return;

   /*  Print vector list  */
   LOOP (ilist, NumList)
      if (!UseSelect || IS_SELECT2(Select,ilist))
         {
         fprintf (OutputFile, "%4i", Type[ilist]+1);
         fprintf (OutputFile, " %10.4e %10.4e %10.4e", 
               CM *VectorList[ilist]->Origin[X],
               CM *VectorList[ilist]->Origin[Y],
               CM *VectorList[ilist]->Origin[Z]
					);
         if (VectorList[ilist]==NULL)
            fprintf (OutputFile, " %10.4e %10.4e %10.4e\n",
					0.0, 0.0, 0.0);
         else
            fprintf (OutputFile, " %10.4e %10.4e %10.4e\n",
               VectorList[ilist]->Vector[X],
               VectorList[ilist]->Vector[Y],
               VectorList[ilist]->Vector[Z]
               );
         }

   }


void WriteXMOL (FILE *OutputFile, Particle_t *a, BOOLEAN UseSelect)
   {
   int ipart;
   int NumWrite;
   char *NameStr;

   if (UseSelect)
      NumWrite = a->nsel;
   else
      NumWrite = a->np;

   /*  Print number of particles written  */
   fprintf (OutputFile, "%i\n", NumWrite);

   /*  Print run and step number (this will be step "label")  */
   fprintf (OutputFile,
      "Run %li  Step %li  %s\n", a->run, a->step, GetCurrentTimeStr() );

   /*  Write particle coordinates  */
   LOOP (ipart, a->np)
	if (!UseSelect || IS_SELECT(ipart))
      {

      /*  Print type name  */
      NameStr = a->TypeNames[a->type[ipart]];
      if (NameStr == NULL)
         fprintf (OutputFile, " %i ", a->type[ipart]+1);
      else
         fprintf (OutputFile, " %s ", NameStr);

      /*  Print coordinates in angstroms  */
      fprintf
         (
         OutputFile,
         " %e %e %e\n",
         CM * a->cur[NDIR*ipart+X],
         CM * a->cur[NDIR*ipart+Y],
         CM * a->cur[NDIR*ipart+Z]
         );

      }

   }

void WritePDB (FILE *OutputFile, Particle_t *a, BOOLEAN UseSelect)
   {
   int ipart;
   int NumWrite;
   char *NameStr;

   /*  Find number of atoms  */
   if (UseSelect)
      NumWrite = a->nsel;
   else
      NumWrite = a->np;

   /*  Print run and step number (this will be step "label")  */
   fprintf (OutputFile, "REMARK Output from XMD\n");
   fprintf (OutputFile, "REMARK (c) 1996, Center for Materials Simulation,");
   fprintf (OutputFile, " University of Connecticut\n");
   fprintf (OutputFile,
      "REMARK Run %li  Step %li  %s\n", a->run, a->step, GetCurrentTimeStr() );
   fprintf (OutputFile, "REMARK  Number of Atoms = %i\n", NumWrite);


   /*  Write particle coordinates  */
   NumWrite = 0;
   LOOP (ipart, a->np)
	if (!UseSelect || IS_SELECT(ipart))
      {

      /*  Count number of particles written  */
      NumWrite++;

      /*  Print tag type and atom number  */
      fprintf (OutputFile, "HETATM%6i", NumWrite);

      /*  Print type name  */
      NameStr = a->TypeNames[a->type[ipart]];
      if (NameStr == NULL)
         fprintf (OutputFile, "%-3i", a->type[ipart]+1);
      else
         fprintf (OutputFile, "%-3s", NameStr);

      /*  Print stuff  */
#if 0
      fprintf (OutputFile, "  AUT     1    ");
#endif
      fprintf (OutputFile, "          1    ");

      /*  Print coordinates in angstroms  */
      fprintf
         (
         OutputFile,
#if 0
         "%8.3lf%8.3lf%8.3lf  0.00  0.00\n",
#endif
         "%8.3lf%8.3lf%8.3lf\n",
         CM * a->cur[NDIR*ipart+X],
         CM * a->cur[NDIR*ipart+Y],
         CM * a->cur[NDIR*ipart+Z]
         );

      }

   /*  Print end tag  */
   fprintf (OutputFile, "END\n");

   }


/*  Calculate average, minimum and maximum values of an array  */
void CalcAvgMinMax
(
double *Array,
int    NumValue,
int    NumDim,
double *Avg,
double *Min,
double *Max,
BOOLEAN UseSelect
)
	{
	int idim;
	int ival;
	int NumAvg;
	double  Value;
	BOOLEAN IsMinMaxInitialized = FALSE;


	/*  Empty array, return  */
	if (NumValue==0 || NumDim==0 || Array==NULL)
		return;
	
	ASSERT(Avg!=NULL)
	ASSERT(Min!=NULL)
	ASSERT(Max!=NULL)

	/*  Initialize output to zero  */
	LOOP (idim, NumDim)
		{
		Avg[idim] = 0.0;
		}


	/*  Loop through array  */
	NumAvg=0;
	LOOP (ival, NumValue)
	if   (!UseSelect || IS_SELECT(ival))
		{

		/*  Treat each dimension of current array group  */
		LOOP (idim, NumDim  )
			{
	
			/*  Get value of this array element */
			Value = Array[ival*NumDim + idim];

			/*  Add to average sum  */
			Avg[idim] += Value;

			/*  If Min,Max uninitialized, then intialize them to current value */
			if (!IsMinMaxInitialized)
				{
				Min[idim] = Value;
				Max[idim] = Value;
				}

			/*  Compare current value to Min,Max  */
			else
				{
				if (Value < Min[idim])
					{
					Min[idim] = Value;
					}
				if (Value > Max[idim])
					{
					Max[idim] = Value;
					}
				}
			}

		/*  Increment number of average sums  */
		NumAvg++;

		/*  Min and Max have been initialized  */
		IsMinMaxInitialized = TRUE;
		}

	/*  Done if no particles selected  */
	if (NumAvg==0)
		return;

	/*  Take average  */
	LOOP (idim, NumDim)
		{
		Avg[idim] /= NumAvg;
		}
	}


void CalcRdf
(
WORD ipart,
WORD *NeighList,
double *NeighRadSqr,
WORD NumNeigh
)
	{
	int ineigh;
	int jpart;
	int Index;
	BOOLEAN TypesOK;
	double  Radius;

	ASSERT (NeighRadSqr!=NULL)
	ASSERT (NeighList!=NULL)

	TypesOK = TRUE;

	LOOP (ineigh, NumNeigh)
		{

		/*
		Using only pairs of specific type
		*/
		if (UseTypes_m)
			{
			jpart = NeighList[ineigh];
			TypesOK =
				(Type_m[ipart]==Type1_m && Type_m[jpart]==Type2_m) ||
				(Type_m[jpart]==Type1_m && Type_m[ipart]==Type2_m);
		
			}

		if (TypesOK)
			{
			Radius = sqrt(NeighRadSqr[ineigh]);
			Index  = RdfTableFactor_m*(Radius-MinRadius_m);
			if (Index>=0 && Index<NumTable_m)
				RdfTable_m[Index]++;
			}
		}
	}



/*  
Convert symbol (t,m,x,y,z,vx,vy,vz,ep)  
to corresponding default print format
*/
char *get_pfmt (char *symbol) {

	if (!strcmp("i",symbol)) return "%d";
	if (!strcmp("g",symbol)) return "%d";
	if (!strcmp("t",symbol)) return "%d";
	if (!strcmp("m",symbol)) return "%10.4lf";
	if (!strcmp("x",symbol)) return "%10.4lf";
	if (!strcmp("y",symbol)) return "%10.4lf";
	if (!strcmp("z",symbol)) return "%10.4lf";
	if (!strcmp("vx",symbol)) return "%12.4e";
	if (!strcmp("vy",symbol)) return "%12.4e";
	if (!strcmp("vz",symbol)) return "%12.4e";
	if (!strcmp("ep",symbol)) return "%13.6e";
	/*  Default format  */
	return "%13.6e";
}


/*  
Convert symbol (t,m,x,y,z,vx,vy,vz,ep), index and format into
a pointer static formatted string.
*/
char *get_pptr (Particle_t *a, char *symbol, int index, char *fmt) {
	static char buffer[NSTR];

	/*  Index  */
	if (!strcmp("i",symbol)) return  sprintf(buffer,fmt,index+1),buffer;

	/*  Type  */
	if (!strcmp("t",symbol)) return  PFMT(a->type, a->type[index]+1);

	/*  Tag   */
	if (!strcmp("g",symbol)) return  PFMT(a->tag, a->tag[index]+1);

	/*  Mass  */
	if (!strcmp("m",symbol)) return  PFMT(a->mass, AVAG*a->mass[index]);

	/*  Position  */
	if (!strcmp("x",symbol)) return  PFMT(a->cur, 1e8*a->cur[NDIR*index+X]);
	if (!strcmp("y",symbol)) return  PFMT(a->cur, 1e8*a->cur[NDIR*index+Y]);
	if (!strcmp("z",symbol)) return  PFMT(a->cur, 1e8*a->cur[NDIR*index+Z]);

	/*  Velocity */
	if (!strcmp("vx",symbol)) return PFMT(a->v, a->v[NDIR*index+X]/s->dtime);
	if (!strcmp("vy",symbol)) return PFMT(a->v, a->v[NDIR*index+Y]/s->dtime);
	if (!strcmp("vz",symbol)) return PFMT(a->v, a->v[NDIR*index+Z]/s->dtime);

	/*  Eatom  */
	if (!strcmp("ep",symbol)) return PFMT(a->eatom, s->eunit*a->eatom[index]);

	/*  Ekin  */
	if (!strcmp("ek",symbol))  {
		double vx,vy,vz;
		if (a->v==NULL) return "0";
		if (a->mass==NULL) return "0";
		vx=a->v[NDIR*index+X];
		vy=a->v[NDIR*index+Y];
		vz=a->v[NDIR*index+Z];
		return PFMT(a->mass, s->eunit*0.5*a->mass[index]*(vx*vx+vy*vy+vz*vz)/(s->dtime*s->dtime));
		return buffer;
	}

	/*  Temperature  */
	if (!strcmp("tm",symbol)) {
		double vx,vy,vz,df;
		if (a->v==NULL || IS_FIX(index) || a->mass==NULL)  return "0";
		vx=a->v[NDIR*index+X];
		vy=a->v[NDIR*index+Y];
		vz=a->v[NDIR*index+Z];
		df = (a->cons && a->cons[index]) ? 
				((a->cons[index]->type==CONSTR_LINE) ? 1 : 2) : 3;
		return PFMT(a->mass, a->mass[index]*(vx*vx+vy*vy+vz*vz)/(s->dtime*s->dtime)/(df*BK));
		return buffer;
	}

	/*  Degree of freedom  */
	/*
	if (!strcmp("df",symbol))  {
		return PFMT(1, IS_FIX(index): 0 ? (( (a->cons && a->cons[index]) ?  ((a->cons[index]->type==0) ? 1 : 2) : 3));
	}
	*/
	if (!strcmp("df",symbol))  {
		return PFMT(1, 
			IS_FIX(index) ? 0 : 
			a->cons && a->cons[index] ? (a->cons[index]->type==0 ? 1 : 2) : 
			3
			);
	}

	/*  Time step  */
	if (!strcmp("ti",symbol))  {
		return PFMT(1, a->time);
	}

	return "0";
}