cdboun.c

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

#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];
	}