cdcsi.c

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

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