neigh.c

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

/*
************************************************************************
History
************************************************************************
*/
/*
31 Jan 1997
   Make new version of sort-style neighbor search
*/

/*
Function CalcNeighborPerAtom()

This function accepts a list of particles and a box, and outputs
a list of particle indices whose separation lies in range of
a cutoff.

It can calculate all the neighbors for each atom, resulting in
duplication of atom pairs, or only unique atoms pairs.  This
is controlled by the parameter CalcUniquePairs.

If UseSelect is TRUE, then only neighobrs between pairs of 
selected atoms are returned.

*/


/*
************************************************************************
Includes
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*  
File particle.h defines SEL_T, IS_SELECT2 and APPLY_SELECT2
*/
#include "particle.h"
#include "cdhouse.h"


/*
************************************************************************
Defines
************************************************************************
*/

/*
Allow compilation of either float or double
Default is double
*/
#ifndef FLOAT
#define FLOAT   double
#endif

/*  Square root of 1, 2 and 3 (rounded up)  */
#define SQRT1  1.0
#define SQRT2  1.414215
#define SQRT3  1.732052


/*
************************************************************************
Compiler Switches
************************************************************************
*/

/*
Compile code to test drive the function CalcNeighborPerAtom()
*/
#define TEST_DRIVER
#undef  TEST_DRIVER


/*
************************************************************************
Module-wide variables
************************************************************************
*/
static WORD    *_x;
static FLOAT   _xd, _td;
static FLOAT  (*_d)[NDIR];
static int      _j, _dir;


/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void sortdx3 (FLOAT (*c)[NDIR], int Dir, WORD  *Ix, int N);


