cdsubs.c

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

/*  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.
*/

/*
************************************************************************
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
************************************************************************
*/
/*  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;
   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();
      }