plotc.c

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

int ver
)
   {
   float xpos, ypos, srad, textsize;
   char tstr[64];
   float Disp [NDIR];
   float Image[NDIR];

   if (pl->symbol[i] || pl->UseDisplacement)
      {
      xpos = xstart+scale*a->cur[hor + 3*i];
      ypos = ystart+scale*a->cur[ver + 3*i];
      srad =        scale*pl->radius[i];
      grdev_moveto (xpos, ypos);
      grdev_setfill (pl->fill[i]);
		grdev_setcolor (pl->ColorModel[i], pl->ColorTable[i]);
      switch (pl->symbol[i])
         {
			case CIR : DrawCircle (i, srad, pl); break;
#if 0
         case CIR : grdev_circle     (srad); break;
#endif
         case CRO : grdev_cross      (srad); break;
         case TRI : grdev_triangle   (srad); break;
         case ITR : grdev_itriangle  (srad); break;
         case DIA : grdev_diamond    (srad); break;
         case AST : grdev_asterisk   (srad); break;
         case SQU : grdev_square     (srad); break;
         case NUM : textsize = grdev_gettextsize();
                    grdev_settextsize (srad);
                    sprintf (tstr, "%i", i+1);
                    grdev_writetext  (tstr);
                    grdev_settextsize (textsize);
                    break;
         }
      if (pl->UseDisplacement)
         {
         /*  Find displacment (include wrapping)  */
         CalcDisp (Disp, a, i);

         /*  Scale displacement  */
         Disp[hor] *= pl->dispscale;
         Disp[ver] *= pl->dispscale;

         /*  Calculate image position  */
         Image[hor] = a->cur[hor + 3*i] + Disp[hor];
         Image[ver] = a->cur[ver + 3*i] + Disp[ver];

         /*  Clip image position if repeating boundary  */
         ClipAtom (&a->cur[3*i], Image, a->bcur, hor, ver);

         /*  Draw tail to image position  */
         xpos = xstart + scale*Image[hor];
         ypos = ystart + scale*Image[ver];
         grdev_drawto (xpos, ypos);
         }
      }
   }


