xmd Code

/*
xmd - molecular dynamics for metals and ceramics
By Jonathan Rifkin <jon.rifkin@uconn.edu>
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
************************************************************************
*/
/*
2004-11-11

***EXPERIMENTAL***

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

/*
Function CalcNeighborPerAtom_Cell()

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.

*/
/*
2010-06-05  In CalcNeighborPerAtom_Cell(), call to CallBackFunction()
            was not done if NumNeigh==0, but that meant that 
				value of a->NeigborListIndex[] and a->NeighborListLength[]
				were unset and thus kept their previous value, which had
				resulting in an atom being treated as its own neighbor,
				which would then give r=0, and then cause problems
				with force calculation.
*/


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

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

#ifdef TEST_DRIVER
int CheckMem_g = 0;
int Debug_g = 1;
void CleanAfterError(void) {exit (1);}
#endif



/*
************************************************************************
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"
#include "hash.h"


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

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

#define NRSEQ 256

#define APPROX_PARTS_PER_CELL 16

/*  Number of first neighbor cells  */
#define NNEIGHCELLS 26

/*  Set CellSpan[] to (1+ERROR_FAC)*CellSpan[] to correct 
 *  for round-off */
#define ERROR_FAC 0.01

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

/* Hash functions  */
void init_rand_sequence(void);
double mrand (WORD a);
WORD makehash(WORD a);
WORD foldkey  (BYTE *key, int keylength);


