/* xmd - molecular dynamics for metals and ceramics By Jonathan Rifkin 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. */ /* ************************************************************************ History ************************************************************************ */ /* 10 Feb 1998 Stillinger-Weber Si Potential Converted from cdcsi.c Tersoff potential */ /* ************************************************************************ Include Files ************************************************************************ */ #include #include #include #include "strsub.h" #include "parse.h" #include "iomngr.h" #include "cdhouse.h" #include "cdsubs.h" #include "cdsearch.h" #include "cdpairf.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 /* Maximum number of neighbors for any one atom */ #define INIT_MAX_BOND_NEIGHBOR 32 /* Size of input string buffer */ #define NBUFF 256 /* ************************************************************************ Macros ************************************************************************ */ #define MAG_SQR(V)\ ((V)[X]*(V)[X] + (V)[Y]*(V)[Y] + (V)[Z]*(V)[Z]) #define DOT(A,B) \ ((A)[X]*(B)[X] + (A)[Y]*(B)[Y] + (A)[Z]*(B)[Z]) /* ************************************************************************ Global Variables ************************************************************************ */ extern Simulation_t *s; extern LIST *inlist; extern FILE *out; /* ************************************************************************ Module-Scope Variables ************************************************************************ */ /* Potential parameters from Stillinger */ static double A_m = 7.049556277; static double B_m = 0.6022245584; static double p_m = 4; static double a_m = 1.80; static double lamda_m = 21.0; static double gamma_m = 1.20; static double sigma_m = 2.0951e-8; /* cm */ static double epsilon_m = 3.4723e-12; /* erg/atom-pair*/ /* Values obtained from above parameters */ static double Cutoff_m; static double PairFac1_m; static double PairFac2_m; static double PairFac3_m; static double AngleFac1_m; static double AngleFac2_m; static double HalfBox_m [NDIR]; static double Box_m [NDIR]; static BOOLEAN BoundRep_m [NDIR]; static BOOLEAN UsePairPotential_m = FALSE; static int MaxBondNeighbor_m = INIT_MAX_BOND_NEIGHBOR; /* ************************************************************************ Local Function Prototypes ************************************************************************ */ static void InitParameters(void); static void ScaleDistance(double); static void ScaleEnergy (double); static double GetEnergy (Particle_t *a); static BOOLEAN CalcSeparation (double *Coord1, double *Coord2, double *Vector, double Cutoff); static double GetSeparationDir(double Separation, int idir); /* ************************************************************************ Exported Functions ************************************************************************ */ #pragma argsused void stil_energy_list (Particle_t *a, Simulation_t *s, int UseSelect) { a->epot = GetEnergy (a); } #pragma argsused void STILL_Parse (char *instr) { int FUNC_ScaleType; int ifunc; int nfunc; int ErrorCode; char *tokstr = NULL; double Scale; Function_t **FuncListPtr = NULL; if (IsFirstToken("MAXBONDCNT",instr)) { /* Remove leading MAXBONDCNT token */ strhed (&instr); /* Get next number */ MaxBondNeighbor_m = intstrf (&instr); /* Sanity check */ if (MaxBondNeighbor_m<4) { ERROR_PREFIX printf ("Maximum number of bond neighbors must be 4 or more.\n"); CleanAfterError(); } } /* Scale potential energy or distance */ else if (IsFirstToken("SCALE",instr)) { /* Remove leading "SCALE" token */ strhed(&instr); /* Get next token */ tokstr = strhed (&instr); if (!strcmpi(tokstr, "DISTANCE")) { Scale = dblstrf (&instr); ScaleDistance (Scale); FUNC_ScaleType = FUNC_INPUT; } else if (!strcmpi(tokstr, "ENERGY")) { Scale = dblstrf (&instr); ScaleEnergy (Scale); FUNC_ScaleType = FUNC_OUTPUT; } /* Unknown SCALE option, return with Warning */ else { printf("WARNING: Unknown option to Still-Weber Si"); printf(" potential SCALE command (%s).\n", tokstr); IncrementNumberOfWarnings(); return; } /* unknown SCALE option */ /* Apply scale to auxillary Pair potentials */ if (UsePairPotential_m) { FuncListPtr = PAIR_GetFuncListPtr(); nfunc = PAIR_GetNumFunction(); LOOP (ifunc, nfunc) { FUNC_Scale (FuncListPtr[ifunc], FUNC_ScaleType, Scale); } } return; } /* SCALE potion */ /* Pass option onto PAIR routine */ else { if (UsePairPotential_m) { ErrorCode = PAIR_Parse (instr); } if (!UsePairPotential_m || ErrorCode==PAIR_COMMAND_ERROR) { printf ("WARNING: Unknown option to Still-Weber Si potential(%s).\n", tokstr); IncrementNumberOfWarnings(); } } } /* Initialze potential */ void STILL_Init(char *instr) { char *tokstr = NULL; /* Fill out paramters tables */ InitParameters(); /* Scale distance from angstroms to cm */ ScaleDistance (sigma_m); /* Scale energy from eV to ergs */ ScaleEnergy (epsilon_m); /* Check for auxillary pair potentials */ tokstr = strhed (&instr); if (!strcmpi("PAIR",tokstr)) { UsePairPotential_m = TRUE; PAIR_Init (instr); } /* Print info */ printf ("\n"); printf (" Stillinger-Weber silicon model\n"); printf (" reference:\n"); printf (" \"Computer simulation of local order in"); printf (" condensed phases of silicon\".\n"); printf (" F. H. Stllinger, T. A. Weber,"); printf (" Phys. Rev. B 31 (8) 5262 (1985)\n"); printf ("\n"); } /* Force routine called externally */ #pragma argsused void stil_calcforc (Particle_t *a, Simulation_t *s) { BOOLEAN UseStress; int iatom; int jatom; int katom; 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 (*Force)[NDIR]; double *Eatom; double CutoffSqr; double Radius; double RadiusSqr; double InvRadius; double PairEnergy; double PairTerm1; double AtomBondEnergy; double CosTerm; double CosTerm3; double Radius4; double PairCutoff; double EnergyFac1; double BondEnergy; double ForceFac1; double ForceFac2; double ForceFac1a; double ForceFac2a; double PairForce; double ForceComponent[NDIR]; double TotalEnergy; int idir; int NumAtom; /* Test for atoms */ NumAtom = a->np; if (NumAtom==0) return; /* 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; Force = (double (*)[NDIR]) a->f; 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]; } /* Set Eatom[] and Force[] to 0 */ LOOP (iatom, NumAtom) { Eatom[iatom] = 0.0; Force[iatom][X] = 0.0; Force[iatom][Y] = 0.0; Force[iatom][Z] = 0.0; } /* Initialize stresses to zero */ UseStress = (s->StressStepOutput.Save || a->BoxMotionAlgorithm!=BMA_NONE); if (UseStress) { 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; } /* Set cutoff for call to search_all() */ s->cutoff = Cutoff_m; /* Include pair cutoff if using auxillary pair routines */ if (UsePairPotential_m && Cutoff_mcutoff = 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) { /* Use only Silicon atoms (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; /* Convert BondVector[] from x,y,z to x/r,y/r,z/r */ InvRadius = 1.0 / Radius; BondVector[NumBondNeighbors][X] *= InvRadius; BondVector[NumBondNeighbors][Y] *= InvRadius; BondVector[NumBondNeighbors][Z] *= InvRadius; /* If we have gotten here, then we are including jneigh as a neighbor of iatom which forms a 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: "); printf ("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 (iatomStress[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_mcutoff = 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=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; }