cdfill.c

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

         /*  Scale point to fit within margin  */
         MatrixVectorMult (ScaledPoint, ScaleMatrix, Point);

         LOOP (idir, NDIR)
            ScaledPoint[idir] += CenterPoint[idir];

         /*  Test point to see if in boundary  */
         if (BoundaryType_m==FILL_BOX)
            {
            UseThisPoint = IsWithinMarginBOX(ScaledPoint);
            }

         else if (BoundaryType_m==FILL_CYLINDER)
            {
            UseThisPoint = IsWithinMarginCYLINDER(ScaledPoint);
            }

         else if (BoundaryType_m==FILL_SPHERE)
            {
            UseThisPoint = IsWithinMarginSPHERE(ScaledPoint);
            }

         if (UseThisPoint)
            {
            AddParticle (a, FillType_m[ipart], ScaledPoint);
				NumNewPart++;
            }
         }
      }

   FREE (RealFillPart)

	return NumNewPart;

   }



static void AddParticle (Particle_t *a, int Type, Coord_t Point)
   {
   int ipart;

   ipart = a->np;
   ReallocateParticle (a, a->np+1);
   a->np++;
   a->cur [ipart*NDIR + X] = Point[X];
   a->cur [ipart*NDIR + Y] = Point[Y];
   a->cur [ipart*NDIR + Z] = Point[Z];
   a->type[ipart]          = (BYTE) Type;
	if (a->tag!=NULL)
		a->tag [ipart]          = (BYTE) 0;

   /*  Unspecified required fields set to 0  */
   a->Selection[ipart]     = 0;

	/*  Select this particle  */
	APPLY_SELECT2(a->Selection, ipart)

   /*  Mark coordinates as initialized  */
   a->IsInitializedCoord   = TRUE;

   }



static void GetLimits(int Limits[NVECT][2], BOUNDARY_TYPE BoundaryType)
	{
	int ivect;
	int icorner;
   int NumCorners;
   double   FloatLimits[NVECT][2];
   Coord_t  Decompose;
   Coord_t *Corners = NULL;
   BOOLEAN  Uninitialized = TRUE;


   if (BoundaryType==FILL_BOX)
      {
      GetCornersBOX(&Corners, &NumCorners);
		}
   else if (BoundaryType==FILL_CYLINDER)
      {
      GetCornersCYLINDER(&Corners, &NumCorners);
      }
   else if (BoundaryType==FILL_SPHERE)
      {
      GetCornersSPHERE(&Corners, &NumCorners);
      }

   ASSERT (Corners!=NULL)


   LOOP (icorner, NumCorners)
      {

      /*  Express box corner as multiples of cell vectors  */
      MatrixVectorMult
         (Decompose, RealToOrientCell_m, Corners[icorner]);

      /*  Save min and max multiples of vectors  */
      if (Uninitialized)
         {
         Uninitialized = FALSE;
			FloatLimits[U][0] = FloatLimits[U][1] = Decompose[U];
         FloatLimits[V][0] = FloatLimits[V][1] = Decompose[V];
         FloatLimits[W][0] = FloatLimits[W][1] = Decompose[W];
         }
      else
         {
         LOOP (ivect, NVECT)
            {
            if (Decompose[ivect] < FloatLimits[ivect][0])
               FloatLimits[ivect][0] = Decompose[ivect];
            if (Decompose[ivect] > FloatLimits[ivect][1])
               FloatLimits[ivect][1] = Decompose[ivect];
            }
         }
      }

   /*
   Convert float limits to integer limits
   (Add/subtract 2 to compensate for
    - roundoff error
    - alignment within cell
   )
   */
   LOOP (ivect, NVECT)
      {
      Limits[ivect][0] = FloatLimits[ivect][0]-2;
      Limits[ivect][1] = FloatLimits[ivect][1]+2;
      }


   /*  Free corners  */
   FREE (Corners)
	}




/*
------------------------------------------------------------------------
Box Specific routines
------------------------------------------------------------------------
*/
static BOOLEAN IsWithinMarginBOX (Coord_t Point)
   {
   return
      Point[X] >= Box_m[0][X]+Margin_m  &&
      Point[X] <  Box_m[1][X]-Margin_m  &&
      Point[Y] >= Box_m[0][Y]+Margin_m  &&
      Point[Y] <  Box_m[1][Y]-Margin_m  &&
      Point[Z] >= Box_m[0][Z]+Margin_m  &&
      Point[Z] <  Box_m[1][Z]-Margin_m;
   }