void plot_coord (char *outname, Particle_t *a, PlotInfo_t *pl, int sflag)
   {
   time_t tstruct;

   FILE  *fout;
   int   *layerbot = NULL;
	int   *layertop = NULL;
	unsigned *ix = NULL;
	int      *vect = NULL;
   float xlimit, ylimit, ximage, yimage, xframe, yframe;
   float xgap, ygap , xorg, yorg;
   float xstart, ystart, scale, ratio;
   float tempbox, t;
   int   ipart, jpart, ilayerbox, ilayerplot;
   int   nx, ny, nbox;
   int   curpagelayer, npagelayer, curpage;
   int   iline, iind, jind;

   int    hor  = pl->orient[0];
   int    ver  = pl->orient[1];
   int    dep  = pl->orient[2];

   double (*c)[3] = (double (*)[3]) a->cur;

   /*
   TEST DATA INTEGRITY
   */

   /*  TEST BOX  */
   if (a->bcur[0]<=0 || a->bcur[1]<=0 || a->bcur[2]<=0)
      {
		ERROR_PREFIX
      printf ("Box dimensions less or equal to zero.\n");
		CleanAfterError();
      }

   /*  TEST FOR INITIALIZE OF PlotInfo_t  */
   if (pl->npl<a->np)
      {
		ERROR_PREFIX
      printf ("Cannot plot, not symbols specified.\n");
		CleanAfterError();
      }

   /*  WHEN PLOTTING DISPLACEMENTS, TEST FOR REFERENCE PARTICLES  */
   if (pl->UseDisplacement && a->ref==NULL)
         {
			ERROR_PREFIX
      	printf ("Cannot plot displacements\n");
      	printf ("   because no reference particle have been specified.\n");
			CleanAfterError();
         }

	/*  Wrap particles  */
	WrapParticles (a);

   /* set directions on plotted page */

   /*  room on page for drawing - in inches  */
   xlimit = pl->PageSize[X] - MLEFT - MRIGHT;  
   ylimit = pl->PageSize[Y] - MBOT  - MTOP  ;

	/* size of image in simulation units   */
   ximage = a->bcur[hor];       
   yimage = a->bcur[ver];

   /* Calculate maximum number of layers drawn on page */
   npagelayer = (pl->nlayerplot-1) / pl->page + 1;

   /*
     If all layers can't fit in one column then rescale xFrame and yFrame
     so that they fit into the least number of multiple columns
   */
   if (npagelayer==1)
      nx = ny = 1;
   else
      {
      ratio = (ylimit/yimage) / (xlimit/ximage);
      if (ratio>=1.0)
         {
		    ny = sqrt (npagelayer * ratio) + 1;
   		 nx = (npagelayer-1) / ny + 1;
         }
      else
         {
         nx = sqrt (npagelayer / ratio) + 1;
         ny = (npagelayer-1) / nx + 1;
         }
      }


   /* conversion factor s.u. -> in.*/
   t      = (xlimit-2*nx*MARGIN)/ximage/nx;
   scale  = (ylimit-2*ny*MARGIN)/yimage/ny;
   if (t<scale)
      scale=t;
   xframe = scale*ximage + 2*MARGIN;  /* size of frame around float layer  */
   yframe = scale*yimage + 2*MARGIN;
   /* gap between frame and allotted space (inches) */
   xgap   = (xlimit/nx - xframe);
   ygap   = (ylimit/ny - yframe);
   /*x, y  coordiantes of first layer drawing origin (inches) */
   xorg   = MLEFT + xgap/2.0;
   yorg   = MBOT  + ygap/2.0;

         /* order atoms by layer */

	ALLOCATE(ix,   unsigned int, a->np)
	ALLOCATE(vect, int, a->np)
   tempbox = a->bcur[dep];
   if (tempbox==0)
      {
      printf ("INTERNAL ERROR (plot_coord): Box in depth direction is zero.\n");
		CleanAfterError();
      }

	/*  Wrap particles  */
	LOOP (ipart, a->np)
      {
      t = fabs(c[ipart][dep])/tempbox;
      if (t>=1)
         t -= (long) t;
      if (c[ipart][dep]<0)   /* RECOVER SIGN  */
         t  = -t;
      c   [ipart][dep] = tempbox * t;
      vect[ipart]       = MAXUNS * t;
      ix  [ipart]       = ipart;
      }


   sortix (vect, ix, a->np);
	FREE(vect)


   /*  ASSIGN SPACE FOR LAYER TOP AND BOX INDICES  */
	ALLOCATE(layertop, int, pl->nlayerbox)
	ALLOCATE(layerbot, int, pl->nlayerbox)

	LOOP (ilayerbox,  pl->nlayerbox)
      {
      layertop[ilayerbox] = 1;
      layerbot[ilayerbox] = 0;
      }
   tempbox = a->bcur[dep];
   jind    = a->np-1;
	LOOP (iind, a->np)
      {
      ipart  = ix[iind];
      jpart  = ix[jind];
      ilayerbox =    pl->nlayerbox*c[ipart][dep]/tempbox;
      layerbot[ilayerbox] = iind;
      ilayerbox =    pl->nlayerbox*c[jpart][dep]/tempbox;
      layertop[ilayerbox] = jind;
      jind--;
      }


   /*  OPEN OUTPUT FILE  */
   if (*outname=='+')
      fout = fopen (outname+1, "at");
   else
      fout = fopen (outname  , "wt");

   /*  TEST FOR FILE OK  */
   if (fout==NULL)
      {
		ERROR_PREFIX
      printf ("Cannot open output file\n");
		CleanAfterError();
      }

   /*  SET DEVICE  */
   if (pl->DEVflag==1)
      grdev_canon();
   else if (pl->DEVflag==0)
      grdev_ps();

   /*  ASSIGN FILE TO PLOTTER  */
   grdev_assignfile (fout);

   /*  INITIALIZE PAGE  */
   grdev_init ();
   grdev_initfile();
   grdev_settextsize(0.11111111);  /*  1/9  */

          /*  loop over layers  */

   nbox = nx * ny;
	LOOP (ilayerplot, pl->nlayerplot)
      {
      ilayerbox = pl->layerplot[ilayerplot]-1;
      /*   determine current page */
      curpage      = ilayerplot / npagelayer;
      /*   determine position of this layer on current page */
      curpagelayer = ilayerplot % npagelayer;

      /*   eject current page if starting new page  */
      if ( (curpagelayer==0) && (curpage>0)  )
         grdev_writepage();

      /*   plot  identification   */
      if (curpagelayer==0)
         {
         grdev_moveto (MLEFT, MBOT);
			LOOP (iline, 8)
         if (a->title[iline] && *(a->title[iline]) )
            {
            grdev_moveto (MLEFT, MBOT-(iline+1)*HLETT);
            grdev_writetext (a->title[iline]);
            }

         grdev_moveto (MLEFT, MBOT-9*HLETT);
         if (a->run>=0)
            grdev_printf    ("RUN  %6li", a->run);
         grdev_moverel (1.25, 0);
         if (a->step>=0)
            grdev_printf    ("STEP %6li", a->step);
         grdev_moverel (1.25, 0);
         grdev_printf
            ("ORIENT %c/%c", OrientSymbol_m[ver], OrientSymbol_m[hor]);
         grdev_moverel (1.25, 0);
         time(&tstruct);
         grdev_writetext (ctime(&tstruct));
         grdev_moverel (1.75, 0);
         grdev_printf    ("Page %i of %i", curpage+1, pl->page);
         }

      xstart = xorg +  (     curpagelayer  ) / ny * (xframe + xgap);
      ystart = yorg +  (nbox-curpagelayer-1) % ny * (yframe + ygap);
      /* draw frame around current layer */
      grdev_moveto ( xstart,  ystart);
      grdev_drawrel( xframe,  0     );
      grdev_drawrel( 0,       yframe);
      grdev_drawrel(-xframe,  0     );
      grdev_drawrel( 0,      -yframe);
      /* write layer number */
      grdev_moveto (xstart+0.25*MARGIN, ystart+0.25*MARGIN);
      grdev_printf ("%i of %i", ilayerbox+1, pl->nlayerbox);

      xstart = xstart + MARGIN;
      ystart = ystart + MARGIN;

          /*  loop over particles in current layer  */

      for (iind=layertop[ilayerbox];iind<=layerbot[ilayerbox];iind++)
         {
         ipart= ix[iind];
			if ( !sflag || IS_SELECT(ipart))
            draw_part (a, ipart, pl, xstart, ystart, scale, hor, ver);
         }

		/*  Reset color to black after drawing particles  */
		grdev_setcolor (COLOR_RGB, BlackTable_m);

      /*  DRAW BONDS IN CURRENT LAYER */
      if (pl->UseBond)
         {
         float rsq,dt,dx,dy,dz, x1,y1,z1;
         float bond1,bond2,bond3,blosq,bhisq, xpos,ypos;

         bond1 = pl->bondhi;
         bond2 = 1.414215*bond1;  /*  sqrt(2)  */
         bond3 = 1.732052*bond1;  /*  sqrt(3)  */
         blosq = pl->bondlo * pl->bondlo;
         bhisq = pl->bondhi * pl->bondhi;
         for (iind=layertop[ilayerbox];iind<=layerbot[ilayerbox];iind++)
            {
            ipart = ix[iind];
				if ( !sflag || IS_SELECT(ipart))
               {
               x1 = c[ipart][hor];
               y1 = c[ipart][ver];
               z1 = c[ipart][dep];
               for (jind=iind+1;jind<=layerbot[ilayerbox];jind++)
                  {
                  jpart = ix[jind];
						if ( !sflag || IS_SELECT(jpart))
                     {
                     dx = fabs(x1-c[jpart][hor]);
                     if (dx<bond1)
                        {
                        dy = fabs(y1-c[jpart][ver]);
                        dt = dx+dy;
                        if (dt<bond2)
                           {
                           dz = fabs(z1-c[jpart][dep]);
                           if (dt+dz<bond3)
                              {
                              rsq = dx*dx + dy*dy + dz*dz;
                              /*  DRAW CURRENT BOND  */
                              if (rsq>=blosq && rsq<=bhisq)
                                 {
                                 xpos = xstart+scale*c[ipart][hor];
                                 ypos = ystart+scale*c[ipart][ver];
                                 grdev_moveto (xpos, ypos);
                                 xpos = xstart+scale*c[jpart][hor];
                                 ypos = ystart+scale*c[jpart][ver];
                                 grdev_drawto (xpos, ypos);
                                 }
                              }
                           }
                        }
                     }
                  }
               }
            }
         }
      }

   /*  PRINT FINAL PAGE ??  */
   grdev_writepage();


   /*  CLOSE FILE  */
   fclose (fout);

   /*  FREE STORAGE  */
	FREE(ix)
	FREE(layerbot)
	FREE(layertop)
   }