/*
************************************************************************
Exported Functions
************************************************************************
*/
void CalcNeighborPerAtom
(
FLOAT    (*Coord)[NDIR],
SEL_T   *Select,
int      NumPart,
BOOLEAN  UseSelect,
FLOAT    MinRadius,
FLOAT    MaxRadius,
FLOAT    Box[NDIR],
BOOLEAN  Repeat[NDIR],
BOOLEAN  CalcUniquePairs,
void     (*CallBackFunction) (WORD, WORD *, double *, WORD)
)
   {
	WORD  *NeighList = NULL;
   WORD  *IxSortToPart = NULL;
	double *NeighRadSqr = NULL;
   WORD  NumNeigh;
   int   idir;
   int   ipart;
	int   jpart;
	int   isort;
   int   jsort;
	BOOLEAN NeighborInRange;
	BOOLEAN OutOfRange;
	FLOAT   MaxRadiusSqr;
	FLOAT   MinRadiusSqr;
	FLOAT   MaxRange;
	FLOAT   MinRange;
	FLOAT   MaxRadius2;
	FLOAT   MaxRadius3;
	FLOAT   MaxCoord[NDIR];
	FLOAT   MinCoord[NDIR];
	FLOAT   RangeCoord[NDIR];
   FLOAT   Radius;
	FLOAT   BoxLimit[NDIR];
   FLOAT TempCoord;
   int   MinDir;
   int   MaxDir;
   int   Hor;
	int   Ver;
   int   Dep;
	FLOAT MinLimit;
	int   IntervalStart;
	FLOAT xs;
	FLOAT ys;
	FLOAT zs;


	/*  Check Data Integrity  */
	if (NumPart==0  ||  MinRadius >= MaxRadius)
		return;

	if (Coord==NULL)
      {
      printf
      ("INTERNAL ERROR (CalcNeighborPerAtom):  Coordinate arrays not allocated.\n");
      exit(1);
		}

	if (UseSelect && Select==NULL)
		{
      printf
      ("INTERNAL ERROR (CalcNeighborPerAtom):  Select array not allocated.\n");
      exit(1);
      }

#if 0
	if (CalcUniquePairs && UseSelect)
		{
      printf ("INTERNAL ERROR (CalcNeighborPerAtom):  ");
      printf ("Cannot set both UseSelect and CalcUniquePairs simultaneously.\n");
      exit   (1);
      }
#endif


   /*   Find dimension with maximum range  */
	LOOP (idir, NDIR)
      {
		MaxCoord[idir] = MinCoord[idir] = Coord[0][idir];
      }

	for (ipart=1; ipart<NumPart; ipart++)
      {
		LOOP (idir, NDIR)
         {
         TempCoord = Coord[ipart][idir];
         if      (TempCoord < MinCoord[idir])
            MinCoord[idir] = TempCoord;
         else if (TempCoord > MaxCoord[idir])
				MaxCoord[idir] = TempCoord;
			}
      }

	LOOP (idir,NDIR)
      {
      RangeCoord[idir] = MaxCoord[idir] - MinCoord[idir];
      }

	LOOP (idir, NDIR)
		{
		BoxLimit[idir] = Box[idir] - MaxRadius;
		}

	/*  Initialize variations on cutoff radius  */
	MaxRadiusSqr = MaxRadius*MaxRadius;
	MinRadiusSqr = MinRadius*MinRadius;
	MaxRadius2   = sqrt(2.0) * MaxRadius;
	MaxRadius3   = sqrt(3.0) * MaxRadius;

	MinRange = RangeCoord[X];
   MaxRange = RangeCoord[X];
   MinDir   = X;
   MaxDir   = X;
   for (idir=Y; idir<=Z; idir++)
      {
      if (RangeCoord[idir]<MinRange)
         {
			MinDir = idir;
         MinRange = RangeCoord[idir];
         }
      else if (RangeCoord[idir]>MaxRange)
         {
			MaxDir = idir;
         MaxRange = RangeCoord[idir];
         }
      }

   /* Order direction Hor, Ver and Dep by size of range  */
   if (MinDir==MaxDir)
      {
      Hor = X;
      Ver = Y;
      Dep = Z;
		}
	else
      {
      Hor = MaxDir;
		Ver = X + Y + Z  - MinDir - MaxDir;
		Dep = MinDir;
      }


	ALLOCATE(IxSortToPart, WORD, NumPart)

	LOOP (ipart, NumPart)
		{
		IxSortToPart[ipart] = ipart;
		}


	/*  SORT BY CURRENT X DIRECTION  */
	sortdx3 (Coord, Hor, IxSortToPart, NumPart);


	/*  Search backward in sorted direction for furthest neighbor form particle 0  */
	if (Repeat[Hor] && 2*MaxRadius<Box[Hor])
		{
		MinLimit = Box[Hor] + Coord[IxSortToPart[0]][Hor] - MaxRadius;
		IntervalStart = NumPart-1;
		while
		(
		IntervalStart > 0 &&
		Coord[IxSortToPart[IntervalStart]][Hor] > MinLimit
		)
			{
			IntervalStart--;
			}
		}

	/*  Do not search, use first  particle  */
	else
		IntervalStart = 0;



	/*  Allocate storage for neighbors  */
	ALLOCATE (NeighList, WORD, NumPart)
	ALLOCATE (NeighRadSqr,  double, NumPart)

	/*  FIND NEIGHBORS  */
	LOOP (isort, NumPart)

	/*  Do neighbor for current particle if selected  */
	if  ( !UseSelect || IS_SELECT2(Select,IxSortToPart[isort]) )
		{

		/*  Get particle index  of isort  */
		ipart = IxSortToPart[isort];


		/*  Move start of interval forward until in range   */
		OutOfRange = TRUE;
		while (OutOfRange)
			{

			/*  Calculate x separation between current and test particle  */
			xs = fabs (Coord[ipart][Hor] - Coord[IxSortToPart[IntervalStart]][Hor]);
			if (Repeat[Hor] && xs>BoxLimit[Hor])
				xs = Box[Hor] - xs;

			/*  Is current starting particle in range  */
			OutOfRange =  (xs > MaxRadius  &&  IntervalStart!=isort);

			/*  Move start of interval forward if out of range  */
			if (OutOfRange)
				{
				IntervalStart++;
				if (IntervalStart >= NumPart)
					IntervalStart = 0;
				}

			}


      NumNeigh = 0;
      NeighborInRange = TRUE;

      jsort = IntervalStart;
		jpart = IxSortToPart[jsort];

      while (NeighborInRange)
         {

         /*
         For Half Search, end search if particles
         current and test particles are identical.
         */
         if (CalcUniquePairs)
				{

            /*
            If current particle and test particle identical,
            end current search
				*/
            NeighborInRange = (isort!=jsort);
            if (!NeighborInRange)
               continue;

            /*  Calculate X separation   */
            xs = fabs(Coord[ipart][Hor]-Coord[jpart][Hor]);
            if (Repeat[Hor] && xs>BoxLimit[Hor])
					xs = Box[Hor]-xs;
            }

			/*
         For Full search, end search if x separation
         greater than cutoff.
         */
         else
            {

            /*  Calculate X separation   */
				xs = fabs(Coord[ipart][Hor]-Coord[jpart][Hor]);
            if (Repeat[Hor] && xs>BoxLimit[Hor])
               xs = Box[Hor]-xs;

				/*
            If X separation out of range,
            then end neighbor search for current particle
            */
				if (xs>MaxRadius)
					{
               NeighborInRange = FALSE;
               continue;
               }

            }

         /*  
			Continue with this neighbor 
				- if particles not identical
				- if not using selection jsort selected
			*/
			if (isort!=jsort && (!UseSelect || IS_SELECT2(Select,jpart)))
            {

				/*  Calculate Y Separation  */
				ys = fabs(Coord[ipart][Ver]-Coord[jpart][Ver]);
				if (Repeat[Ver] && ys>BoxLimit[Ver])
					ys = Box[Ver]-ys;

				/*  Test Y Separation  */
				if (ys>=MaxRadius || xs+ys>=MaxRadius2)
					goto NextLoop;


				/*  Calculate Z Separation  */
            zs = fabs(Coord[ipart][Dep]-Coord[jpart][Dep]);
            if (Repeat[Dep] && zs>BoxLimit[Dep])
					zs = Box[Dep]-zs;

            /*  Test Z Separation  */
            if (zs>=MaxRadius || xs+ys+zs>=MaxRadius3)
					goto NextLoop;

            /*  Calculate radius of separation  */
				Radius = xs*xs + ys*ys + zs*zs;

            /*  Test radius of separation  */
            if (Radius>=MinRadiusSqr && Radius<MaxRadiusSqr)
               {

					/*  Add to neighbor list of current particle  */
               NeighList  [NumNeigh] = jpart;
					NeighRadSqr[NumNeigh] = Radius;
               NumNeigh++;

               /*  Test for internal consistency  */
               if (NumNeigh > NumPart)
                  {
                  printf ("INTERNAL ERROR (CalcNeighborPerAtom): ");
                  printf ("Too Many neighbors.\n");
                  exit(1);
						}
               }
            }

			/*
			Label: Jump here if current pair not elgible
			*/
			NextLoop:

         /*  Step to next X coordinate in sequence  */
			jsort++;
         if (jsort==NumPart)
            jsort=0;
			jpart = IxSortToPart[jsort];

         /*  Make sure jsort has gone round the whole list  */
         NeighborInRange &= (jsort!=IntervalStart);
         }


		/*  Call the Callback function with neighbor list  */
		CallBackFunction (ipart, NeighList, NeighRadSqr, NumNeigh);
		}

	/*  RETURN  STORAGE  */
	FREE (NeighRadSqr)
	FREE (NeighList)
	FREE (IxSortToPart)
	}


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



