cdstil.c

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

/*
************************************************************************
History
************************************************************************
*/
/*
10 Feb 1998
Stillinger-Weber Si Potential
Converted from cdcsi.c Tersoff potential
*/






/*
************************************************************************
Include Files
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "strsub.h"
#include "parse.h"
#include "iomngr.h"
#include "cdhouse.h"
#include "cdsubs.h"
#include "cdsearch.h"
#include "cdpairf.h"


/*
************************************************************************
Compiler Switches
************************************************************************
*/
/*  COMPILER OPTION FOR SPEED  */
#ifdef __TURBOC__
#pragma option -G
#endif

#define PAIR_ONLY
#undef  PAIR_ONLY

#define CORRECTION
#undef CORRECTION

#define TEST_ZETA
#undef  TEST_ZETA

/*
************************************************************************
Defines
************************************************************************
*/
/*  DEFINES THAT SHOULD HAVE COME FROM math.h  */
#ifndef M_PI
#define M_PI       3.1415926535897931160E0  /*Hex  2^ 1 * 1.921FB54442D18 */
#endif


/*  Maximum number of neighbors for any one atom  */
#define INIT_MAX_BOND_NEIGHBOR 32

/*  Size of input string buffer  */
#define NBUFF 256





/*
************************************************************************
Macros
************************************************************************
*/
#define MAG_SQR(V)\
   ((V)[X]*(V)[X] + (V)[Y]*(V)[Y] + (V)[Z]*(V)[Z])


#define DOT(A,B) \
   ((A)[X]*(B)[X] + (A)[Y]*(B)[Y] + (A)[Z]*(B)[Z])


/*
************************************************************************
Global Variables
************************************************************************
*/
extern Simulation_t *s;
extern LIST         *inlist;
extern FILE         *out;




/*
************************************************************************
Module-Scope Variables
************************************************************************
*/
/*  Potential parameters from Stillinger  */
static double A_m       =  7.049556277;
static double B_m       =  0.6022245584;
static double p_m       =  4;
static double a_m       =  1.80;
static double lamda_m   = 21.0;
static double gamma_m   =  1.20;
static double sigma_m   =  2.0951e-8;    /* cm  */
static double epsilon_m =  3.4723e-12;   /* erg/atom-pair*/

/*  Values obtained from above parameters  */
static double Cutoff_m;
static double PairFac1_m;
static double PairFac2_m;
static double PairFac3_m;
static double AngleFac1_m;
static double AngleFac2_m;

static double  HalfBox_m   [NDIR];
static double  Box_m       [NDIR];
static BOOLEAN BoundRep_m  [NDIR];

static BOOLEAN UsePairPotential_m = FALSE;

static int     MaxBondNeighbor_m = INIT_MAX_BOND_NEIGHBOR;


/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
static void   InitParameters(void);
static void   ScaleDistance(double);
static void   ScaleEnergy (double);
static double GetEnergy (Particle_t *a);

static BOOLEAN CalcSeparation
(double *Coord1, double *Coord2, double *Vector, double Cutoff);
static double GetSeparationDir(double Separation, int idir);

/*
************************************************************************
Exported Functions
************************************************************************
*/
#pragma argsused
void stil_energy_list  (Particle_t *a, Simulation_t *s, int UseSelect)
   {
   a->epot   = GetEnergy (a);
   }


