cdcsi.c

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

/*
************************************************************************
Comments
************************************************************************
*/

/*
TERSOFF Si-C POTENTIAL
*/

/*
************************************************************************
History
************************************************************************
*/

/*

6 Nov 1996 - Got Forces to run correctly.  Problem was caused by
             macro SQR() in cdhouse.h.  It was defined as (X)*(X),
             should have been ((X)*(X)).  The old definition was
             a problem when calculating GtermD in csi_calcforc().
             The statement  ?/SQR(GtermDeninominator) was not doing
             what had been intended.
           - There is another problem with forces.  The pure pair
             forces are fine, but an examination of the sum of KE
             and PE under CLAMP OFF shows that in order for sum to
             be constant, fluctuations in KE need to be reduced by
             half.  So those forces due to the angular terms are
             reduced by a factor of 2.
22 Nov 1996
   VERSION 2.2.1
          - Previos problems were fixed a week ago.
   VERSION 2.2.2
          - Add stress calculation to Tersoff force function

16 Jul 1997
   Version 2.2.5
          - Add separate pair potential
          - Add routine CalcTotalCutoff() to calculate the maximum
            cutoff per species for both Tersoff and auxillary pair potential
          - command  POTENTIAL PAIR ncoeff cutoff
                     a1 a2 a3 a4
            resulting energy P(r) = a1*(r-rc)^2 + a2*(r-rc)^3 + ..
          - set s->cutoff=GetMaxCutoff() before each call  to
            search_all()

22 Aug 1997
   Version 2.4.8
          - When implementing separate pair potential above, forgot
            to initilize AtomBondEnergy=0.0 before calcualting each
            atom


*/




/*
************************************************************************
Include Files
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "strsub.h"
#include "parse.h"
#include "iomngr.h"
#include "cdhouse.h"
#include "cdsubs.h"
#include "cdsearch.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


/*  Number of Tersoff atom species:  Carbon and Silicon  */
#define NSPECIES 2

/*
Maximum number of atom types -
types other than 1 and 2 are pair only
*/
#define MAX_TYPE 5
#define MAX_PAIR MAX_TYPE*(MAX_TYPE+1)/2


/*  Maximum number of neighbors for any one atom  */
#define MAX_BOND_NEIGHBOR 24

/*  Size of input string buffer  */
#define NBUFF 256





/*
************************************************************************
Macros
************************************************************************
*/



#define ARITH_AVERAGE(A,I,J) \
   A[I][J] = 0.5* (A[I][I] + A[J][J]);

#define GEOM_AVERAGE(A,I,J) \
   A[I][J] = sqrt (A[I][I] * A[J][J]);

#define SQR_RAD(A,B) \
   ( \
      ((A)[X]-(B)[X])*((A)[X]-(B)[X]) + \
      ((A)[Y]-(B)[Y])*((A)[Y]-(B)[Y]) + \
      ((A)[Z]-(B)[Z])*((A)[Z]-(B)[Z])   \
   )

#define MAG_SQR(A) \
   ((A)[X]*(A)[X] + (A)[Y]*(A)[Y] + (A)[Z]*(A)[Z])

#define DOT(A,B) \
   ((A)[X]*(B)[X] + (A)[Y]*(B)[Y] + (A)[Z]*(B)[Z])


#define FA(I,J,R,F,F1) \
   (F ) = -B[I][J] * exp(-Mu[I][J]*(R)); \
   (F1) = -Lamda[I][J] * (F);

#define FR(I,J,R,F,F1) \
   (F ) = A[I][J] * exp(-Lamda[I][J]*(R)); \
   (F1) = -Lamda[I][J] * (F);

#define IS_TERSOFF(I,J) \
   ((I)<NSPECIES && (J)<NSPECIES)


/*
************************************************************************
Global Variables
************************************************************************
*/
extern Simulation_t *s;
extern LIST         *inlist;
extern FILE         *out;

#ifdef TEST_ZETA
   extern double MinZeta_g;
   extern int    MinStep_g;
   extern int    MinI_g;
   extern int    MinJ_g;
#endif



/*
************************************************************************
Module-Scope Variables
************************************************************************
*/
static double  HalfBox_m  [NDIR];
static double  Box_m      [NDIR];
static BOOLEAN BoundRep_m [NDIR];


/*
This parameters depend on species of both atoms in pair
The cross terms must be set by the initialization routine
*/


/*  Cross terms are arithmetic average of pure terms  */
static double Lamda            [NSPECIES][NSPECIES] = {{3.4879, 0},{0,        2.4799}};
static double Mu               [NSPECIES][NSPECIES] = {{2.2119, 0},{0,        1.7322}};