/*  INDEXED SORT - FLOAT[NDIR]  */

void _qsortdx3rec (int lo, int hi)
	{
	int i;
	 i=lo;
	_j=hi;
	_xd=_d[_x[(lo+hi)/2]][_dir];
	while (i<_j)
		{
		while (_d[_x[i]][_dir]<_xd)
			i++;
		while (_xd<_d[_x[_j]][_dir])
			_j--;
		if (i<=_j)
			{
			_td    = _x[ i];
			_x[ i] = _x[_j];
			_x[_j] = _td;
			i++;
			_j--;
			}
		}
	if (lo<_j ) _qsortdx3rec (lo, _j );
	if (i < hi) _qsortdx3rec (i , hi);
	}

void sortdx3 (FLOAT (*z)[NDIR], int dir, WORD *x, int n)
	{
	_d = z;
	_x = x;
	_dir = dir;
	_qsortdx3rec (0, n-1);
	}



/*
************************************************************************
End of routines
************************************************************************
*/







#ifdef TEST_DRIVER

/*
************************************************************************
                           TEST DRIVER
************************************************************************
*/




/*
************************************************************************
Defines
************************************************************************
*/
#define MAXPART 128


/*
************************************************************************
Module Variables
************************************************************************
*/
FLOAT  Coord_m[MAXPART][NDIR];
int    TotalNumNeigh_m = 0;


