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 Files
************************************************************************
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "iomngr.h"
#include "strsub.h"
#include "parse.h"
#include "sortsub.h"
#include "particle.h"
#include "cdhouse.h"
#include "cdsubs.h"


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


/*
************************************************************************
Global Variables
************************************************************************
*/
extern Particle_t *a;
extern Simulation_t *s;
extern FILE 	 *out;
extern LIST 	 *inlist;


/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void nearneighbor (
		double, 
		double, 
		double, 
		BYTE *, 
		int, 
		double, 
		double, 
		BOOLEAN, 
		BOOLEAN, 
		BOOLEAN);

								

/*
************************************************************************
Exported Functions
************************************************************************
*/
void read_select	 (char *instr)
	{
	typedef enum { _AND, _OR, _NOT, _XOR, _ONLY } loptype;
	loptype lop;
	char	  *tokstr;
	BYTE    *tsel = NULL;
	BYTE     PrevSel;
	double	x1,y1,z1;
	double   x2,y2,z2;
	double   xc,yc,zc;
	double   r1=0, r2=0;
	double   xdel, ydel, zdel;
	double   nx,ny,nz;
	double   umin,umax,utmp, dbox;
	double   dot1, dot2, t;
	double   lo;
	double   hi;
	double   Value;
	int      idir;
	int		i,n,n1,n2,nskip,u;
	int		ntag, ntype, iset;
	double  (*c)[NDIR] = (double (*)[NDIR]) a->cur;
	int	  np = a->np;
	BOOLEAN AtomFound;
	BOOLEAN NoNeighbors;
	BOOLEAN OmitPoint;
	BOOLEAN UseAbsoluteValue;
	BOOLEAN UseMagnitude;
	BOOLEAN UseRadius = FALSE;
	int     index;
	BOOLEAN UseAllNeighbors;

	/*  GET FIRST TOKEN	*/
	tokstr = strhed (&instr);

	/*  Check for KEEP option  */
	if (!strcmpi("KEEP",tokstr))
		{
		tokstr = strhed (&instr);
		if       (!strcmpi("ON",tokstr))
			a->selkeep = TRUE;
		else if (!strcmpi("OFF",tokstr))
			a->selkeep = FALSE;
		else
			{
			IncrementNumberOfWarnings();
			printf 
				("WARNING:  Unknown option (%s) to SELECT KEEP\n",
				tokstr);
			}
		return;
		}
	
	/*  Check for non-zero number of particles  */
	if (a->np==0)
		return;

	/*  Allocate arrays  */
	ALLOCATE (tsel, BYTE, a->np)

	CHECK_HEAP

	/*  INITIALIZE tsel	*/
	LOOP (i, a->np)
		tsel[i] = FALSE;

	CHECK_HEAP


	/*  PARSE LOGICAL OPERATION  */
	if 	  (!strcmpi(tokstr,"AND" )) lop =  _AND ;
	else if (!strcmpi(tokstr,"OR"  )) lop =  _OR  ;
	else if (!strcmpi(tokstr,"NOT" )) lop =  _NOT ;
	else if (!strcmpi(tokstr,"XOR" )) lop =  _XOR ;
	else										 lop =  _ONLY;

	CHECK_HEAP

	/*  IF NOT DEFAULT lop PARSE COMMAND OPTION	*/
	if (lop!=_ONLY)  tokstr = strhed (&instr);

	CHECK_HEAP

	/*  SELECT ALL  */
	if (!strcmpi(tokstr,"ALL"))
		{
		LOOP (i, a->np)
			tsel[i] = TRUE;
		}

	/*  SELECT NONE  */
	else if (!strcmpi(tokstr,"NONE"))
		{
		LOOP (i, a->np)
			tsel[i] = FALSE;
		}

	/*  SELECT BOX  */
	else if (!strcmpi(tokstr,"BOX"))
		{
		x1 = ANGS*dblstrf (&instr);
		y1 = ANGS*dblstrf (&instr);
		z1 = ANGS*dblstrf (&instr);
		x2 = ANGS*dblstrf (&instr);
		y2 = ANGS*dblstrf (&instr);
		z2 = ANGS*dblstrf (&instr);
		LOOP (i,np)
			tsel[i] =
				(
				x1<=c[i][0] && c[i][0]<x2 &&
				y1<=c[i][1] && c[i][1]<y2 &&
				z1<=c[i][2] && c[i][2]<z2
				);
		}

	/*  SELECT EATOM  */
	else if (!strcmpi(tokstr,"EATOM"))
		{
		lo = dblstrf (&instr)/s->eunit;
		hi = dblstrf (&instr)/s->eunit;
	
		if (a->eatom!=NULL)
			{
			LOOP (i, a->np)
				{
				tsel[i] = (a->eatom[i]>=lo && a->eatom[i]<=hi);
				}
			}
		}

	/*  SELECT ELLIPSE  */
	else if (!strcmpi(tokstr,"ELLIPSE"))
		{
		x1 = ANGS*dblstrf (&instr);
		y1 = ANGS*dblstrf (&instr);
		z1 = ANGS*dblstrf (&instr);
		x2 = ANGS*dblstrf (&instr);
		y2 = ANGS*dblstrf (&instr);
		z2 = ANGS*dblstrf (&instr);
		LOOP (i, a->np)
			{
			xdel = (c[i][X]-x1)/x2;
			ydel = (c[i][Y]-y1)/y2;
			zdel = (c[i][Z]-z1)/z2;
			tsel[i] =  (xdel*xdel + ydel*ydel + zdel*zdel) < 1.0;
			}
		}

	/*  SELECT INDEX  */
	else if (!strcmpi(tokstr,"INDEX"))
		{
		n1 	= intstrf (&instr);
		n2 	= intstrf (&instr);
		nskip = intstrf (&instr);
		if (n1 < 1		)
			n1 = 1;
		if (n2 > np)
			n2=np;
		if (n2 == 0)
			n2 = n1;
		if (nskip==0)
			nskip = 1;
		LOOP (i, a->np)
			tsel[i] = FALSE;
		for (i=n1-1; i<n2; i+=nskip) tsel[i] = TRUE;
		}

	/*  NEW SELECT OPTIONS (28 Mar 1992)  */

	/*  SELECT ILIST  */
	else if (!strcmpi(tokstr,"ILIST"))
		{
		char str[256], *s;
		int Index;

		n = intstrf(&instr);
		while (n && m_gets_list (str, 256, inlist))
			{
			s = str;
			while (*s)
				{
				/*  Test for index in range  */
				Index = intstrf (&s) - 1;
				if (Index>=a->np || Index<0)
					{
					printf
						("ERROR (read_select): Index in ILIST command (%i)",
						Index);
					printf (" is out of range [0..%i]\n", a->np-1);
					CleanAfterError();
					}
				tsel[Index] = 1;
				n--;
				}
			}
		}

	/*  SELECT LAYER  */
	else if (!strcmpi(tokstr,"LAYER"))
		{
			/*  find layer parameters	*/
			tokstr = strhed (&instr);
			if 	  (!strcmpi(tokstr,"X"))
				u = 0;
			else if (!strcmpi(tokstr,"Y"))
				u = 1;
			else if (!strcmpi(tokstr,"Z"))
				u = 2;
			else {
				printf ("WARNING:  select layer direction unrecognized.");
				IncrementNumberOfWarnings();
				FREE (tsel)
				return;
			}
			/*  test if box is needed	*/
			if (!a->surf[u] && (!a->IsInitializedBox))
				{
				printf ("WARNING: cannot use select layer ");
				printf ("without first specifying box in %c direction.\n",'X'+u);
				IncrementNumberOfWarnings();
				FREE (tsel)	
				return;
				}

			/*  read number of layers	*/
			n = intstrf (&instr);
			if (n<=0)
				{
				printf 
					("WARNING: select layer must have positive number of layers.");
				IncrementNumberOfWarnings();
				FREE (tsel)
				return;
				}

			/*  test for existance of indivual layer specification  */
			if (!instr)
				{
				printf ("WARNING: select layer command did not specify an layers.");
				IncrementNumberOfWarnings();
				FREE (tsel)
				return;
				}

			/*  find umin and umax: 													 */
			/* 	  if repeating boundary (umin,umax)==(0,box)                            */
			/* 	  else                                  (umin,umax) == min,max coordinate  */
			if (!a->surf[u])
				{
				umin = 0;
				umax = a->bcur[u];
				}
			else 
				{
				umin = c[0][u];
				umax = umin;
				LOOP (i,np)
					{
					utmp = c[i][u];
					if (utmp<umin)
						umin = utmp;
					else if (utmp>umax)
						umax = utmp;
					}
				}
			dbox = (umax-umin) / n;
			while (*instr)
				{
				i	= intstrf(&instr);
				x1 = (i-1) * dbox + umin;
				x2 =	i	  * dbox + umin;
				LOOP(i,np)
					tsel[i] = tsel[i] ||
								(x1 <= c[i][u]) && (c[i][u] <= x2);
				}
			}

	/*  SELECT NEAR  */
	else if (!strcmpi(tokstr,"NEAR"))
		{

		/*  Initilize flag for no neighbors  */
		NoNeighbors = FALSE;

		/*  Test first token   */
		tokstr = strhed (&instr);

		/*  Identify first token, if RADIUS then not using n */
		if (!strcmpi(tokstr,"RADIUS")) {
			UseAllNeighbors = TRUE;
		/*  First token is not RADIUS, so read number of neighbors  */
		} else {
			UseAllNeighbors = FALSE;
			n      = intstrf(&tokstr);
			tokstr = strhed (&instr);
		}

		/*  Parse RADIUS options if present  */
		if (!strcmpi(tokstr,"RADIUS")) {
			r1 = ANGS*dblstrf(&instr);
			r2 = ANGS*dblstrf(&instr);
			UseRadius = TRUE;
			tokstr    = strhed (&instr);
		} else {
			UseRadius = FALSE;
		}

		/* Next, which atom's neighborhood do we search?  */
		
		/*  Find neighbors of indexed atom  */
		if (!strcmpi(tokstr,"INDEX"))
			{
			OmitPoint = TRUE;
			index = intstrf (&instr)-1;
			xc    = c[index][X];
			yc    = c[index][Y];
			zc    = c[index][Z];
			}

		/*  Find neighbors of first point in set  */
		else if (!strcmpi(tokstr,"SET"))
			{
			OmitPoint = TRUE;
			AtomFound = FALSE;
			iset   = intstrf (&instr);
			LOOP (i,np)
				{
				if (IS_SET(i,iset))
					{
					xc = c[i][X];
					yc = c[i][Y];
					zc = c[i][Z];
					AtomFound = TRUE;
					break;
					}
				}
			if (!AtomFound)
				{
				NoNeighbors = TRUE;
				}
			}

		/*  Find neighbors of point  */
		else if (!strcmpi(tokstr,"POINT"))
			{
			OmitPoint = FALSE;
			xc = dblstrf(&instr)*1e-8;
			yc = dblstrf(&instr)*1e-8;
			zc = dblstrf(&instr)*1e-8;
			}

		/*  Find neighbors of first selected atom  */
		else if ( !strcmpi(tokstr,"SEL") || IsBlank(tokstr) )
			{
			OmitPoint = TRUE;
			AtomFound = FALSE;
			LOOP (i,np)
				{
				if (IS_SELECT(i))
					{
					xc = c[i][X];
					yc = c[i][Y];
					zc = c[i][Z];
					AtomFound = TRUE;
					break;
					}
				}
			if (!AtomFound)
				{
				NoNeighbors = TRUE;
				}
			}

		/*  Unrecognized option  */
		else
			{
			printf ("WARNING:  Unrecognized SELECT NEAR option.\n");
			IncrementNumberOfWarnings();
			NoNeighbors = TRUE;
			}

		/*  Call near neighbor routine  */
		if (!NoNeighbors)
			nearneighbor (xc,yc,zc, tsel, n, r1, r2, 
					OmitPoint, UseRadius, UseAllNeighbors);
		}


	/*  SELECT PLANE  */
	else if (!strcmpi(tokstr,"PLANE"))
		{
		nx = dblstrf(&instr);
		ny = dblstrf(&instr);
		nz = dblstrf(&instr);
		x1 = 1e-8*dblstrf(&instr);
		y1 = 1e-8*dblstrf(&instr);
		z1 = 1e-8*dblstrf(&instr);
		x2 = 1e-8*dblstrf(&instr);
		y2 = 1e-8*dblstrf(&instr);
		z2 = 1e-8*dblstrf(&instr);
		if (nx==0 && ny==0 && nz==0)
			{
			printf ( "ERROR (read_select): plane normal is zero.");
			FREE (tsel)
			return;
			}
		dot1 = nx*x1 + ny*y1 + nz*z1;
		dot2 = nx*x2 + ny*y2 + nz*z2;
		if (dot1>dot2)
			{
			t	  = dot1;
			dot1 = dot2;
			dot2 = t;
			}
		LOOP (i,np)
			{
			t = nx*c[i][0] + ny*c[i][1] + nz*c[i][2];
			tsel[i] = (dot1<=t) && (t<=dot2);
			}
		}

	/*  SELECT SET  */
	else if (!strcmpi(tokstr,"SET"))
		{

		/*  Read set number  */
		iset = intstrf(&instr);
		if (iset<1 || iset>MAX_SET)
			{
			printf("ERROR (read_select):  Set");
			printf(" (%i) out of range [1,%i].", iset, MAX_SET);
			FREE (tsel)
			return;
			}

		/*  Test set  */
		LOOP (i, np)
			{
		   tsel[i] = IS_SET(i,iset);
			}
		}

	/*  SELECT TAG  */
	else if (!strcmpi(tokstr,"TAG"))
		{
		ntag = intstrf(&instr);
		if (a->tag!=NULL)
		LOOP (i, np)
			tsel[i] = (a->tag[i]==ntag);
		}

	/*  SELECT TYPE  */
	else if (!strcmpi(tokstr,"TYPE"))
		{
		ntype = intstrf(&instr)-1;
		LOOP(i, np)
			tsel[i] = (a->type[i]==ntype);
		}

	/*  SELECT VELOCITY  */
	else if (!strcmpi(tokstr,"VELOCITY"))
		{
		tokstr = strhed (&instr);

		/*  Check Absolute Value Flag  */
		UseAbsoluteValue = FALSE;

		if (!strcmpi("ABS",tokstr))
			{
			UseAbsoluteValue = TRUE;
			tokstr = strhed (&instr);
			}

		/*  Read direction / magnitude flag  */
		UseMagnitude = FALSE;
		if (!strcmpi("X",tokstr))
			{
			idir = X;
			}
		else if (!strcmpi("Y",tokstr))
			{
			idir = Y;
			}
		else if (!strcmpi("Z",tokstr))
			{
			idir = Z;
			}
		else if (!strcmpi("MAG",tokstr))
			{
			UseMagnitude = TRUE;
			}

		/*  
		Read velocity limits in units of cm/sec,
		convert to cm/time step
		*/
		lo = dblstrf (&instr)*s->dtime;
		hi = dblstrf (&instr)*s->dtime;
	
		if (a->v!=NULL)
			{
			LOOP (i, a->np)
				{
				if (UseMagnitude)
					{
					Value = sqrt
						(
						a->v[i*NDIR+X]*a->v[i*NDIR+X] +
						a->v[i*NDIR+Y]*a->v[i*NDIR+Y] +
						a->v[i*NDIR+Z]*a->v[i*NDIR+Z]
						);
					}
				else
					{
					Value = a->v[i*NDIR+idir];
					}
				if (UseAbsoluteValue)
					Value = fabs(Value);
				tsel[i] = (Value >= lo && Value <= hi);
				}
			}
		}

	/*   Unrecognized SELECT option  */
	else
		{
		ERROR_PREFIX
		printf ("Unrecognized SELECT option (%s).\n", tokstr);
		CleanAfterError();
		}


		CHECK_HEAP

	/*  COMBINE SELECTION WITH PREVIOUS SELECTION  */
	a->nsel = 0;
	LOOP (i, np)
		{

		/*  Store previous value of selection  */
		PrevSel = IS_SELECT(i);

		/*  Intialize to unselected  */
		REMOVE_SELECT(i)

		/*  Combine with current value of select  */
		switch (lop)
			{
			case _ONLY:
				if (tsel[i])
					APPLY_SELECT(i)
				break;
			case _AND:
				if (tsel[i] && PrevSel)
					APPLY_SELECT(i);
				break;
			case _OR:
				if (tsel[i] || PrevSel)
					APPLY_SELECT(i)
				break;
			case _NOT:
				if (!tsel[i])
					APPLY_SELECT(i)
				break;
			case _XOR:
				if (	(tsel[i] || PrevSel) &&
									!(tsel[i] && PrevSel) 	)
					APPLY_SELECT(i)
			
				break;
			}

		/*  Sum selected atoms  */
		if (IS_SELECT(i))
			a->nsel++;
		}

	CHECK_HEAP

	/*  Print number selected  */
	if (s->PrintInfo)
		printf ("***  NUMBER SELECTED %i\n", a->nsel);

	CHECK_HEAP

	/*  FREE STORAGE	*/
	FREE (tsel)
	}

	/*  END OF read_select	*/


