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.
*/
/*
************************************************************************
History
************************************************************************
*/
/*
4 Nov 1992 - in SEARCH2 - changed index[NPL] to *index = calloc (NPL, ...)
             This corrected (?) "Internal Stack Overflow" error
21 Jan 1997 - VERSION 2.3.3
    Corrected WrapParticle(), was not wrapping particles in the negative
    direction correctly.
26 Feb 2006 - Add GetMinMaxCoord()
13 Jun 2010 - In GetMinMaxCoord(), rename UseBox to DontUseBox.
*/

/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "iomngr.h"
#include "parse.h"
#include "strsub.h"
#include "sortsub.h"
#include "rcvio.h"
#include "plotc.h"
#include "cdalloc.h"
#include "cdsubs.h"
#include "cditemp.h"
#include "cdhouse.h"
#include "cdcor.h"
#include "cdboun.h"
#include "cdstep.h"
#include "cdvect.h"
#include "cdpairf.h"

/*
************************************************************************
Defines
************************************************************************
*/
#define NULL_MACRO_TEST
/*  Logical constants  */
#define BOOLEAN int
#define TRUE 1
#define FALSE 0

/*  Dimension constants  */
#define NDIR 3
#define X 0
#define Y 1
#define Z 2

#ifndef M_PI
#define M_PI        3.14159265358979323846264338327950288
#endif


#define NSTR 256


#if 0
#define MAX_NUM_FORMULA 16
#endif

#define MAX_BUFFER 1024

#define PUSH_FORWARD  0
#define PUSH_BACKWARD 1

/*
************************************************************************
Compile Switches
************************************************************************
*/


/*
************************************************************************
Global Variables
************************************************************************
*/
Particle_t   *a = NULL;
Particle_t   *b = NULL;
Simulation_t *s = NULL;

FILE         *out    = NULL;
LIST         *inlist = NULL;

/*
************************************************************************
Module-wide variables
************************************************************************
*/


/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void ZeroDerivatives (int, Particle_t *);
void RotateCoord (double[NDIR], double[NDIR][NDIR], double[NDIR]);
void RotateCoordinates
   (double [NDIR][NDIR], double[NDIR], BOOLEAN, double *, int);
void RotateVectorList
   (
   double [NDIR][NDIR], double[NDIR], BOOLEAN,
   ExternalVector_t **, int
   );
void ReadStepOutput (StepOutput_t *StepOutput, char *InputString);

#if 0
void StoreCalcCoord (double Coord[NDIR]);
#endif

void PushArray (double Value, double *Array, int NumArray, int Direction);

void POT_Free(void);


/*
************************************************************************
Functions
************************************************************************
*/

/*  Initialize Random Number Generator  */
void initrand (unsigned seed)
{
   srand(seed);
}

/*  Random Number Generator - returns between [-1..1]  */

double randfloat (void) {
   int   rand_half = RAND_MAX >> 1;
   float rfac      = 1.0/rand_half;
   return (rfac * (rand() - rand_half) );
}


void init()
   {

   /*  ALLOCATE SPACE FOR COODINATES  */
   a = palloc(4,4);
   b = NULL;
   s = NULL;
	ALLOCATE (s, Simulation_t, 1) ;

   /*  INITIALIZE SIMULTION PARAMETERS  */
   s->uratio = 1.1;
	s->nsearch = NEIGH_SEARCH_SORT;
   strcpy (s->eulabel, "K");
   s->eunit = 1.0/BK;
   s->nstep = 0;
   s->nrej  = 0;
   s->nacc  = 0;
   s->forcepnt = (void *) qu_calcforc;
   s->energylistpnt = (void *) qu_energy_list;
   InitStepOutput ( &(s->EnergyStepOutput) );
   InitStepOutput ( &(s->BoxStepOutput) );
   InitStepOutput ( &(s->StressStepOutput) );
   InitStepOutput ( &(s->TrajectoryStepOutput) );
   s->UseCompare  = FALSE;
   s->tempc  = -1;
   s->quench =  0;
   s->cstep  = 33;
   s->dtime  = 1.0;
   s->echo   = 1;
   s->UseSelectStress = FALSE;
	s->UseVerboseOutputFile=TRUE;
	s->PrintInfo=TRUE;
	s->FreePotPtr = (void *) POT_Free;
   }

        /*****************************************************/
        /**                                                 **/
        /**             CALCULATION SUBROUTINES             **/
        /**                                                 **/
        /**  There is one such subroutine for each command  **/
        /*****************************************************/



void RotateCoord
(
double Coord[NDIR],
double RotMatrix[NDIR][NDIR],
double Center[NDIR]
)
   {
   int idir;
   double TempCoord[NDIR];

   /*  Rotate coordinates, store in temporary space  */
   LOOP (idir,NDIR)
      {
      TempCoord[idir] =
         RotMatrix[idir][X] * (Coord[X] - Center[X]) +
         RotMatrix[idir][Y] * (Coord[Y] - Center[Y]) +
         RotMatrix[idir][Z] * (Coord[Z] - Center[Z]);
      }

   /*  Store rotated coordinate  */
   LOOP (idir, NDIR)
      Coord[idir] = TempCoord[idir] + Center[idir];

   }



void RotateCoordinates
(
double  RotMatrix[NDIR][NDIR],
double  Center[NDIR],
BOOLEAN UseSelect,
double  *InputCoord,
int     NumCoord
)
   {
   int ipart;
   double (*Coord)[NDIR] = (double (*)[NDIR]) InputCoord;


   /*  Return if null list  */
   if (InputCoord==NULL)
      return;

   /*  Loop over coordinates  */
   LOOP (ipart, NumCoord)
   if (!UseSelect || IS_SELECT(ipart) )
      {
      RotateCoord (Coord[ipart], RotMatrix, Center);
      }
   }



void RotateVectorList
(
double RotMatrix[NDIR][NDIR],
double Center[NDIR],
BOOLEAN UseSelect,
ExternalVector_t **List,
int NumPart
)
   {
   int ipart;


   /*  Return if null list  */
   if (List == NULL)
      return;

   /*  Loop over particles  */
   LOOP (ipart, NumPart)
   if (List[ipart] != NULL)
   if (!UseSelect || IS_SELECT(ipart) )
      {

      /*  Rotate origin  */
      RotateCoord (List[ipart]->Origin, RotMatrix, Center);

      /*  Rotate Vector  */
      RotateCoord (List[ipart]->Vector, RotMatrix, Center);

      }
   }


#pragma argsused
void energy_all (Particle_t *a, Simulation_t *s, int sflag)
    {
   /*  CHECK FOR VALID BOX  */
   if (    !a->surf[0] && a->bcur[0]==0
        || !a->surf[1] && a->bcur[1]==0
        || !a->surf[2] && a->bcur[2]==0
      )
      {
      printf ("ERROR (energy_all):  Cannot have repeating boundary without\n");
      printf ("                     specifying box size.\n");
      return;
      }
	/*  ENERGY  */
   ( (PotEnergy_f *) s->energylistpnt) (a, s, 0);
	}


/*
************************************************************************
Command Fuctions - for handling user commands
************************************************************************
*/

                       /*  ALLOCATE SPACE FOR ARRAYS  */

#pragma argsused
void read_allocate (char *instr)
   {
   printf ("WARNING:  Allocate command is now obsolete.\n");
   IncrementNumberOfWarnings();
   }

                            /*  READ BOX  */

void read_box      (char *instr)
   {
   char   *tokstr = NULL;
   BOOLEAN ScaleOption;
   double *box;
   double scale;
   int    idir;
   int    ipart;

   /*  TEST FOR SCALE OPTION  */
   tokstr = strhed(&instr);
   ScaleOption = !strcmpi(tokstr,"SCALE");

   /*  PARSE BOX SIZE  (x) */
   if (ScaleOption)
      a->bcur[X] = dblstrf (&instr);
   else
      a->bcur[X] = dblstrf (&tokstr);

   /*  PARSE BOX SIZE  (y & z) */
   a->bcur[Y] = dblstrf (&instr);
   a->bcur[Z] = dblstrf (&instr);

   /*  CONVERT UNITS FOR SCALE BOX  */
   LOOP (idir, NDIR)
      a->bcur[idir] *= ANGS;

   /*  Check for box zero  */
   if (a->bcur[X]==0.0 || a->bcur[Y]==0.0 || a->bcur[Z]==0.0)
      {
      printf ("ERROR:  One or more box dimensions is zero.\n");
      CleanAfterError();
      }

   if (ScaleOption)
   if (a->IsInitializedCoord)
      {
      box = a->bcur;
      LOOP (idir, NDIR)
         {
         scale = a->bcur[idir] / box[idir];
         LOOP (ipart, a->np)
            a->cur[3*ipart + idir] *= scale;
         }
      }

   else printf ("%s\n",
    "ERROR (read_box):  Cannot scale particles without initializing them.");

	/*  If first call to box, and no call to surface, set surface on  */
	if (a->IsInitializedBox==FALSE && a->IsInitializedSurface==FALSE)
		{
		a->surf[X] = a->surf[Y] = a->surf[Z] = FALSE;
		}

   /*  Set flag indicating box has been set  */
   a->IsInitializedBox = TRUE;

   /*  Set flag indicating current neighbor list is invalid  */
   a->InvalidNeighborList = TRUE;

   }


                           /*  READ BSAVE FLAG  */

void read_bsave  (char *instr)
   {
   ReadStepOutput ( &(s->BoxStepOutput), instr);
   }


                              /*  READ CALC   */

void read_calc  (char *instr)
   {
   int err;
   evalform (&instr, &err);
   if (err)
      printf ("ERROR: %s\n", parsemsg(err));
   }


                              /*  READ CLAMP  */


