cdfill.c

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

/*
------------------------------------------------------------------------
Compile Switches
------------------------------------------------------------------------
*/



/*
------------------------------------------------------------------------
File Includes
------------------------------------------------------------------------
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "cdalloc.h"
#include "cdvect.h"
#include "iomngr.h"

#include "particle.h"
#include "cdfill.h"
#include  "cdsubs.h"
#include "cdhouse.h"
#include "strsub.h"
#include "parse.h"



/*
------------------------------------------------------------------------
Defines
------------------------------------------------------------------------
*/
#define BOUNDARY_TYPE int
#define FILL_NONE 0
#define FILL_BOX  1
#define FILL_PARALLELPIPED 2
#define FILL_CYLINDER 3
#define FILL_SPHERE  4


/*  Constants representing the three monoclinic directions  */
#define U 0
#define V 1
#define W 2
#define NVECT 3

#ifndef NSTR
#define NSTR 1024
#endif

#ifndef ANGS
#define ANGS 1e-8
#endif



/*
------------------------------------------------------------------------
Macros
------------------------------------------------------------------------
*/

#define CLEAN_AFTER_ERROR CleanAfterError();

/*
Replacement macros when testing program alone from XMD
*/
#ifdef TEST_DRIVER

#define INCREMENT_WARNING
#define GETS(STR,N,FILE) fgets(STR,N,stdin)


#undef WRITEMSG
#undef WRITEVAR
#define WRITEMSG \
   { \
   printf ("In file %s at line %i.\n", __FILE__, __LINE__); \
   }