static void GetCornersBOX (Coord_t **Corners, int *NumCorners)
   {
   int i;
   int j;
   int k;
   int icorner;

   *NumCorners = 8;
   if (*Corners!=NULL)
      FREE (*Corners);
   ALLOCATE (*Corners, Coord_t, *NumCorners)

   /*  Loop over all 8 box corners  */
   icorner = 0;
   LOOP (i, 2)
   LOOP (j, 2)
   LOOP (k, 2)
      {

      /*  Current box corner  */
      (*Corners)[icorner][X] = Box_m[i][X];
      (*Corners)[icorner][Y] = Box_m[j][Y];
      (*Corners)[icorner][Z] = Box_m[k][Z];
      icorner++;
      }

   }


/*
Find appropriate scale matrix and scale center
*/
static void GetMarginScaleBOX (Matrix_t Scale, Coord_t Center)
   {
   int idir;
   int jdir;
   double Length;
   BOOLEAN IsError = FALSE;

   /*  Initialize off diagonal terms to zero  */
   LOOP (idir, NDIR)
   LOOP (jdir, NDIR)
      Scale[idir][jdir] = 0.0;

   LOOP (idir, NDIR)
      {
      Center[idir] = 0.5 * (Box_m[0][idir] + Box_m[1][idir]);
      Length = (Box_m[1][idir] - Box_m[0][idir]);
      if (Length-2*Margin_m <= 0)
         {
         Scale[idir][idir] = 1.0;
         IsError = TRUE;
         }
      else
         {
         Scale[idir][idir] = (Length-2*Margin_m)/Length;
         if (Scale[idir][idir]<0.5)
            {
            Scale[idir][idir] = 1.0;
            IsError = TRUE;
            }
         }
      }

   /*  Print warning  */
   if (IsError)
      {
      printf
         ("WARNING:  Boundary margin is too large.  Will ignore it.\n");
      INCREMENT_WARNING
      }

   }




/*
------------------------------------------------------------------------
Cylinder Specific routines
------------------------------------------------------------------------
*/
static BOOLEAN IsWithinMarginCYLINDER (Coord_t Point)
   {
   Coord_t CenteredPoint;
   Coord_t OrientedPoint;
   double  Radius;

   /*
   Translate point to cylindrical coordinates
   (Cylinder axis tranlates to Z)
   */
   CenteredPoint[X] = Point[X] - CylCenter_m[X];
   CenteredPoint[Y] = Point[Y] - CylCenter_m[Y];
   CenteredPoint[Z] = Point[Z] - CylCenter_m[Z];
   MatrixVectorMult (OrientedPoint, CylRealToOrient_m, CenteredPoint);

   /* Test point radius  */
   Radius =
      sqrt
      (
      OrientedPoint[X]*OrientedPoint[X] +
      OrientedPoint[Y]*OrientedPoint[Y]
      );

   if (Radius>=CylRadius_m-Margin_m)
      return FALSE;

   /*  Test z projection  */
   if
      (
      OrientedPoint[Z] <  -0.5*CylLength_m+Margin_m ||
      OrientedPoint[Z] >=  0.5*CylLength_m-Margin_m
      )
      return FALSE;

   /*  Point is inside  */
   return TRUE;
   }


static void GetCornersCYLINDER (Coord_t **Corners, int *NumCorners)
   {
   int icorners;
   double a;
   double b;
   double c;
   Coord_t TempCorn[8];

   *NumCorners = 8;

   /*  Free previous allocation  */
   if (*Corners!=NULL)
      FREE (*Corners)

   /*  Allocate local variable  */
   ALLOCATE (*Corners, Coord_t, *NumCorners)

   /*  Set corners in cylinderical coordinates  */
   a =     CylRadius_m;
   b =     CylRadius_m;
   c = 0.5*CylLength_m;
   TempCorn[0][X] = -a;   TempCorn[0][Y] = -b;   TempCorn[0][Z] = -c;
   TempCorn[1][X] =  a;   TempCorn[1][Y] = -b;   TempCorn[1][Z] = -c;
   TempCorn[2][X] = -a;   TempCorn[2][Y] =  b;   TempCorn[2][Z] = -c;
   TempCorn[3][X] =  a;   TempCorn[3][Y] =  b;   TempCorn[3][Z] = -c;
   TempCorn[4][X] = -a;   TempCorn[4][Y] = -b;   TempCorn[4][Z] =  c;
   TempCorn[5][X] =  a;   TempCorn[5][Y] = -b;   TempCorn[5][Z] =  c;
   TempCorn[6][X] = -a;   TempCorn[6][Y] =  b;   TempCorn[6][Z] =  c;
   TempCorn[7][X] =  a;   TempCorn[7][Y] =  b;   TempCorn[7][Z] =  c;

   /*  Translate from cylinder orientation about cylinder center  */
   LOOP (icorners, *NumCorners)
      {
      MatrixVectorMult
         ( (*Corners)[icorners], CylOrientToReal_m, TempCorn[icorners]);


      /*  Translate corners relative to origin (add cylinder center)  */
      (*Corners)[icorners][X] += CylCenter_m[X];
      (*Corners)[icorners][Y] += CylCenter_m[Y];
      (*Corners)[icorners][Z] += CylCenter_m[Z];
      }


   }


