cdpairf.c

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

				/*  Sum internal stresses   */
            if
            (
            s->StressStepOutput.Save ||
            a->BoxMotionAlgorithm!=BMA_NONE
            )
               {
               a->Stress[XX] -= tfx*Disp[X];
               a->Stress[YY] -= tfy*Disp[Y];
               a->Stress[ZZ] -= tfz*Disp[Z];
               a->Stress[YZ] -= tfy*Disp[Z];
               a->Stress[ZX] -= tfz*Disp[X];
               a->Stress[XY] -= tfx*Disp[Y];
               }
            }
            /*  if (RadiusSqr < CutoffSqr_m[tc]) */
         }
         /*  LOOP (jatom)  */
      }
      /* LOOP(iatom)  */


   }


/*
Reads subcommands from POTENTIAL statement
   (POTENTIAL ...) INIT npair
   (POTENTIAL ...) EVAL PAIR t1 t2 deriv x
   (POTENTIAL ...) PAIR t1 t2 (.. input to FUNC_Read )

*/
int PAIR_Parse (char *instr)
	{
	int   MaxType;
	int   PotType;
	int   tc;
	int   Deriv;
	int   BoundaryType;
	int   FuncCommand;
	char *tokstr = NULL;
	double x;
	double Value;


	/*  Command format '.. EVAL PAIR t1 t2 nderv x' */
	tokstr = strhed (&instr);
	if (!strcmpi("EVAL",tokstr))
      {
		/*  Get next token  */
      tokstr = strhed (&instr);

		/*  Determin which elements this pair potential is for  */
		ParsePot (&PotType, &tc, tokstr, &instr);
      Deriv = intstrf (&instr);
      x     = dblstrf (&instr) * 1e-8;
		if (POT_PAIR==PotType  &&  NULL!=PairFunction_m[tc])
			{
			if (Deriv==0)
				Value = F0(tc,x);
			else if (Deriv==1)
				Value = F1(tc,x);
			else
				Value = 0.0;
         }
		else
         Value = 0.0;
      printf ("        %g\n", Value);
		return PAIR_COMMAND_EVAL;
      }


   /*
	'Standard' potential command
	Command format .. PAIR t1 t2 (FUNC_command...)
	*/
	ParsePot (&PotType, &tc, tokstr, &instr);
	if (PotType==0)
		{
		ERROR_PREFIX
      printf ("Unrecognized potential type %s.\n", tokstr);
		CleanAfterError();
		/* 
		NOTE:  Next command not executed because previous commands
		halts execution.  But next command included to avoid compiler
		error.
		*/
		return PAIR_COMMAND_ERROR;
		}

   /*
   Pass function to FUNC_ routines to parse
	*/
	if (PotType==POT_PAIR)
		{

		/*  Set boundary value flags  */
		PairFunction_m[tc]->BoundaryType = 
         MAKE_BOUND_FLAG
				(
            FUNC_EXTRAP_CONST_SLOPE, FUNC_BOUND_CONST_SLOPE,
            FUNC_EXTRAP_ZERO,        FUNC_BOUND_ZERO
            );

      /*  Read in function   */
		FuncCommand = FUNC_Parse (PairFunction_m[tc], instr);

		switch (FuncCommand)
			{
			case FUNC_COMMAND_ERROR:
				ERROR_PREFIX
				printf ("Error returned by FUNC routine.\n");
				CleanAfterError();
				break;
			case FUNC_COMMAND_READ:
      		/*
				Scale function from units (ang,1/s->eunit) to (cm,ergs)
				*/
				FUNC_Scale (PairFunction_m[tc], 1e-8,         FUNC_INPUT );
				FUNC_Scale (PairFunction_m[tc], 1.0/s->eunit, FUNC_OUTPUT);
				break;
			case FUNC_COMMAND_SCALE:
				break;
			default:
				ERROR_PREFIX
				printf ("INTERNAL ERROR:  Unknown FUNC_COMMAND_ type.\n");
				CleanAfterError();
				break;
			}

      /*  Set cutoff arrays  */
      CalcCutoff();

		return PAIR_COMMAND_READ;
		}


	/*  Program should not reach here  */
	return PAIR_COMMAND_ERROR;
	}


/*
Release storage for pair potential
*/
void PAIR_Free (void)
	{
	int ipair;

	/*  Free previous allocations  */
	if (Cutoff_m   !=NULL) FREE(Cutoff_m   )
	if (CutoffSqr_m!=NULL) FREE(CutoffSqr_m)

	LOOP (ipair, MaxTypePair_m)
		{
		FUNC_Free (&(PairFunction_m[ipair]));
		}
	FREE(PairFunction_m)
	}


