cdfunc.c

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

/*
************************************************************************
History
************************************************************************
*/

/*
First pass at general function routine
*/


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

/*  LOCAL:  Free standing compilation, does not use XMD routines  */
/*
#define LOCAL
*/




/*
************************************************************************
Include Files
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "cdfunc.h"
#include "cdhouse.h"

#ifndef LOCAL
#include "iomngr.h"
#include "strsub.h"
#include "parse.h"
#endif



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

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

#define NBUFF 1024


#define DATA_NULL  0
#define DATA_TABLE 1
#define DATA_POLY  2


#define DERIV_ZERO 0
#define DERIV_ONE  1


#ifdef LOCAL
#define LIST FILE
#endif


/*
************************************************************************
Macros
************************************************************************
*/
#ifdef LOCAL

#define FGETS(STRING,SIZE,FILE) \
   fgets (STRING,SIZE,FILE)
#define DBLSTRF(STRING)  dblstr(STRING)
#define INTSTRF(STRING)  intstr(STRING)
#define STRCMPI(S1,S2)   strcasecmp(S1,S2)
#define STRNCMPI(S1,S2,N)  strncasecmp(S1,S2,N)

#else

#define FGETS(STRING,SIZE,FILE) \
   m_gets_list (STRING,SIZE,FILE)
#define DBLSTRF(STRING)  dblstrf(STRING)
#define INTSTRF(STRING)  intstrf(STRING)
#define STRCMPI(S1,S2)   strcmpi(S1,S2)
#define STRNCMPI(S1,S2)  strncmpi(S1,S2)

#endif

/*
************************************************************************
Type Definitions
************************************************************************
*/


/*
************************************************************************
Global Variables
************************************************************************
*/

#ifdef LOCAL
int Debug_g = TRUE;

FILE *inlist = stdin;
BOOLEAN CheckMem_g = FALSE;

#else

extern LIST         *inlist;

#endif

extern BOOLEAN UseCheckFunc_g;



/*
************************************************************************
Module-Wide Variables
************************************************************************
*/




/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
static void ReadTableFunction      (Function_t *, char *, LIST *);
static void ReadPolyFunction       (Function_t *, char *, LIST *);
static void ScaleTableFunction     (Function_t *, double, int);
static void ScalePolyFunction      (Function_t *, double, int);
#if 0
static void CalcArrayTableFunction (Function_t *, char *, double *);
#endif

static void InterpolateTable       (Function_t *, double *);

static double GetTableValue (Function_t *, double);
static double GetTableDeriv (Function_t *, double);
static double GetNullValue (Function_t *F, double x);
static double GetPolyValue  (Function_t *, double);
static double GetPolyDeriv  (Function_t *, double);
static double OutOfRangeValue (Function_t *F, double x, int Deriv);

static double poly (double, int, double *);


#ifdef LOCAL
BOOLEAN IsWhiteSpace(char);
BOOLEAN IsBlank (char *);
char    *strhed (char **);
double  dblstr(char **);
double  intstr(char **);
void    CleanAfterError(void);
#endif



/*
************************************************************************
Exported Routines
************************************************************************
*/

int FUNC_Parse (Function_t *F, char *InputStr)
	{
	char *HeadStr = NULL;

	
	/*  Check commands other than input commands  */
	if (IsFirstToken("SCALE",InputStr))
		{
		strhed (&InputStr);
		FUNC_ParseScale (F, InputStr);
		return FUNC_COMMAND_SCALE;
		}

	/*  
	Pass on to FUNC_ParseRead
	*/
	else
		{
		FUNC_ParseRead (F, InputStr);
		return FUNC_COMMAND_READ;
		}
	}



void FUNC_ParseScale (Function_t *F, char *InputStr)
	{
	char *HeadStr = NULL;
	double Scale;

	HeadStr = strhed (&InputStr);
	Scale   = dblstrf(&InputStr);

	if (!STRCMPI("INPUT",HeadStr))
		{
		FUNC_Scale (F, Scale, FUNC_INPUT);
		}
	else if (!STRCMPI("OUTPUT",HeadStr))
		{
		FUNC_Scale (F, Scale, FUNC_OUTPUT);
		}
	else
		{
		ERROR_PREFIX
		printf ("Unrecognized option for SCALE command.\n");
		CleanAfterError();
		}
	}


