cdcsi.c

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

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


         /*