cdsearch.c

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

/*
************************************************************************
Compile Switches
************************************************************************
*/
#define DEBUG
#undef  DEBUG

/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include "neigh.h"
#include "cdhouse.h"
#include "cdsubs.h"
#include "cdsearch.h"
#include "cdalloc.h"
#include "cdboun.h"

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

/*
************************************************************************
Defines
************************************************************************
*/
/*  Logical constants  */
#define BOOLEAN int
#define TRUE 1
#define FALSE 0

/*  Dimension constants  */
#define NDIR 3
#define X 0
#define Y 1
#define Z 2

/*  Constants for passing to neighbor routine CalcNeighborPerAtom  */
#define DO_NOT_USE_SELECT       FALSE
#define CALCULATE_UNIQUE_PAIRS  TRUE
#define MIN_RADIUS              0.0



/*
************************************************************************
Global Variables
************************************************************************
*/
extern Particle_t   *a;
extern Simulation_t *s;



/*
************************************************************************
Module-wide variables
************************************************************************
*/
static BOOLEAN TooManyNeighbors_m;
static long    NumNeighborSearchCalls_m = 0;
static time_t  ElapsedTime_m = 0;
static int     CurNumNeigh_m;


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

double GetMaxDisp (Particle_t *a);
void   TestCutoff (Particle_t *a);
void   NeighborCallBackFunction
(
WORD ipart,
WORD *NeighList,
double *NeighRadSqr,
WORD SizeList
);


/*
************************************************************************
Exported Functions
************************************************************************
*/

/*
Test if neighbor search needed, 
if so then 
	Allocate neighbor list if needed
	call UpdateNeighborList().
called by
  energy and force functions.
*/
void search_all  (Particle_t *a)
   {
   double   MaxAllowedDisp;
   BOOLEAN  MakeNeighborList;

   /*  DETERMINE IF CUTOFF IS TOO LARGE FOR BOX  */
   TestCutoff(a);

   /*  Start by assuming neighbor list calculation not needed  */
   MakeNeighborList = FALSE;

   /*  Check conditions which require calculating neighbor list  */

   /*  Determine if previous neighbor coordinate array exists  */
   if (a->ngh == NULL)
      {
      MakeNeighborList = TRUE;
      ALLOCATE (a->ngh, double, NDIR*a->NumPartAlloc)
      }

   /*  Check neighborlist flag  */
   else if (a->InvalidNeighborList)
      MakeNeighborList = TRUE;

   /*  See if particle displacements since last neighbor list exceeds limit  */
   else
      {
      MaxAllowedDisp = s->cutoff * (s->uratio-1)/2;
      if (GetMaxDisp(a) > MaxAllowedDisp)
         MakeNeighborList = TRUE;
      }

   /*  Neighbor List Needed  */
   if (MakeNeighborList)
		UpdateNeighborList(a);

	}

/*
Update Neighbor List, 
called by 
	cdsearch.c search_all()
	cdwrite.c  read_write()   after it calls writestate().
*/
void UpdateNeighborList(Particle_t *a)
	{
   double   Cutoff;
   time_t   InitialTime;
   time_t   FinalTime;
   int      icoord;
   BOOLEAN  RepeatingBoundary[NDIR];

	/*  Return if no neighbor list was allocated  */
	if (a->ngh==NULL)
		return;

	/*  Wrap particle positions  */
	WrapParticles (a);

	/*  Save current positions (for subsequent displacment testing  */
	LOOP (icoord, NDIR*a->np)
		a->ngh[icoord] = a->cur[icoord];

	/*   Calcualate effective cutoff  */
	Cutoff = s->cutoff * s->uratio;

	/*  Convert surface-on to repeating-on  */
	RepeatingBoundary[X] = !a->surf[X];
	RepeatingBoundary[Y] = !a->surf[Y];
	RepeatingBoundary[Z] = !a->surf[Z];

	/*  Set module-wide variables  */
	TooManyNeighbors_m = FALSE;
	CurNumNeigh_m = 0;
	NumNeighborSearchCalls_m++;
	time (&InitialTime);

	a->TotalNumNeighbors = 0;

	CalcNeighborPerAtom
	(
	(double (*)[NDIR]) a->cur,
	a->Selection,
	a->np,
	DO_NOT_USE_SELECT,
	MIN_RADIUS,
	Cutoff,
	a->bcur,
	RepeatingBoundary,
	a->CalcUniqueNeighbors,
	NeighborCallBackFunction
	);

	/*  Get elapsed time for neighbor calculation  */
	time (&FinalTime);
	ElapsedTime_m += difftime (FinalTime, InitialTime);


	/*  Save state of neighbor list  */
	a->InvalidNeighborList = FALSE;
	a->ng = CurNumNeigh_m;

	CHECK_HEAP
	}