/*
Read function table or polynomial from input
  - Input is throught macro FGETS(STRING,SIZE,FILE)
  - Received string of the form(s)
      [ TABLE ] n InteriorCutoff ExteriorCutoff
           x0 x1 ... xn-1
      POLY      n0 nn x0 ExteriorCutoff
           an0 an1 .. ann
      f(x) = an0*(x-xc)^n0 + ... ann*(x-xc)^nn
*/
void FUNC_ParseRead (Function_t *F, char *InputStr)
   {

   /*
   Test for valid function boundary
   */
   if (FUNC_EXTRAP_UNDEFINED==GET_INT_EXTRAP(F->BoundaryType))
      {
      printf ("INTERNAL ERROR: File %s Line %i\n",
         __FILE__, __LINE__);
      printf ("Invalid Interior Extrap value\n");
      CleanAfterError();
      }

   if (FUNC_EXTRAP_UNDEFINED==GET_EXT_EXTRAP(F->BoundaryType))
      {
      printf ("INTERNAL ERROR: File %s Line %i\n",
         __FILE__, __LINE__);
      printf ("Invalid Exterior Extrap value\n");
      CleanAfterError();
		}


   /*
   Polynomial Data Type
   */
	if (IsFirstToken("POLY",InputStr))
      {
		strhed (&InputStr);
      F->DataType = DATA_POLY;
      ReadPolyFunction (F, InputStr, inlist);
      }

   /*
   Table Data type (default)
   */
   else
      {
      F->DataType = DATA_TABLE;

      /*
      TABLE token found, remove it
      */
		if (IsFirstToken("TABLE",InputStr))
         {
			strhed (&InputStr);
         }

      ReadTableFunction (F, InputStr, inlist);
      }

   }


Function_t *FUNC_CreateNullFunction(void)
	{
	Function_t *F = NULL;

	/*  Allocate storage  */
	ALLOCATE (F, Function_t, 1)

	/*  Set function pointer to 'null' function  */
	F->f0 = (FunctionCall_t*) GetNullValue;
	F->f1 = (FunctionCall_t*) GetNullValue;
   F->DataType = DATA_NULL;

	return F;
	
	}

void FUNC_Free (Function_t **F)
	{
	if (*F!=NULL)
		{
		FREE ((*F)->Coeff)
		FREE ((*F))
		}
   }



/*  Scale function  */
void FUNC_Scale (Function_t *F, double Scale, int ScaleType)
   {

   /*  Check scale  */
      if (Scale==0 && FUNC_INPUT==ScaleType)
         {
         printf ("ERROR:  Cannot scale function input by 0.\n");
         CleanAfterError();
         }

   /*  Scale value common to both types  */
   if (FUNC_INPUT==ScaleType)
      {
      F->InteriorCutoff *= Scale;
      F->ExteriorCutoff *= Scale;
      F->InteriorSlope  /= Scale;
      F->ExteriorSlope  /= Scale;
      }
   else if (FUNC_OUTPUT==ScaleType)
      {
      F->InteriorValue *= Scale;
      F->ExteriorValue *= Scale;
      F->InteriorSlope *= Scale;
      F->ExteriorSlope *= Scale;
      }

   /*  Scale function specific values  */
   if (DATA_POLY==F->DataType)
      ScalePolyFunction(F, Scale, ScaleType);
   else
      ScaleTableFunction(F, Scale, ScaleType);
   }


double FUNC_GetValue (Function_t *F, double x)
   {
   if      (DATA_POLY ==F->DataType)
      return GetPolyValue (F, x);
	else if (DATA_TABLE==F->DataType)
      return GetTableValue (F, x);
	else
		return GetNullValue  (F, x);
   }



double FUNC_GetDeriv (Function_t *F, double x)
   {
   if (DATA_POLY==F->DataType)
      return GetPolyDeriv (F, x);
   else
      return GetTableDeriv (F, x);
   }




double FUNC_GetCutoff(Function_t *F)
   {
   return F->ExteriorCutoff;
   }