/*
************************************************************************
Call Back Function
************************************************************************
*/


void CallBackFunction_m 
(
WORD ipart, 
WORD *NeighList, 
double *NeighRadSqr, 
WORD NumNeigh
)
   {
   WORD jlist;
   FLOAT  r;
   FLOAT  dx;
   FLOAT  dy;
   FLOAT  dz;
   int   jpart;

	if (NumNeigh>0 && NeighList==NULL)
		{
		printf ("ERROR (CallBackFunction): NeighList is null.\n");
		exit(1);
		}

	TotalNumNeigh_m += NumNeigh;

	for (jlist=0; jlist<NumNeigh; jlist++)
		{

      jpart = NeighList[jlist];

      dx = fabs(Coord_m[ipart][X] - Coord_m[jpart][X]);
      dy = fabs(Coord_m[ipart][Y] - Coord_m[jpart][Y]);
      dz = fabs(Coord_m[ipart][Z] - Coord_m[jpart][Z]);

		if (dx>1)
			dx = 2-dx;
		if (dy>1)
			dy = 2-dy;
		if (dz>1)
			dz = 2-dz;

      r  = sqrt(dx*dx + dy*dy + dz*dz);

      printf ("%u  %u  %f\n", ipart, NeighList[jlist], r);

		}
	}

/*
************************************************************************
Main Routine
************************************************************************
*/


int main (int argc, char *argv[])
   {
   SEL_T   *Select = NULL;
   BOOLEAN UseSelect;
   FLOAT   MinRadius;
   FLOAT   MaxRadius;
   FLOAT   Box[NDIR];
   BOOLEAN Repeat[NDIR];
   int     NumPart;
   int     i,j,k;
   BOOLEAN CalcUniquePairs;

   CalcUniquePairs = FALSE;

	if (argc>1 && argv[1][0]=='h')
			CalcUniquePairs = TRUE;
	CalcUniquePairs = TRUE;


   /*  Initialize coordinates  */
	NumPart = 0;
	for (i=0; i<2; i++)
	for (j=0; j<2; j++)
	for (k=0; k<2; k++)
		{
		Coord_m[NumPart][X] = 0.25 + i;
		Coord_m[NumPart][Y] = 0.25 + j;
		Coord_m[NumPart][Z] = 0.25 + k;
		NumPart++;
		Coord_m[NumPart][X] = 0.75 + i;
		Coord_m[NumPart][Y] = 0.75 + j;
		Coord_m[NumPart][Z] = 0.75 + k;
		NumPart++;
		}

	if (NumPart>MAXPART)
		{
		printf ("ERROR (main):  Too many particles.\n");
		exit(1);
		}

	Box[X]    = Box[Y]    = Box[Z]    = 2.0;
	Repeat[X] = Repeat[Y] = Repeat[Z] = TRUE;
	MinRadius = 0.0;
	MaxRadius = 0.9;
	ALLOCATE (Select, SEL_T, NumPart)

	APPLY_SELECT2(Select,16)
	a->nsel++;
   UseSelect  = FALSE;


	/*  Call routine  */
	TotalNumNeigh_m = 0;
	CalcNeighborPerAtom
		(
		Coord_m,
		Select,
		NumPart,
		UseSelect,
		MinRadius,
		MaxRadius,
		Box,
      Repeat,
		CalcUniquePairs,
      CallBackFunction_m
      );

	printf ("TotalNumber of Neighbors = %i\n", TotalNumNeigh_m);
	printf ("Number of Atoms = %i\n", NumPart);
	/*  End  */
	return 0;
	}

#endif
/*
neigh2.c: In function `CalcNeighborPerAtom':
neigh2.c:297: parse error before `;'
neigh2.c: In function `Sort':
neigh2.c:441: array subscript is not an integer
*/