void read_tag (char *instr) 
	{
	int i, itag;

	/*  DO IF (1) Argument string  AND (2) selected particles  */
	/*  Assign tag values if particle selected and tag value entered	*/
	if (a->nsel != 0	 &&	instr[0] != '\0' )
		{

		/*  Allocate array new	*/
		if (a->tag == NULL)
			{

			/*  Allocate space for tag array  */
			ALLOCATE (a->tag, BYTE, a->np)
			}

		/*  Read and assign tag value  */
		itag = intstrf (&instr);
		LOOP (i, a->np)
			if (IS_SELECT(i))
				a->tag[i] = itag;

		}
	}


void read_set (char *instr) 
	{
	unsigned int i;
	unsigned int iset;
	BOOLEAN addflag;
	BOOLEAN subflag;
	BOOLEAN clrflag;
	char *tokstr;

	tokstr = strhed (&instr);

	/*  PARSE OPTION	*/
	addflag = subflag = clrflag = FALSE;
	if 	  (!strcmpi(tokstr,"ADD"))
		addflag = TRUE;
	else if (!strcmpi(tokstr,"SUB"))
		subflag = TRUE;
	else if (!strcmpi(tokstr,"CLEAR"))
		clrflag = TRUE;
	else 
		{
		printf ("ERROR (read_set):  unknown option.\n");
		return;
		}

	/*  RETURN IF NO PARTICLES SELECTED  */
	CheckForNoneSelected();
	if (a->nsel==0)
		return;

	/*  RETURN IF NO SET ARRAY AND EITHER SUB OR CLEAR  */
	if (a->Selection==NULL && (subflag || clrflag))
		return;

	/*  SET BITS  */
	iset = intstrf(&instr);

	/*  Is set in allowed range?  */
	if (iset<1 || iset>MAX_SET)
		{
		printf ("ERROR (read_set):  set (%i) out of range [1,%i].\n",
				iset, MAX_SET);
		return;
		}

	/*  apply set to particles  */
	LOOP (i, a->np)
		{
		if (IS_SELECT(i)	|| clrflag) 
			{

			/*  Add particle to set  */
			if (addflag)
				{
				APPLY_SET(i, iset)
				}

			/*  Remove particle from set  */
			else if (subflag)
				REMOVE_SET(i, iset)

			/*  Clear set entirely  */
			else if (clrflag)
				REMOVE_SET(i, iset)
			}
		}
	}



