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


/*
************************************************************************
Comments
************************************************************************
*/

/*
TERSOFF Si-C POTENTIAL
*/

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

/*

6 Nov 1996 - Got Forces to run correctly.  Problem was caused by
             macro SQR() in cdhouse.h.  It was defined as (X)*(X),
             should have been ((X)*(X)).  The old definition was
             a problem when calculating GtermD in csi_calcforc().
             The statement  ?/SQR(GtermDeninominator) was not doing
             what had been intended.
           - There is another problem with forces.  The pure pair
             forces are fine, but an examination of the sum of KE
             and PE under CLAMP OFF shows that in order for sum to
             be constant, fluctuations in KE need to be reduced by
             half.  So those forces due to the angular terms are
             reduced by a factor of 2.
22 Nov 1996
   VERSION 2.2.1
          - Previos problems were fixed a week ago.
   VERSION 2.2.2
          - Add stress calculation to Tersoff force function

16 Jul 1997
   Version 2.2.5
          - Add separate pair potential
          - Add routine CalcTotalCutoff() to calculate the maximum
            cutoff per species for both Tersoff and auxillary pair potential
          - command  POTENTIAL PAIR ncoeff cutoff
                     a1 a2 a3 a4
            resulting energy P(r) = a1*(r-rc)^2 + a2*(r-rc)^3 + ..
          - set s->cutoff=GetMaxCutoff() before each call  to
            search_all()

22 Aug 1997
   Version 2.4.8
          - When implementing separate pair potential above, forgot
            to initilize AtomBondEnergy=0.0 before calcualting each
            atom


*/




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


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


/*  Number of Tersoff atom species:  Carbon and Silicon  */
#define NSPECIES 2

/*
Maximum number of atom types -
types other than 1 and 2 are pair only
*/
#define MAX_TYPE 5
#define MAX_PAIR MAX_TYPE*(MAX_TYPE+1)/2


/*  Maximum number of neighbors for any one atom  */
#define MAX_BOND_NEIGHBOR 24

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





/*
************************************************************************
Macros
************************************************************************
*/



#define ARITH_AVERAGE(A,I,J) \
   A[I][J] = 0.5* (A[I][I] + A[J][J]);

#define GEOM_AVERAGE(A,I,J) \
   A[I][J] = sqrt (A[I][I] * A[J][J]);

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

#define MAG_SQR(A) \
   ((A)[X]*(A)[X] + (A)[Y]*(A)[Y] + (A)[Z]*(A)[Z])

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


#define FA(I,J,R,F,F1) \
   (F ) = -B[I][J] * exp(-Mu[I][J]*(R)); \
   (F1) = -Lamda[I][J] * (F);

#define FR(I,J,R,F,F1) \
   (F ) = A[I][J] * exp(-Lamda[I][J]*(R)); \
   (F1) = -Lamda[I][J] * (F);

#define IS_TERSOFF(I,J) \
   ((I)<NSPECIES && (J)<NSPECIES)


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

#ifdef TEST_ZETA
   extern double MinZeta_g;
   extern int    MinStep_g;
   extern int    MinI_g;
   extern int    MinJ_g;
#endif



/*
************************************************************************
Module-Scope Variables
************************************************************************
*/
static double  HalfBox_m  [NDIR];
static double  Box_m      [NDIR];
static BOOLEAN BoundRep_m [NDIR];


/*
This parameters depend on species of both atoms in pair
The cross terms must be set by the initialization routine
*/


/*  Cross terms are arithmetic average of pure terms  */
static double Lamda            [NSPECIES][NSPECIES] = {{3.4879, 0},{0,        2.4799}};
static double Mu               [NSPECIES][NSPECIES] = {{2.2119, 0},{0,        1.7322}};

/*  Cross terms are geometric  average of pure terms  */
static double A                [NSPECIES][NSPECIES] = {{ 1.3936e3, 0},{0,  1.8308e3}};
static double B                [NSPECIES][NSPECIES] = {{ 3.467e2,  0},{0,  4.7118e2}};
static double InteriorCutoff   [NSPECIES][NSPECIES] = {{1.8, 0},{0,  2.7}};
static double ExteriorCutoff   [NSPECIES][NSPECIES] = {{2.1, 0},{0,  3.0}};
static double CutoffInterval   [NSPECIES][NSPECIES] = {{0.3, 0},{0,  0.3}};
static double ExteriorCutoffSqr[NSPECIES][NSPECIES];

/*  This paramter cross terms not an average of pure terms  */
static double Chi              [NSPECIES][NSPECIES] = {{1.0, 0.9776}, {0.9776, 1.0}};

/*  These parameters depend on upon species of one atom  */
static double Beta             [NSPECIES]           = {1.5724e-7,    1.1000e-6};
static double n                [NSPECIES]           = {7.2751e-1,    7.8734e-1};
static double c                [NSPECIES]           = {3.8049e4,     1.0039e5};
static double d                [NSPECIES]           = {4.384,        1.6217e1};
static double h                [NSPECIES]           = {-5.7058e-1,  -5.9825e-1};

/*  The values of these parameters depend on previous ones  */
static double c2  [NSPECIES];
static double d2  [NSPECIES];
static double cd2 [NSPECIES];


/*  These values control the auxiallary pair potential  */
static BOOLEAN UsePairPotential_m = FALSE;
static BOOLEAN IsPairPotentialInitialized_m = FALSE;
static int     NumPairCoeff_m[MAX_PAIR];
static double  PairCutoff_m  [MAX_PAIR];
static double *PairCoeff_m   [MAX_PAIR];

/*  These array hold temporary values for Tersoff calculation  */
static double BondZetaD_m [MAX_BOND_NEIGHBOR] [MAX_BOND_NEIGHBOR] [NDIR];
static double BondVector_m[MAX_BOND_NEIGHBOR] [NDIR];
static double TotalCutoff_m    [MAX_PAIR] [MAX_PAIR];
static double TotalCutoffSqr_m [MAX_PAIR] [MAX_PAIR];