/*
Command syntax:

CLAMP { OFF | temp [cstep]  }

Set all particles to clamp value.
   (1) set all particles to tag 0.
   (2) set tag 0 properties

CLAMP SEL { OFF | temp cstep }
Set selected particles to clamp value.
   (1) Count number of particle remaining in each group

*/
void read_clamp (char *instr)
   {
   double  InputTemperature;
   double  InputCstep;
   int     ipart;
   int     itag;
   int     NumClamp;
   int     TagValue;
   int     ClampTagCount[NUM_CLAMP_TAG];
	BOOLEAN UseSelect = FALSE;
	BOOLEAN UseClamp  = FALSE;
	BOOLEAN ConditionFound = FALSE;
   BOOLEAN OpenTagFound   = FALSE;
   char   *tokstr = NULL;



   /*  Get first token  */
   tokstr = strhed (&instr);

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

   /*  Initialize tag count  */
   LOOP (itag, NUM_CLAMP_TAG)
      ClampTagCount[itag] = 0;

   /*  Count number of atoms in each tag group  */
	LOOP (ipart, a->np)
		{
		/*  Don't include selected atoms  */
		if (!UseSelect || !IS_SELECT(ipart))
			ClampTagCount[GET_CLAMP_TAG(ipart)]++;
		}

	/*  Test for "INFO" option  */
	if (!strcmpi(tokstr,"INFO"))
      {

		/*  Print info  */
		LOOP (itag, NUM_CLAMP_TAG)
			{
			/*  Do not print info if no atoms in group  */
			if (ClampTagCount[itag]==0)
				continue;

			printf ("CLAMP TAG %i\n",  itag+1);
			printf ("CLAMP NP  %i\n", ClampTagCount[itag]);
			if (!a->UseClamp[itag])
				{
				printf ("CLAMP OFF\n");
				}
			else
				{
				printf ("CLAMP TEMPERATURE %e\n", a->ClampTemp[itag]);
				printf ("CLAMP STEP        %e\n", a->ClampStep[itag]);
				}
			}
		return;
		}


   /*  Read clamp parameters  */
   if (!strcmpi(tokstr,"OFF"))
      {
      UseClamp = FALSE;
      }

   else
      {

      /*  Set clamp value  */
      UseClamp = TRUE;

      /*  Read temperature  */
      InputTemperature = dblstrf(&tokstr);

      /*
      To maintain compatability, negative temperature means clamp off
      */
      if (InputTemperature<0)
         UseClamp = FALSE;

      /*  Read cstep if present - ignore value of 0 or less  */
      InputCstep = dblstrf(&instr);

      /*  If no Cstep specified, use 33  */
      if (InputCstep==0)
         InputCstep = 33;
      }


   /*
   **************************
   Determine tag value to use
   **************************
   */

   /*  Test, start with flagged TagValue  */
   TagValue = -1;

   /*  Set up flags  */
   ConditionFound = FALSE;
   OpenTagFound   = FALSE;

   /*  First, if not using SEL, assign all atoms to tag 0  */
   if (!UseSelect)
		{
      TagValue = 0;
		OpenTagFound = TRUE;
		}

   /*  Second, if using SEL, is this clamp condition already in use?  */
   if (UseSelect)
		{
		for (itag=0; (itag<NUM_CLAMP_TAG) & (!ConditionFound); itag++)
         {
         if
         (
         UseClamp==a->UseClamp[itag] &&
            (
            UseClamp==FALSE ||
               (
               InputTemperature==a->ClampTemp[itag] &&
               InputCstep      ==a->ClampStep[itag]
               )
            )
         )
            {
            ConditionFound = TRUE;
            TagValue       = itag;
            }
         }
      }


   /*  If clamp condition not found, test for open tag value  */
   if (UseSelect && !ConditionFound)
      {
      for (itag=0; itag<NUM_CLAMP_TAG && !OpenTagFound; itag++)
         {
         OpenTagFound = (ClampTagCount[itag]==0);
         if (OpenTagFound)
            {
            TagValue = itag;
            }
         }
      }


   /*  If condition not yet in use, and no open tag found, error  */
   if (UseSelect && !ConditionFound && !OpenTagFound)
      {
      printf ("ERROR: ");
      printf ("Cannot maintain more than %i separate clamp conditions\n",
         NUM_CLAMP_TAG);
      CleanAfterError();
      }



   ASSERT(TagValue!=-1)
   ASSERT(TagValue>=0&&TagValue<NUM_CLAMP_TAG)

   /*  Set condition for this tag group  */
   if (OpenTagFound)
      {
      a->UseClamp [TagValue] = UseClamp;
      a->ClampTemp[TagValue] = InputTemperature;
      a->ClampStep[TagValue] = InputCstep;
      }

   /*  Place selected atom in tag group  */
   NumClamp = 0;
   LOOP (ipart, a->np)
      {

      /*  This particle is selected, set to itag  */
      if (!UseSelect || IS_SELECT(ipart))
         {
			APPLY_CLAMP_TAG(ipart,TagValue)
         }

      /*  This particle is in current clamp set, count it  */
      if (IS_CLAMP_TAG(ipart,TagValue))
         {
         NumClamp++;
         }
      }


   }


/*  Read COR file  */
void read_cor (char *instr)
   {
   long step;
   FILE *fin;
   char *tokstr = strhed (&instr);

   /*  OPEN FILE  */
   fin = fopen (tokstr, "rb");
   if (fin==NULL)
      {
      printf ("ERROR :  Cannot open input file %s.\n", tokstr);
      CleanAfterError();
      }

   /*  READ JUST FIRST STEP  */
   if (*instr==0)
      ReadCorFile (fin, a);
   else
      {
      step = intstrf(&instr);
      do
         {
         ReadCorFile (fin, a);
         }
         while (!feof(fin) && a->step!=step);
      }
   fclose (fin);

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }


                              /* READ CMD  */

void read_cmd (char *instr)
   {
   int nstep = intstrf (&instr);

	/*     Check for needed initializations  */
	if (!a->IsInitializedBox)
		{
		printf
			("ERROR:  Cannot peform CMD command without first setting BOX.\n");
		CleanAfterError();
		}

	if (!a->IsInitializedCoord)
		{
      printf ("ERROR:  Cannot peform CMD command without");
      printf (" first setting particle coordinates.\n");
      CleanAfterError();
      }

	if (!a->IsInitializedMass)
		{
		printf ("ERROR:  Cannot peform CMD command without");
		printf (" first setting all atoms' masses.\n");
		CleanAfterError();
		}

   if (!a->IsInitializedPotential)
      {
      printf ("WARNING:  Potential was not set");
      printf (" before performing molecular dynamics\n");
		IncrementNumberOfWarnings();
		}

	/*  Perform molecular dynamics  */
	runcmd (a, s, nstep);

	/*  Write out ending step  */
	if (s->PrintInfo)
		printf ("***   Current step is %li\n", a->step);

   }

                              /*  READ DTIME  */

void read_damp  (char *instr)
   {
   int ipart;
   BOOLEAN IsDampZero;
   char  *HeadStr    = NULL;
   double Coord[NDIR];

   HeadStr = strhed (&instr);

   if (!strcmpi(HeadStr,"ON"))
      a->UseDamp = TRUE;
   else if (!strcmpi(HeadStr, "OFF"))
      a->UseDamp = FALSE;
   else
      {
      printf ("WARNING:  Unrecognized option to DAMP command.\n");
      IncrementNumberOfWarnings();
      return;
      }

   /*  If damping is turned off, then return now  */
   if (a->UseDamp==FALSE)
      return;

   /*
   If damping requested,
   but no current damping or previous damping values,
   set damping off
   */
   HeadStr = strhed (&instr);
   if (a->Damp==NULL && HeadStr[0]=='\0')
      {
      a->UseDamp = FALSE;
      return;
      }

   /*
   If damping selected, and value specified, but no storage allocated,
   allocate now
   */
   if (a->Damp==NULL)
      {
      ALLOCATE(a->Damp, double, a->NumPartAlloc)
      }


   /*  Loop through particles  */
	CheckForNoneSelected();
   IsDampZero = FALSE;
   LOOP (ipart, a->np)
      {
      if (IS_SELECT(ipart) )
         {

         /*  Set values for X, Y, Z in formula  */
         Coord[X] = 1e8*a->cur[NDIR*ipart+X];
         Coord[Y] = 1e8*a->cur[NDIR*ipart+Y];
         Coord[Z] = 1e8*a->cur[NDIR*ipart+Z];
         EvaluateFormula
            (HeadStr, Coord, 1e+8, &a->Damp[ipart], 1.0, 1);
         }

      /*  Check if all damping coefficients are zero  */
      IsDampZero = ( IsDampZero || (a->Damp[ipart]==0.0) );
      }

   /*  All damping coeffiecients are zero  */
   if (IsDampZero)
      {

      /*  Set damp flag to false  */
      a->UseDamp = FALSE;

      /*  Free damp array  */
      FREE (a->Damp)
      }

   }
                              /*  READ DTIME  */