void CheckForNoneSelected(void)
	{
	if (a->nsel==0)
		{
		IncrementNumberOfWarnings();
		printf ("WARNING:  No atoms are selected\n");
		}
	}



/*
************************************************************************
Local Functions
************************************************************************
*/
void nearneighbor 
(
double xc, 
double yc, 
double zc,
unsigned char *tsel, 
int nnear,
double r1,
double r2,
BOOLEAN OmitPoint,
BOOLEAN UseRadius,
BOOLEAN UseAllNeighbors
) 
	{
	unsigned i;
	unsigned u;
	unsigned *index = NULL;
   double   *rlist = NULL;
	double   rsq;
	double   r[NDIR];
	double   rd;
	double  (*c)[NDIR] = (double (*)[NDIR]) a->cur;
	double   r1sq;
	double   r2sq;
	int     NumAtom;
	int     ifirst;
	int     ncount;

	ALLOCATE (index, unsigned, a->np)
	ALLOCATE (rlist, double,   a->np)

	r[0] = xc;
	r[1] = yc;
	r[2] = zc;
	NumAtom = 0;
	LOOP (i,a->np)
		{
		rsq = 0;
		LOOP(u,NDIR)
			{
			rd = fabs(r[u] - c[i][u]);
			if (!a->surf[u] && 2*rd>a->bcur[u])
				rd = a->bcur[u] - rd;
			rsq = rsq + rd*rd;
			}
		if (!OmitPoint || rsq>0.0)
			{
			index[i] = i;
			rlist[i] = rsq;
			NumAtom++;
			}
		}

	/*  sort rlist  */
	sortdx (rlist, index, NumAtom);

	/*  Select nearest atoms   */

	/*  Find first neighbor particle. 
	 *  If using RADIUS command, first neighbor is at radius=r1, 
	 *  otherwise first neighbor is at radius=0
	 */
	if (UseRadius) 
		{
		/*  Unless we find otherwise, assume all particle within r1  */
		ifirst = NumAtom;
		r1sq = r1*r1;
		LOOP(i, NumAtom) 
			{
			if (rlist[index[i]] >= r1sq) 
				{
				ifirst = i;
				break;
				}
			} 
		}
	else 
		ifirst = 0;
	
	/*  Select particles until exceed count nnear or neighbor distance r2  */
	if (UseRadius) r2sq = r2*r2;
	ncount = 0;
	for (i=ifirst; i<NumAtom; i++) 
		{
		/*  End collection, radius exceeded  */
		if (UseRadius && rlist[index[i]]>r2sq)
			break;
		/*  End collection, desired count reached  */
		if (! UseAllNeighbors && ncount>=nnear)
			break;
		ncount++;
		tsel[index[i]] = 1;
		}


	/*  release storage	*/
	FREE (rlist)
	FREE (index)
	}