cdpairf.c

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

/*
************************************************************************
History
************************************************************************
*/
/*
Pre-history
  cpair.c  Pair potential

12 June 1998
  Use cdfunc.c to make provide an auxillary pair potential function
  for other potential models and to be used internally by this model

TODO
  (1)
  Make INIT function to intialize number of types and storage for
  functions and cutoffs.

   (2)  Calling function must inialize eatom[] and f[]
*/


/*
************************************************************************
Compile Switches
************************************************************************
*/


/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "cdhouse.h"
#include "strsub.h"
#include "parse.h"
#include "cdsubs.h"
#include "cdsearch.h"
#include "cdpairf.h"
#include "cdfunc.h"
#include "iomngr.h"

#ifdef __TURBOC__
#include <alloc.h>
#endif

/*
************************************************************************
Defines
************************************************************************
*/
#define POT_NONE 0
#define POT_PAIR 1


/*
************************************************************************
Macros
************************************************************************
*/
#define MAG_SQR(V)\
	((V)[X]*(V)[X] + (V)[Y]*(V)[Y] + (V)[Z]*(V)[Z])

#define IS_TYPE_OK(t) ((t)>=0 && (t)<MaxType_m)

/*  Calls evaluation functions for function types  */
#define F0(t,x) (PairFunction_m[t])->f0(PairFunction_m[t],x)
#define F1(t,x) (PairFunction_m[t])->f1(PairFunction_m[t],x)



/*
************************************************************************
Type definitions
************************************************************************
*/


/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
static BOOLEAN CalcSeparation
	(double *Coord1, double *Coord2, double *Vector, double Cutoff);
static double GetSeparationDir(double Separation, int idir);
static void   CalcCutoff  (void);
void ParsePot (int *, int *, char *, char **);



/*
************************************************************************
External Variables (Define in cdsubs.c)
************************************************************************
*/
extern Simulation_t *s;
extern LIST       *inlist;

BOOLEAN UseCheckFunc_g;

/*
************************************************************************
Module Variables
************************************************************************
*/
static int MaxType_m = 0;
static int MaxTypePair_m = 0;
static double      *Cutoff_m       = NULL;
static double      *CutoffSqr_m    = NULL;
static Function_t **PairFunction_m = NULL;

static BOOLEAN BoundRep_m[NDIR];
static double  Box_m     [NDIR];
static double  HalfBox_m [NDIR];


/*
************************************************************************
Exported Functions:  Interface to XMD Potential Calls
************************************************************************
*/



#pragma argsused
void pr_energy_list  (Particle_t *a, Simulation_t *s, int sflag)
	{
   int iatom;

   /*  NEIGHBOR LIST FOR ALL PARTICLES  */
	s->cutoff = PAIR_GetMaxCutoff();
   a->CalcUniqueNeighbors = TRUE;
   search_all (a);

   /*  Initialize sums  */
   a->epot = 0;
   LOOP (iatom, a->np)
      a->eatom[iatom] = 0.0;

   /*  Call Auxillary Pair Potential  */
   PAIR_CalcEnergy(a);

   /*  Sum pair potential  */
	LOOP (iatom, a->np)
      {
      a->epot  += a->eatom[iatom];
      }

   }


/*
Force function
*/

#pragma argsused
void pr_calcforc  (Particle_t *a, Simulation_t *s )
   {
   int iatom;
   int istress;

   /*  RECALCULATE NEIGHBORS  */
   s->cutoff = PAIR_GetMaxCutoff();
   a->CalcUniqueNeighbors = TRUE;
   search_all (a);

   /*  Initialize sums  */
   a->epot = 0;
   LOOP (iatom, a->np)
      {
      a->eatom[iatom] = 0.0;
      a->f[NDIR*iatom+X] = 0.0;
      a->f[NDIR*iatom+Y] = 0.0;
      a->f[NDIR*iatom+Z] = 0.0;
      }

   /*  Zero Internal Stresses  */
   LOOP (istress, NVOIGHT)
      a->Stress[istress] = 0.0;

   /*  Call Auxillary Pair Potential  */
   PAIR_CalcForce(a);

   /*  Sum pair potential  */
	LOOP (iatom, a->np)
      {
      a->epot  += a->eatom[iatom];
      }

   }



/*
************************************************************************
Exported Functions:  Generic Pair Potential Routines
************************************************************************
*/