/*
************************************************************************
Local Function Prototypes
************************************************************************
*/

static void    InitParameters   (void);
static void    InitPairPotential(void);
static void    InitializeCutoff (void);
static double  GetEnergy        (Particle_t *);
static void    ScaleDistance    (double);
static void    ScaleEnergy      (double);
static double  GetMaxCutoff     (void);
static BOOLEAN CalcSeparation   (double *, double *, double *, double);
static double  GetSeparationDir (double, int);

static double  fr (int,int,double);
static double  fa (int,int,double);
static double  fc (int,int,double);
static double  fr1(int,int,double);
static double  fa1(int,int,double);
static double  fc1(int,int,double);
static double  PairPotentialEnergy (int, int, double);
static double  PairPotentialForce  (int, int, double);

void CalcTotalCutoff
   (double [MAX_PAIR][MAX_PAIR], double [MAX_PAIR][MAX_PAIR]);


double ***Allocate3D(int, int, int);
void      Free3D    (double ***, int, int);
double  **Allocate2D(int, int);
void      Free2D    (double ** , int);

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

#pragma argsused
void csi_read_potential (char *tokstr, char *instr)
   {
   double Scale;
   int    Type;
   int    iread;
   BOOLEAN EndOfFile;
   double InputInteriorCutoff;
   double InputExteriorCutoff;
   char   InputBuffer[NBUFF];
   int    Type1;
   int    Type2;
   int    CombinedType;


   /*  Initialize potential  */
   if (!strcmpi(tokstr, "INIT"))
      {

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

      /*  Initialize Pair Potential  */
      InitPairPotential();

      /*  Scale distance from angstroms to cm  */
      ScaleDistance (1.0e-8);

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

      /*  Print info  */
      printf ("\n");
      printf ("   Tersoff's Silicon-Carbon potential\n");
      printf ("   reference:\n");
      printf ("   \"Modeling solid-state chemistry:  Interatomic potentials\n");
      printf ("   for multicomponent systems\",\n");
      printf ("   J. Tersoff, Phys Rev. B 39 (8), 5566-5568 (15 March 1989)\n");
#ifdef PAIR_ONLY
      printf ("\nWARNING:  ");
      printf ("This version of potential uses PAIR POTENTIAL ONLY.\n");
      printf ("No angular forces are used.\n");
#endif
      printf ("\n");

      return;
      }


   /*  Scale potential energy or distance  */
   if (!strcmpi(tokstr, "SCALE"))
      {

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

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

      else
         {
         printf ("WARNING:  Unknown option to Tersoff C-Si potential SCALE command (%s).\n",
                 tokstr);
         IncrementNumberOfWarnings();
         return;
         }
      }



   /*  Set potential cutoff  */
   if (!strcmpi(tokstr, "CUTOFF"))
      {
      Type = intstrf (&instr);
      if (Type!=1 && Type!=2)
         {
         printf ("WARNING:  Type is out of range (1,2) for Tersoff C-Si Cutoff command.\n");
         IncrementNumberOfWarnings();
         return;
         }

      InputInteriorCutoff = dblstrf (&instr) * 1e-8;
      InputExteriorCutoff = dblstrf (&instr) * 1e-8;
      if (InputInteriorCutoff>InputExteriorCutoff)
         {
         printf ("WARNING:  Interior Cutoff exceeds Exterior Cutoff.  Will ignore CUTOFF command.\n");
         IncrementNumberOfWarnings();
         return;
         }

      /*  Set up cutoffs  */
      Type--;
      InteriorCutoff[Type][Type] = InputInteriorCutoff;
      ExteriorCutoff[Type][Type] = InputExteriorCutoff;

      /*  Initialize cutoff  */
      InitializeCutoff();

      /*  return  */
      return;
      }

   /*  Use Auxillary Pair Potential  */
   if (!strcmpi(tokstr, "PAIR"))
      {
      tokstr = strhed (&instr);

      /*  Turn auxillary pair potential off  */
      if (!strcmpi(tokstr, "OFF"))
         {
         UsePairPotential_m = FALSE;
         }

      else
         {

         /*  Skip ON token  */
         if (!strcmpi(tokstr,"ON"))
            {
            tokstr= strhed (&instr);
            UsePairPotential_m = TRUE;
            }

         /*  Read parameters if not blank  */
         if (!IsBlank(tokstr))
            {

            /*  Read types  */
            Type1 = intstrf(&tokstr)-1;
            Type2 = intstrf(&instr )-1;

            /*  Check types  */
            if (Type1<0 || Type1>=MAX_TYPE)
               {
               printf ("ERROR:  Type 1 (%i) is out of range [1..%i]\n",
                  Type1+1, MAX_TYPE);
               CleanAfterError();
               }

            if (Type2<0 || Type2>=MAX_TYPE)
               {
               printf ("ERROR:  Type 2 (%i) is out of range [1..%i]\n",
                  Type2+1, MAX_TYPE);
               CleanAfterError();
               }

            CombinedType = COMBINED_TYPE(Type1,Type2);
            NumPairCoeff_m[CombinedType] =      dblstrf (&instr);
            PairCutoff_m  [CombinedType] = 1e-8*dblstrf (&instr);

            /*  Allocate coefficients  */
            FREE (PairCoeff_m[CombinedType])
            ALLOCATE
               (
               PairCoeff_m[CombinedType],
               double,
               NumPairCoeff_m[CombinedType]
               )

            iread = 0;
            EndOfFile      = FALSE;
            while (iread<NumPairCoeff_m[CombinedType] && !EndOfFile)
               {

               /*  Read new line if current line is blank  */
               if (IsBlank(instr))
                  {
                  EndOfFile =
                     (NULL==m_gets_list(InputBuffer, NBUFF, inlist));
                  instr = InputBuffer;
                  }

               /*  Premature end-of-file, end program  */
               if (EndOfFile)
                  {
                  printf ("ERROR: File ends before coefficients were found\n");
                  CleanAfterError();
                  }

               PairCoeff_m[CombinedType][iread] =
                  dblstrf (&instr)/s->eunit;
               iread++;
               }

            UsePairPotential_m = TRUE;
            IsPairPotentialInitialized_m = TRUE;
            }

         /*  Make sure that auxillary pair potential initialized  */
         if (!IsPairPotentialInitialized_m && UsePairPotential_m)
            {
            printf
               ("%s %s\n",
               "ERROR:  Cannot use C-Si auxillary pair potential",
               "because the necessary parameters have not been read."
               );
            CleanAfterError();
            }
         }
      }

   else
      {
      printf ("WARNING:  Unknown option to Tersoff Si-C potential(%s).\n",
              tokstr);
      IncrementNumberOfWarnings();
      return;
      }
   }