#define WRITEVAR(VAR_NAME,VAR_TYPE) \
      { \
      printf ("FILE %s LINE %i :", __FILE__, __LINE__); \
      printf ("%s = ", #VAR_NAME); \
      printf (#VAR_TYPE, (VAR_NAME) ); \
      printf ("\n"); \
      }


#else

#define INCREMENT_WARNING IncrementNumberOfWarnings();
#define GETS(STR,N,FILE) m_gets_list(STR,N,FILE)

#endif


/*
------------------------------------------------------------------------
Defines
------------------------------------------------------------------------
*/
typedef double Matrix_t[NDIR][NDIR];


/*
------------------------------------------------------------------------
External Variables
------------------------------------------------------------------------
*/

#ifdef TEST_DRIVER
Particle_t *a;
LIST       *inlist;

#else
extern Particle_t *a;
extern LIST       *inlist;
#endif



/*
------------------------------------------------------------------------
Module wide variables
------------------------------------------------------------------------
*/
static int      NumPart_m   = 0;
static int     *FillType_m = NULL;
static Coord_t *FillPart_m = NULL;
static BOUNDARY_TYPE BoundaryType_m = FILL_NONE;

/*
Vector sets describing orientation,
repeating cell and re-oriented repeating cell
*/
static Coord_t Orient_m    [NDIR] = { {1,0,0}, {0,1,0}, {0,0,1} };
static Coord_t Cell_m      [NDIR] = { {1e-8,0,0}, {0,1e-8,0}, {0,0,1e-8} };



/*
Matrix for transforming between
original ("real") coordinates and re-oriented cell coordinates

Note:  This matrix does necessarily represent rotations.
This matrix has units of cm or 1/cm (since they convert between
coordinates (cm) and "cell vectors" (unitless).
*/
static Matrix_t RealToOrientCell_m;

/*  Alignment of cell origin within space  */
static Coord_t  Align_m        = {0.0, 0.0, 0.0};

/*
Margin to leave between boundary and particles.
This space created by scaling the coordinates after they
have been included within boundary
*/
static double   Margin_m = 0.0;

static BOOLEAN  CellRead_m = FALSE;


/*  BOX boundary  */
static Coord_t  Box_m[2];

/*  CYLINDER boundary  */
static double   CylLength_m;
static double   CylRadius_m;
static Coord_t  CylCenter_m;
static Coord_t  CylAxis_m;
static Matrix_t CylRealToOrient_m;
static Matrix_t CylOrientToReal_m;


/*  SPHERE boundary  */
static double   SphRadius_m;
static double   SphRadiusSqr_m;
static Coord_t  SphCenter_m;



/*
------------------------------------------------------------------------
Local Function Prototypes
------------------------------------------------------------------------
*/
static int  FillBoundary     (void);
static void ReadCellAngs     (Coord_t [NDIR]);
static void AddParticle      (Particle_t *, int, Coord_t);
static void GetLimits        (int [NVECT][2], BOUNDARY_TYPE);

/*  Box routines  */
static void    GetCornersBOX     (Coord_t **, int *);
static BOOLEAN IsWithinMarginBOX (Coord_t );
static void    GetMarginScaleBOX (Matrix_t, Coord_t );

/*  Cylinder routines  */
static void    GetCornersCYLINDER     (Coord_t **, int *);
static BOOLEAN IsWithinMarginCYLINDER (Coord_t );
static void    GetMarginScaleCYLINDER (Matrix_t, Coord_t );

/*  Sphere routines  */
static void    GetCornersSPHERE     (Coord_t **, int *);
static BOOLEAN IsWithinMarginSPHERE (Coord_t );
static void    GetMarginScaleSPHERE (Matrix_t, Coord_t );


#ifdef TEST_DRIVER
void PrintVector(Coord_t , char *);
void PrintMatrix (Matrix_t m, char *);
void PrintCoord (Particle_t *);
#endif


/*
------------------------------------------------------------------------
Exported Functions
------------------------------------------------------------------------
*/


                             /*   READ FILL COMMAND */

/*
Fill specified boundary with particles using defined orientation
*/
void read_fill (char *instr)
   {
   int  idir;
   int  ivector;
   int  jvector;
   int  ipart;
	int  NumNewPart;
   char *BufferPtr;
   char BufferStr[NSTR];
   char *tokstr = NULL;
   double swap;
   static Coord_t Zaxis = {0.0, 0.0, 1.0};


   tokstr = strhed (&instr);

   /*  Align option  */
   if (!strcmpi(tokstr,"ALIGN" ))
      {
      Align_m[X] = ANGS*dblstrf (&instr);
      Align_m[Y] = ANGS*dblstrf (&instr);
      Align_m[Z] = ANGS*dblstrf (&instr);
      }

   /*  Read boundary  */
   else if (!strcmpi(tokstr,"BOUNDARY" ))
      {

      tokstr = strhed (&instr);

      if (!strcmpi(tokstr, "BOX"))
         {
         BoundaryType_m = FILL_BOX;
         Box_m[0][X] = ANGS*dblstrf(&instr);
         Box_m[0][Y] = ANGS*dblstrf(&instr);
         Box_m[0][Z] = ANGS*dblstrf(&instr);
         Box_m[1][X] = ANGS*dblstrf(&instr);
         Box_m[1][Y] = ANGS*dblstrf(&instr);
         Box_m[1][Z] = ANGS*dblstrf(&instr);
         LOOP (idir, NDIR)
            {
            if (Box_m[0][idir] > Box_m[1][idir])
               {
               swap           = Box_m[0][idir];
               Box_m[0][idir] = Box_m[1][idir];
               Box_m[1][idir] = swap;
               }
            }
         }

      else if (!strcmpi(tokstr, "CYLINDER"))
         {
         BoundaryType_m = FILL_CYLINDER;
         CylRadius_m    = ANGS*dblstrf(&instr);
         CylLength_m    = ANGS*dblstrf(&instr);
         CylCenter_m[X] = ANGS*dblstrf(&instr);
         CylCenter_m[Y] = ANGS*dblstrf(&instr);
         CylCenter_m[Z] = ANGS*dblstrf(&instr);
         CylAxis_m  [X] = dblstrf(&instr);
         CylAxis_m  [Y] = dblstrf(&instr);
         CylAxis_m  [Z] = dblstrf(&instr);

         /*  Orientation matrix rotates cylinder axis into z  */
         MakeRotationMatrix(CylRealToOrient_m, Zaxis, CylAxis_m);
         InvertMatrix      (CylOrientToReal_m, CylRealToOrient_m);
         }

      else if (!strcmpi(tokstr, "SPHERE"))
         {
         BoundaryType_m = FILL_SPHERE;
         SphRadius_m    = ANGS*dblstrf(&instr);
         SphCenter_m[X] = ANGS*dblstrf(&instr);
         SphCenter_m[Y] = ANGS*dblstrf(&instr);
         SphCenter_m[Z] = ANGS*dblstrf(&instr);
         SphRadiusSqr_m = SphRadius_m * SphRadius_m;
         }

      else
         {
         printf ("WARNING:  Unrecognized BOUNDARY option.\n");
         INCREMENT_WARNING
         }
      }


   /*  Cell option  */
   else if      (!strcmpi(tokstr,"CELL"))
      {
      ReadCellAngs (Cell_m);
      CellRead_m = TRUE;
      }


   /*  go option  */
   else if (!strcmpi(tokstr,"GO" ))
      {

      /*  If no boundary specified, use repeating box  */
      if (FILL_NONE==BoundaryType_m)
         {
         BoundaryType_m = FILL_BOX;
         Box_m[0][X] = 0.0;
         Box_m[0][Y] = 0.0;
         Box_m[0][Z] = 0.0;
         Box_m[1][X] = a->bcur[X];
         Box_m[1][Y] = a->bcur[Y];
         Box_m[1][Z] = a->bcur[Z];
         }

      /*  Test for fill prepared  */
      if
      (
      0==NumPart_m ||
      FILL_NONE==BoundaryType_m
      )
         {
         printf ("ERROR:  Cannot perform fill without first calling");
         printf (" FILL PARTICLE.\n");
         CLEAN_AFTER_ERROR
         }

      /*  Do actual fill  */
      NumNewPart = 
			FillBoundary();
		printf ("***   Number new particles  %i\n", NumNewPart);

		/*  Make neighbor list as invalid  */
		a->InvalidNeighborList = TRUE;
      }

	/*  Margin option  */
   else if (!strcmpi(tokstr,"MARGIN" ))
		{
		Margin_m = ANGS*dblstrf(&instr);
		if (Margin_m < 0.0)
			{
			printf ("ERROR:  FILL MARGIN (%le) is less than zero.\n",
				Margin_m/ANGS);
			CleanAfterError();
			}
		}

   /*  Orient option  */
   else if (!strcmpi(tokstr,"ORIENT" ))
      {

      /*  Read vector  */
      LOOP (ivector, NVECT)
      LOOP (idir,  NDIR )
         Orient_m[ivector][idir] = dblstrf(&instr);

      /*  Normalized orientation vectors  */
      Normalize(Orient_m[U]);
      Normalize(Orient_m[V]);
      Normalize(Orient_m[W]);

      /*  Insure that orientation vectors are orthogonal  */
      LOOP (ivector, NDIR)
      LOOP (jvector, ivector)
         {
         if (fabs(GetDot(Orient_m[ivector], Orient_m[jvector])) > 1e-6)
            {
            printf
               (
               "ERROR:  Orientation vectors %i and %i are not orthogonal\n",
               jvector+1,
               ivector+1
               );
            CLEAN_AFTER_ERROR
            }
         }

      }


   /*  particle option  */
   else if (!strcmpi(tokstr,"PARTICLE" ))
      {

      /*  Free previous allocation  */
      FREE(FillType_m)
      FREE(FillPart_m)

      /*  Allocate particles  */
      NumPart_m = intstrf (&instr);
      ALLOCATE (FillType_m, int,     NumPart_m)
      ALLOCATE (FillPart_m, Coord_t, NumPart_m);

      /*  Read particles  */
      LOOP (ipart, NumPart_m)
         {
         GETS ( BufferStr, NSTR, inlist);
         BufferPtr = BufferStr;
         FillType_m[ipart]    = intstrf (&BufferPtr)-1;
         FillPart_m[ipart][X] = ANGS*dblstrf (&BufferPtr);
         FillPart_m[ipart][Y] = ANGS*dblstrf (&BufferPtr);
         FillPart_m[ipart][Z] = ANGS*dblstrf (&BufferPtr);
         }
      }

   else
      {
      printf ("Warning:  Unrecognized option to FILL command.\n");
      INCREMENT_WARNING
      }
   }


/*
Read and store cell vectors in format
  x1 y1 z1
  x2 y2 z2
  x3 y3 z3
*/
static void ReadCellAngs (Matrix_t Matrix)
   {
   char  BufferStr[NSTR];
   char *BufferPtr = NULL;
   int   ivector;

   LOOP (ivector, NDIR)
      {
      GETS (BufferStr, NSTR, inlist);
      BufferPtr = BufferStr;
      Matrix[ivector][X] = ANGS*dblstrf (&BufferPtr);
      Matrix[ivector][Y] = ANGS*dblstrf (&BufferPtr);
      Matrix[ivector][Z] = ANGS*dblstrf (&BufferPtr);
      }
   }


/*
Fill boundary with particles.
Set all new particles to SELECTED.
Unselect all previous particles
*/
static int FillBoundary (void)
   {
   int ipart;
   int idir;
   int ivect;
   int jvect;
   int kvect;
   int Limits[NVECT][2];
	int NumNewPart = 0;
   BOOLEAN UseThisPoint = FALSE;
   Coord_t Point;
   Coord_t  Decompose;
   Matrix_t ScaleMatrix;
	int      CellMultiple;
   Coord_t  CenterPoint;
   Coord_t  ScaledPoint;
	Matrix_t OrientToReal;
	Matrix_t RealToOrient;
	Coord_t  RealCell[NDIR];
   Coord_t  *RealFillPart = NULL;

	/*  Remove selection from all current particles  */
	LOOP (ipart, a->np)
		REMOVE_SELECT(ipart)

   /*
   Calculate matrix for converting new orientation to real coordinates
   NOTE: Since Orient_m is unitary matrix (rotation matrix), can
         invert it simply by transposing it.
   */
   CopyMatrix       (OrientToReal,  Orient_m);
   CopyMatrix       (RealToOrient,  OrientToReal);
   TransposeMatrix  (RealToOrient);


   /*  Calculate re-oriented cell directions  */
	MatrixVectorMultN(RealCell, OrientToReal, Cell_m, 3);

   /*  Calculate particles in real space */
   ALLOCATE (RealFillPart, Coord_t, NumPart_m)
	MatrixVectorMultN
		(RealFillPart, OrientToReal, FillPart_m, NumPart_m);

	/*
	Calculate matrix that converts real space coordinate to 
   coordinates expressed as combination of real space cell vectors
	*/
	InvertMatrix    (RealToOrientCell_m, RealCell);
	TransposeMatrix (RealToOrientCell_m);
	

   /*
   Linearly translate alignment point by cell vectors
   so that it is within a repeating cell centered at
   the origin.
      The resulting displacement from the origin
   will be added to all atoms which fill the boundary.
   */
   MatrixVectorMult (Decompose, RealToOrientCell_m, Align_m);
   LOOP (ivect, NVECT)
      {
      CellMultiple = fabs(Decompose[ivect])+0.5;
      if (Decompose[ivect] < 0)
         CellMultiple = -CellMultiple;
      LOOP (idir, NDIR)
         {
         Align_m[idir] -= CellMultiple * RealCell[ivect][idir];
         }
      }


   /*  Get lower and upper cell vector counts  */
   GetLimits (Limits, BoundaryType_m);

   /*
   Get values for scaling atoms withing margin distance of boundary surface
   */
   if (BoundaryType_m==FILL_BOX)
      {
      GetMarginScaleBOX (ScaleMatrix, CenterPoint);
      }

   else if (BoundaryType_m==FILL_CYLINDER)
      {
      GetMarginScaleCYLINDER (ScaleMatrix, CenterPoint);
      }

   else if (BoundaryType_m==FILL_SPHERE)
      {
      GetMarginScaleSPHERE (ScaleMatrix, CenterPoint);
      }


   /*  Sum over vector limits  */
	NumNewPart = 0;
   for (ivect=Limits[U][0]; ivect <=Limits[U][1]; ivect++)
   for (jvect=Limits[V][0]; jvect <=Limits[V][1]; jvect++)
   for (kvect=Limits[W][0]; kvect <=Limits[W][1]; kvect++)
      {

      /*  Sum over individual particles  */
      LOOP (ipart, NumPart_m)
         {

         /*  Loop over directions  */
         LOOP (idir, NDIR)
            {
            Point[idir] =
               RealFillPart[ipart][idir] +
               ivect * RealCell[U][idir] +
               jvect * RealCell[V][idir] +
               kvect * RealCell[W][idir] +
               Align_m[idir] -
               CenterPoint[idir];
            }