/*
Create storage for new pair potential
*/
void PAIR_Init (char *instr)
	{
	int ipair;
	int NewMaxType;

	NewMaxType = intstrf (&instr);

	/*  Test for valid value  */
	if (NewMaxType<=0)
		{
		ERROR_PREFIX
		printf ("Number of types must be greater than zero.\n");
		CleanAfterError();
		}

	/*  Free previous allocations  */
	PAIR_Free();

	/*  Calculate new Allocation  */
	MaxType_m     = NewMaxType;
	MaxTypePair_m = NewMaxType*(NewMaxType+1)/2;

	/*  Allocate storage  */
	ALLOCATE(Cutoff_m   ,    double,       MaxTypePair_m)
	ALLOCATE(CutoffSqr_m,    double,       MaxTypePair_m)
	ALLOCATE(PairFunction_m, Function_t *, MaxTypePair_m)

	/* 
	Create Null Functions  
	This prevents null pointer being called if a paricular pair
   of elements has no assigned potential between them.
	*/
	LOOP (ipair, MaxTypePair_m)
		{
		PairFunction_m[ipair] = FUNC_CreateNullFunction();
		}
	}


/*  Find maximum cutoff for pair potentials  */
double PAIR_GetMaxCutoff(void)
	{
	int ti;
	int tj;
	int tc;
	double TempCutoff = 0.0;
	double MaxCutoff  = 0.0;

	LOOP (ti, MaxType_m)
	LOOP (tj, ti+1     )
		{
		tc = COMBINED_TYPE(ti,tj);
      if (NULL!=PairFunction_m[tc])
         TempCutoff = PairFunction_m[tc]->ExteriorCutoff;
      if (TempCutoff>MaxCutoff)
         MaxCutoff=TempCutoff;
      }

   return MaxCutoff;
   }


/*
Return pointer to array of Function Pointers
so calling program can call Pair functions directly.
This is provided because the calling program
may way to use the PAIR_Parse() function to handle function
input, but not use the PAIR_CalcEnergy() and PAIR_CalcForce()
functions, but rather use the functions to calculate these
quanties.
*/
Function_t **PAIR_GetFuncListPtr(void)
	{
	return PairFunction_m;
	}

int PAIR_GetNumFunction(void)
	{
	return MaxTypePair_m;
	}



/*
************************************************************************
Local Functions
************************************************************************
*/



                           /*  LOCAL SUBROUTINES  */

void ParsePot (int *ptype, int *type, char *tokstr, char **instr)
	{
	int t1;
	int t2;

	if       (!strcmpi(tokstr, "PAIR"))
		{
      *ptype = POT_PAIR;
      t1 = intstrf (instr)-1;
		t2 = intstrf (instr)-1;
		if (!IS_TYPE_OK(t1) || !IS_TYPE_OK(t2))
			{
			printf ("ERROR (ParsePot):  Types out of range\n");
			CleanAfterError();
			}
		*type = COMBINED_TYPE(t1,t2);
		}

   else
		*ptype = POT_NONE;
   }



/*  Find Maximum cutoff for each pair of types  */
void CalcCutoff(void)
	{
   int ti;
   int tj;
   int tc;
   double TempCutoff;

	LOOP (ti, MaxType_m)
   LOOP (tj, ti+1     )
      {
      tc = COMBINED_TYPE(ti,tj);
      if (NULL!=PairFunction_m[tc])
         TempCutoff = PairFunction_m[tc]->ExteriorCutoff;
      Cutoff_m   [tc] = TempCutoff;
		CutoffSqr_m[tc] = TempCutoff*TempCutoff;
      }
   }



/*
Calculate separation vector if within cutoff.
Return TRUE if in cutoff,
Return FALSE otherwise
*/
static BOOLEAN CalcSeparation
(double *Coord1, double *Coord2, double *Vector, double Cutoff)
   {

   /*  Returns false if atoms out of range  */
   Vector[X] = GetSeparationDir(Coord2[X]-Coord1[X],X);
   if (Vector[X] >= Cutoff || Vector[X] <= -Cutoff)
      return FALSE;
   Vector[Y] = GetSeparationDir(Coord2[Y]-Coord1[Y],Y);
   if (Vector[Y] >= Cutoff || Vector[Y] <= -Cutoff)
      return FALSE;
   Vector[Z] = GetSeparationDir(Coord2[Z]-Coord1[Z],Z);
	if (Vector[Z] >= Cutoff || Vector[Z] <= -Cutoff)
      return FALSE;

   /*  Got thru three tests, atoms are within cutoff  */
   return TRUE;
   }



static double GetSeparationDir(double Separation, int idir)
   {
   if (BoundRep_m[idir])
      {
		if (Separation > HalfBox_m[idir])
         Separation -= Box_m[idir];
      else if (Separation < -HalfBox_m[idir])
         Separation += Box_m[idir];
      }
   return Separation;
   }



/*
************************************************************************

************************************************************************
*/

/*
Assembly language routines for Borland C on Intel 386 and above
*/
#ifdef ASM_CODE
#pragma inline
double __pr_poly (double x, int n, double *p)
{
   asm {
		.386
      /*  Load x  */
		fld qword ptr x;
      /*  Load n  */
      mov bx, n;
      shl bx, 3;
      /*  Load and multiply coeff  */
      push ds;
      lds si, p;
      fld qword ptr [si+bx];
	}
   loop:
   asm {
		.386
      or bx,bx;
		jz  done;
      fmul st[0],st[1];
      sub bx, 8;
      fadd qword ptr [si+bx];
      jmp loop;
   }
   done:
   asm {
		.386
      pop ds;
      ffree st[1];
	}
}


double __pr_fabs (double x)
{
   asm {
      .386
      /*  Load x  */
		fld qword ptr x;
		fabs;
	}
}
#endif