/*
Force routine called externally
*/
void csi_calcforc (Particle_t *a, Simulation_t *s)
   {
   int idir;
   int iatom;
   int jatom;
   int katom;
   int jneigh;
   int jbond;
   int kbond;
   int TypeI;
   int TypeJ;
   int *Neighbors = NULL;
   int NumNeighbors;
   double TotalBondEnergy;
   double TotalPairEnergy;
   double AtomBondEnergy;
   int    NumBondNeighbors;
   double BondRadius        [MAX_BOND_NEIGHBOR];
   double BondRadiusSqr     [MAX_BOND_NEIGHBOR];
   double BondZeta          [MAX_BOND_NEIGHBOR];
   double Bondfc            [MAX_BOND_NEIGHBOR];
   double Bondfc1           [MAX_BOND_NEIGHBOR];
   int    BondNeighborIndex [MAX_BOND_NEIGHBOR];
   double BondStrength;
   double BondStrengthD;
   double DotNormalization;
   double Gterm;
   double GtermD;
   double GtermDenominator;
   double CosTerm;
   double DCosJ;
   double DCosK;
   double Bondfa;
   double Bondfr;
   double Bondfa1;
   double Bondfr1;
   double BetaZetaPowN;
   double ForceX;
   double ForceY;
   double ForceZ;
   double ForceR;
   double TempFactor;
   double (*Coord)[NDIR];
   double (*Force)[NDIR];
   BYTE    *Type;
   int    NumAtom;
   int    CombinedType;


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

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

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

   /*  Setup cutoff info  */
   CalcTotalCutoff (TotalCutoff_m, TotalCutoffSqr_m);

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

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


   /*  Intialize all forces to zero  */
   LOOP (iatom, NumAtom)
   LOOP (idir,  NDIR)
      {
      Force[iatom][idir] = 0.0;
      }

   /*  Initialize stresses to zero  */
   if (s->StressStepOutput.Save || a->BoxMotionAlgorithm!=BMA_NONE)
      {
      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;
      }


   /*
   *******************************
   Calculate energy for every atom
   *******************************
   */
   TotalPairEnergy = 0.0;
   TotalBondEnergy = 0.0;
   LOOP (iatom, NumAtom)
      {

      /*
      Store information about bond neighbors
      */
      NumBondNeighbors = 0;


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

      LOOP (jneigh, NumNeighbors)
         {

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

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

         /*  Get atom types  */
         TypeI = Type[iatom];
         TypeJ = Type[jatom];

          if
         (
         !CalcSeparation
            (
            Coord[iatom],
            Coord[jatom],
            BondVector_m[NumBondNeighbors],
            TotalCutoff_m[TypeI][TypeJ]
            )
         )
            continue;

         /*  Calculate BondRadius  */
         BondRadiusSqr[NumBondNeighbors] =
            MAG_SQR(BondVector_m[NumBondNeighbors]);

         /*  Test cutoff  */
         if (
            BondRadiusSqr[NumBondNeighbors] >=
            TotalCutoffSqr_m[TypeI][TypeJ]
            )
            continue;

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

         /*  Calculate pair potential  */
         if ( UsePairPotential_m )
            {

            /*  Calculate combined type  */
            CombinedType = COMBINED_TYPE(TypeI,TypeJ);

            /*  Test if particle within pair potential range  */
            if (BondRadius[NumBondNeighbors]<PairCutoff_m[CombinedType])

            /*  Calculate only iatom<jatom (to prevent double counting)  */
            if (iatom<jatom)

               {
               /*  Calculate energy  */
               TotalPairEnergy +=
                  PairPotentialEnergy(TypeI,TypeJ,BondRadius[NumBondNeighbors]);

               /*  Force along radius  */
               ForceR  =
                  PairPotentialForce(TypeI,TypeJ,BondRadius[NumBondNeighbors]);

               /*  Force in X, Y, Z directions  */
               ForceR /= BondRadius[NumBondNeighbors];
               ForceX = ForceR*BondVector_m[NumBondNeighbors][X];
               ForceY = ForceR*BondVector_m[NumBondNeighbors][Y];
               ForceZ = ForceR*BondVector_m[NumBondNeighbors][Z];

               /*  Force on atoms  */
               /*
               NOTE ON SIGNS:
                Theory:
                  Force[iatom][X] = - dE(r_ij)/dx_i
                                  = - dE(r_ij)/dr_ij (x_ij/r_ij)
                Implementation:
                   PairPotendialForce() returns  -dE(r)/dr
                   BondVector_m[][X] contains x_ji = (x_j - x_i) = - x_ij

                Thus Force[iatom][X] has - sign in summation below, and
                  Force[jatom] has + sign

               */
               Force[jatom][X] += ForceX;
               Force[jatom][Y] += ForceY;
               Force[jatom][Z] += ForceZ;

               Force[iatom][X] -= ForceX;
               Force[iatom][Y] -= ForceY;
               Force[iatom][Z] -= ForceZ;

               }

            }
            /*  Finished pair interaction with neighbors of iatom  */

         /*
         If both particles are not either Si or C,
         then skip to next pair of iatom,jneigh
         */
         if (!IS_TERSOFF(TypeI,TypeJ))
            continue;

         /*
         If particle outside Tersoff range
         then skip to next pair of iatom, jneigh
         */
         if ( BondRadius[NumBondNeighbors] >= ExteriorCutoff[TypeI][TypeJ] )
            continue;


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

         /*  Save value of cutoff function between atom i and j  */
         Bondfc[NumBondNeighbors] =
            fc(TypeI,TypeJ,BondRadius[NumBondNeighbors]);

         Bondfc1[NumBondNeighbors] =
            fc1(TypeI, TypeJ, BondRadius[NumBondNeighbors]);


         /*  Increment number of temporary neighbors  */
         NumBondNeighbors++;
         if (NumBondNeighbors>=MAX_BOND_NEIGHBOR)
            {
            printf ("ERROR (FILE %s Line %i): %s",
               __FILE__, __LINE__,
               "Too many temporary neighbors needed for C-Si calculation.\n");
            printf ("Number of neighbors exceeds limit %i.\n",
               MAX_BOND_NEIGHBOR);
            printf ("Perhaps your atoms are scaled too small?\n");
            exit   (1);
            }

         }
         /*  Finished first pass at neighbors pairs with iatom  */


      /*  Skip Tersoff calculation if not C or Si  */
      if (!IS_TERSOFF(TypeI,TypeI))
         continue;
      ASSERT(TypeI<NSPECIES)


      /*
      Following calculations on iatom are
      exclusively for Tersoff potential
      */

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

      /*  Initialize Zeta sums to zero  */
      LOOP (jbond, NumBondNeighbors)
         {
         BondZeta[jbond] = 0.0;
         }

      LOOP (jbond, NumBondNeighbors)
      LOOP (kbond, NumBondNeighbors)
      LOOP (idir,  NDIR)
         {
         BondZetaD_m[jbond][kbond][idir] = 0.0;
         }

      /*
      Sum over bonds of iatom
      */
      LOOP (jbond, NumBondNeighbors)
      LOOP (kbond, jbond           )
         {

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


         GtermDenominator = d2[TypeI] + SQR(h[TypeI] - CosTerm);
         Gterm   =
            1 + cd2[TypeI] -
            c2 [TypeI] / GtermDenominator;


         /*  Derivative of Gterm with respect to CosTerm  */
         GtermD =
            -2.0 * c2[TypeI] * (h[TypeI] - CosTerm)
            / SQR(GtermDenominator);


         /*  Contribution to Zeta for this bond  */
         BondZeta[jbond] += Bondfc[kbond] * Gterm;
         BondZeta[kbond] += Bondfc[jbond] * Gterm;


         /*  Calculate force terms  */
         LOOP (idir, NDIR)
            {

            /*  Derivative of Cosine term between j,k wrt j  */
            DCosJ =
               BondVector_m[kbond][idir]*DotNormalization -
               BondVector_m[jbond][idir]*CosTerm/BondRadiusSqr[jbond];

            /*
            Derivative of Cosine term between j,k wrt k
            (NOTE: dCosTerm/dx.k = DCosK * x.k )
            */
            DCosK =
               BondVector_m[jbond][idir]*DotNormalization -
               BondVector_m[kbond][idir]*CosTerm/BondRadiusSqr[kbond];


            /*  Derivative wrt to atoms j,i of Zeta for bond i-j  */
            BondZetaD_m[jbond][jbond][idir] +=
               Bondfc[kbond]*GtermD*DCosJ;


            /*  Derivative wrt to atoms k,i of Zeta for bond i-k  */
            BondZetaD_m[kbond][kbond][idir] +=
               Bondfc[jbond]*GtermD*DCosK;


            /*  Derivative wrt to atoms k,i of Zeta for bond i-j  */
            BondZetaD_m[jbond][kbond][idir] +=
               Bondfc1[kbond] * Gterm *
               (BondVector_m[kbond][idir]/BondRadius[kbond]) +
               Bondfc[kbond]*GtermD*DCosK;


            /*  Derivative wrt to atoms j,i of Zeta for bond i-k  */
            BondZetaD_m[kbond][jbond][idir] +=
               Bondfc1[jbond] * Gterm *
               (BondVector_m[jbond][idir]/BondRadius[jbond]) +
               Bondfc[jbond]*GtermD*DCosJ;

            }

         }
         /*
         Finished with storage bond-two-body terms
         (bond interactions)
         */




      /*  Sum up terms for each bond iatom-jatom  */
      AtomBondEnergy = 0.0;
      LOOP (jbond, NumBondNeighbors)
         {

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

         /*  Get atom type  */
         TypeJ = Type[jatom];
         ASSERT(TypeJ<NSPECIES)

         BetaZetaPowN = pow (Beta[TypeI]*BondZeta[jbond], (double) n[TypeI]);
         BondStrength =
            Chi[TypeI][TypeJ] *
            pow
            (
            1 + BetaZetaPowN,
            -1.0/(2.0*n[TypeI])
            );

#ifdef PAIR_ONLY
BondStrength = 1.0;
#endif

         Bondfa  = fa (TypeI, TypeJ, BondRadius[jbond]);
         Bondfr  = fr (TypeI, TypeJ, BondRadius[jbond]);
         Bondfa1 = fa1(TypeI, TypeJ, BondRadius[jbond]);
         Bondfr1 = fr1(TypeI, TypeJ, BondRadius[jbond]);
#if 0
         FA(TypeI,TypeJ,BondRadius[jbond],Bondfa,Bondfa1)
         FR(TypeI,TypeJ,BondRadius[jbond],Bondfr,Bondfr1)
#endif

         AtomBondEnergy +=
            0.5 * Bondfc[jbond] * (Bondfr + BondStrength * Bondfa);


         /*
         ********************************************
         Forces due to interactions between two atoms
         (Two-body forces)
         ********************************************
         */


         /*
         Forces due to bond i-j,
         independent of interaction with neigboring bonds
         (independent of the BondStrength term)
         */

         /*
         Save computation time by ignoring cutoff-derivative
         if before cutoff
         */
         if (BondRadius[jbond]<InteriorCutoff[TypeI][TypeJ])
            {
            TempFactor =
               0.5* Bondfc [jbond] * (Bondfr1 + BondStrength * Bondfa1)
               / BondRadius[jbond];
            }

         /*  Include cutoff derivative  */
         else
            {
            TempFactor =
               0.5 *
               (
               Bondfc1[jbond] * (Bondfr  + BondStrength * Bondfa ) +
               Bondfc [jbond] * (Bondfr1 + BondStrength * Bondfa1)
               )
               / BondRadius[jbond];
            }


         /*
         NOTE:  Recall that Vector[jbond] depends
         on both iatom and jatom
         */
         ForceX = TempFactor * BondVector_m[jbond][X];
         ForceY = TempFactor * BondVector_m[jbond][Y];
         ForceZ = TempFactor * BondVector_m[jbond][Z];

         Force[jatom][X] -= ForceX;
         Force[jatom][Y] -= ForceY;
         Force[jatom][Z] -= ForceZ;

         Force[iatom][X] += ForceX;
         Force[iatom][Y] += ForceY;
         Force[iatom][Z] += ForceZ;

         /*  Calculate system-wide stresses  */
         if (s->StressStepOutput.Save || a->BoxMotionAlgorithm!=BMA_NONE)
            {
            a->Stress[XX] -= ForceX * BondVector_m[jbond][X];
            a->Stress[YY] -= ForceY * BondVector_m[jbond][Y];
            a->Stress[ZZ] -= ForceZ * BondVector_m[jbond][Z];
            a->Stress[YZ] -= ForceY * BondVector_m[jbond][Z];
            a->Stress[ZX] -= ForceZ * BondVector_m[jbond][X];
            a->Stress[XY] -= ForceX * BondVector_m[jbond][Y];
            }

         /*
         ********************************************
         Forces due to interactions between two bonds
         (or three atoms:  three-body forces)
         ********************************************
         */

         /*
         If iatom does not have two or more neighbors,
         then there are no bond interactions
         */
         if (NumBondNeighbors<=1)
            continue;

#ifdef TEST_ZETA
         if (MinStep_g==0 || MinZeta_g>BondZeta[jbond])
            {
            MinStep_g = a->step;
            MinZeta_g = BondZeta[jbond];
            MinI_g    = iatom;
            MinJ_g    = jatom;
            }
#endif

         /*  Calculate derivative of BondStrength wrt Zeta  */
         BondStrengthD =
            -0.25 *
            Bondfc[jbond] * Bondfa * Chi[TypeI][TypeJ] *
            pow( 1 + BetaZetaPowN, -0.5/n[TypeI] - 1.0) *
            BetaZetaPowN / BondZeta[jbond];

         /*
         For currently unknown reasons BondStrengthD is off
         by a factor of 2, correct it.
         */
#ifdef CORRECTION
         BondStrengthD *= 0.5;
#endif

#ifdef PAIR_ONLY
BondStrengthD = 0.0;
#endif

         /*
         Sum contribution to forces from BondStrength term.
         This term contains interactions *between* bonds
         */

         /*
         Derivative wrt j on Energy due to bond between i and j
         - Contribution to force on atom j
         */
         ForceX = BondStrengthD * BondZetaD_m[jbond][jbond][X];
         ForceY = BondStrengthD * BondZetaD_m[jbond][jbond][Y];
         ForceZ = BondStrengthD * BondZetaD_m[jbond][jbond][Z];
         Force[jatom][X] -= ForceX;
         Force[jatom][Y] -= ForceY;
         Force[jatom][Z] -= ForceZ;

         /*
         Derivative wrt j on Energy due to bond between i and j
         - Contribution to force on atom i
         */
         Force[iatom][X] += ForceX;
         Force[iatom][Y] += ForceY;
         Force[iatom][Z] += ForceZ;

         /*  Calculate system-wide stresses  */
         if (s->StressStepOutput.Save || a->BoxMotionAlgorithm!=BMA_NONE)
            {
            a->Stress[XX] -= ForceX * BondVector_m[jbond][X];
            a->Stress[YY] -= ForceY * BondVector_m[jbond][Y];
            a->Stress[ZZ] -= ForceZ * BondVector_m[jbond][Z];
            a->Stress[YZ] -= ForceY * BondVector_m[jbond][Z];
            a->Stress[ZX] -= ForceZ * BondVector_m[jbond][X];
            a->Stress[XY] -= ForceX * BondVector_m[jbond][Y];
            }


         /*
         Sum forces on atom k (and i) due to change of energy
         of bond i-k or i-j due to derivative wrt k (and i)
         */
         LOOP (kbond, NumBondNeighbors)
            {

            /*  Insure that jbond and kbond are different  */
            if (kbond==jbond)
               continue;


            katom = BondNeighborIndex[kbond];

            /*
            Derivative wrt k on Energy due to bond between i and j
            - Contribution to force on atom k
            */
            ForceX = BondStrengthD * BondZetaD_m[jbond][kbond][X];
            ForceY = BondStrengthD * BondZetaD_m[jbond][kbond][Y];
            ForceZ = BondStrengthD * BondZetaD_m[jbond][kbond][Z];
            Force[katom][X] -= ForceX;
            Force[katom][Y] -= ForceY;
            Force[katom][Z] -= ForceZ;


            /*
            Derivative wrt k on Energy due to bond between i and j
            - Contribution to force on atom i
            */
            Force[iatom][X] += ForceX;
            Force[iatom][Y] += ForceY;
            Force[iatom][Z] += ForceZ;

            /*  Calculate system-wide stresses  */
            if (s->StressStepOutput.Save || a->BoxMotionAlgorithm!=BMA_NONE)
               {
               a->Stress[XX] -= ForceX * BondVector_m[kbond][X];
               a->Stress[YY] -= ForceY * BondVector_m[kbond][Y];
               a->Stress[ZZ] -= ForceZ * BondVector_m[kbond][Z];
               a->Stress[YZ] -= ForceY * BondVector_m[kbond][Z];
               a->Stress[ZX] -= ForceZ * BondVector_m[kbond][X];
               a->Stress[XY] -= ForceX * BondVector_m[kbond][Y];
               }

            }

          }
         /*  Finished with bond "one-body" contributions  */


      /*  Sum total energy  */
      TotalBondEnergy += AtomBondEnergy;
      }
      /*   Finished with each atom (iatom)  */

   a->epot = TotalBondEnergy + TotalPairEnergy;

   }





