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


#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"

#ifndef FLOAT
#define FLOAT double
#endif

#define CELL_MAX 5
#define MAX_NEIGHBORS 200

void CalcNeighborPerAtom_Cell2
(
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)
)
	{
// Despite the use of NDIR, this routine is not set up for
// calculations in less than 3 dimensions.
	
	WORD  *NeighList = NULL;
	WORD  *IxSortToPart = NULL;
	FLOAT *NeighRadSqr = NULL;
	FLOAT NDIR_INV = 1.0 / (1.0*NDIR);
	WORD  NumNeigh;
	int   idir;
	int   ipart, jpart;
	FLOAT   MaxRadiusSqr;
	FLOAT   MinRadiusSqr;
	FLOAT   MaxCoord[NDIR];
	FLOAT   MinCoord[NDIR];
	FLOAT   RangeCoord[NDIR];
	FLOAT xs, xsw;
	FLOAT TempCoord;
	int d1, d2, i; 
	int Nd[NDIR], nd[NDIR], dn[NDIR];
	int hi[NDIR], lo[NDIR];
	/*  Repeat is global; Wrap is particular to one particle--
	is it in a position where it could be wrapped? */
	BOOLEAN Wrap[NDIR]; 
	/*  
	For some reason, gcc on AIX 5.2 didn't like the variable name 'hz'
	(?!?!?), gave non-sensical compile error
		cell2.c: In function `CalcNeighborPerAtom_Cell2':
		cell2.c:76: parse error before `100'
	Changing hz->hzz (and hx,hy accordingly just for aestetics)
	fixed it (?!?!?!).

	Jon 2005-01-07
	*/
	int gx, gy, gz, hxx, hyy, hzz;
	FLOAT rsqr;
	FLOAT *r1, *r2;
	int m, mm;
	int CubeRootN;
	int * S;
	int * SC;
	int * GP;
// int itest = 5323;
// int iDummy = 0;

	BOOLEAN bNoPBC = (Repeat[0] == 0) && (Repeat[1] == 0) && (Repeat[2] == 0);

// clock_t eam_nl_s3, eam_nl_f3;
// clock_t eam_nl_s4, eam_nl_f4;
// clock_t eam_nl_s5, eam_nl_f5;

// eam_nl_s3 = clock();
	//  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=0; 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];
		}


	//  Initialize variations on cutoff radius
	MaxRadiusSqr = MaxRadius*MaxRadius;
	MinRadiusSqr = MinRadius*MinRadius;

	CubeRootN = pow((FLOAT)(NumPart), NDIR_INV);

// Calculate size of the 3d grid; get ~ NAtoms grid cells.
	for (idir = 0; idir < NDIR; idir++)  {
		d1 = idir + 1;
		if (d1 > 2) d1 = d1 - 3;
		d2 = idir + 2;
		if (d2 > 2) d2 = d2 - 3;
		Nd[idir] = 1 + CubeRootN * (int)(0.5 + pow(RangeCoord[idir]*RangeCoord[idir]/RangeCoord[d1]/RangeCoord[d2], 0.333333) );
	}

	S =  (int *)calloc(Nd[0] * Nd[1] * Nd[2] * CELL_MAX, sizeof(int));
	SC = (int *)calloc(Nd[0] * Nd[1] * Nd[2], sizeof(int)); // Should be zero at allocation
	GP = (int *)calloc(3 * NumPart, sizeof(int));

// Assign each atom to a 3-d grid location
	for (ipart = 0; ipart < NumPart; ipart++) {

//    if (ipart == itest) {
//       iDummy++;
//    }
		
		for (idir = 0; idir < NDIR; idir++)
			nd[idir] = (int)((FLOAT)(Nd[idir] - 1) * (Coord[ipart][idir] - MinCoord[idir]) / RangeCoord[idir] );
	
		m = nd[0] + Nd[0] * nd[1] + Nd[0]*Nd[1] * nd[2];
		mm = m * CELL_MAX + SC[m];

		S[mm] = ipart;
		SC[m]++;
		if (SC[m] > CELL_MAX) {
			printf ("INTERNAL ERROR (CalcNeighborPerAtom):  ");
			printf ("Insufficient storage for particles in cell; increase CELL_MAX.\n");
			exit   (1);
		}

		// Store cell location for each atom...could have recalculated it from position.
		for (idir = 0; idir < NDIR; idir++)
			GP[3*ipart + idir] = nd[idir];
	
	} // end for ipart

	for (idir = 0; idir < NDIR; idir++)  {
		dn[idir] = (int)((FLOAT)(Nd[idir]) * MaxRadius / RangeCoord[idir] ) + 1;
	} // end for idir