void PAIR_CalcEnergy (Particle_t *a)
	{
	int idir;
	int iatom;
	int jatom;
	int ineigh;
	int ti;
	int tj;
	int tc;
	int StartInterval;
	double Radius;
	double RadiusSqr;
	double Disp[NDIR];
	double (*Coord)[NDIR];
	double Etemp;

   /*  Set pointer to coordinates  */
   Coord = (double (*)[NDIR]) a->cur;

   /*  Set box values  */
	LOOP (idir, NDIR)
      {
      BoundRep_m[idir]  = !a->surf[idir];
      Box_m[idir]     = a->bcur[idir];
      HalfBox_m[idir] = Box_m[idir]/2;
      }

	/*  Sum over atom pairs  */
   LOOP (iatom, a->np)
      {
      ti    = a->type[iatom];

      /* Skip this atom if type is out of range  */
      if (ti>=MaxType_m)
         continue;

      StartInterval = a->NeighborListIndex[iatom];
      LOOP (ineigh, a->NeighborListLength[iatom])
         {
         jatom = a->NeighborList[StartInterval+ineigh];
			/*
         If neighbor list has duplicate neighbors entries,
         (iatom<jatom) and (iatom>jatom)
         ignore the second (iatom>jatom)
         */
         if (!a->CalcUniqueNeighbors && iatom>jatom)
            continue;

         tj    = a->type[jatom];

         /* Skip this atom if type is out of range  */
         if (tj>=MaxType_m)
            continue;

         /*  Get unique type for pair of atom types  */
         tc = COMBINED_TYPE(ti,tj);

         ASSERT(tc<MaxTypePair_m)

         /*  Calculate separation between two atoms
         while accounting for repeating boundary conditions (if present)
         returns false if atoms are out of range
         */
         if
         (
         !CalcSeparation
				(
            Coord[iatom],
            Coord[jatom],
            Disp,
            Cutoff_m[tc]
            )
         )
            continue;

			RadiusSqr = MAG_SQR(Disp);

			if (RadiusSqr < CutoffSqr_m[tc])
            {

				Radius = sqrt(RadiusSqr);
				Etemp  = 0.5*F0(tc,Radius);
				a->eatom[iatom] += Etemp;
            a->eatom[jatom] += Etemp;
            }
         }
      }
   }


void PAIR_CalcForce (Particle_t *a)
	{
	int idir;
   int iatom;
   int jatom;
   int ineigh;
   int ti;
   int tj;
	int tc;
	int StartInterval;
	double Radius;
	double RadiusSqr;
	double Disp[NDIR];
	double (*Coord)[NDIR];
	double (*Force)[NDIR];
	double Epair;
	double forc;
	double tfx;
   double tfy;
   double tfz;


	/*  Set pointer to coordinates  */
   Coord = (double (*)[NDIR]) a->cur;
   Force = (double (*)[NDIR]) a->f;

   /*  Set box values  */
	LOOP (idir, NDIR)
      {
      BoundRep_m[idir]  = !a->surf[idir];
      Box_m[idir]     = a->bcur[idir];
      HalfBox_m[idir] = Box_m[idir]/2;
      }

   /*  Sum over atom pairs  */
   LOOP (iatom, a->np)
      {
      ti    = a->type[iatom];

      /* Skip this atom if type is out of range  */
      if (ti>=MaxType_m)
         continue;

      StartInterval = a->NeighborListIndex[iatom];
      LOOP (ineigh, a->NeighborListLength[iatom])
         {
         jatom = a->NeighborList[StartInterval+ineigh];
         /*
         If neighbor list has duplicate neighbors entries,
         (iatom<jatom) and (iatom>jatom)
         ignore the second (iatom>jatom)
         */
         if (!a->CalcUniqueNeighbors && iatom>jatom)
            continue;

         tj    = a->type[jatom];

         /* Skip this atom if type is out of range  */
         if (tj>=MaxType_m)
            continue;

         /*  Get unique type for pair of atom types  */
         tc = COMBINED_TYPE(ti,tj);

         ASSERT(tc<MaxTypePair_m)

         /*  Calculate separation between two atoms
         while accounting for repeating boundary conditions (if present)
         returns false if atoms are out of range
         */
         if
         (
         !CalcSeparation
            (
            Coord[iatom],
            Coord[jatom],
            Disp,
            Cutoff_m[tc]
            )
         )
            continue;

			RadiusSqr = MAG_SQR(Disp);

			if (RadiusSqr < CutoffSqr_m[tc])
            {

				Radius = sqrt(RadiusSqr);
				Epair  = 0.5*F0(tc,Radius);
				forc   = F1(tc,Radius)/Radius;
				a->eatom[iatom] += Epair;
				a->eatom[jatom] += Epair;
				tfx = forc*Disp[X];
            tfy = forc*Disp[Y];
            tfz = forc*Disp[Z];
            Force[jatom][X] -= tfx;
            Force[jatom][Y] -= tfy;
            Force[jatom][Z] -= tfz;
            Force[iatom][X] += tfx;
            Force[iatom][Y] += tfy;
            Force[iatom][Z] += tfz;