/*  Cross terms are geometric  average of pure terms  */
static double A                [NSPECIES][NSPECIES] = {{ 1.3936e3, 0},{0,  1.8308e3}};
static double B                [NSPECIES][NSPECIES] = {{ 3.467e2,  0},{0,  4.7118e2}};
static double InteriorCutoff   [NSPECIES][NSPECIES] = {{1.8, 0},{0,  2.7}};
static double ExteriorCutoff   [NSPECIES][NSPECIES] = {{2.1, 0},{0,  3.0}};
static double CutoffInterval   [NSPECIES][NSPECIES] = {{0.3, 0},{0,  0.3}};
static double ExteriorCutoffSqr[NSPECIES][NSPECIES];

/*  This paramter cross terms not an average of pure terms  */
static double Chi              [NSPECIES][NSPECIES] = {{1.0, 0.9776}, {0.9776, 1.0}};

/*  These parameters depend on upon species of one atom  */
static double Beta             [NSPECIES]           = {1.5724e-7,    1.1000e-6};
static double n                [NSPECIES]           = {7.2751e-1,    7.8734e-1};
static double c                [NSPECIES]           = {3.8049e4,     1.0039e5};
static double d                [NSPECIES]           = {4.384,        1.6217e1};
static double h                [NSPECIES]           = {-5.7058e-1,  -5.9825e-1};

/*  The values of these parameters depend on previous ones  */
static double c2  [NSPECIES];
static double d2  [NSPECIES];
static double cd2 [NSPECIES];


/*  These values control the auxiallary pair potential  */
static BOOLEAN UsePairPotential_m = FALSE;
static BOOLEAN IsPairPotentialInitialized_m = FALSE;
static int     NumPairCoeff_m[MAX_PAIR];
static double  PairCutoff_m  [MAX_PAIR];
static double *PairCoeff_m   [MAX_PAIR];

/*  These array hold temporary values for Tersoff calculation  */
static double BondZetaD_m [MAX_BOND_NEIGHBOR] [MAX_BOND_NEIGHBOR] [NDIR];
static double BondVector_m[MAX_BOND_NEIGHBOR] [NDIR];
static double TotalCutoff_m    [MAX_PAIR] [MAX_PAIR];
static double TotalCutoffSqr_m [MAX_PAIR] [MAX_PAIR];


/*
************************************************************************
Local Function Prototypes
************************************************************************
*/

static void    InitParameters   (void);
static void    InitPairPotential(void);
static void    InitializeCutoff (void);
static double  GetEnergy        (Particle_t *);
static void    ScaleDistance    (double);
static void    ScaleEnergy      (double);
static double  GetMaxCutoff     (void);
static BOOLEAN CalcSeparation   (double *, double *, double *, double);
static double  GetSeparationDir (double, int);

static double  fr (int,int,double);
static double  fa (int,int,double);
static double  fc (int,int,double);
static double  fr1(int,int,double);
static double  fa1(int,int,double);
static double  fc1(int,int,double);
static double  PairPotentialEnergy (int, int, double);
static double  PairPotentialForce  (int, int, double);

void CalcTotalCutoff
   (double [MAX_PAIR][MAX_PAIR], double [MAX_PAIR][MAX_PAIR]);


double ***Allocate3D(int, int, int);
void      Free3D    (double ***, int, int);
double  **Allocate2D(int, int);
void      Free2D    (double ** , int);

/*
************************************************************************
Exported Functions
************************************************************************
*/
#pragma argsused
void csi_energy_list  (Particle_t *a, Simulation_t *s, int UseSelect)
   {
   a->epot   = GetEnergy (a);
   }

