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


/*
************************************************************************
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;
            Eatom[jatom] += PairEnergy;

            LOOP (idir, NDIR)
               {
               ForceComponent[idir] =
                  PairForce * BondVector[NumBondNeighbors][idir];
               Force[iatom][idir] += ForceComponent[idir];
               Force[jatom][idir] -= ForceComponent[idir];
               }

            /*  Calculate system-wide stresses  */
            if (UseStress)
               {
               ForceComponent[X] *= BondRadius[NumBondNeighbors];
               ForceComponent[Y] *= BondRadius[NumBondNeighbors];
               ForceComponent[Z] *= BondRadius[NumBondNeighbors];
               a->Stress[XX] -=
                  ForceComponent[X] * BondVector[NumBondNeighbors][X];
               a->Stress[YY] -=
                  ForceComponent[Y] * BondVector[NumBondNeighbors][Y];
               a->Stress[ZZ] -=
                  ForceComponent[Z] * BondVector[NumBondNeighbors][Z];
               a->Stress[YZ] -=
                  ForceComponent[Y] * BondVector[NumBondNeighbors][Z];
               a->Stress[ZX] -=
                  ForceComponent[Z] * BondVector[NumBondNeighbors][X];
               a->Stress[XY] -=
                  ForceComponent[X] * BondVector[NumBondNeighbors][Y];
               }


            }
				/* Done with pair potential terms for atoms  i-j  */


         /*  Increment number of temporary neighbors  */
         NumBondNeighbors++;
         if (NumBondNeighbors>=MaxBondNeighbor_m)
            {
            printf ("ERROR (FILE %s Line %i): %s",
               __FILE__, __LINE__,
               "Too many temporary neighbors needed for Stillinger-Weber calculation.\n");
            printf ("Number of neighbors exceeds limit %i.\n",
               MaxBondNeighbor_m);
            printf ("Perhaps your atoms are scaled too small?\n");
				CleanAfterError();
            }
         }
         /*  Finished first pass at neighbors pairs with iatom  */



      /*
      Now, sum over bonds of iatom
      */

      /*  Sum bond interactions (i-j) and (i-k) for each bond  */
      AtomBondEnergy = 0.0;
      LOOP (jbond, NumBondNeighbors)
      LOOP (kbond, jbond           )
         {
         jatom = BondNeighborIndex[jbond];
         katom = BondNeighborIndex[kbond];

         /*
         Find cos() of angle between bonds i-j and i-k
         */
         CosTerm  = DOT(BondVector[jbond],BondVector[kbond]);
         CosTerm3 = CosTerm + 1./3.;


         /*  Calculate bond energy  */
         EnergyFac1 =
            AngleFac1_m *
            BondCutoffFactor[jbond] *
            BondCutoffFactor[kbond] *
            CosTerm3;
         BondEnergy      = EnergyFac1 * CosTerm3;
         AtomBondEnergy += BondEnergy;


         /*  Calculate bond force  */
         ForceFac1 = BondEnergy *  AngleFac2_m;
         ForceFac2 =  2 * EnergyFac1;


         /*  Forces due to atom j  */
         ForceFac2a  =  ForceFac2  / BondRadius[jbond];
         ForceFac1a  = -ForceFac1  / SQR(BondRadius[jbond]-Cutoff_m)
                       - ForceFac2a * CosTerm;

         LOOP (idir, NDIR)
            {
            ForceComponent[idir] =
               ForceFac1a * BondVector[jbond][idir] +
               ForceFac2a * BondVector[kbond][idir];
            Force[jatom][idir] -= ForceComponent[idir];
            Force[iatom][idir] += ForceComponent[idir];
            }


         /*  
			Contribution to system-stress from 
			force between atoms i and j
			*/
         if (UseStress)
            {
            ForceComponent[X] *= BondRadius[jbond];
            ForceComponent[Y] *= BondRadius[jbond];
            ForceComponent[Z] *= BondRadius[jbond];
            a->Stress[XX] -= ForceComponent[X] * BondVector[jbond][X];
            a->Stress[YY] -= ForceComponent[Y] * BondVector[jbond][Y];
            a->Stress[ZZ] -= ForceComponent[Z] * BondVector[jbond][Z];
            a->Stress[YZ] -= ForceComponent[Y] * BondVector[jbond][Z];
            a->Stress[ZX] -= ForceComponent[Z] * BondVector[jbond][X];
            a->Stress[XY] -= ForceComponent[X] * BondVector[jbond][Y];
            }


         /*  Forces due to atom k  */
         ForceFac2a  = ForceFac2  / BondRadius[kbond];
         ForceFac1a  = -ForceFac1  / SQR(BondRadius[kbond]-Cutoff_m)
                       - ForceFac2a * CosTerm;

         LOOP (idir, NDIR)
            {
            ForceComponent[idir] =
               ForceFac1a * BondVector[kbond][idir] +
               ForceFac2a * BondVector[jbond][idir];
            Force[katom][idir] -= ForceComponent[idir];
            Force[iatom][idir] += ForceComponent[idir];
            }

         /*  
			Contribution to system-stress from 
			force between atoms i and k
			*/
         if (UseStress)
            {
            ForceComponent[X] *= BondRadius[kbond];
            ForceComponent[Y] *= BondRadius[kbond];
            ForceComponent[Z] *= BondRadius[kbond];
            a->Stress[XX] -= ForceComponent[X] * BondVector[kbond][X];
            a->Stress[YY] -= ForceComponent[Y] * BondVector[kbond][Y];
            a->Stress[ZZ] -= ForceComponent[Z] * BondVector[kbond][Z];
            a->Stress[YZ] -= ForceComponent[Y] * BondVector[kbond][Z];
            a->Stress[ZX] -= ForceComponent[Z] * BondVector[kbond][X];
            a->Stress[XY] -= ForceComponent[X] * BondVector[kbond][Y];
            }
         }


      /*
      Sum bond energy to atom energy
      (factor of 2 to  correct for later  pair normalization)
      */
      Eatom[iatom] += AtomBondEnergy;
      }
      /*  End of loop over iatom  */

	/*  Include pair forces  */
	if (UsePairPotential_m)
		{
		PAIR_CalcForce(a);
		}


   /*  Normalize pair energy  */
   TotalEnergy = 0.0;
   LOOP (iatom, NumAtom)
      {
      TotalEnergy += Eatom[iatom];
      }

   a->epot = TotalEnergy;

	/*  Free storage  */
	FREE (BondVector)
	FREE (BondRadius)
	FREE (BondCutoffFactor)
	FREE (BondNeighborIndex)


   }