#pragma argsused
void STILL_Parse (char *instr)
   {
	int    FUNC_ScaleType;
	int    ifunc;
	int    nfunc;
	int    ErrorCode;
	char  *tokstr = NULL;
   double Scale;
	Function_t **FuncListPtr = NULL;

	if      (IsFirstToken("MAXBONDCNT",instr))
		{
		
		/*  Remove leading MAXBONDCNT token  */
		strhed (&instr);

		/*  Get next number  */
		MaxBondNeighbor_m = intstrf (&instr);

		/*  Sanity check  */
		if (MaxBondNeighbor_m<4)
			{
			ERROR_PREFIX
			printf ("Maximum number of bond neighbors must be 4 or more.\n");
			CleanAfterError();
			}
		}

   /*  Scale potential energy or distance  */
	else if (IsFirstToken("SCALE",instr))
      {

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

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

      if (!strcmpi(tokstr, "DISTANCE"))
      	{   
         Scale = dblstrf (&instr);
         ScaleDistance (Scale);
			FUNC_ScaleType = FUNC_INPUT;
         }

      else if (!strcmpi(tokstr, "ENERGY"))
         {
         Scale = dblstrf (&instr);
         ScaleEnergy (Scale);
			FUNC_ScaleType = FUNC_OUTPUT;
         }

		/*  Unknown SCALE option, return with Warning  */
      else
         {
         printf("WARNING:  Unknown option to Still-Weber Si");
			printf(" potential SCALE command (%s).\n", tokstr);
         IncrementNumberOfWarnings();
         return;
         }
			/* unknown SCALE option  */

		/*  Apply scale to auxillary Pair potentials  */
		if (UsePairPotential_m)
			{
			FuncListPtr = PAIR_GetFuncListPtr();
			nfunc = PAIR_GetNumFunction();
			LOOP (ifunc, nfunc)
				{
				FUNC_Scale (FuncListPtr[ifunc], FUNC_ScaleType, Scale);
				}
			}
			return;
      }  
		/*  SCALE potion  */


	/*  Pass option onto PAIR routine  */
	else 
		{
		if (UsePairPotential_m)
			{
			ErrorCode = PAIR_Parse (instr);
			}

		if (!UsePairPotential_m || ErrorCode==PAIR_COMMAND_ERROR)
			{
	      printf ("WARNING:  Unknown option to Still-Weber Si potential(%s).\n",
   	           tokstr);
     		IncrementNumberOfWarnings();
			}
		}

   }



/*  Initialze potential  */
void STILL_Init(char *instr)
	{
	char *tokstr = NULL;

   /*  Fill out paramters tables  */
   InitParameters();

   /*  Scale distance from angstroms to cm  */
   ScaleDistance (sigma_m);

   /*  Scale energy from eV to ergs  */
   ScaleEnergy (epsilon_m);

   /*  Check for auxillary pair potentials  */
	tokstr = strhed (&instr);
	if (!strcmpi("PAIR",tokstr))
		{
		UsePairPotential_m = TRUE;
		PAIR_Init (instr);			
		}

      /*  Print info  */
   printf ("\n");
   printf ("   Stillinger-Weber silicon model\n");
   printf ("   reference:\n");
   printf ("   \"Computer simulation of local order in");
	printf (" condensed phases of silicon\".\n");
   printf ("   F. H. Stllinger, T. A. Weber,");
	printf (" Phys. Rev. B 31 (8) 5262 (1985)\n");
   printf ("\n");
	}