/*
************************************************************************
Local Functions
************************************************************************
*/
static void InitParameters (void)
   {
   int i;
   int j;

   /*  Initialize just off-diagonal terms  */
   LOOP (i,NSPECIES)
   LOOP (j,NSPECIES)
      {
      if (i==j)
         continue;
      ARITH_AVERAGE(Lamda, i, j)
      ARITH_AVERAGE(Mu,    i, j)
      GEOM_AVERAGE (A,     i, j)
      GEOM_AVERAGE (B,     i, j)
      }

   /*  Initialize cutoffs  */
   InitializeCutoff();


   LOOP (i, NSPECIES)
      {
      c2 [i] = SQR(c[i]);
      d2 [i] = SQR(d[i]);
      if (d2[i]==0.0)
         {
         printf ("ERROR:  Value of parameter d[%i] cannot be zero.\n", i);
         exit   (1);
         }
      cd2[i] = c2[i]/d2[i];
      }
   }



static void InitPairPotential(void)
   {
   int i;
   LOOP (i, MAX_PAIR)
      {
      PairCoeff_m [i] = NULL;
      PairCutoff_m[i] = 0.0;
      NumPairCoeff_m[i] = 0;
      }
   }



void InitializeCutoff(void)
   {
   int i;
   int j;


   /*  Intialize off diagnol terms  */
   LOOP (i,NSPECIES)
   LOOP (j,NSPECIES)
      {
      if (i==j)
         continue;
      GEOM_AVERAGE (InteriorCutoff,     i, j)
      GEOM_AVERAGE (ExteriorCutoff,     i, j)
      }

   /*  Intialize entire arrays  */
   LOOP (i,NSPECIES)
   LOOP (j,NSPECIES)
      {
      ExteriorCutoffSqr[i][j] = SQR(ExteriorCutoff[i][j]);
      CutoffInterval   [i][j] = ExteriorCutoff[i][j] - InteriorCutoff[i][j];
      }

   }