/*
************************************************************************
Local Functions
************************************************************************
*/
static void ReadTableFunction
(
Function_t *F,
char *InputStr,
LIST *InFile
)
	{
   int idata;
   int nread;
   int ErrorCode;
   char ReadBuffer[NBUFF];
   char TempBuffer[NBUFF];
	char *TokenStr = NULL;
	char *FormPtr  = NULL;
	char *TempPtr  = NULL;
	double *a = NULL;
	double *atmp = NULL;
	double DelR;
	double r;

	F->NumCoeff = INTSTRF(&InputStr) - 1;
	F->InteriorCutoff = DBLSTRF (&InputStr);
	F->ExteriorCutoff = DBLSTRF (&InputStr);
	F->f0 = (FunctionCall_t *) GetTableValue;
	F->f1 = (FunctionCall_t *) GetTableDeriv;

	/*  Sanity Check  */

   if (F->NumCoeff <  2)
      {
      printf
         ("ERROR:  ");
      printf
         ("There must be at least 2 data points in table potenial.\n");
		CleanAfterError();
      }

   if (F->ExteriorCutoff <= F->InteriorCutoff)
		{
      printf ( "ERROR:  Exterior Cutoff (%e)",
         F->ExteriorCutoff);
		printf (" is less than or equal to interior cutoff (%e).\n",
			F->InteriorCutoff);
      CleanAfterError();
      }

   F->TableFactor =
      F->NumCoeff / ( F->ExteriorCutoff - F->InteriorCutoff );

   /*  Make temporary space for input values  */
	ALLOCATE (atmp, double, F->NumCoeff+5);
   a     = atmp+2;

   /*  Make permanant space for coefficients  */
	ALLOCATE (F->Coeff, double, 6*F->NumCoeff)

	/*  Get next token  */
	TokenStr = strhed (&InputStr);

	/*  Check for valid token  */
	if (!IsBlank(TokenStr) && STRCMPI("FORMULA",TokenStr))
		{
		ERROR_PREFIX
		printf ("Unknown option \"%s\"\n", TokenStr);
		CleanAfterError();
		}

   /*  Read formula  */
	if (!STRCMPI("FORMULA",TokenStr))
      {

		FGETS(ReadBuffer,NBUFF,InFile);
		DelR = (F->ExteriorCutoff - F->InteriorCutoff)/F->NumCoeff;
		/*  
		Loop to NumCoeff+1 because, in addition to values at 
		start of all intervales, also include value at end of 
		of last interval.
		*/
      LOOP (idata, F->NumCoeff+1)
         {

         /*  Write R (in angstroms) to parsing variable R  */
         r = F->InteriorCutoff + idata * DelR;
         sprintf (TempBuffer, "R=%le", r);
         TempPtr = TempBuffer;
         dblstrf (&TempPtr);

         /*  Evaluate formula (will use value of R passed via dblstrf()  */
         FormPtr = ReadBuffer;
			a[idata] = evalform (&FormPtr, &ErrorCode);

         /*  Error checking  */
			if (ErrorCode)
            {
				printf
					(
					"ERROR (pr_read_potential):  Cannot parse formula\n"
					);
				printf
					(
					"   Parser returns following message: %s\n",
					parsemsg(ErrorCode)
					);
				CleanAfterError();
				}
			}
		}

	/*  Read input table  */
	else
		{
      nread = 0;
		while (nread<F->NumCoeff+1 && NULL!=FGETS(ReadBuffer,NBUFF,InFile))
			{
			TokenStr = ReadBuffer;
			while (*TokenStr && nread < F->NumCoeff+1)
				{
				a[nread] = DBLSTRF (&TokenStr);
				nread++;
            }
         }

      if (nread<F->NumCoeff+1)
         {
         printf ("ERROR: Potential table too short.\n");
			CleanAfterError();
         }
      }

   ASSERT (FUNC_EXTRAP_UNDEFINED!=GET_INT_EXTRAP(F->BoundaryType))
   ASSERT (FUNC_EXTRAP_UNDEFINED!=GET_EXT_EXTRAP(F->BoundaryType))

	/*  With data array in hand, call table interpolation routine  */
	InterpolateTable(F, a);

	/*  Free unneeded allocation  */
	FREE (atmp)
	}


#if 0
/*
Like ReadTableFunction() except data table passed as an array.
NOTE:  Array starts at a[-2] and runs to a[n+2]
Data stored in a[0] to a[n-1].
*/
static void CalcArrayTableFunction
(
Function_t *F,
char *InputStr,
double *a
)
	{

   F->NumCoeff = INTSTRF(&InputStr) - 1;
   F->InteriorCutoff = DBLSTRF (&InputStr);
   F->ExteriorCutoff = DBLSTRF (&InputStr);
   F->f0 = (FunctionCall_t *) GetTableValue;
   F->f1 = (FunctionCall_t *) GetTableDeriv;

   /*  Sanity Check  */

   if (F->NumCoeff <  2)
		{
		printf
         ("ERROR:  ");
      printf
         ("There must be at least 2 data points in table potenial.\n");
      CleanAfterError();
      }

   if (F->ExteriorCutoff <= F->InteriorCutoff)
      {
      printf ( "ERROR:  Exterior Cutoff (%e)",
         F->ExteriorCutoff);
		printf (" is less than or equal to interior cutoff (%e).\n",
			F->InteriorCutoff);
      CleanAfterError();
      }

   F->TableFactor =
      F->NumCoeff / ( F->ExteriorCutoff - F->InteriorCutoff );

   /*  Make permanant space for coefficients  */
   ALLOCATE (F->Coeff, double, 6*F->NumCoeff)

   /*  With data array in hand, call table interpolation routine  */
	InterpolateTable(F, a);

	}
#endif


void InterpolateTable (Function_t *F, double *a)
	{
	double *p;
	int i,j,k,l,m,n;

	/*  Condition temporary array according to boundary conditions  */
	switch (GET_INT_EXTRAP(F->BoundaryType))
		{
		case FUNC_EXTRAP_CONST_SLOPE:
			a[ -1] = 2.0*a[  0] - a[1];