void CalcDisp (float Disp[NDIR], Particle_t *a, int i)
   {
   int idir;

   LOOP (idir, NDIR)
      {
      Disp[idir] = a->ref[idir + 3*i] - a->cur[idir + 3*i];
      if (!a->surf[idir])
         {
         if      (2*Disp[idir] < -a->bcur[idir])
            Disp[idir] += a->bcur[idir];
         else if (2*Disp[idir] >  a->bcur[idir])
            Disp[idir] -= a->bcur[idir];
         }
      }
   }


/*
Assume position 1 is inside box which runs from 0 to Box[]
Move   position 2 inside box along line joining 1 and 2
*/
void ClipAtom
(
double  Pos1 [NDIR],
float   Pos2 [NDIR],
double  Box [NDIR],
int     hor,
int     ver
)
   {
   double hDel;
   double vDel;
   double Slope;

	/*  
	If both atoms are outside the box, don't clip  
	This is a quick-dirty hack to correct for the
	fact that plot_coord() does not handle free
   surface simulations well when atoms are outside
   the nominal box.
	*/
	if 
	(
		(
		Pos1[hor] < 0 || Pos1[hor]>=Box[hor] ||
		Pos1[ver] < 0 || Pos1[ver]>=Box[ver]
		) 
		&&
		(
		Pos2[hor] < 0 || Pos2[hor]>=Box[hor] ||
		Pos2[ver] < 0 || Pos2[ver]>=Box[ver]
		)
	)
		return;
	

   /*  Displacement between two positions  */
   hDel = Pos2[hor] - Pos1[hor];
   vDel = Pos2[ver] - Pos1[ver];

   /*  Two points are coincident - no need to clip  */
   if (hDel==0 && vDel==0)
      return;

   /*  Horizontal Line  */
   if (vDel==0.0)
      {
      if (Pos2[hor] < 0.0)
         Pos2[hor] = 0.0;
      else if (Pos2[hor] > Box[hor])
         Pos2[hor] = Box[hor];
      return ;
      }

   /*  Vertical Line  */
   if (hDel==0)
      {
      if (Pos2[ver] < 0.0)
         Pos2[ver] = 0.0;
      else if (Pos2[ver] > Box[ver])
         Pos2[ver] = Box[ver];
      return;
      }


   /*  General case  */
   Slope = vDel / hDel;

   if (Pos2[hor] < 0.0)
      {
      Pos2[ver] -= Slope * Pos2[hor];
      Pos2[hor]  = 0.0;
      }
   else if (Pos2[hor] > Box[hor])
      {
      Pos2[ver] += Slope * (Box[hor] - Pos2[hor]);
      Pos2[hor]  = Box[hor];
      }

   if (Pos2[ver] < 0.0)
      {
      Pos2[hor] -= Pos2[ver] / Slope;
      Pos2[ver]  = 0.0;
      }
   else if (Pos2[ver] > Box[ver])
      {
      Pos2[hor] += (Box[ver] - Pos2[ver]) / Slope;
      Pos2[ver]  = Box[ver];
      }

   }


void DrawCircle (int ipart, float radius, PlotInfo_t *p)
	{
	/*  Filled circle - draw black outline  */
	if (p->fill[ipart])
		{
		/*  Draw fill  */
		grdev_setfill  (TRUE);
		grdev_setcolor (p->ColorModel[ipart], p->ColorTable[ipart]);
		grdev_circle   (radius);

		/* Draw outline  */
      grdev_setfill  (FALSE);
		grdev_setcolor (COLOR_RGB, BlackTable_m);
		grdev_circle   (radius);
		}

	/*  Unfilled circle - draw circle in color  */
	else
		{
		grdev_setfill  (FALSE);
		grdev_setcolor (COLOR_RGB, BlackTable_m);
		grdev_circle   (radius);
		}
	
	}