cdcsi.c

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

         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