static void ScaleDistance (double Scale)
   {
   int ispecies;
   int jspecies;
   int icoeff;
   int icombinedtype;


   /*  Scale original Tersoff potential  */
   LOOP(ispecies, NSPECIES)
   LOOP(jspecies, NSPECIES)
      {
      Mu                [ispecies][jspecies]  /=  Scale;
      Lamda             [ispecies][jspecies]  /=  Scale;
      InteriorCutoff    [ispecies][jspecies]  *=  Scale;
      ExteriorCutoff    [ispecies][jspecies]  *=  Scale;
      CutoffInterval    [ispecies][jspecies]  *=  Scale;
      ExteriorCutoffSqr [ispecies][jspecies]  *=  (Scale*Scale);
      }


   /*  If no pair potential then return  */
   if (!UsePairPotential_m)
      return;


   /*  Scale auxillary pair potential  */
   LOOP (icombinedtype, MAX_PAIR)
      {
      LOOP (icoeff, NumPairCoeff_m[icombinedtype])
         {
         ASSERT (PairCoeff_m[icombinedtype]!=NULL)
         PairCoeff_m[icombinedtype][icoeff]
            /= pow(Scale, (double) (icoeff+2) );
         }
      PairCutoff_m[icombinedtype] *= Scale;
      }

   }