/*
Find appropriate scale matrix and scale center
*/
static void GetMarginScaleCYLINDER (Matrix_t Scale, Coord_t Center)
	{
	int idir;
	int jdir;
	Matrix_t TempScale;
	BOOLEAN IsError = FALSE;

   /*  Initialize off diagonal terms to zero  */
   LOOP (idir, NDIR)
   LOOP (jdir, NDIR)
      Scale[idir][jdir] = 0.0;

   /*  Copy cylinder center  */
   LOOP (idir, NDIR)
      {
      Center[idir] = CylCenter_m[idir];
      }

   Scale[0][0] = (CylRadius_m-1.0*Margin_m)/CylRadius_m;
   Scale[1][1] = Scale[0][0];
   Scale[2][2] = (CylLength_m-2.0*Margin_m)/CylLength_m;
   if (Scale[0][0]<=0.0 || Scale[1][1]<=0.0 || Scale[2][2]<=0.0)
      {
      IsError=TRUE;
      }
   if (Scale[0][0]<=0.5 || Scale[1][1]<=0.5 || Scale[2][2]<=0.5)
      {
      IsError=TRUE;
      }

   /*  Print warning  */
   if (IsError)
      {
      Scale[0][0] = Scale[1][1] = Scale[2][2] = 1.0;
      printf
         ("WARNING:  Boundary margin is too large.  Will ignore it\n");
      INCREMENT_WARNING
      return;
      }

   /*  Re-orient scale to work in cylinder orientation  */
   MatrixMatrixMult (TempScale, Scale, CylRealToOrient_m);
   MatrixMatrixMult (Scale, CylOrientToReal_m, TempScale);

	}


/*
------------------------------------------------------------------------
Sphere Specific routines
------------------------------------------------------------------------
*/
static BOOLEAN IsWithinMarginSPHERE (Coord_t Point)
   {
   Coord_t CenteredPoint;
   double  RadiusSqr;

   /*  Translate point to sphere coordinates  */
   CenteredPoint[X] = Point[X] - SphCenter_m[X];
   CenteredPoint[Y] = Point[Y] - SphCenter_m[Y];
   CenteredPoint[Z] = Point[Z] - SphCenter_m[Z];
   RadiusSqr =
      CenteredPoint[X]*CenteredPoint[X] +
      CenteredPoint[Y]*CenteredPoint[Y] +
      CenteredPoint[Z]*CenteredPoint[Z];

   /*  Test radius against original sphere radius */
	if (RadiusSqr >= SphRadiusSqr_m)
		return FALSE;

	/*  Test radius against sphere radius modified by margin  */
	if (sqrt(RadiusSqr) >= SphRadius_m - Margin_m)
		return FALSE;

   /*  Point is inside  */
   return TRUE;
   }


