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 "cdboun.h"
#include "cdhouse.h"

/*  Place particles within Repeating Boundaries  */
int WrapParticles (Particle_t *a)
	{
	int idir;
	int ipart;
	int NumRepeat;
	int CntWrap;
	BOOLEAN IsParticleWrapped;
	double (*Position  )[NDIR] = (double (*)[NDIR]) a->cur;
	double  *RepeatDist			=							a->bcur;

	/*  Cannot wrap particles if no box was specified  */
	if (!a->IsInitializedBox)
		return 0;

	/*  Initialize number  of particles wrapped  */
	CntWrap = 0;
	LOOP (ipart, a->np)
		{

		/*  Set wrap flag  */
		IsParticleWrapped = FALSE;


		/*  Loop over directions  */
		LOOP (idir, NDIR)
			{

			/*  Test for repeating boundary conditions  */
			if ( ! a->surf [idir] )
				{


				/*  Translate ipart backward by repeat distance  */
				if (Position[ipart][idir] > RepeatDist[idir])
					{
					NumRepeat = Position[ipart][idir] / RepeatDist[idir];
					Position[ipart][idir] -= NumRepeat*RepeatDist[idir];
					if (a->BoxMotionAlgorithm!=BMA_NONE)
						CorrectDeriv (a, ipart, idir, (double) -NumRepeat);
					IsParticleWrapped = TRUE;
					}


				/*  Translate ipart forward by repeat distance	*/
				if (Position[ipart][idir] < 0.0)
					{
					/*  Corrected version 2.3.3 (was missing 1)  */
					NumRepeat = 1 - Position[ipart][idir] / RepeatDist[idir];
					Position[ipart][idir] += NumRepeat*RepeatDist[idir];
					if (a->BoxMotionAlgorithm!=BMA_NONE)
						CorrectDeriv (a, ipart, idir, (double) NumRepeat);
					IsParticleWrapped = TRUE;
					}

				}

			}

		/*  Count if particle was wrapped  */
		if (IsParticleWrapped)
			{
			CntWrap++;
			}

		}

	return CntWrap;
	}



/*  Correct particle derivatives for box motion  */
void CorrectDeriv (Particle_t *a, int ipart, int idir, double Factor)
	{
	int i;

	i = NDIR * ipart + idir;
	a->v [i] += Factor * a->BoxMotion[idir][1];
	a->c2[i] += Factor * a->BoxMotion[idir][2];
	a->c3[i] += Factor * a->BoxMotion[idir][3];
	a->c4[i] += Factor * a->BoxMotion[idir][4];
	a->c5[i] += Factor * a->BoxMotion[idir][5];
	}