static void ScaleEnergy (double Scale)
   {
   int ispecies;
   int jspecies;
   int icoeff;
   int icombinedtype;


   /*  Scale original Tersoff potential  */
   LOOP(ispecies, NSPECIES)
   LOOP(jspecies, NSPECIES)
      {
      A[ispecies][jspecies] *= Scale;
      B[ispecies][jspecies] *= Scale;
      }

   /*  If no pair potential then return  */
   if (!UsePairPotential_m)
      return;


   /*  Scale auxillary pair potential  */
   LOOP (icombinedtype, MAX_PAIR)
   LOOP (icoeff, NumPairCoeff_m[icombinedtype])
      {
      ASSERT (PairCoeff_m[icombinedtype]!=NULL)
      PairCoeff_m[icombinedtype][icoeff] *= Scale;
      }

   }



static double GetEnergy (Particle_t *a)
   {
   int iatom;
   int jatom;
   int jbond;
   int kbond;
   int TypeI;
   int TypeJ;
   int *Neighbors = NULL;
   int    NumNeighbors;
   double PairEnergy;
   double AtomBondEnergy;
   double TotalPairEnergy;
   double TotalBondEnergy;
   int    NumBondNeighbors;
   double BondRadius        [MAX_BOND_NEIGHBOR];
   double BondRadiusSqr     [MAX_BOND_NEIGHBOR];
   double BondZeta          [MAX_BOND_NEIGHBOR];
   double Bondfc            [MAX_BOND_NEIGHBOR];
   int    BondNeighborIndex [MAX_BOND_NEIGHBOR];
   double BondStrength;
   double Gterm;
   double CosTerm;
   double (*Coord)[NDIR];
   double  *Eatom;
   BYTE    *Type;
   int    idir;
   int    NumAtom;
   int    CombinedType;


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


   /*  Intialize pointers  */
   Coord = (double (*)[NDIR]) a->cur;
   Type  =                    a->type;
   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];
      }


   /*  Setup cutoff info  */
   CalcTotalCutoff (TotalCutoff_m, TotalCutoffSqr_m);

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

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

   /*  Set Eatom[] to 0  */
   LOOP (iatom, NumAtom)
      Eatom[iatom] = 0.0;

   /*  Calculate energy for every atom  */
   TotalPairEnergy = 0.0;
   TotalBondEnergy = 0.0;
   LOOP (iatom, NumAtom)
      {


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

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

         /*  Get atom types  */
         TypeI = Type[iatom];
         TypeJ = Type[jatom];

         if
         (
         !CalcSeparation
            (
            Coord[iatom],
            Coord[jatom],
            BondVector_m[NumBondNeighbors],
            TotalCutoff_m[TypeI][TypeJ]
            )
         )
            continue;

         /*  Calculate BondRadius  */
         BondRadiusSqr[NumBondNeighbors] =
            MAG_SQR(BondVector_m[NumBondNeighbors]);

         /*  Test cutoff  */
         if (
            BondRadiusSqr[NumBondNeighbors] >
            TotalCutoffSqr_m[TypeI][TypeJ]
            )
            continue;

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

         /*  Calculate pair potential  */
         if ( UsePairPotential_m )
            {

            /*  Calcualte combined type  */
            CombinedType = COMBINED_TYPE(TypeI,TypeJ);


            /*  Test if particle within pair potential range  */
            if (BondRadius[NumBondNeighbors]<PairCutoff_m[CombinedType])

            /*  Calculate only iatom<jatom (to prevent double counting)  */
            if (iatom<jatom)

               {
               /*  ADD ENERGY CALCULATE HERE  */
               PairEnergy =
                  PairPotentialEnergy(TypeI,TypeJ,BondRadius[NumBondNeighbors]);
               TotalPairEnergy += PairEnergy;
               Eatom[iatom] += 0.5*PairEnergy;
               Eatom[jatom] += 0.5*PairEnergy;
               }

            }
            /*  Finished pair interaction with neighbors of iatom  */


         /*
         If both particles are not either Si or C,
         then skip to next pair of iatom,jneigh
         */
         if (!IS_TERSOFF(TypeI,TypeJ))
            continue;

         /*
         If particle outside Tersoff range then skip to next neighbor
         */
         if ( BondRadius[NumBondNeighbors] >= ExteriorCutoff[TypeI][TypeJ] )
            continue;

         /*
         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 (BondRadius[NumBondNeighbors]==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  */
         Bondfc[NumBondNeighbors] =
            fc(TypeI,TypeJ,BondRadius[NumBondNeighbors]);

         /*  Initialize Zeta sum to zero  */
         BondZeta[NumBondNeighbors] = 0.0;

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



      /*  Skip Tersoff calculation if not C or Si  */
      if (!IS_TERSOFF(TypeI,TypeI))
         continue;
      ASSERT(TypeI<NSPECIES)

      /*
      Following calculations on iatom are
      exclusively for Tersoff potential
      */

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


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

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


         Gterm   =
           1 +
           cd2[TypeI] -
           c2 [TypeI] / (d2[TypeI] + SQR(h[TypeI] - CosTerm));


         BondZeta[jbond] += Bondfc[kbond] * Gterm;
         BondZeta[kbond] += Bondfc[jbond] * Gterm;
         }
         /*
         Finished with storage bond-two-body terms
         (bond interactions)
         */



      /*  Intiailize energy for bonds to iatom  */
      AtomBondEnergy = 0.0;

      /*  Sum up terms for each bond iatom-jatom  */
      LOOP (jbond, NumBondNeighbors)
         {

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

         /*  Get atom type  */
         TypeJ = Type[jatom];
         ASSERT(TypeJ<NSPECIES)

         BondStrength =
            Chi[TypeI][TypeJ] *
            pow
            (
            1 + pow ( Beta[TypeI]*BondZeta[jbond], (double) n[TypeI] ),
            -1.0/(2.0*n[TypeI])
            );

#ifdef PAIR_ONLY
BondStrength = 1.0;
#endif


         AtomBondEnergy +=
            0.5*
            Bondfc[jbond] *
            (
            fr(TypeI,TypeJ,BondRadius[jbond])  +
            BondStrength * fa(TypeI,TypeJ,BondRadius[jbond])
            );

          }


      /*  Sum total energy  */
      Eatom[iatom]    += AtomBondEnergy;
      TotalBondEnergy += AtomBondEnergy;
      }

   /* Internal test  */
   {
   double TempEnergy=0;
   LOOP (iatom, NumAtom)
      {
      TempEnergy += Eatom[iatom];
      }
   ASSERT(fabs(TempEnergy-(TotalBondEnergy+TotalPairEnergy))<=1e-6*fabs(TempEnergy))
   }


   return TotalBondEnergy + TotalPairEnergy;
   }