/*
************************************************************************
Exported Functions
************************************************************************
*/
void CalcNeighborPerAtom_Cell
(
Coord_t  *Coord,
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)
) {
	FLOAT TempCoord;
	WORD  *NeighList = NULL;
	WORD  *NeighCandidates = NULL;
	double *NeighRadSqr = NULL;
	int    ncell = 0;
   WORD  NumNeigh;
	WORD  NumHash;
   int   idir;
	int   ihash, jhash;
   int   ipart, jpart;
	int   isort, jsort;
	int   icand, jcand;
	int   ix,iy,iz;
	int   dx,dy,dz;
	int   n_this_cell, n_all_cells;
	int   last;
	BOOLEAN NeighborInRange;
	BOOLEAN OutOfRange;
	FLOAT   MaxRadiusSqr;
	FLOAT   MinRadiusSqr;
	FLOAT   MaxCoord[NDIR];
	FLOAT   MinCoord[NDIR];
	FLOAT   CellSpan[NDIR];
   FLOAT   Radius;
	FLOAT   BoxLimit[NDIR];
	FLOAT xs;
	FLOAT ys;
	FLOAT zs;
	WORD  *thiscell;
	WORD  NumCell[NDIR];
	WORD  icells[NDIR];
	WORD  *PartToHash=NULL;
	WORD  *HashToSort=NULL;
	WORD  *PartIndex=NULL;
	BOOLEAN is_both_origin;
	WORD  NeighCell[26][NDIR];
	static BOOLEAN is_rand_init = FALSE;
	htable_t *hcell = NULL;
	helem_t  *t = NULL;
	BOOLEAN  is_new, is_found;
	int      *iptr;
	int      isize;

	/*  Initialize random sequence used by hash function  */
	if (!is_rand_init) init_rand_sequence();


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

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

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

	/*  Test box size - each repeating dimension must be at least 3*MaxRadius  */
	LOOP (idir,NDIR) {
		if (Repeat[idir] && 3.0*MaxRadius>=Box[idir])  {
			printf ("ERROR (CalcNeighborPerAtom_Cell):  When using the Cell based neighbor search\n");
			printf ("   the repeating dimensions must be at least 3 times the maximum cutoff distance.\n");
			printf ("   This distance is determined by the potential cutoff and the\n");
			printf ("   SearchCutoffRatio (see the NRANGE command).\n");
			CleanAfterError();
		}
	}

	/*  Variations on Min/Max Radius for quick comparisons  */
	MaxRadiusSqr = MaxRadius*MaxRadius;
	MinRadiusSqr = MinRadius*MinRadius;

	LOOP (idir,NDIR) 
	if (Repeat[idir])
		BoxLimit[idir] = Box[idir] - MaxRadius;

   /*   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;
		}
   }

	/*  Set MinCoord,MaxCoord if repeating boundary  */
	LOOP (idir,NDIR)
	if (Repeat[idir]) {
		MinCoord[idir] = 0.0;
		MaxCoord[idir] = Box[idir];
	}

	/*  Find X,Y,Z dimension of neighbor cells  */
	LOOP (idir,NDIR) {
      NumCell[idir] = (MaxCoord[idir] - MinCoord[idir])/((1 + ERROR_FAC)*MaxRadius);
		if (NumCell[idir]==0) NumCell[idir] = 1;
		CellSpan[idir] = (1 + ERROR_FAC) * (MaxCoord[idir] - MinCoord[idir])/NumCell[idir];
		ASSERT(CellSpan[idir]>0);
	}

	/*  Initialize hash tables  */
	ALLOCATE(PartToHash, WORD, NumPart);       /*  Gives cell hash value for particle      */
	ALLOCATE(PartIndex,  WORD, NumPart);       /*  Gives atom indexes sorted by cell hash  */

	/*  Allocation hash table  */
	/*  Hash table size can be smaller than number of items stored, items
	 *  will just stack-up within a hash table slot  */
	hcell = ht_init(NumPart/APPROX_PARTS_PER_CELL,0);

	/*  Find to which cell each atom belongs  */
	NumHash = 0;
	isort   = 0;
	LOOP (ipart,NumPart) {
		LOOP (idir,NDIR) 
			icells[idir] = (Coord[ipart][idir]-MinCoord[idir])/CellSpan[idir];
		/*  Find hash index for this cell  */
		is_new = ! ht_findkey(hcell,(BYTE *)&icells, sizeof(icells), (BYTE **) &iptr, (int *) &isize);
		/*  Give new cell a new hash index  */
		if (is_new) {
			ht_storekey (hcell, (BYTE *) icells, sizeof(icells), (BYTE *) &NumHash, sizeof(NumHash));
			ihash = NumHash;
			NumHash++;
		/*  This cell already assigned a hash index  */
		} else {
			ihash = *iptr;
		}
		/*  Record cell index for this particle  */
		PartToHash[ipart] = ihash;
		PartIndex [isort] = ipart;
		isort++;
	}

	/*  Make index to point to particles belonging to hash index  */
	sortux (PartToHash, PartIndex, NumPart);
	ALLOCATE(HashToSort, WORD, NumHash+1)
	for (isort=0;isort<NumPart;isort++) {
		ipart = PartIndex[isort];
		ihash = PartToHash[ipart];
		HashToSort[ihash+1] = isort+1;
	}

	/*  Now ready to find neighbors  */
	ALLOCATE (NeighList, WORD, NumPart)
	ALLOCATE (NeighCandidates, WORD, NumPart)
	ALLOCATE (NeighRadSqr, double, NumPart)

	/*  Loop over cells stored in hash  */
   ht_initwalk (hcell);
	while ((t=ht_getnext(hcell))) {
		/*  Store this cell's particles in temporary array  */
		ihash = * (int *) t->data;
		thiscell = (WORD *) t->key;
		n_this_cell = 0;
		for (isort=HashToSort[ihash]; isort<HashToSort[ihash+1]; isort++) {
			ipart = PartIndex[isort];
			NeighCandidates[n_this_cell] = ipart;
			n_this_cell++;
		}
		n_all_cells = n_this_cell;
		/*  Test all neighbors to thiscell[]  */
		for (dx=-1;dx<=1;dx++) {
			ix = thiscell[X] + dx;
			if (Repeat[X]) {
				if (ix<0) ix=NumCell[X]-1;
				else if (ix>=NumCell[X]) ix=0;
			} else {
				if (ix<0 || ix>=NumCell[X]) continue;
			}
		for (dy=-1;dy<=1;dy++) {
			iy = thiscell[Y] + dy;
			if (Repeat[Y]) {
				if (iy<0) iy=NumCell[Y]-1;
				else if (iy>=NumCell[Y]) iy=0;
			} else {
				if (iy<0 || iy>=NumCell[Y]) continue;
			}
		for (dz=-1;dz<=1;dz++) {
			/*  Skip origin cell  */
			if (dx==0 && dy==0 && dz==0) continue;
			iz = thiscell[Z] + dz;
			if (Repeat[Z]) {
				if (iz<0) iz=NumCell[Z]-1;
				else if (iz>=NumCell[Z]) iz=0;
			} else {
				if (iz<0 || iz>=NumCell[Z]) continue;
			}
			icells[X] = ix;
			icells[Y] = iy;
			icells[Z] = iz;
			/*  Skip unoccupied cell  */
			is_found = ht_findkey(hcell,(BYTE *)&icells, sizeof(icells), (BYTE **) &iptr, (int *) &isize);
			if (!is_found) continue;
			jhash = *iptr;
			/*  Two cells cannot give the same hash  */
			ASSERT(ihash!=jhash)
			/*  If using Unique pairs, don't include those cells whose
			 *  hash is less than origin cell */
			if (CalcUniquePairs && ihash>jhash) continue;
			for (jsort=HashToSort[jhash]; jsort<HashToSort[jhash+1]; jsort++) {
				jpart = PartIndex[jsort];
				ASSERT(n_all_cells<NumPart);
				NeighCandidates[n_all_cells] = jpart;
				n_all_cells++;
			}
		}}}

		/*  Build neighbor lists for origin cell particles  */
		LOOP(icand,n_this_cell) {
			/*  Only make neighbor list for selected atoms  */
			ipart = NeighCandidates[icand];
			if (UseSelect && !IS_SELECT2(Select,ipart)) continue;
			NumNeigh = 0;
			LOOP(jcand,n_all_cells) {
				/*  Don't do particle with itself  */
				if (icand==jcand) continue;
				jpart = NeighCandidates[jcand];
				/*  Only make neighbor list for selected atoms  */
				if (UseSelect && !IS_SELECT2(Select,jpart)) continue;
				/*  If using Unique pairs, don't include particles pairs
				 *  from origin cell where jpart<ipart  */
				is_both_origin = (jcand<n_this_cell);
				if (is_both_origin && CalcUniquePairs && ipart>jpart) continue;
				/*  Calculate X Separation  */
				xs = fabs(Coord[ipart][X]-Coord[jpart][X]);
				if (Repeat[X] && xs>BoxLimit[X]) xs = Box[X]-xs;
				/*  Test Y Separation  */
				if (xs>=MaxRadius) continue;

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

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

            /*  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_Cell): ");
                  printf ("Too Many neighbors.\n");
                  exit(1);
					}
				}
			}
			/*  Invoke callback function  */
			CallBackFunction (ipart, NeighList, NeighRadSqr, NumNeigh);
		}
	}

	/*  RETURN  STORAGE  */
	FREE (NeighRadSqr)
	FREE (NeighList)
	FREE (PartIndex)
	FREE (HashToSort)
	FREE (PartToHash)
	ht_free(hcell);
}


/*
************************************************************************
End of Functions
************************************************************************
*/







#ifdef TEST_DRIVER

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




/*
************************************************************************
Defines
************************************************************************
*/
#define MAXPART 256


/*
************************************************************************
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;

	/*  Print this particle set  */
#if 0
	printf ("ipart (%d) (%e,%e,%e)\n", ipart, Coord_m[ipart][X], Coord_m[ipart][Y], Coord_m[ipart][Z]);
	for (jlist=0; jlist<NumNeigh; jlist++) {
      jpart = NeighList[jlist];
		printf ("   jpart (%d) (%e,%e,%e)\n", jpart, Coord_m[jpart][X], Coord_m[jpart][Y], Coord_m[jpart][Z]);
	}
#endif	
	for (jlist=0; jlist<NumNeigh; jlist++) {
		jpart = NeighList[jlist];
		printf ("ipart jpart %4d %4d\n", ipart, jpart);
	}
#if 0
	printf ("ipart NumNeigh %4d %4d (%e,%e,%e)\n", 
			ipart, NumNeigh, 
			Coord_m[ipart][X], Coord_m[ipart][Y], Coord_m[ipart][Z]);
#endif


}