/*
************************************************************************
Local Functions
************************************************************************
*/

/*
   Returns maximum displacements of particles
   relative to their positions prior to previous search
*/
double GetMaxDisp (Particle_t *a)
   {
   int i;
   double MaxDisp;
   double CurrentDisp;
   double (*c)[NDIR] = (double (*)[NDIR]) a->cur;
   double (*n)[NDIR] = (double (*)[NDIR]) a->ngh;
   double xdel;
   double ydel;
   double zdel;

   /*  Test for a->ngh allocated  */
   if (a->ngh == NULL)
      {
      printf ("INTENAL ERROR (max_disp):  a->ngh not allocated.\n");
      CleanAfterError();
      }

   MaxDisp = 0;
   LOOP (i, a->np)
      {
      xdel = c[i][X] - n[i][X];
      ydel = c[i][Y] - n[i][Y];
      zdel = c[i][Z] - n[i][Z];
      CurrentDisp =  xdel*xdel + ydel*ydel + zdel*zdel;
      if (CurrentDisp>MaxDisp)
         MaxDisp = CurrentDisp;
      }
   return sqrt(MaxDisp);
   }


/*  IF CUTOFF TOO LARGE WRITE MESSAGE AND RETURN 1  */
void TestCutoff (Particle_t *a)
   {
   int    i, repeat;
   double mindim=0.0;

   repeat = 0;
   for (i=0;i<3;i++)
      if (!a->surf[i])
         {
         if (repeat)
            {
            if (a->bcur[i]<mindim)
               mindim = a->bcur[i];
            }
         else
            {
            mindim = a->bcur[i];
            repeat = 1;
            }
         }

   if (repeat && mindim<2*s->cutoff*s->uratio)
      {
      printf ("\n");
      printf ("ERROR (TestCutoff): Potential cutoff distance (%le angs) exceeds half \n",
                  1e8*s->cutoff*s->uratio);
      printf ("                    the minimum repeat boundary length (%f angs).\n",
                  1e8*mindim);
      printf ("\n");
      CleanAfterError();
      }
   else
      return;
   }


void NeighborCallBackFunction
(
WORD    ipart,
WORD   *NeighList,
double *NeighRadSqr,
WORD    SizeList
)
   {
   int OldNumNeigh;
   int ilist;
   int ineigh;

   /*  Update number of current neighors  */
   OldNumNeigh   = CurNumNeigh_m;
   CurNumNeigh_m = OldNumNeigh + SizeList;

   /*  New allocation routine  */
   if (CurNumNeigh_m > a->NumNeighAlloc)
      ReallocateNeighbors (a, CurNumNeigh_m);


   /*  Add current neighbors to list  */
   a->NeighborListIndex [ipart]  = OldNumNeigh;
   a->NeighborListLength[ipart]  = SizeList;
   a->TotalNumNeighbors         += SizeList;
   for (ilist=0; ilist<SizeList; ilist++)
      {
      ineigh = OldNumNeigh + ilist;
      a->NeighborList[ineigh] = NeighList[ilist];
/*  Print out neighbors of 0  */
#ifdef DEBUG
      if (ipart==0 || NeighList[ilist]==0)
         {
         int jpart;
         double (*Coord)[NDIR];

         Coord = (double (*)[NDIR]) a->cur;
         jpart = NeighList[ilist];
         printf ("ipart jpart  %i  %i", ipart, jpart);
         if (ipart==0)
            printf ("  %12.3e  %12.3e  %12.3e\n",
               Coord[jpart][X], Coord[jpart][Y], Coord[jpart][Z]);
         else
            printf ("  %12.3e  %12.3e  %12.3e\n",
               Coord[ipart][X], Coord[ipart][Y], Coord[ipart][Z]);
         }
#endif
      }

   }

long GetNumNeighborSearchCalls (void)
   {
   return NumNeighborSearchCalls_m;
   }

double GetElapsedNeighborSearchTime (void)
   {
   return ElapsedTime_m;
   }