cdstil.c

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

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