/*
************************************************************************
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;


   /*  Initialize coordinates  */
	NumPart = 0;
	for (i=0; i<6; i++)
	for (j=0; j<6; j++)
	for (k=0; k<6; k++)
		{
		Coord_m[NumPart][X] = 1.0e-8 * (0.50 + i);
		Coord_m[NumPart][Y] = 1.0e-8 * (0.50 + j);
		Coord_m[NumPart][Z] = 1.0e-8 * (0.50 + k);
		NumPart++;
#if 0
		Coord_m[NumPart][X] = 1.0e-8 * (0.25 + i);
		Coord_m[NumPart][Y] = 1.0e-8 * (0.25 + j);
		Coord_m[NumPart][Z] = 1.0e-8 * (0.25 + k);
		NumPart++;
		Coord_m[NumPart][X] = 1.0e-8 * (0.75 + i);
		Coord_m[NumPart][Y] = 1.0e-8 * (0.75 + j);
		Coord_m[NumPart][Z] = 1.0e-8 * (0.75 + k);
		NumPart++;
#endif
		}

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

	Box[X]    = Box[Y]    = Box[Z]    = 6.0e-8;
	Repeat[X] = Repeat[Y] = Repeat[Z] = FALSE;
	MinRadius = 0.0e-8;
	MaxRadius = 1.1e-8;
	ALLOCATE (Select, SEL_T, NumPart)

	APPLY_SELECT2(Select,16)
#if 0
	a->nsel++;
#endif
   UseSelect  = FALSE;



	/*  Call routine  */
	TotalNumNeigh_m = 0;
	CalcNeighborPerAtom_Cell
		(
		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