/*
Force routine called externally
*/
#pragma argsused
void stil_calcforc (Particle_t *a, Simulation_t *s)
   {
   BOOLEAN UseStress;
   int iatom;
   int jatom;
   int katom;
   int jbond;
   int kbond;
   int *Neighbors = NULL;
   int    NumNeighbors;
   int    NumBondNeighbors;
	Coord_t *BondVector        = NULL;
   double  *BondRadius        = NULL;
   double  *BondCutoffFactor  = NULL;
   int     *BondNeighborIndex = NULL;

   double (*Coord)[NDIR];
   double (*Force)[NDIR];
   double  *Eatom;
   double CutoffSqr;
   double Radius;
   double RadiusSqr;
   double InvRadius;
   double PairEnergy;
   double PairTerm1;
   double AtomBondEnergy;
   double CosTerm;
   double CosTerm3;
   double Radius4;
   double PairCutoff;
   double EnergyFac1;
   double BondEnergy;
   double ForceFac1;
   double ForceFac2;
   double ForceFac1a;
   double ForceFac2a;
   double PairForce;
   double ForceComponent[NDIR];
   double TotalEnergy;
   int    idir;
   int    NumAtom;

   /*  Test for atoms  */
   NumAtom = a->np;
   if (NumAtom==0)
      return;

	/*  Allocate storage  */
	ALLOCATE (BondVector,        Coord_t, MaxBondNeighbor_m)
	ALLOCATE (BondRadius,        double,  MaxBondNeighbor_m)
	ALLOCATE (BondCutoffFactor,  double,  MaxBondNeighbor_m)
	ALLOCATE (BondNeighborIndex, int,     MaxBondNeighbor_m)

   /*  Intialize pointers  */
   Coord = (double (*)[NDIR]) a->cur;
   Force = (double (*)[NDIR]) a->f;
   Eatom =                    a->eatom;

   /*  Set up box  */
   LOOP (idir,NDIR)
      {
      Box_m     [idir] = a->bcur[idir];
      HalfBox_m [idir] = 0.5 * Box_m[idir];
      BoundRep_m[idir] = !a->surf[idir];
      }


   /*  Set Eatom[] and Force[] to 0  */
   LOOP (iatom, NumAtom)
      {
      Eatom[iatom] = 0.0;
      Force[iatom][X] = 0.0;
      Force[iatom][Y] = 0.0;
      Force[iatom][Z] = 0.0;
      }

   /*  Initialize stresses to zero  */
   UseStress =
      (s->StressStepOutput.Save || a->BoxMotionAlgorithm!=BMA_NONE);
   if (UseStress)
      {
      a->Stress[XX] = 0.0;
      a->Stress[YY] = 0.0;
      a->Stress[ZZ] = 0.0;
      a->Stress[YZ] = 0.0;
      a->Stress[ZX] = 0.0;
      a->Stress[XY] = 0.0;
      }

   /*  Set cutoff for call to search_all()  */
   s->cutoff = Cutoff_m;

	/*  Include pair cutoff if using auxillary pair routines  */
	if (UsePairPotential_m && Cutoff_m<PAIR_GetMaxCutoff())
		s->cutoff = PAIR_GetMaxCutoff();

   /*  Get neighbors (get all neighbors for each atom)  */
   a->CalcUniqueNeighbors = FALSE;
   search_all (a);

   /*  Calculate energy for every atom  */
   CutoffSqr = Cutoff_m * Cutoff_m;
   LOOP (iatom, NumAtom)
      {

		/*  Use only Silicon atoms (type 1)  */
		if (a->type[iatom]!=0)
			continue;

      /*  Neighbors[] points to start of neighbors of iatom  */
      Neighbors    = &(a->NeighborList      [a->NeighborListIndex[iatom]]);
      NumNeighbors =   a->NeighborListLength[                     iatom ];


      /*  Store information about bond neighbors  */
      NumBondNeighbors = 0;
      LOOP (jbond, NumNeighbors)
         {

         /*  Get indice of neighbor atom  */
         jatom = Neighbors[jbond];

			/*  Skip atom if not type 1  */
			if (a->type[jatom]!=0)
				continue;

         /*  Cannot use same atom twice  */
         ASSERT (iatom!=jatom)


         if
         (
         !CalcSeparation
            (
            Coord[iatom],
            Coord[jatom],
            BondVector[NumBondNeighbors],
            Cutoff_m
            )
         )
            continue;

         /*  Calculate BondRadius  */
         RadiusSqr =
            MAG_SQR(BondVector[NumBondNeighbors]);

         /*  Test cutoff  */
         if (RadiusSqr > CutoffSqr)
            continue;

         /*  We will need the radius, so calculate it now  */
         Radius = sqrt(RadiusSqr);
         BondRadius[NumBondNeighbors] = Radius;

         /*  Convert BondVector[] from x,y,z to x/r,y/r,z/r  */
         InvRadius = 1.0 / Radius;
         BondVector[NumBondNeighbors][X] *= InvRadius;
         BondVector[NumBondNeighbors][Y] *= InvRadius;
         BondVector[NumBondNeighbors][Z] *= InvRadius;

         /*
         If we have gotten here, then we are including jneigh
         as a neighbor of iatom which forms a bond
         */

         /*  Store index of this neighbor for Tersoff calculation  */
         BondNeighborIndex[NumBondNeighbors] = jatom;


         /*  Test for particle pair with zero separation  */
         if (Radius==0.0)
            {
				printf ("ERROR: ");
            printf ("Particles %i and %i have the exact same position.\n",
               iatom+1, jatom+1);
            CleanAfterError();
            }

         /*  Save value of cutoff function between atom i and j  */
         BondCutoffFactor[NumBondNeighbors] =
            exp(AngleFac2_m/(Radius-Cutoff_m));

         /*  Sum pair potential term  */
         if (iatom<jatom)
            {

            /*  Temporary variables  */
            Radius4    = RadiusSqr*RadiusSqr;
            InvRadius  = 1.0/(Radius-Cutoff_m);
            PairTerm1  = PairFac1_m / Radius4;
            PairCutoff = exp(PairFac3_m*InvRadius);

            /*  Potential Energy  */
            PairEnergy = 
               (PairTerm1 - PairFac2_m) * PairCutoff;

            /*  Forces  */
            PairForce =
               -PairEnergy * PairFac3_m / SQR(Radius-Cutoff_m) -
               4 * PairTerm1 * PairCutoff / Radius;

				/*  Sum pair energy  */
				PairEnergy *= 0.5;
            Eatom[iatom] += PairEnergy;