static void GetCornersSPHERE (Coord_t **Corners, int *NumCorners)
   {
   int icorners;
   double a;
   double b;
   double c;
   Coord_t TempCorn[8];

   *NumCorners = 8;

   /*  Free previous allocation  */
   if (*Corners!=NULL)
      FREE (*Corners)

   /*  Allocate local variable  */
   ALLOCATE (*Corners, Coord_t, *NumCorners)

   /*  Set corners in cylinderical coordinates  */
   a =     SphRadius_m;
   b =     SphRadius_m;
   c =     SphRadius_m;
   TempCorn[0][X] = -a;   TempCorn[0][Y] = -b;   TempCorn[0][Z] = -c;
   TempCorn[1][X] =  a;   TempCorn[1][Y] = -b;   TempCorn[1][Z] = -c;
   TempCorn[2][X] = -a;   TempCorn[2][Y] =  b;   TempCorn[2][Z] = -c;
   TempCorn[3][X] =  a;   TempCorn[3][Y] =  b;   TempCorn[3][Z] = -c;
   TempCorn[4][X] = -a;   TempCorn[4][Y] = -b;   TempCorn[4][Z] =  c;
   TempCorn[5][X] =  a;   TempCorn[5][Y] = -b;   TempCorn[5][Z] =  c;
   TempCorn[6][X] = -a;   TempCorn[6][Y] =  b;   TempCorn[6][Z] =  c;
   TempCorn[7][X] =  a;   TempCorn[7][Y] =  b;   TempCorn[7][Z] =  c;

   /*  Translate from cylinder orientation about cylinder center  */
   LOOP (icorners, *NumCorners)
      {

      /*  Translate corners relative to origin (add cylinder center)  */
      (*Corners)[icorners][X] = TempCorn[icorners][X] + SphCenter_m[X];
      (*Corners)[icorners][Y] = TempCorn[icorners][Y] + SphCenter_m[Y];
      (*Corners)[icorners][Z] = TempCorn[icorners][Z] + SphCenter_m[Z];
      }


   }


/*
Find appropriate scale matrix and scale center
*/
static void GetMarginScaleSPHERE (Matrix_t Scale, Coord_t Center)
   {
   int idir;
   int jdir;
   BOOLEAN IsError = FALSE;

   /*  Copy cylinder center  */
   LOOP (idir, NDIR)
      {
      Center[idir] =SphCenter_m[idir];
      }


   /*  Initialize off diagonal terms to zero  */
   LOOP (idir, NDIR)
   LOOP (jdir, NDIR)
      Scale[idir][jdir] = 0.0;
   Scale[0][0] = (SphRadius_m-Margin_m)/SphRadius_m;
   Scale[1][1] = (SphRadius_m-Margin_m)/SphRadius_m;
   Scale[2][2] = (SphRadius_m-Margin_m)/SphRadius_m;

   if (Scale[0][0]<=0.0 || Scale[1][1]<=0.0 || Scale[2][2]<=0.0)
      {
      IsError=TRUE;
      }
   if (Scale[0][0]<=0.5 || Scale[1][1]<=0.5 || Scale[2][2]<=0.5)
      {
      IsError=TRUE;
      }

   /*  Print warning  */
   if (IsError)
      {
      Scale[0][0] = Scale[1][1] = Scale[2][2] = 1.0;
      printf
         ("WARNING:  Boundary margin is too large.  Will ignore it\n");
      INCREMENT_WARNING
      return;
      }

   }




/*
------------------------------------------------------------------------
Test Driver
------------------------------------------------------------------------
*/
#ifdef TEST_DRIVER
int main()
   {
   Coord_t Point;
   FILE *infile;
   char InputStr[NSTR];
   char *InputPtr;
   char *TokenPtr;

   infile = stdin;
   a = palloc (8192,16);

   printf ("np = %i\n", a->np);
   Point[X] = 4.0*1e-8;
   Point[Y] = 4.0*1e-8;
   Point[Z] = 8.0*1e-8;
#if 0
   AddParticle (a, 0, Point);
   printf ("np = %i\n", a->np);
#endif

   /*  Read data  */
   while (!feof(infile))
      {
      fgets (InputStr, NSTR, infile);
      if (feof(infile))
         break;
      InputPtr = InputStr;
      TokenPtr = strhed (&InputPtr);
      if (!strcmpi(TokenPtr ,"FILL"))
         read_fill (InputPtr);
      }

   printf ("np = %i\n", a->np);
   PrintCoord (a);

   return 0;
   }


void CleanAfterError(void)
   {
   exit(1);
   }

void IncrementNumberOfWarnings()
   {
   }

void PrintVector(Coord_t v, char *String)
   {
   printf ("%s: %e %e %e\n", String, v[0],v[1],v[2]);
   }

void PrintMatrix (Matrix_t m, char *String)
   {
   PrintVector(m[0], String);
   PrintVector(m[1], String);
   PrintVector(m[2], String);
   }

void PrintCoord (Particle_t *a)
   {
   int i;
   LOOP (i, a->np)
      {
      printf
         (
         "%4i  %10.4f  %10.4f  %10.4f\n",
         a->type[i]+1,
         1e8*a->cur[i*NDIR+X],
         1e8*a->cur[i*NDIR+Y],
         1e8*a->cur[i*NDIR+Z]
         );
      }
   }
#endif