static double fa (int Type1, int Type2, double Radius)
   {
   return -B[Type1][Type2] * exp( -Mu  [Type1][Type2] * Radius);
   }



static double fa1 (int Type1, int Type2, double Radius)
   {
   return
      B[Type1][Type2] * Mu[Type1][Type2] *
      exp( -Mu  [Type1][Type2] * Radius);
   }



static double fr (int Type1, int Type2, double Radius)
   {
   return  A[Type1][Type2] * exp(-Lamda[Type1][Type2] * Radius);
   }



static double fr1 (int Type1, int Type2, double Radius)
   {
   return
      -A[Type1][Type2] * Lamda[Type1][Type2] *
      exp(-Lamda[Type1][Type2] * Radius);
   }



static double fc (int Type1, int Type2, double Radius)
   {
   if (Radius < InteriorCutoff [Type1][Type2])
      return 1.0;
   if (Radius > ExteriorCutoff[Type1][Type2])
      return 0.0;
   return 0.5 *
      (
      1.0 +
      cos
      (
      M_PI *
      (Radius-InteriorCutoff[Type1][Type2]) /
      CutoffInterval[Type1][Type2]
      )
      );
   }



static double fc1 (int Type1, int Type2, double Radius)
   {
   if (Radius < InteriorCutoff [Type1][Type2])
      return 0.0;
   if (Radius > ExteriorCutoff[Type1][Type2])
      return 0.0;
   return -(0.5 * M_PI / CutoffInterval[Type1][Type2]) *
      sin
      (
      M_PI *
      (Radius-InteriorCutoff[Type1][Type2]) /
      CutoffInterval[Type1][Type2]
      );
   }




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