// Find neighbors of each atom
	ALLOCATE (NeighList, WORD, MAX_NEIGHBORS);
	ALLOCATE (NeighRadSqr,  FLOAT, MAX_NEIGHBORS);

	for (ipart = 0; ipart < NumPart; ipart++) {

//    if (ipart == itest) {
//       iDummy++;
//    }
		
		// Default to non-periodic boundary conditions
		Wrap[0] = Wrap[1] = Wrap[2] = 0;
		// Now test for periodic boundary conditions
		if (!bNoPBC) {
			for (idir = 0; idir < 3; idir++)  {

				if (!Repeat[idir]) continue; // Can't wrap. Not periodic in this idir if it doesn't repeat

			// Test for wrap
				Wrap[idir] = ( Coord[ipart][idir] - MaxRadius < 0.) || (Coord[ipart][idir] + MaxRadius > Box[idir]);
			} // end for idir
		}

		r1 = Coord[ipart];
		NumNeigh = 0;

		for (idir = 0; idir < NDIR; idir++)  {
		// Locate this atom in the grid
			nd[idir] = GP[3 * ipart + idir];
		// Find region within grid to test for neighbors
			lo[idir] = nd[idir] - dn[idir];
			if (lo[idir] < 0)  {
				if (!Wrap[idir])
					lo[idir] = 0;
			}
			if (CalcUniquePairs) hi[idir] = nd[idir];
			else  {
				hi[idir] = nd[idir] + dn[idir];
				if (hi[idir] > Nd[idir]-1)  {
					if (!Wrap[idir])
						hi[idir] = Nd[idir]-1;
						}
				}
		} // end for idir


		for (gx = lo[0]; gx <= hi[0]; gx ++)  {
			hxx = gx;
			if (Wrap[0]) {
				if (hxx < 0)       hxx += Nd[0];
				if (hxx > Nd[0] -1)   hxx -= Nd[0];
			}
			for (gy = lo[1]; gy <= hi[1]; gy ++)  {
				hyy = gy;
				if (Wrap[1]) {
					if (hyy < 0)       hyy += Nd[1];
					if (hyy > Nd[1] -1)   hyy -= Nd[1];
				}
				for (gz = lo[2]; gz <= hi[2]; gz ++)  {
					hzz = gz;
					if (Wrap[2]) {
						if (hzz < 0)       hzz += Nd[2];
						if (hzz > Nd[2] -1)   hzz -= Nd[2];
					}
					m = hxx + Nd[0] * hyy + Nd[0]*Nd[1] * hzz;
					mm = m * CELL_MAX;

//             if (hxx == 0 && hyy == 0 && hzz ==1)  {
//                iDummy++;
//             }
					
					for (i = 0; i < SC[m]; i++) {
						jpart = S[mm+i];
						if (jpart == ipart) continue;
						if (CalcUniquePairs && jpart > ipart) continue;

						// Now censor to get those corner lurkers out.
						r2 = Coord[jpart];
			
						rsqr = 0.;
						for (idir = 0; idir < NDIR; idir++)  {
							xs = r2[idir] - r1[idir];
							if (Wrap[idir]) { // Means, might be a wrapped pair.
								if (r2[idir] > r1[idir]) xsw = xs - Box[idir];
								else xsw = -xs - Box[idir];
								if (fabs(xsw) < fabs(xs) ) xs = xsw;
							}
							rsqr += xs * xs;
						} // end for idir
						if (rsqr >= MaxRadiusSqr) continue;
			
						NeighList[NumNeigh] = jpart;
						NeighRadSqr[NumNeigh] = rsqr;
						NumNeigh++;

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

					} // end for i

				} // end for gz
			} // end for gy
		} // end for gx

		CallBackFunction (ipart, NeighList, NeighRadSqr, NumNeigh);

	} // end for ipart
	
	free (S);
	free (SC);
	free (GP);

// eam_nl_f5 = clock();
// eam_nl_d5 += double(eam_nl_f5 - eam_nl_s5) / CLOCKS_PER_SEC;
	//  RETURN  STORAGE
	FREE (NeighRadSqr)
	FREE (NeighList)
	}