/*
************************************************************************
Local Functions
************************************************************************
*/
static void InitParameters (void)
   {

   PairFac1_m  = A_m*B_m;
   PairFac2_m  = A_m;
   PairFac3_m  = 1.0;
   Cutoff_m    = a_m;
   AngleFac1_m = lamda_m;
   AngleFac2_m = gamma_m;
   }



/*
Alter potential constants so that
for instance Scale=2 means cutoff is doubled
*/
static void ScaleDistance (double Scale)
   {
   PairFac1_m  *= pow(Scale, p_m);
   PairFac3_m  *= Scale;
   AngleFac2_m *= Scale;
   Cutoff_m    *= Scale;
   }


static void ScaleEnergy (double Scale)
   {
   PairFac1_m  *= Scale;
   PairFac2_m  *= Scale;
   AngleFac1_m *= Scale;
   }


static double GetEnergy (Particle_t *a)
   {
   int iatom;
   int jatom;
   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  *Eatom;
   double CutoffSqr;
   double Radius;
   double RadiusSqr;
   double PairTerm;
   double AtomBondEnergy;
   double CosTerm;
   double TotalEnergy;
   int    idir;
   int    NumAtom;

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

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

   /*  Initialize Eatom  */
   LOOP (iatom, NumAtom)
      Eatom[iatom] = 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)
      {

		/*  Skip atom if not 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;

         /*
         If we have gotten here, then we are including jneigh
         as a neighbor of iatom which forms a Tersoff 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:  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)
            {
            PairTerm = 0.5*
               (PairFac1_m/(RadiusSqr*RadiusSqr) - PairFac2_m)
               * exp(PairFac3_m/(Radius-Cutoff_m));
            Eatom[iatom] += PairTerm;
            Eatom[jatom] += PairTerm;
            }

         /*  Increment number of temporary neighbors  */
         NumBondNeighbors++;
         if (NumBondNeighbors>=MaxBondNeighbor_m)
            {
				ERROR_PREFIX
				printf ("Too many temporary neighbors needed for ");
				printf ("Stillinger-Weber calculation.\n");
            printf ("Number of neighbors exceeds limit %i.\n",
               MaxBondNeighbor_m);
            printf ("Perhaps your atoms are scaled too small?\n");
				CleanAfterError();
            }
         }
         /*  Finished first pass at neighbors pairs with iatom  */



      /*
      Now, prepare bond information for sum over bonds of iatom
      */

      /*  Sum bond interactions (i-j) and (i-k) for each bond  */
      AtomBondEnergy = 0.0;
      LOOP (jbond, NumBondNeighbors)
      LOOP (kbond, jbond           )
         {

         /*
         Find cos() of angle between bonds i-j and i-k
         */
         CosTerm =
            DOT(BondVector[jbond],BondVector[kbond]) /
            (BondRadius[jbond]*BondRadius[kbond]);

         AtomBondEnergy +=
            AngleFac1_m  *
            BondCutoffFactor[jbond] *
            BondCutoffFactor[kbond] *
            SQR(CosTerm+1./3.);


         }


      /*
      Sum bond energy to atom energy
      (factor of 2 to  correct for later  pair normalization)
      */
      Eatom[iatom] += AtomBondEnergy;
      }
      /*  End of loop over iatom  */


	/*  Call pair energy routine  */
	if (UsePairPotential_m)
		{
		PAIR_CalcEnergy(a);
		}

   /*  Normalize pair energy  */
   TotalEnergy = 0.0;
   LOOP (iatom, NumAtom)
      {
      TotalEnergy += Eatom[iatom];
      }

	/*  Free storage  */
	FREE (BondVector)
	FREE (BondRadius)
	FREE (BondCutoffFactor)
	FREE (BondNeighborIndex)


   return TotalEnergy;
   }




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