void read_dtime (char *instr)
   {
   double NewTimeStep;
   double OldTimeStep;
   double TimeScale;
   double CumulativeTimeScale;
   int    i;
   int    j;

   /*  No tim step entered - leave tim step unchanged  */
   if (*instr==0)
      return;

   /*  Read new time step  */
   NewTimeStep = dblstrf (&instr);

   /*  Check time step  */
   if (NewTimeStep<=0)
      {
      printf ("ERROR: DTIME (%f) <= 0\n", NewTimeStep);
      return;
      }

   /*  Save old time step  */
   OldTimeStep    = s->dtime;

   /*  Store new time step  */
   s->dtime = NewTimeStep;

   /*  RESCALE TIME DERIVATIVES  */
   if (OldTimeStep != NewTimeStep)
      {

      /*  Scaling Factor  */
      TimeScale  = NewTimeStep/OldTimeStep;

      /*  Rescale velocities  */
      CumulativeTimeScale = TimeScale;
      if (a->v != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->v[i + 3*j] *= CumulativeTimeScale;

      /*  Rescale accelerations  */
      CumulativeTimeScale *= TimeScale;
      if (a->c2 != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->c2[i + 3*j] *= CumulativeTimeScale;

      /*  Rescale 3rd time derivatives  */
      CumulativeTimeScale *= TimeScale;
      if (a->c3 != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->c3[i + 3*j] *= CumulativeTimeScale;

      /*  Rescale 4th time derivatives  */
      CumulativeTimeScale *= TimeScale;
      if (a->c4 != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->c4[i + 3*j] *= CumulativeTimeScale;

      /*  Rescale 5th time derivatives  */
      CumulativeTimeScale *= TimeScale;
      if (a->c5 != NULL)
         for (i=0;i<3;i++)
            for (j=0;j<a->np;j++)
               a->c5[i + 3*j] *= CumulativeTimeScale;
      }
   }

                              /*  READ DEBUG  */

void read_debug (char *instr)
   {
   char *HeadStr;

   HeadStr = strhed (&instr);

   if      (!strcmpi(HeadStr, "ON" ))
      s->debug = TRUE;
   else if (!strcmpi(HeadStr, "OFF"))
      s->debug = FALSE;

   else
      printf ("ERROR:  Unrecognized option to DEBUG command.\n");
   Debug_g = s->debug;
   }

                            /*  READ DUP  */

void read_dup      (char *instr)
   {
   int OldNumPart;
   int NewNumPart;
   int idir;
   int idup;
   int iold;
   int inew;
   int indexnew;
   int indexold;
   int ndup;
   int NumWrap;
   double dx;
   double dy;
   double dz;
   double tx;
   double ty;
   double tz;
   double (*Coord)[NDIR] = NULL;

   ndup = intstrf (&instr);
   dx = ANGS*dblstrf (&instr);
   dy = ANGS*dblstrf (&instr);
   dz = ANGS*dblstrf (&instr);
   tx = ty = tz = 0;

   /*  RETURN IF NO DUPLICATIONS SPECIFIED  */
   if (ndup==0)
      return;

   OldNumPart = a->np;
   NewNumPart = a->np + ndup * a->nsel;

   /*  REALLOCATE Particles  */
   ReallocateParticle (a, NewNumPart);

   /*  Set pointer to coordinates  */
   Coord = (double (*)[NDIR]) a->cur;

   /*  DUPLICATE PARTICLES  */
   inew = OldNumPart;
   LOOP(idup,ndup)
      {
      tx += dx;
      ty += dy;
      tz += dz;
      LOOP(iold,OldNumPart)
      if (IS_SELECT(iold) )
         {
         Coord[inew][X] = Coord[iold][X] + tx;
         Coord[inew][Y] = Coord[iold][Y] + ty;
         Coord[inew][Z] = Coord[iold][Z] + tz;
         a->type[inew] = a->type [iold];
         if (a->mass!=NULL) a->mass[inew] = a->mass[iold];

         /*  Copy velocity and higher derivatives  */
         LOOP (idir, NDIR)
            {
            indexnew = inew * NDIR + idir;
            indexold = iold * NDIR + idir;
            if (a->v  != NULL)  a->v [indexnew] = a->v [indexold];
            if (a->c2 != NULL)  a->c2[indexnew] = a->c2[indexold];
            if (a->c3 != NULL)  a->c3[indexnew] = a->c3[indexold];
            if (a->c4 != NULL)  a->c4[indexnew] = a->c4[indexold];
            if (a->c5 != NULL)  a->c5[indexnew] = a->c5[indexold];
            }
         APPLY_SELECT(inew)
         a->nsel++;
         inew ++;
         }
      }

   /*  Set new number of particles  */
   ASSERT(inew==NewNumPart)
   a->np = NewNumPart;

   /*  Wrap particles   */
   NumWrap = WrapParticles(a);

   if (NumWrap>0)
      PrintWrapWarning (NumWrap);

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }
                           /*  READ ENERGY UNITS  */

void read_eunit  (char *instr)
   {
   double tunit;
   char *tokstr = strhed (&instr);

   /*  NO UNITS SPECIFIED - PRINT OUT CURRENT UNIT  */
   if (!*tokstr)
      {
      printf ("   EUNIT %s %e\n", s->eulabel, s->eunit);
      return;
      }

   /* KELVIN */
   if (!strcmpi(tokstr,"K"))
      {
      strcpy (s->eulabel, "K");
      s->eunit = 1/BK;
      return;
      }

   /* ERGS   */
   if (!strcmpi(tokstr,"ERG"))
      {
      strcpy (s->eulabel, "erg");
      s->eunit = 1;
      return;
      }

   /* JOULES  */
   if (!strcmpi(tokstr,"JOULE"))
      {
      strcpy (s->eulabel, "joule");
      s->eunit = JOULE;
      return;
      }

   /* ELECTRON VOLTS  */
   if (!strcmpi(tokstr,"EV"))
      {
      strcpy (s->eulabel, "eV");
      s->eunit = EV;
      return;
      }

   /* ELECTRON VOLTS  */
   if (!strcmpi(tokstr,"HARTREE"))
      {
      strcpy (s->eulabel, "Hartree");
      s->eunit = HARTREE;
      return;
      }

   /*  USER SPECIFIED UNITS  */
   tunit = dblstrf (&instr);
   if (tunit==0)
      printf ("ERROR (read_eunit): energy unit 0 or unreadable.");
   else
      {
      strcpy (s->eulabel, tokstr);
      s->eunit   = tunit;
      }
   }

                          /*  READ EXTERNAL FORCE  */
void read_extforce (char *InputString)
   {
   double InputForce[NDIR];
   double (*Coord)[NDIR];
   int ipart;
   int idir;
   int NumActiveForces;

   if (!strcmpi(InputString,"CLEAR"))
      {
      /*  Free force constant and origin storage for individual particles  */
      if (a->ForcePtrList != NULL)
         LOOP (ipart, a->np)
            FREE (a->ForcePtrList[ipart])

      /*  Free storage pointer storage of entire list  */
      FREE (a->ForcePtrList);

      return;
      }

#if 0
   /*  Read input force components  */
   LOOP (idir, NDIR)
      InputForce[idir] = dblstrf (&InputString) / s->eunit;

   /*  Set flag if no force  */
   ZeroForce =
      (
      InputForce[X]==0.0 &&
      InputForce[Y]==0.0 &&
      InputForce[Z]==0.0
      );

   /*  If zero force and currently no list, than don't create one  */
   if (ZeroForce && a->ForcePtrList==NULL)
      return;
#endif

   /*  ALLOCATE ARRAY  */
   if (a->ForcePtrList==NULL)
      ALLOCATE (a->ForcePtrList, ExternalVector_t *, a->NumPartAlloc)

   /*
   For each selected particle, allocate and store for force
   (NOTE:  If all remaining force constants are set to zero,
   then list will be removed from memory below)
   */
	CheckForNoneSelected();
   Coord = (double (*)[NDIR]) a->cur;
   NumActiveForces = 0;
   LOOP (ipart, a->np)
      {
      if (IS_SELECT(ipart) )
         {

         EvaluateFormula
            (
            InputString,
            Coord[ipart],
            1e8,
            InputForce,
            1.0,
            NDIR
            );

         /*  Free Storage if applied force constant is zero  */
         if (InputForce[X]==0.0 && InputForce[Y]==0.0 && InputForce[Z]==0.0)
            {
            if (a->ForcePtrList[ipart] != NULL)
               FREE (a->ForcePtrList[ipart]);
            }

         else
            {
            /*  Allocate new storage if needed  */
            if (a->ForcePtrList[ipart] == NULL)
					ALLOCATE (a->ForcePtrList[ipart], ExternalVector_t, 1);

            /*  Copy origin and force constant  */
            LOOP (idir, NDIR)
               {
               a->ForcePtrList[ipart]->Origin   [idir] =
                  Coord[ipart][idir];
               a->ForcePtrList[ipart]->Vector  [idir] =
                  InputForce[idir];
               }
            }
         }

      /*  Count number of active forces  */
      if (a->ForcePtrList[ipart] != NULL)
         NumActiveForces++;
      }

   /*  If no active forces remain, free storage for list  */
   if (NumActiveForces==0)
      {
      FREE (a->ForcePtrList)
      }

   }


                          /*  READ EXTERNAL SPRING  */
void read_extspring (char *InputString)
   {
   double  (*Coord)[NDIR];
   double InputSpringConstant[NDIR];
   int NumActiveSprings;
   int ipart;
   int idir;

   /*  Remove storage for external spring  */
   if (!strcmpi(InputString,"CLEAR"))
      {

      /*  Free spring constant and origin storage for individual particles  */
      if (a->SpringPtrList != NULL)
         LOOP (ipart, a->np)
            FREE (a->SpringPtrList[ipart])

      /*  Free storage pointer storage of entire list  */
      FREE (a->SpringPtrList);

      /*  Return to calling function  */
      return;
      }

#if 0
   /*  Read input spring constant  */
   LOOP (idir, NDIR)
      InputSpringConstant[idir] = dblstrf (&InputString) / s->eunit;

   /*  Set flag if no spring  */
   ZeroSpring =
      (
      InputSpringConstant[X]==0.0 &&
      InputSpringConstant[Y]==0.0 &&
      InputSpringConstant[Z]==0.0
      );

   /*  If zero input spring and currently no list, than don't create one  */
   if (ZeroSpring && a->SpringPtrList==NULL)
      return;
#endif


   /*  ALLOCATE ARRAY  */
   if (a->SpringPtrList==NULL)
      ALLOCATE (a->SpringPtrList, ExternalVector_t *, a->NumPartAlloc)

   /*
   For each selected particle, allocate and store for spring
   (NOTE:  If all remaining spring constants are set to zero,
   then list will be removed from memory below)
   */
	CheckForNoneSelected();
   Coord = (double (*)[NDIR]) a->cur;
   NumActiveSprings = 0;
   LOOP (ipart, a->np)
      {
      if (IS_SELECT(ipart) )
         {

         EvaluateFormula
            (
            InputString,
            Coord[ipart],
            1e8,
            InputSpringConstant,
            1.0,
            NDIR
            );

         /*  Free Storage if applied spring constant is zero  */
         if (
            InputSpringConstant[X]==0.0 &&
            InputSpringConstant[Y]==0.0 &&
            InputSpringConstant[Z]==0.0
            )
            {
            if (a->SpringPtrList[ipart] != NULL)
               FREE (a->SpringPtrList[ipart]);
            }

         else
            {
            /*  Allocate new storage if needed  */
            if (a->SpringPtrList[ipart] == NULL)
					ALLOCATE (a->SpringPtrList[ipart], ExternalVector_t, 1);

            /*  Copy origin and spring constant  */
            LOOP (idir, NDIR)
               {
               a->SpringPtrList[ipart]->Origin  [idir] =
                  Coord[ipart][idir];
               a->SpringPtrList[ipart]->Vector  [idir] =
                  InputSpringConstant[idir];
               }

            }

         }

      /*  Count number of active springs  */
      if (a->SpringPtrList[ipart] != NULL)
         NumActiveSprings++;
      }

   /*  If no active springs remain, free storage for list  */
   if (NumActiveSprings==0)
      {
      FREE (a->SpringPtrList)
      }

   }

                           /*  READ FIX FLAG    */
void read_fix (char *instr)
   {
   int i;
   char *tokstr = NULL;

   tokstr = strhed (&instr);

   if (!strcmpi(tokstr,"ON"))
      {

      /*  
		Make sure that fix is allocated  
		(Am pretty sure this ALLOCATION unnecessary, because
      Selection is allocated during program inialization May 6, 1998)
		*/
      if (a->Selection==NULL)
         ALLOCATE(a->Selection,SEL_T,a->NumPartAlloc)

      /*  Set selected particles as fixed  */
		CheckForNoneSelected();
      LOOP(i,a->np)
         {

         /*  If selected ..  */
         if (IS_SELECT(i) )
            {

            /*  Set fix flag  */
            APPLY_FIX(i)

            /*  Zero derivatives  */
            ZeroDerivatives (i, a);
            }
         }
      }
   else if (!strcmpi(tokstr,"OFF"))
      {
		CheckForNoneSelected();
      if (a->Selection!=NULL)
         {
         LOOP (i, a->np)
            {
            if (IS_SELECT(i) )
               REMOVE_FIX(i)
            }
         }
      }
   else
      {
      printf ("WARNING:  Unrecognized  option to FIX command.\n");
      IncrementNumberOfWarnings();
      }

   /*  Count number of fixed particles */
   a->nfix = 0;
   if (a->Selection!=NULL)
      {
      LOOP (i, a->np)
         {
         if (IS_FIX(i))
            a->nfix++;
         }
      }

	/*  Update number of degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }
                           /*  READ ECHO FLAG   */
void read_echo (char *instr)
   {
   char *tokstr = NULL;

   tokstr = strhed (&instr);

   if      (!strcmpi(tokstr,"ON" )) s->echo = TRUE;
   else if (!strcmpi(tokstr,"OFF")) s->echo = FALSE;
   }


/*  Erase file  */
void read_erase (char *instr)
   {
   int result;

   result = remove(instr);
   if (result!=0)
      {
      printf ("WARNING:  Could not erase file %s\n", instr);
      IncrementNumberOfWarnings();
      }

   }

                           /*  READ ESAVE FLAG  */

void read_esave  (char *instr)
   {
   ReadStepOutput ( &(s->EnergyStepOutput), instr);
   }

                               /*  READ ITEMP  */

/*
NOTE:
To intialize Temp for a subset of particle,
use FIX command prior to calling ITEMP.
*/
void read_itemp (char *instr)
   {
   double Temp;
   char *tokstr = NULL;
   BOOLEAN  xflag;
   BOOLEAN  yflag;
   BOOLEAN  zflag;
   BOOLEAN  OptionRead;
   BOOLEAN  UseSelect      = FALSE;
   BOOLEAN  UseOldVelocity = FALSE;
   int  ndim;


   /*  Check for SEL and OLD options  */
   do
      {
      OptionRead = FALSE;
      tokstr = strhed (&instr);
      if (!strcmpi(tokstr, "SEL"))
         {
         UseSelect = TRUE;
         OptionRead = TRUE;
			CheckForNoneSelected();
         }
      if (!strcmpi(tokstr, "OLD"))
         {
         UseOldVelocity = TRUE;
         OptionRead = TRUE;
         }
      }
      while (OptionRead);

   /*  Read input temperature  */
   Temp = dblstrf (&tokstr);

   /*  Find Dimensions to Initialize  */
   if (*instr==0)
      {
      xflag = TRUE;
      yflag = TRUE;
      zflag = TRUE;
      ndim  = 3;
      }

   else
      {
      xflag = FALSE;
      yflag = FALSE;
      zflag = FALSE;
      tokstr = strhed (&instr);
      while (*tokstr!=0)
         {
         if      (!strcmpi(tokstr, "X"))   xflag = TRUE;
         else if (!strcmpi(tokstr, "Y"))   yflag = TRUE;
         else if (!strcmpi(tokstr, "Z"))   zflag = TRUE;
         else
            {
            printf ("ERROR IN TEMPERATURE INITIALIZATION.\n");
            printf ("    Invalid dimension specified \"%s\".\n", tokstr);
            }
         tokstr = strhed (&instr);
         }
      ndim = 0;
      if (xflag) ndim++;
      if (yflag) ndim++;
      if (zflag) ndim++;
      if (ndim==0)
         {
         printf ("ERROR (read_itemp) IN TEMPERATURE INITIALIZATION.\n");
         printf ("    No valid dimension specified.\n");
         CleanAfterError();
         }
      }

   /*  Call routine to initialize velocities  */
   InitialTemp (a, s, Temp, UseSelect, UseOldVelocity, xflag,yflag,zflag);
   }


                               /*  READ LABEL  */
void read_label (char *instr)
{
   int i,ntitle;
   char bufstr[NTITLESTR];

   ntitle = intstrf(&instr);
   i = 0;
   while (   i<ntitle
          && NULL!=m_gets_list(bufstr,NTITLESTR,inlist))
   {
      if (i<8)
         strcpy (a->title[i], bufstr);
      i++;
   }
}


                           /*  READ MASS  */

void read_mass (char *instr)
   {
   int ipart;
   double Mass;

   Mass = dblstrf (&instr) / AVAG;
   a->IsInitializedMass = TRUE;

	CheckForNoneSelected();
   LOOP (ipart, a->np)
      {

      /*  Assign mass to selected particles  */
      if (IS_SELECT(ipart) )
         a->mass[ipart] = Mass;

      /*  Test if all masses initialized  */
      a->IsInitializedMass =
         (a->IsInitializedMass && a->mass[ipart]!=0.0);
      }

   }

                           /*  READ MC  */

void read_mc (char *instr) {
   int nstep;
   double temp;
   nstep = intstrf (&instr);
   temp  = dblstrf (&instr);
   runmc (a, s, nstep, temp);
}

                           /*  READ MOVE  */

void read_move     (char *instr)
   {
   int    ipart;
   double (*Coord)[NDIR] = (double (*)[NDIR]) a->cur;
   double Values[NDIR];


	CheckForNoneSelected();
   LOOP (ipart, a->np)
   if (IS_SELECT(ipart) )
      {

      /*  Map formula string into Values[] array  */
      EvaluateFormula (instr, Coord[ipart], 1e+8, Values, 1e-8, NDIR);

/*test*/
#if 0
{
char  buffer[160];
char* tokstr=buffer;

sprintf(tokstr, "X=%le Y=%le Z=%le  X Y Z", 0.1, 0.1, 0.1);
dblstrf(&tokstr);
dblstrf(&tokstr);
dblstrf(&tokstr);
Values[X] = dblstrf(&tokstr)*1e-8;
Values[Y] = dblstrf(&tokstr)*1e-8;
Values[Z] = dblstrf(&tokstr)*1e-8;
}
#endif
/*end test*/

      /*  Increment position by Values[]  */
      Coord[ipart][X] += Values[X];
      Coord[ipart][Y] += Values[Y];
      Coord[ipart][Z] += Values[Z];

      }

   /*  Wrap displaced particles throught repeating boundary conditions  */
   WrapParticles (a);
   }


                           /*  READ NRANGE */

void read_nrange   (char *instr)
   {
   double SearchCutoffRatio;

   /*  Read ratio of neighbor search cutoff to potential cutoff  */
   SearchCutoffRatio = dblstrf (&instr);


   /*  Check value of ratio  */
   if (SearchCutoffRatio < 1)
      {
      printf ("ERROR (read_nrange):  NRANGE must be 1 or greater.\n");
      CleanAfterError();
      }

   if (SearchCutoffRatio > 2)
      {
      printf ("WARNING:  Are you sure you want NRANGE greater than 2?\n");
      printf ("  This may result in a large and expensive neighbor list.\n");
      IncrementNumberOfWarnings();
      }

   /*  Save value  */
   s->uratio  = SearchCutoffRatio;
   }


                           /*  READ NSEARCH */

void read_nsearch  (char *instr) {
	if (!strcmpi(instr,"SORT")) {
		s->nsearch = NEIGH_SEARCH_SORT;
	} else if (!strcmpi(instr,"CELL")) {
		s->nsearch = NEIGH_SEARCH_CELL;
	} else if (!strcmpi(instr,"CELL2")) {
		s->nsearch = NEIGH_SEARCH_CELL2;
	} else {
		printf ("ERROR:  Unknown value for NSEARCH.\n");
		CleanAfterError();
	}
}


                         /*  READ PARTICLE  */

void read_particle (char *instr)
   {
   char *tokstr=NULL;
   char bufstr[NSTR];
   int  LastNewPart;
   int  LastOldPart;
   int  FirstNewPart;
   int  NumRead;
   int  ipart;
   int  Type;
   int  NumWrap;
   BOOLEAN addflag;
   int    i;
   double x,y,z;
   double (*Coord)[NDIR];

   addflag = FALSE;

   /*  Parse first token  */
   tokstr = strhed (&instr);

   /*  ADD option  */
   if (!strcmpi(tokstr,"ADD"))
      {
      addflag = TRUE;
      tokstr = strhed (&instr);
      }

   /*  No particle number entered  */
   if (tokstr==NULL)
      {
      printf ("ERROR (read_particle):  No number specified for PARTICLE.");
      CleanAfterError();
      }

   /*  Number of particles to input  */
   NumRead = intstrf (&tokstr);

   /*  Invalid number  */
   if (NumRead==0)
      {
      printf ("ERROR (read_particle): ");
      printf ("Number for PARTICLE is 0 or unrecognizable.\n");
      CleanAfterError();
      }

   /*  Save old number of particles  */
   LastOldPart = a->np;

   /*  Clear selected particles  */
   LOOP (i, a->np)
		REMOVE_SELECT(i)

   /*  READ COORDINATES  */
   if (addflag)
      {
      FirstNewPart = a->np;
      }
   else
      {
      FirstNewPart = 0;
      }

	/*  Last particle after new particles are read  */
   LastNewPart  = FirstNewPart + NumRead;

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

   /*  Set pointer to particles  */
   Coord = (double (*)[NDIR]) a->cur;

	
   /*  Save only new particles  */
   ipart = FirstNewPart;
   while (   ipart < LastNewPart
          && NULL!=m_gets_list(bufstr,NSTR,inlist))
      {
      tokstr = bufstr;
      Type = intstrf (&tokstr)-1;
      x = dblstrf (&tokstr);
      y = dblstrf (&tokstr);
      z = dblstrf (&tokstr);
		CLEAR_ALL_FLAGS(ipart)
      APPLY_SELECT(ipart)
		if (a->tag != NULL)
			a->tag [ipart] = (BYTE) 0;
		a->type[ipart] = (BYTE) Type;
      Coord[ipart][X] = x * ANGS;
      Coord[ipart][Y] = y * ANGS;
      Coord[ipart][Z] = z * ANGS;
      ipart++;
      }

   if (ipart != LastNewPart)
      {
      printf ("ERROR (read_particle):  Could not find all particles.");
      printf ("Number of particles expected is %d.\n", LastNewPart);
      printf ("Number of particles read     is %d.\n",
               ipart-LastOldPart);
      CleanAfterError();
      }

   /*  Set number of particles  */
   a->np   = LastNewPart;
   a->nsel = NumRead;

   /*  Set Coordinate flag  */
   a->IsInitializedCoord = TRUE;


   /*  Wrap particles   */
   NumWrap = WrapParticles(a);


   if (NumWrap>0)
      PrintWrapWarning (NumWrap);

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }


                         /*  READ POSVEL */

void read_posvel (char *instr)
   {
   char *tokstr=NULL;
   char bufstr[NSTR];
   int  LastNewPart;
   int  LastOldPart;
   int  FirstNewPart;
   int  NumRead;
   int  ipart;
   int  Type;
   int  NumWrap;
   BOOLEAN addflag;
   int    i;
   double x,y,z;
   double vx,vy,vz;
   double (*Coord)[NDIR];
   double (*Vel  )[NDIR];

   addflag = FALSE;

   /*  Parse first token  */
   tokstr = strhed (&instr);

   /*  ADD option  */
   if (!strcmpi(tokstr,"ADD"))
      {
      addflag = TRUE;
      tokstr = strhed (&instr);
      }

   /*  No particle number entered  */
   if (tokstr==NULL)
      {
      printf ("ERROR (read_posvel):  No number specified for POSVEL.");
      CleanAfterError();
      }

   /*  Number of particles to input  */
   NumRead = intstrf (&tokstr);

   /*  Invalid number  */
   if (NumRead==0)
      {
      printf ("ERROR (read_posvel): ");
      printf ("Number for POSVEL is 0 or unrecognizable.\n");
      CleanAfterError();
      }

   /*  Save old number of particles  */
   LastOldPart = a->np;

   /*  Clear previous particle selection  */
   LOOP (i, a->np)
		REMOVE_SELECT(i)

   /*  READ COORDINATES  */
   if (addflag)
      {
      FirstNewPart = a->np;
      }
   else
      {
      FirstNewPart = 0;
      }

	/*  Last particle after new particles are read  */
   LastNewPart  = FirstNewPart + NumRead;

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

	/*  Allocate velocities if not already present  */
	if (a->v==NULL)
		ALLOCATE (a->v, double, NDIR*a->NumPartAlloc)

   /*  Set pointer to particles  */
   Coord = (double (*)[NDIR]) a->cur;
   Vel   = (double (*)[NDIR]) a->v;

	
   /*  Save only new particles  */
   ipart = FirstNewPart;
   while (   ipart < LastNewPart
          && NULL!=m_gets_list(bufstr,NSTR,inlist))
      {
      tokstr = bufstr;
      Type = intstrf (&tokstr)-1;
      x  = dblstrf (&tokstr);
      y  = dblstrf (&tokstr);
      z  = dblstrf (&tokstr);
      vx = dblstrf (&tokstr);
      vy = dblstrf (&tokstr);
      vz = dblstrf (&tokstr);
		CLEAR_ALL_FLAGS(ipart)
      APPLY_SELECT(ipart)
		a->type[ipart] = Type;
      Coord[ipart][X] = x  * ANGS;
      Coord[ipart][Y] = y  * ANGS;
      Coord[ipart][Z] = z  * ANGS;
      Vel  [ipart][X] = vx * s->dtime;
      Vel  [ipart][Y] = vy * s->dtime;
      Vel  [ipart][Z] = vz * s->dtime;
      ipart++;
      }

   if (ipart != LastNewPart)
      {
      printf ("ERROR (read_particle):  Could not find all particles.");
      printf ("Number of particles expected is %d.\n", LastNewPart);
      printf ("Number of particles read     is %d.\n",
               ipart-LastOldPart);
      CleanAfterError();
      }

   /*  Set number of particles  */
   a->np   = LastNewPart;
   a->nsel = NumRead;

   /*  Set Coordinate flag  */
   a->IsInitializedCoord = TRUE;

   /*  Wrap particles   */
   NumWrap = WrapParticles(a);

	/*  Warn user if input particles needed wrapping  */
   if (NumWrap>0)
      PrintWrapWarning (NumWrap);

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }



                           /*  READ POTENTIAL  */

/*  potential { eam | quartic | tersoff }          */

void read_potential (char *instr) {
   char *tokstr = NULL;

   /*  SET POTENTIAL TYPES   */
	if (IsFirstToken("SET",instr)) {

		/*  Record that potential was initialized  */
		a->IsInitializedPotential = TRUE;

		/*  Remove leading "SET" token  */
		strhed (&instr);

		/*  Get Next Token  */
		tokstr = strhed (&instr);

		/*  QUARTIC option  */
		if (!strcmpi(tokstr,"QUARTIC")) {
			strcpy (s->pottype, tokstr);
			s->energylistpnt = (void *) qu_energy_list;
			s->forcepnt      = (void *) qu_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) POT_Free;
			return;
		}

		/*  TERSOFF option   */
		if (!strcmpi(tokstr,"TERSOFF")) {
			strcpy (s->pottype, tokstr);
			s->energylistpnt = (void *) csi_energy_list;
			s->forcepnt      = (void *) csi_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) POT_Free;
			csi_read_potential ("INIT", "");
			return;
		}

		/*  Stillinger Weber  option  */
		if (!strcmpi(tokstr,"STILL")) {
			strcpy (s->pottype, tokstr);
			s->energylistpnt = (void *) stil_energy_list;
			s->forcepnt      = (void *) stil_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) POT_Free;
			STILL_Init (instr);
			return;
		}

		/*  EAM option  */
		if (!strcmpi(tokstr,"EAM")) {
			strcpy (s->pottype, tokstr);
			s->energylistpnt = (void *) em_energy_list;
			s->forcepnt      = (void *) em_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) POT_Free;
			em_read_potential ("INIT", instr);
			return;
		}

		/*  PAIR option  */
		if (!strcmpi(tokstr,"PAIR")) {
			strcpy (s->pottype, tokstr);
			s->energylistpnt = (void *) pr_energy_list;
			s->forcepnt      = (void *) pr_calcforc;
			( (PotFree_f *) s->FreePotPtr) ();
			s->FreePotPtr = (void *) PAIR_Free;
			PAIR_Init(instr);
			return;
		}

		/*  Unrecognized option   */
		printf ("ERROR:  Unrecognized Potential Type (%s)\n", tokstr);
		CleanAfterError();

	}

	/*   Should not be here if potential has not been initalized (SET command)  */
	if (a->IsInitializedPotential==FALSE) {
		printf ("ERROR:  You must first use the POTENTIAL SET command.\n");
		CleanAfterError();
	}
        

	/*  Pass commands to appropriate potential routine  */

	/*  This routine wants a single string  */
	if (!strcmpi(s->pottype,"PAIR")) {
		PAIR_Parse(instr);

	} else if (!strcmpi(s->pottype,"STILL")) {
		STILL_Parse(instr);

	/*  These routines want front token split off  */
	} else {

		tokstr = strhed (&instr);

      if      (!strcmpi(s->pottype,"QUARTIC"))
         qu_read_potential (tokstr, instr);

      else if (!strcmpi(s->pottype,"EAM"))
         em_read_potential (tokstr, instr);

      else if (!strcmpi(s->pottype,"TERSOFF"))
         csi_read_potential (tokstr, instr);
	}

}



/*  read PRESSURE  */

void read_pressure   (char *instr)
   {
   char *LeadingToken;

   while (!IsBlank (instr))
      {

      LeadingToken = strhed (&instr);

      if (!strcmpi("OFF",LeadingToken))
         a->BoxMotionAlgorithm = BMA_NONE;

      else if (!strcmpi("ANDERSEN",LeadingToken))
			{
			/*  Set pressure algorithm to Andersen  */
         a->BoxMotionAlgorithm = BMA_ANDERSEN;

         /*  Read in box masses  */
         a->BoxMass[X] = dblstrf (&instr)/AVAG;
         a->BoxMass[Y] = dblstrf (&instr)/AVAG;
         a->BoxMass[Z] = dblstrf (&instr)/AVAG;

         /*  Quit if X mass is invalid  */
         if (a->BoxMass[X]<=0)
            {
            printf ("ERROR:  Box mass is less or equal to zero.\n");
            CleanAfterError();
            }

         /*  Set Y and Z to X mass if they are invalid  */
         if (a->BoxMass[Y]<=0)
            a->BoxMass[Y] = a->BoxMass[X];
         if (a->BoxMass[Z]<=0)
            a->BoxMass[Z] = a->BoxMass[Y];
         }

		else if (!strcmpi("CLAMP",LeadingToken))
			{

			/*  Set pressure algorithm to Andersen  */
         a->BoxMotionAlgorithm = BMA_CLAMP;

			/*  Convert from MBAR to erg/cm^3  */
			a->BulkModulus = dblstrf (&instr)/MBAR;

			/*  Test Bulkmodulus  */
			if (a->BulkModulus <= 0.0)
				{
				printf
					(
					"ERROR:  Bulk Modulus (%le) equal to or less than zero\n",
					a->BulkModulus*MBAR
					);
				CleanAfterError();
				}

			/*  Read in cstep if present  */
			if (!IsBlank(instr))
				a->VolClampStep = dblstrf(&instr);
			}

		/*  
		Read X,Y,Z external pressures
		if only one pressure,  then X,Y,Z pressure equal
		if two      pressures, then Y=Z
		*/
      else if (!strcmpi("EXTERNAL",LeadingToken))
			{
         a->Pressure[X] = dblstrf (&instr)/MBAR;
			if (IsBlank(instr))
				{
				a->Pressure[Y] = a->Pressure[Z] = a->Pressure[X];
				}
			else
				{
				a->Pressure[Y] = dblstrf (&instr)/MBAR;
				if (IsBlank(instr))
					{
					a->Pressure[Z] = a->Pressure[Y];
					}
				else
					{
					a->Pressure[Z] = dblstrf (&instr)/MBAR;
					}
				}
			}

		/*  Box dimensions expands/contracts uniformly  */
      else if (!strcmpi("ISOTROPIC",LeadingToken))
         a->BoxMotionGeometry = BMG_ISOTROPIC;

		/*  Box dimensions expand/contract independently  */
      else if 
			(
			!strcmpi("ORTHONORMAL",LeadingToken) ||
			!strcmpi("ORTHORHOMBIC",LeadingToken)
			)
         a->BoxMotionGeometry = BMG_ORTHORHOMBIC;

      else
         {
         printf ("ERROR:  Unrecognized option to PRESSURE command.\n");
         CleanAfterError();
         }

      }

   }



                           /*  READ PSTRAIN  */

void read_pstrain (char *instr)
   {
   int  i;
   double dot, eps, dx,dy,dz, nx,ny,nz;
   double (*Coord)[NDIR] = (double (*)[NDIR]) a->cur;

   eps  = dblstrf (&instr);

   dx   = dblstrf (&instr);
   dy   = dblstrf (&instr);
   dz   = dblstrf (&instr);

   nx   = dblstrf (&instr);
   ny   = dblstrf (&instr);
   nz   = dblstrf (&instr);

	CheckForNoneSelected();
   LOOP (i, a->np)
   if (IS_SELECT(i) )
      {
      dot = eps * (nx*Coord[i][0] + ny*Coord[i][1] + nz*Coord[i][2]);
      Coord[i][0] += dot*dx;
      Coord[i][1] += dot*dy;
      Coord[i][2] += dot*dz;
      }
   WrapParticles (a);
   }

                           /*  READ QUENCH  */
void read_quench (char *instr)
	{
   int nstep = intstrf (&instr);


	/*     Check for needed initializations  */
	if (!a->IsInitializedBox)
		{
		printf
			("ERROR:  Cannot peform QUENCH command without first setting BOX.\n");
		CleanAfterError();
		}

	if (!a->IsInitializedCoord)
		{
      printf ("ERROR:  Cannot peform QUENCH command without");
      printf (" first setting particle coordinates.\n");
      CleanAfterError();
      }

	if (!a->IsInitializedMass)
		{
		printf ("ERROR:  Cannot peform QUENCH command without");
		printf (" first setting all atoms' masses.\n");
		CleanAfterError();
		}

   if (!a->IsInitializedPotential)
      {
      printf ("WARNING:  Potential was not set");
      printf (" before performing molecular dynamics\n");
		IncrementNumberOfWarnings();
		}

   /*  TURN QUENCH FLAG ON  */
   if (*instr==0)
      s->quench=1;   

	/*  Read non-standard quench value  */
   else
      s->quench= intstrf (&instr);

	/*  Perform dynamics  */
   runcmd (a, s, nstep);

	/*  TURN QUENCH FLAG OFF */
   s->quench=0;   

	/*  Write out ending step  */
	if (s->PrintInfo)
		printf ("***   Current step is %li\n", a->step);

	}

                           /*  READ RCV FILE  */

/*  ADD FILE TO STACK in[]  */
void read_rcv (char *instr)  {
   long step;
   FILE *fin;
   char *tokstr = strhed (&instr);

   /*  OPEN FILE  */
   fin = fopen (tokstr, "rt");
   if (fin==NULL) {
      printf ("ERROR (read_rcv):  Cannot open input file %s.\n", tokstr);
      return;
   }
   /*  READ JUST FIRST STEP  */
   if (*instr==0)
      readrcv (fin, a);
   else {
      step = intstrf(&instr);
      do {
         readrcv (fin, a);
      } while (!feof(fin) && a->step!=step);
   }
   fclose (fin);

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
}


                           /*  READ FILE  */

/*  ADD FILE TO STACK in[]  */
void read_read (char *InputStr)
   {
   FILE *InputFile;
   char *LeadingTokenStr;
   char  BufferStr[NSTR];
   BOOLEAN QuitFlag;

   /*  Get leading token  */
   LeadingTokenStr = strhed (&InputStr);


   /*  Open file  */
   InputFile = fopen (LeadingTokenStr, "rt");
   if (InputFile==NULL)
      {
      printf ("ERROR:  Cannot open file %s.\n", LeadingTokenStr);
      CleanAfterError();
      }

   /*  Add file to input list  */
   m_add_list (&inlist, InputFile, "f");


   /*  Read file until end  */
   QuitFlag = FALSE;
   while (!QuitFlag && m_gets_list(BufferStr,NSTR,inlist)!=NULL)
      {
      read_command (BufferStr, &QuitFlag);
      }

   /*  Remove file from input list (also closes file)   */
   m_del_list(&inlist);

   }

                              /*  READ REFSTEP  */
void read_refstep (char *instr)
   {
   double *cp;
   double (*c)[3], (*r)[3];
   double b[3];
   int   i,j;
   unsigned char *ti;
   char *tokstr = NULL;

   tokstr = strhed (&instr);

   /*  TEST DATA INTEGRITY  */
   if (a->ref == NULL)
      {
      if      (!strcmpi(tokstr,"SWAP"))
         {
         printf ("ERROR (read_refstep):  Cannot swap particles because\n");
         printf ("   there are no reference particles.\n");
         CleanAfterError();
         }
      else if (!strcmpi(tokstr,"COPY"))
         {
         printf ("ERROR (read_refstep):  Cannot copy particles because\n");
         printf ("   there are no reference particles.\n");
         CleanAfterError();
         }
      }


   /*  Swap current and reference step  */
   if      (!strcmpi(tokstr,"SWAP"))
      {

      /*  Save current particle info  */
      cp   = a->cur;
      b[0] = a->bcur[0];
      b[1] = a->bcur[1];
      b[2] = a->bcur[2];
      ti   = a->type;

      /*  Copy reference info into current step  */
      a->cur     = a->ref;
      a->bcur[0] = a->bref[0];
      a->bcur[1] = a->bref[1];
      a->bcur[2] = a->bref[2];
      a->type    = a->rtype;

      /*  Copy saved current info into reference step  */
      a->ref     = cp;
      a->bref[0] = b[0];
      a->bref[1] = b[1];
      a->bref[2] = b[2];
      a->rtype   = ti;
      }

   /*  COPY REFERENCE PARTICLES TO CURRENT STEP  */
   else if (!strcmpi(tokstr,"COPY"))
      {
      /*  Copy box  */
      for (i=0;i<3;i++)
         a->bcur[i] = a->bref[i];

      /*  Copy types  */
      for (j=0;j<a->np;j++)
         a->type[j] = a->rtype[j];

      /*  Copy coordinates  */
      c = (double (*)[3]) a->cur;
      r = (double (*)[3]) a->ref;
      for (i=0;i<a->np;i++)
         for (j=0;j<3;j++)
            c[i][j] = r[i][j];
      }

   /*  CLEAR REFERENCE STEP  */
   else if (!strcmpi(tokstr,"CLEAR"))
      {
      FREE (a->ref)
      FREE (a->rtype)
      }

   /*  COPY CURRENT PARTICLES TO REFERENCE STEP  */
   else if (*instr==0)
      {

      /*  Allocate reference particle array  */
      if (a->ref == NULL)
			{
			ALLOCATE (a->ref, double, NDIR*a->NumPartAlloc)
			}

      /*  Allocate reference type array  */
      if (a->rtype == NULL)
         ALLOCATE (a->rtype, BYTE, a->NumPartAlloc)

      /*  Copy Particles  */
      c = (double (*)[3]) a->cur;
      r = (double (*)[3]) a->ref;
      for (i=0;i<3;i++)
         {
         a->bref[i] = a->bcur[i];
         for (j=0;j<a->np;j++)
            r[j][i] = c[j][i];
         }

      /*  Copy reference types  */
      for (j=0;j<a->np;j++)
         a->rtype[j] = a->type[j];
      }

   /*  Unrecognized Option  */
   else
      {
      printf ("WARNING:  Unrecognized option to REFSTEP command.\n");
      IncrementNumberOfWarnings();
      }

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }

                              /*  REMOVE ATOMS  */

#pragma argsused
void read_remove (char *instr)
   {
   int ipart;
   int isave;

   /*
   Copy saved particles (the unselected particles)
   to cover up "hole" that is removed
   */
   isave = 0;
	CheckForNoneSelected();
   LOOP (ipart, a->np)
      {
      if (!IS_SELECT(ipart) )
         {
         if (isave != ipart)
            {
            CopyParticle (a, isave, ipart);
            }
         isave++;
         }
      }

   /*  Set new number of particles  */
   a->np = isave;

   /*  Set number of selected particles to zero  */
   a->nsel = 0;

   /*  Realloate particle info  */
   ReallocateParticle (a, isave);

   /*  Mark neighbor list as invalid  */
   a->InvalidNeighborList = TRUE;

   /*  Print message  */
   printf ("***   NUMBER REMAINING %i\n", isave);

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }


                              /*  READ ROTATE  */
void read_rotate (char *InputString)
   {
   double Axis   [NDIR];
   double Center [NDIR];
   double Cosine;
   double Sine;
   double Rotation[NDIR][NDIR];
   double Asym    [NDIR][NDIR];
   double Magnitude;
   double Angle;
   int    idir;
   int    jdir;
	BOOLEAN UseSelect;


   /*  Read paramters  */
   Axis  [X] = dblstrf (&InputString);
   Axis  [Y] = dblstrf (&InputString);
   Axis  [Z] = dblstrf (&InputString);
   Center[X] = 1e-8 * dblstrf (&InputString);
   Center[Y] = 1e-8 * dblstrf (&InputString);
   Center[Z] = 1e-8 * dblstrf (&InputString);
   Angle     = M_PI * dblstrf (&InputString)/180.0;

   /*  No rotation is angle is zero  */
   if (Angle==0.0)
      return;

   /*  Normalize axis  */
   Magnitude = sqrt (Axis[X]*Axis[X] + Axis[Y]*Axis[Y] + Axis[Z]*Axis[Z]);
   if (Magnitude==0.0)
      {
      printf ("ERROR (read_rotate):  Null axis.\n");
      CleanAfterError();
      }
   Axis[X] /= Magnitude;
   Axis[Y] /= Magnitude;
   Axis[Z] /= Magnitude;

   /*  SET UP ASYMMETRIC MATRIX (REPRESENTS CROSS PRODUCT)  */
   Asym[X][X] = Asym[Y][Y] = Asym[Z][Z] = 0;
   Asym[X][Y] = -Axis[Z];
   Asym[Y][Z] = -Axis[X];
   Asym[Z][X] = -Axis[Y];
   Asym[Y][X] = -Asym[X][Y];
   Asym[Z][Y] = -Asym[Y][Z];
   Asym[X][Z] = -Asym[Z][X];

   /*  SET UP MATRIX  */
   Sine   = sin(Angle);
   Cosine = cos(Angle);
   for (idir=0; idir<NDIR; idir++)
   for (jdir=0; jdir<NDIR; jdir++)
      if (idir==jdir)
         {
         Rotation[idir][jdir] =
            Axis[idir]*Axis[jdir] +
            Cosine * (1 - Axis[idir]*Axis[jdir]) +
            Sine   * Asym[idir][jdir];
         }
      else
         Rotation[idir][jdir] =
            Axis[idir]*Axis[jdir] * (1 - Cosine) +
            Sine * Asym[idir][jdir];

   /*  TRANSFORM COORDINATES  */
	UseSelect = TRUE;
	CheckForNoneSelected();
   RotateCoordinates (Rotation, Center, UseSelect, a->cur, a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->v , a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->c2, a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->c3, a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->c4, a->np);
   RotateCoordinates (Rotation, Center, UseSelect, a->c5, a->np);
   RotateVectorList
      (Rotation, Center, UseSelect, a->ForcePtrList,  a->np);
   RotateVectorList
      (Rotation, Center, UseSelect, a->SpringPtrList, a->np);
   }


                           /*  READ RUN  */

void read_run (char *instr)  {
   a->run = lngstrf (&instr);
}


                            /*  READ SCALE  */

void read_scale    (char *instr)
   {
   double scale[3];
   int   idir, ipart;
   double (*Coord)[3] = (double (*)[3]) a->cur;


   /*  PARSE SCALE  */
   LOOP (idir, NDIR)
       scale[idir] = dblstrf (&instr);

   /*  Test for zero scale  */
   if (scale[X]==0.0)
      {
      printf ("ERROR:  Cannot scale by zero.\n");
      CleanAfterError();
      }

   if (scale[Y]==0.0)
       scale[Y] = scale[X];
   if (scale[Z]==0.0)
       scale[Z] = scale[Y];

   /*  Scale Box  */
   LOOP (idir, NDIR)
      {
      a->bcur[idir] *= scale[idir];
      }

   /*  Scale cavity contraint center and axis  */
   LOOP (idir, NDIR)
      {
      a->CavityCenter[idir] *= scale[idir];
      a->CavityAxis  [idir] *= scale[idir];
      }


   /*  Test for particles present  */
   if (!a->IsInitializedCoord)
      {
      printf ("WARNING:  Cannot scale particle without initializing them.\n");
      IncrementNumberOfWarnings();
      return;
      }

   /*  Scale particles  */
   LOOP (idir, NDIR)
      {
      LOOP (ipart, a->np)
         Coord[ipart][idir] *= scale[idir];
      }

   }





void read_screw    (char *instr)
   {
   int    idir;
   int    ipart;
   double (*Coord)[NDIR] = (double (*)[NDIR]) a->cur;
   double BurgersVector[NDIR];
   double Origin    [NDIR];
   double Axis   [NDIR];
   double Reference [NDIR];
   double Perpendicular[NDIR];
   double Position  [NDIR];
   double Angle;
   double Xpos;
   double Ypos;
   double ScrewFactor;

   /*  Parse Screw Burgers Vector  */
   LOOP (idir, NDIR)
      BurgersVector[idir] = 1e-8*dblstrf (&instr);
   LOOP (idir, NDIR)
      Origin[idir] = 1e-8*dblstrf (&instr);
   LOOP (idir, NDIR)
      Reference[idir] = 1e-8*dblstrf (&instr);


   /*
   If magnitude of burgers vector is zero,
   then no need to move atoms
   */
   if (GetMagnitude(BurgersVector)==0)
      {
      return;
      }

   /*  Set screw axis as normalized Burgers Vector  */
   LOOP (idir, NDIR)
      {
      Axis[idir] = BurgersVector[idir];
      }
   Normalize(Axis);


   /*  If no angle reference vector specified, create one  */
   if (GetMagnitude(Reference)==0.0)
      {
      Reference[X] =  Axis[Y];
      Reference[Y] = -Axis[Z];
      Reference[Z] =  Axis[X];
      }

   /*  Orthogonalize reference vector to axis  */
   Orthogonalize (Reference, Axis);
   Normalize(Reference);
   if (GetMagnitude(Reference)==0.0)
      {
      printf ("INPUT ERROR: Reference vector is parallel to Burgers vector.\n");
      CleanAfterError();
      }

   /*  Find vector perpendicular to Axis[] and Reference[]  */
   CalcCross (Perpendicular, Axis, Reference);


   /*
   For each particle,
      1)  find position vector relative to Origin[]
      2)  orthgonalize position vector relative to Axis[]
      3)  obtain angle between orthoganilized position vector and Reference[]
      4)  displace atom by factor proportional to burgers vector and
          sin() of angle between Reference[] and orthogonalized position vector
   */
	CheckForNoneSelected();
   LOOP (ipart, a->np)
   if (IS_SELECT(ipart))
      {

      /*  Find atom position relative to origin  */
      Position[X] = Coord[ipart][X] - Origin[X];
      Position[Y] = Coord[ipart][Y] - Origin[Y];
      Position[Z] = Coord[ipart][Z] - Origin[Z];

      /*  Orthoginalize position vector to Axis[]  */
      Orthogonalize (Position, Axis);

      /*
      Find X and Y componenent of position given that Axis[] is Z direction
      */
      Xpos = GetDot (Position, Reference);
      Ypos = GetDot (Position, Perpendicular);
      Angle = atan2 (Ypos, Xpos);
      ScrewFactor = Angle/(2.0*M_PI);

      LOOP (idir, NDIR)
         {
         Coord[ipart][idir] += ScrewFactor * BurgersVector[idir];
         }

      }

   /*  Wrap displaced particles throught repeating boundary conditions  */
   WrapParticles (a);
   }


                           /*  READ SEED  */
void read_seed (char *instr)
   {
   long l = lngstrf (&instr);
   l &= 0xffff;
   if (l!=0)
      {
      s->seed = (unsigned) l;
      initrand (s->seed);
      }
   }


                           /*  READ SIZE */

void read_size   (char *instr)
   {
   s->size  = ANGS*dblstrf (&instr);
   }

                           /*  READ SSAVE FLAG  */

void read_ssave  (char *instr)
   {
   ReadStepOutput ( &(s->StressStepOutput), instr);
   }

                           /*  READ STATE */

void read_state   (char *instr)
   {
   readstate (&a, s, instr);

	/*  Update degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }
                           /*  READ STEP */

void read_step   (char *instr)
   {
   char *tokstr = strhed(&instr);

   if (*tokstr=='+')
      {
      tokstr++;
      a->step += intstrf (&tokstr);
      }
   else
      a->step  = intstrf  (&tokstr);
   }


                            /*  READ STRESS */

void read_stress (char *instr) {
	char *tokstr = strhed (&instr);
	if (!strcmpi("THERMAL", tokstr)) {
		tokstr = strhed (&instr);
		if (!strcmpi("ON", tokstr)) {
			a->stress_thermal = TRUE;
		} else if (!strcmpi("OFF", tokstr)) {
			a->stress_thermal = FALSE;
		} else {
			printf ("ERORR:  Unrecognized option to STRESS THERMAL command.\n");
      	CleanAfterError();
		}
	} else {
		printf ("ERORR:  Unrecognized option to STRESS command.\n");
     	CleanAfterError();
	}
}


                           /*  READ SURFACE  */

void read_surface  (char *instr)
   {
   char *tokstr = NULL;
   BOOLEAN OnFlag;
   int  NumWrap;

   /*  GET FIRST "OFF" or "ON"  */
   tokstr = strhed (&instr);
   if      (!strcmpi(tokstr,"ON"))
      OnFlag = TRUE;
   else if (!strcmpi(tokstr,"OFF"))
      OnFlag = FALSE;
   else
      {
      printf ("        ***  ERRROR:  UNKNOWN OPTION");
      return;
      }

   /*  PARSE REMAINING OPTIONS  */
   tokstr = strhed (&instr);
   while (strcmpi(tokstr,""))
      {
      if      (!strcmpi(tokstr,"ON"))
         OnFlag = TRUE;
      else if (!strcmpi(tokstr,"OFF"))
         OnFlag = FALSE;
      else if (!strcmpi(tokstr,"X"))
         a->surf[0] = OnFlag;
      else if (!strcmpi(tokstr,"Y"))
         a->surf[1] = OnFlag;
      else if (!strcmpi(tokstr,"Z"))
         a->surf[2] = OnFlag;
      else
         printf ("        ***  ERRROR:  UNKNOWN OPTION \"%s\"\n",tokstr);
      tokstr = strhed (&instr);
      }

   /*  Wrap particles  */
   NumWrap = WrapParticles (a);

   if (NumWrap>0)
      PrintWrapWarning (NumWrap);

   /*  Mark neighbor list as invalid  */
   a->InvalidNeighborList = TRUE;

	/*  Make surface as initialized  */
	a->IsInitializedSurface = TRUE;

   }



                           /*  READ TSAVE FLAG  */

void read_tsave  (char *instr)
   {
   ReadStepOutput ( &(s->TrajectoryStepOutput), instr);
   }

                           /*  READ TYPE  */

void read_type     (char *instr)
   {
   int i;
   int curtype;

   curtype = intstrf (&instr)-1;
   LOOP (i, a->np)
      if (IS_SELECT(i))
			a->type[i] = curtype;
   }


void read_typelist (char *instr)
   {
   char *tokstr = NULL;
   char bufstr[NSTR];
   int i, nread, ntype, sflag;
   int *typelist = NULL;

   tokstr = strhed (&instr);
   if (!strcmpi(tokstr,"SEL"))
      {
      sflag = TRUE;
      ntype = intstrf (&instr);
      }
   else
      {
      sflag = FALSE;
      ntype = intstrf (&tokstr);
      }

   if (ntype==0)
      printf ("ERROR:  Number of types is 0\n");

      /*  READ TYPES  */
   else
      {
      /*  ALLOCATE ARRAY  */
      ALLOCATE (typelist, int, ntype)

      nread = 0;
      while (   nread<ntype
                && NULL!=m_gets_list(bufstr,NSTR,inlist))
         {
         tokstr = bufstr;
         while (*tokstr)
            typelist[nread++] = intstrf(&tokstr)-1;
         }
      if (nread<ntype)
         {
         printf ("Number of particles expected is %d.\n", ntype);
         printf ("Number of particles read     is %d.\n", nread);
         CleanAfterError();
         }

		/*  Check for no atoms selected  */
		if (sflag)
			CheckForNoneSelected();

      /*  ASSIGN TYPES  */
      nread = 0;
      LOOP (i, a->np)
         {
         if (!sflag || IS_SELECT(i) )
            a->type[i] = typelist[nread++];
         }
      /*  FREE ALLOCATION  */
      FREE (typelist)
      }
   }


/*  Read name string to associate with each type (used to WRITE XMOL command)  */
void read_typename (char *instr)
   {
   int   Type;
   char *NameStr;

   /*  Raad type name pairs on line  */
   while (*instr != '\0')
      {

      /*  Read type and name string   */
      Type    = intstrf (&instr)-1;
      NameStr = strhed  (&instr);

      /*  Insure that type number in range  */
      if (Type>MAX_TYPE_NAMES)
         {
         printf ("ERROR:  Maximum type for TYPENAME command is %i\n",
            MAX_TYPE_NAMES);
         CleanAfterError();
         }

      /*  Store type name string  */
      ALLOCATE (a->TypeNames[Type], char, strlen(NameStr)+1)
      strcpy   (a->TypeNames[Type], NameStr);

      }
   }


void read_verbose  (char *instr)
	{
	char *tokstr = NULL;

	tokstr = strhed (&instr);

	/*  VERBOSE FILE { ON | OFF }  */
	if (!strcmpi("FILE", tokstr))
		{
		tokstr = strhed (&instr);
		if (!strcmpi("ON",tokstr))
			{
			s->UseVerboseOutputFile = TRUE;
			}
		else if (!strcmpi("OFF",tokstr))
			{
			s->UseVerboseOutputFile = FALSE;
			}
		else
			{
   		printf ("WARNING:  Unknown option to VERBOSE FILE.\n");
  			IncrementNumberOfWarnings();
			}
		}
	/*  VERBOSE FILE { ON | OFF }  */
	else if (!strcmpi("INFO", tokstr))
		{
		tokstr = strhed (&instr);
		if (!strcmpi("ON",tokstr))
			{
			s->PrintInfo = TRUE;
			}
		else if (!strcmpi("OFF",tokstr))
			{
			s->PrintInfo = FALSE;
			}
		else
			{
   		printf ("WARNING:  Unknown option to VERBOSE INFO.\n");
  			IncrementNumberOfWarnings();
			}
		}

	/*  Incorrect Verbose command  */
	else
		{
  		printf ("WARNING:  Unknown option to VERBOSE.\n");
		IncrementNumberOfWarnings();
		}
	}


void read_velocity (char *instr)
   {
   char *VelocityTypeStr;
   double Direction[NDIR];
   double Velocity [NDIR];
   double Magnitude;
   double NormalizationMag;
   double MagnitudeChange;


   /*  Chech number of selected atoms  */
	CheckForNoneSelected();
   if (a->nsel==0)
      return;

   /*  Get velocity type (either linear or angular)  */
   VelocityTypeStr = strhed (&instr);

   /*  Get velocity direction (linear) or axis (angular) */
   Direction[X] = dblstrf (&instr);
   Direction[Y] = dblstrf (&instr);
   Direction[Z] = dblstrf (&instr);

   /*  Get velocity magnitude  */
   Magnitude    = dblstrf (&instr);

   /*  Set magnitude of direction  */
   NormalizationMag =
      Direction[X]*Direction[X] +
      Direction[Y]*Direction[Y] +
      Direction[Z]*Direction[Z] ;

   /*  Check for null vector  */
   if (NormalizationMag==0.0)
      {
      printf ("WARNING:  Direction vector is 0.\n");
      IncrementNumberOfWarnings();
      return;
      }

   /*  Normalize Direction  */
   MagnitudeChange = Magnitude / sqrt(NormalizationMag);
   Velocity[X] = Direction[X] * MagnitudeChange;
   Velocity[Y] = Direction[Y] * MagnitudeChange;
   Velocity[Z] = Direction[Z] * MagnitudeChange;

   /*
   Convert velocity from (cm or radians)/sec
   to (cm or radians)/time step
   */
   Velocity[X] *= s->dtime;
   Velocity[Y] *= s->dtime;
   Velocity[Z] *= s->dtime;

   /*  Write error message  */
   if
   (
   strcmpi(VelocityTypeStr,"LINEAR" )!=0 &&
   strcmpi(VelocityTypeStr,"ANGULAR")!=0
   )
      {
      printf ("ERROR:  Velocity type not recognized.\n");
      CleanAfterError();
      }

   /*  Allocate velocity array if not allocated  */
   if (a->v==NULL)
      {
      ALLOCATE (a->v, double, NDIR*a->NumPartAlloc)
      }

   /*  Set velocities  */
   if (!strcmpi(VelocityTypeStr,"LINEAR"))
      {
      SetVelocity
         (
         (double (*)[NDIR]) a->v,
         a->mass,
         a->Selection,
         a->np, Velocity
         );
      }
   else if (!strcmpi(VelocityTypeStr, "ANGULAR"))
      {
      SetAngularVelocity
         (
         (double (*)[NDIR]) a->cur,
         (double (*)[NDIR]) a->v,
         a->mass,
         a->Selection,
         a->np, Velocity
         );
      }
   else
      {
      printf ("ERROR:  Internal error in read_velocity.\n");
      CleanAfterError();
      }

   }



void read_wave (char *InputStr)
   {
   double  (*Coord)[NDIR];
   double  Phase;
   double  WaveVector[NDIR];
   double  Disp[NDIR];
   double  Angle;
   double  Cosine;
   int     icoord;


   /*  Convert coordinate pointer   */
   Coord = (double (*)[NDIR]) a->cur;

   /*  Read paramters  */
   Phase         =        dblstrf (&InputStr);
   Disp      [X] = 1e-8 * dblstrf (&InputStr);
   Disp      [Y] = 1e-8 * dblstrf (&InputStr);
   Disp      [Z] = 1e-8 * dblstrf (&InputStr);
   WaveVector[X] = 1e+8 * dblstrf (&InputStr);
   WaveVector[Y] = 1e+8 * dblstrf (&InputStr);
   WaveVector[Z] = 1e+8 * dblstrf (&InputStr);


	CheckForNoneSelected();
   LOOP (icoord, a->np)
      {
      if (IS_SELECT(icoord) )
         {
         Angle =
            WaveVector[X] * Coord[icoord][X] +
            WaveVector[Y] * Coord[icoord][Y] +
            WaveVector[Z] * Coord[icoord][Z] +
            Phase;
         Cosine = cos(Angle);
         Coord[icoord][X] += Cosine * Disp[X];
         Coord[icoord][Y] += Cosine * Disp[Y];
         Coord[icoord][Z] += Cosine * Disp[Z];
         }
      }
   }



void ZeroDerivatives (int Index, Particle_t *a)
   {
   int MemoryOffset;
   int MemorySize;

   MemoryOffset = 3*Index;
   MemorySize   = 3*sizeof(a->cur[0]);

   if (a->v )  memset (a->v  + MemoryOffset, 0, MemorySize);
   if (a->c2)  memset (a->c2 + MemoryOffset, 0, MemorySize);
   if (a->c3)  memset (a->c3 + MemoryOffset, 0, MemorySize);
   if (a->c4)  memset (a->c4 + MemoryOffset, 0, MemorySize);
   if (a->c5)  memset (a->c5 + MemoryOffset, 0, MemorySize);
   }



void ReadStepOutput (StepOutput_t *StepOutput, char *InputString)
   {
   char *FileName;

   StepOutput->StepIncrement = intstrf (&InputString);
   StepOutput->Save = (StepOutput->StepIncrement!=0);
   StepOutput->CurrentStep  = 0;
   FileName = strhed (&InputString);
   strcpy (StepOutput->FileName, FileName);

   /*  Determine if file opened as write (not append) first time  */
   if (*(StepOutput->FileName)=='+')
      StepOutput->IsFileNew = FALSE;
   else
      StepOutput->IsFileNew = TRUE;
   }


void PrintWrapWarning (int NumWrap)
   {
   if (NumWrap>0)
      {
      printf ("WARNING:  %i input particles have been \"wrapped\"\n", NumWrap);
      printf (" to conform to repeating boundary conditions.\n");
      IncrementNumberOfWarnings();
      }
   }



/*
Evaluate multiple formula in string InputFormula.  Last NumOutput formulas assigned
to OutputValues.  If not enough formulas, missing values assigned to zero.

InputFormula      string containing  multiple formula, separated by one or more blanks
InputValues       are mapped to variables X, Y, Z (maximum of 3 input values)
NumInput          number of input values.
OutputValues      number of output values, assigned values of last NumOutput formulas
NumOutput         number of output values.
*/

#if 0
void EvaluateFormula
(
char   *InputFormula,
double  InputValues[NDIR],
double  InputConversion,
double *OutputValues,
double  OutputConversion,
int     NumOutput
)
   {
   int   NumFormula;
   char  FormulaBuffer[MAX_BUFFER];
   char *FormulaPtr;
   double *TempValue=NULL;

   /*  Assign values to coordinate variables X, Y, Z  */
   AssignVariable("X", InputConversion*InputValue[X]);
   AssignVariable("Y", InputConversion*InputValue[Y]);
   AssignVariable("Z", InputConversion*InputValue[Z]);

   /*  Allocate storage for temporary values  */
   ALLOCATE(TempValue,double,NumOutput)

   NumFormula = 0;
   IndFormula = 0;
   FormulaPtr = InputFormula;
   while (!IsBlank(FormulaPtr))
      {
      TempValue[IndFormula] = dblstrf(&FormulaPtr);
      NumFormula++;
      IndFormula++;
      if (IndFormula==NumOutput)
         IndFormula=0;
      ASSERT(IndFormula<NumOutput)
      }

   /*  Determine number of values than can be used as output  */
   NumRead = (NumFormula<NumOutput ? NumFormula : NumOutput);

   /*  Position index at first usable value  */
   IndFormula -= NumRead;
   if (IndFormula < 0)
      IndFormula += NumOutput;

   /*  Copy temporary values to output array  */
   LOOP (i, NumRead)
      {
      OutputValues[i] = OutputConversion*TempValue[IndFormula];
      IndFormula++;
      if (IndFormula>NumOutput)
         IndFormula -= NumOutput;
      }

   /*  Set unassigned values to last assigned value  */
   if (NumRead < NumOutput)
      {
      for (i==NumRead; i<NumOutput; i++)
         OutputValues[i] = OutputValues[NumRead-1];
      }

   }
#endif




void EvaluateFormula
(
char   *InputFormula,
double *InputValues,
double  InputConversion,
double *OutputValues,
double  OutputConversion,
int     NumOutput
)
	{
	int   NumFormula;
	char  FormulaBuffer[MAX_BUFFER];
   char *FormulaPtr;
   double TempValue;


   /*  Assign values to coordinate variables X, Y, Z  */
   AssignVariable("X", InputConversion*InputValues[X]);
   AssignVariable("Y", InputConversion*InputValues[Y]);
   AssignVariable("Z", InputConversion*InputValues[Z]);


   /*  Copy input formula to buffer   */
	strncpy (FormulaBuffer, InputFormula, MAX_BUFFER);
   FormulaPtr = FormulaBuffer;


   /*
   Evaluate formula,  left to right through string,
   pushing the results into the OuputValues from  its end
   */
   NumFormula = 0;
   while (FormulaPtr[0]!='\0')
      {
      TempValue = OutputConversion*dblstrf (&FormulaPtr);
		PushArray (TempValue, OutputValues, NumOutput, PUSH_BACKWARD);
      NumFormula++;
      }

   /*   Set remaining unspecified output values to zero  */
   while (NumFormula<NumOutput)
      {
      PushArray (0.0, OutputValues, NumOutput, PUSH_BACKWARD);
      NumFormula++;
      }

   }



/*  Push new element into start (end) of array, dropping last (first) element  */
void PushArray (double Value, double *Array, int NumArray, int Direction)
   {
   int i;

   if (Direction==PUSH_FORWARD)
      {
      for (i=NumArray-1; i>0; i--)
         {
         Array[i] = Array[i-1];
         }
      Array[0] = Value;
      }

   else if (Direction==PUSH_BACKWARD)
      {
      LOOP (i, NumArray-1)
         {
         Array[i] = Array[i+1];
         }
      Array[NumArray-1] = Value;
      }

   else
      {
      printf ("INTERNAL ERROR:  Push array called with invalid direction\n");
      CleanAfterError();
      }
   }


/*
Null function called in place of potential model FREE() 
*/
void POT_Free(void)
	{
	}


/*  
Calculate total degrees of freedom for particles, 
DOES NOT including box motion  
*/
int GetTotalDOF (Particle_t *a) {
	int ipart;
	int df = 0;

	LOOP (ipart, a->np) 
		if (!IS_FIX(ipart)) {
			if (a->cons && a->cons[ipart]) {
				df += ( (a->cons[ipart]->type==CONSTR_LINE) ? 1 : 2 );
			} else {
				df += NDIR;
			}
		}
	return df;
}


/*  Calculate Min,Max Coordinates
     * Use only selected coordinates if UseSelect=TRUE
     * Use box limits unless DontUseBox==TRUE and free surface
*/
void GetMinMaxCoord(
Particle_t *a, 
double MinCoord[NDIR], 
double MaxCoord[NDIR], 
BOOLEAN UseSelect, 
BOOLEAN DontUseBox) {
	int ipart, idir, FirstPart;
	double TempCoord;

	/*  Find Min,Max coord in each direction  */
	for (idir=0;idir<NDIR;idir++) {
		/*  Get Min,Max coord for free surface  */
		if (DontUseBox && a->surf[idir]) {
			/*  Find first selected particle  */
			FirstPart = -1;
			if (UseSelect) {
				for (ipart=1;ipart<a->np;ipart++) {
					if (IS_SELECT(ipart)) {
						FirstPart = ipart;
						break;
					}
				}
				if (FirstPart==-1) {
      					printf ("INTERNAL ERROR:  UseSelect is set but there are no selected particles\n");
					CleanAfterError();
				}
			} else {
				FirstPart = 0;
			}
			/*  Intialize Min,Max values  */
			MinCoord[idir] = MaxCoord[idir] = a->cur[FirstPart*NDIR + idir];
			/*  Compare Min,Max value with all (selected) particles  */
			for (ipart=FirstPart+1;ipart<a->np;ipart++) {
				if ( (!UseSelect) || IS_SELECT(ipart) ) {
					TempCoord = a->cur[ipart*NDIR + idir];
					if (MinCoord[idir]>TempCoord)
						MinCoord[idir] = TempCoord;
					else if (MaxCoord[idir]<TempCoord)
						MaxCoord[idir] = TempCoord;
				}
			}
		/*  Min,Max coord taken from box  */
		} else {
			MinCoord[idir] = 0.0;
			MaxCoord[idir] = a->bcur[X];
		}
	}
}