#pragma argsused
void csi_read_potential (char *tokstr, char *instr)
   {
   double Scale;
   int    Type;
   int    iread;
   BOOLEAN EndOfFile;
   double InputInteriorCutoff;
   double InputExteriorCutoff;
   char   InputBuffer[NBUFF];
   int    Type1;
   int    Type2;
   int    CombinedType;


   /*  Initialize potential  */
   if (!strcmpi(tokstr, "INIT"))
      {

      /*  Fill out paramters tables  */
      InitParameters();

      /*  Initialize Pair Potential  */
      InitPairPotential();

      /*  Scale distance from angstroms to cm  */
      ScaleDistance (1.0e-8);

      /*  Scale energy from eV to ergs  */
      ScaleEnergy (1.0/EV);

      /*  Print info  */
      printf ("\n");
      printf ("   Tersoff's Silicon-Carbon potential\n");
      printf ("   reference:\n");
      printf ("   \"Modeling solid-state chemistry:  Interatomic potentials\n");
      printf ("   for multicomponent systems\",\n");
      printf ("   J. Tersoff, Phys Rev. B 39 (8), 5566-5568 (15 March 1989)\n");
#ifdef PAIR_ONLY
      printf ("\nWARNING:  ");
      printf ("This version of potential uses PAIR POTENTIAL ONLY.\n");
      printf ("No angular forces are used.\n");
#endif
      printf ("\n");

      return;
      }


   /*  Scale potential energy or distance  */
   if (!strcmpi(tokstr, "SCALE"))
      {

      tokstr = strhed (&instr);
      if (!strcmpi(tokstr, "DISTANCE"))
         {
         Scale = dblstrf (&instr);
         ScaleDistance (Scale);
         return;
         }

      else if (!strcmpi(tokstr, "ENERGY"))
         {
         Scale = dblstrf (&instr);
         ScaleEnergy (Scale);
         return;
         }

      else
         {
         printf ("WARNING:  Unknown option to Tersoff C-Si potential SCALE command (%s).\n",
                 tokstr);
         IncrementNumberOfWarnings();
         return;
         }
      }



   /*  Set potential cutoff  */
   if (!strcmpi(tokstr, "CUTOFF"))
      {
      Type = intstrf (&instr);
      if (Type!=1 && Type!=2)
         {
         printf ("WARNING:  Type is out of range (1,2) for Tersoff C-Si Cutoff command.\n");
         IncrementNumberOfWarnings();
         return;
         }

      InputInteriorCutoff = dblstrf (&instr) * 1e-8;
      InputExteriorCutoff = dblstrf (&instr) * 1e-8;
      if (InputInteriorCutoff>InputExteriorCutoff)
         {
         printf ("WARNING:  Interior Cutoff exceeds Exterior Cutoff.  Will ignore CUTOFF command.\n");
         IncrementNumberOfWarnings();
         return;
         }

      /*  Set up cutoffs  */
      Type--;
      InteriorCutoff[Type][Type] = InputInteriorCutoff;
      ExteriorCutoff[Type][Type] = InputExteriorCutoff;

      /*  Initialize cutoff  */
      InitializeCutoff();

      /*  return  */
      return;
      }

   /*  Use Auxillary Pair Potential  */
   if (!strcmpi(tokstr, "PAIR"))
      {
      tokstr = strhed (&instr);

      /*  Turn auxillary pair potential off  */
      if (!strcmpi(tokstr, "OFF"))
         {
         UsePairPotential_m = FALSE;
         }

      else
         {

         /*  Skip ON token  */
         if (!strcmpi(tokstr,"ON"))
            {
            tokstr= strhed (&instr);
            UsePairPotential_m = TRUE;
            }

         /*  Read parameters if not blank  */
         if (!IsBlank(tokstr))
            {

            /*  Read types  */
            Type1 = intstrf(&tokstr)-1;
            Type2 = intstrf(&instr )-1;
            CombinedType = COMBINED_TYPE(Type1,Type2);

            /*  Check types  */
            if (Type1<0 || Type1>=MAX_TYPE)
               {
               printf ("ERROR:  Type 1 (%i) is out of range [1..%i]\n",
                  Type1+1, MAX_TYPE);
               CleanAfterError();
               }

            if (Type2<0 || Type2>=MAX_TYPE)
               {
               printf ("ERROR:  Type 2 (%i) is out of range [1..%i]\n",
                  Type2+1, MAX_TYPE);
               CleanAfterError();
               }

            CombinedType = COMBINED_TYPE(Type1,Type2);
            NumPairCoeff_m[CombinedType] =      dblstrf (&instr);
            PairCutoff_m  [CombinedType] = 1e-8*dblstrf (&instr);

            /*  Allocate coefficients  */
            FREE (PairCoeff_m[CombinedType])
            ALLOCATE
               (
               PairCoeff_m[CombinedType],
               double,
               NumPairCoeff_m[CombinedType]
               )

            iread = 0;
            EndOfFile      = FALSE;
            while (iread<NumPairCoeff_m[CombinedType] && !EndOfFile)
               {

               /*  Read new line if current line is blank  */
               if (IsBlank(instr))
                  {
                  EndOfFile =
                     (NULL==m_gets_list(InputBuffer, NBUFF, inlist));
                  instr = InputBuffer;
                  }

               /*  Premature end-of-file, end program  */
               if (EndOfFile)
                  {
                  printf ("ERROR: File ends before coefficients were found\n");
                  CleanAfterError();
                  }

               PairCoeff_m[CombinedType][iread] =
                  dblstrf (&instr)/s->eunit;
               iread++;
               }