static double GetMaxCutoff (void)
   {
   double MaxCutoff = 0;
   int i;
   int j;
   int icombinedtype;

   /*  Get cutoff from all pair potentials */
   LOOP (i, NSPECIES)
   LOOP (j, NSPECIES)
      {
      if (ExteriorCutoff[i][j] > MaxCutoff)
         MaxCutoff = ExteriorCutoff[i][j];
      }

   /*  Include Auxillary pair potential cutoff  */
   if (UsePairPotential_m)
   LOOP (icombinedtype, MAX_PAIR)
      {
      if (PairCutoff_m[icombinedtype] > MaxCutoff)
         MaxCutoff = PairCutoff_m[icombinedtype];
      }
   return MaxCutoff;
   }



static double PairPotentialEnergy (int Type1, int Type2, double r)
   {
   double Sum;
   double u;
   int    icoeff;
   int    CombinedType;


   /*  Get combined type  */
   CombinedType = COMBINED_TYPE(Type1, Type2);

   /*  No coefficients - return 0.0  */
   if (NumPairCoeff_m[CombinedType]==0)
      return 0.0;

   /*  Beyond cutoff, return 0.0  */
   u = PairCutoff_m[CombinedType] - r;
   if (u<0)
      return 0.0;

   /*  Get  a[n-1]*u^(n+1) + a[n-2]*u^n ...  + a[0]*u^2  */
   Sum = PairCoeff_m[CombinedType][NumPairCoeff_m[CombinedType]-1];

   for (icoeff=NumPairCoeff_m[CombinedType]-2; icoeff>=0; icoeff--)
      {
      Sum = Sum*u + PairCoeff_m[CombinedType][icoeff];
      }
   Sum *= u*u;

   return Sum;
   }



static double PairPotentialForce
(
int Type1,
int Type2,
double r
)
   {
   double Sum;
   double u;
   int    icoeff;
   int    CombinedType;
   int    NumCoeff;


   /*  Get combined type  */
   CombinedType = COMBINED_TYPE(Type1, Type2);

   /*  No coefficients - return 0.0  */
   NumCoeff = NumPairCoeff_m[CombinedType];
   if (NumCoeff==0)
      return 0.0;

   /*  Beyond cutoff, return 0.0  */
   u = PairCutoff_m[CombinedType] - r;
   if (u<0)
      return 0.0;

   /*  Get  a[n-1]*u^(n+1) + a[n-2]*u^n ...  + a[0]*u^2  */
   Sum = (NumCoeff+1)*PairCoeff_m[CombinedType][NumCoeff-1];
   for (icoeff=NumCoeff-2; icoeff>=0; icoeff--)
      {
      Sum = Sum*u + (icoeff+2)*PairCoeff_m[CombinedType][icoeff];
      }
   Sum *= u;

   return Sum;

   }



/*
Calculate largest cutoff per species between Tersoff potential and
auxillary pair potential
*/
void CalcTotalCutoff
(
double TotalCutoff[MAX_PAIR][MAX_PAIR],
double TotalCutoffSqr[MAX_PAIR][MAX_PAIR]
)
   {
   int itype;
   int jtype;
   int CombinedType;

   LOOP (itype, MAX_TYPE)
   LOOP (jtype, MAX_TYPE)
      {

      /*  Calculate combined type  */
      CombinedType = COMBINED_TYPE(itype,jtype);

      /*  Get cutoff for Tersoff potential  */
      if (IS_TERSOFF(itype,jtype))
         {
         TotalCutoff[itype][jtype] = ExteriorCutoff[itype][jtype];
         }
      else
         {
         TotalCutoff[itype][jtype] = 0.0;
         }

      /*  Include cutoff for pair potential  */
      if
      (
      UsePairPotential_m &&
      NumPairCoeff_m[CombinedType] > 0 &&
      TotalCutoff_m[itype][jtype]<PairCutoff_m[CombinedType]
      )
         {
         TotalCutoff[itype][jtype] = PairCutoff_m[CombinedType];
         }

      /*  Calculate square of cutoff  */
      TotalCutoffSqr[itype][jtype] = SQR(TotalCutoff_m[itype][jtype]);
      }
   }


/*
Removed 2 Feb 1998, replace dynamic with static arrays to
cut down allocation overhead and to streamline output
for CheckMem
*/
#ifdef NEWALLOC
double ***Allocate3D(int dim1, int dim2, int dim3)
   {
   double ***Array = NULL;
   int i;
   int j;

   ALLOCATE  (Array, double **, dim1)
   LOOP (i, dim1)
      {
      ALLOCATE (Array[i], double *, dim2)
      LOOP (j, dim2)
         {
         ALLOCATE(Array[i][j], double, dim3)
         }
      }
   return Array;
   }

void Free3D(double ***Array,  int dim1, int dim2)
   {
   int i;
   int j;

   LOOP (i, dim1)
      {
      LOOP (j, dim2)
         {
         FREE (Array[i][j])
         }
      FREE (Array[i])
      }
   FREE (Array)
   }

double **Allocate2D(int dim1, int dim2)
   {
   double **Array = NULL;
   int i;

   ALLOCATE  (Array, double *, dim1)
   LOOP (i, dim1)
      {
      ALLOCATE (Array[i], double, dim2)
      }
   return Array;
   }


void Free2D(double **Array,  int dim1)
   {
   int i;

   LOOP (i, dim1)
      {
      FREE (Array[i])
      }
   FREE (Array)
   }
#endif