cdwrite.c

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

      else
         fout = stdout;


      fprintf (fout,"ILIST %i\n", NumPartUsed);
      cnt = 0;
		LOOP (i, a->np)
         {
         /*  If select flag, print only selected particles  */
			if (!UseSelect || IS_SELECT(i))
            {
            fprintf (fout,"%4i", i+1);
            cnt++;

            /*  If end of line, print new line  */
            if (cnt==16)
               {
               cnt=0;
               fprintf (fout,"\n");
               }

            /*  .. else print blank for separator  */
            else
               fprintf (fout, " ");
            }
         }
      if (cnt)
         fprintf (fout,"\n");

      /*  Close file  */
      if (IsFile)
         fclose(fout);
      }

   else if (!strcmpi("TYPELIST", tokstr))
      {
      FILE *fout;
      BOOLEAN IsFile;

      tokstr = strhed (&instr);

      /*  Is output to a file?  */
      IsFile = (tokstr[0]!='\0');

      /*  Open file  */
      if (IsFile)
         {

         /*  Open file  */
         if (*tokstr=='+')
            fout = fopen (tokstr+1, "at");
         else
            fout = fopen (tokstr, "wt");

         /*  Error opening file  */
         if (fout==NULL)
            {
            printf ("ERROR:  Cannot open file %s\n", tokstr);
            return;
            }
         }

      else
         fout = stdout;


      fprintf (fout, "TYPELIST %i\n", NumPartUsed);
      cnt = 0;
		LOOP (i, a->np)
         {
			if (!UseSelect || IS_SELECT(i))
            {
            fprintf (fout," %4i", a->type[i]+1);
            cnt++;
            }
         if (cnt==16)
            {
            cnt=0;
            fprintf (fout,"\n");
            }
         }
      if (cnt)
         fprintf (fout,"\n");

      /*  Close file  */
      if (IsFile)
        fclose(fout);
      }


   /*  Write data for XMOL program  */
   else if (!strcmpi("XMOL", tokstr))
      {
      FILE *fout;

      tokstr = strhed (&instr);
      if (*tokstr=='+')
         fout = fopen (tokstr+1, "at");
      else
         fout = fopen (tokstr, "wt");
      if (fout==NULL)
         {
         printf ("ERROR:  Cannot open file %s\n", tokstr);
         return;
         }

      /*  Write XMOL format  */
      WriteXMOL (fout, a, UseSelect);
      fclose (fout);
      }

   /*  Write VAR format data  */
   else if (!strcmpi("VAR", tokstr)) {
	   fprintf (OutputFile, "VAR %s %8i\n", instr, NumPartUsed);
		LOOP (i, a->np)
			if (!UseSelect || IS_SELECT(i)) {
				/*  Run through list of variables  */
				strncpy (buffer, instr, NSTR);
				bptr = &buffer[0];
				vptr = strhed(&bptr);
				while (*vptr!=0) {
					fprintf (OutputFile, " %s", get_pptr(a,vptr,i,get_pfmt(vptr)) );
				vptr = strhed (&bptr);
				}
			fprintf (OutputFile, "\n");
			}
		}

	/*  Write random number seed  */
	else if (!strcmpi("SEED",tokstr)) {
		fprintf (OutputFile, "SEED %u\n", s->seed);
	}


   /*  TREAT AS NUMBER  */
   else
      {
      char *hedstr = tokstr;
      double x = dblstrf (&tokstr);
      fprintf (OutputFile, "%s %g\n", hedstr, x);
      }

   /*  If file opened then close it  */
   if (UseOutputFile)
      fclose (OutputFile);
	}



/*
************************************************************************
Local functions
************************************************************************
*/
void PrintExternalVectorList
(
FILE             *OutputFile,
ExternalVector_t **VectorList,
BYTE             *Type,
SEL_T            *Select,
int               NumList,
BOOLEAN           UseSelect
)
   {
   int ilist;

   /*  Empty list - print nothing  */
   if (VectorList==NULL)
      return;

   /*  Print vector list  */
   LOOP (ilist, NumList)
      if (!UseSelect || IS_SELECT2(Select,ilist))
         {
         fprintf (OutputFile, "%4i", Type[ilist]+1);
         fprintf (OutputFile, " %10.4e %10.4e %10.4e", 
               CM *VectorList[ilist]->Origin[X],
               CM *VectorList[ilist]->Origin[Y],
               CM *VectorList[ilist]->Origin[Z]
					);
         if (VectorList[ilist]==NULL)
            fprintf (OutputFile, " %10.4e %10.4e %10.4e\n",
					0.0, 0.0, 0.0);
         else
            fprintf (OutputFile, " %10.4e %10.4e %10.4e\n",
               VectorList[ilist]->Vector[X],
               VectorList[ilist]->Vector[Y],
               VectorList[ilist]->Vector[Z]
               );
         }

   }


void WriteXMOL (FILE *OutputFile, Particle_t *a, BOOLEAN UseSelect)
   {
   int ipart;
   int NumWrite;
   char *NameStr;

   if (UseSelect)
      NumWrite = a->nsel;
   else
      NumWrite = a->np;

   /*  Print number of particles written  */
   fprintf (OutputFile, "%i\n", NumWrite);

   /*  Print run and step number (this will be step "label")  */
   fprintf (OutputFile,
      "Run %li  Step %li  %s\n", a->run, a->step, GetCurrentTimeStr() );

   /*  Write particle coordinates  */
   LOOP (ipart, a->np)
	if (!UseSelect || IS_SELECT(ipart))
      {

      /*  Print type name  */
      NameStr = a->TypeNames[a->type[ipart]];
      if (NameStr == NULL)
         fprintf (OutputFile, " %i ", a->type[ipart]+1);
      else
         fprintf (OutputFile, " %s ", NameStr);

      /*  Print coordinates in angstroms  */
      fprintf
         (
         OutputFile,
         " %e %e %e\n",
         CM * a->cur[NDIR*ipart+X],
         CM * a->cur[NDIR*ipart+Y],
         CM * a->cur[NDIR*ipart+Z]
         );

      }

   }

void WritePDB (FILE *OutputFile, Particle_t *a, BOOLEAN UseSelect)
   {
   int ipart;
   int NumWrite;
   char *NameStr;

   /*  Find number of atoms  */
   if (UseSelect)
      NumWrite = a->nsel;
   else
      NumWrite = a->np;

   /*  Print run and step number (this will be step "label")  */
   fprintf (OutputFile, "REMARK Output from XMD\n");
   fprintf (OutputFile, "REMARK (c) 1996, Center for Materials Simulation,");
   fprintf (OutputFile, " University of Connecticut\n");
   fprintf (OutputFile,
      "REMARK Run %li  Step %li  %s\n", a->run, a->step, GetCurrentTimeStr() );
   fprintf (OutputFile, "REMARK  Number of Atoms = %i\n", NumWrite);


   /*  Write particle coordinates  */
   NumWrite = 0;
   LOOP (ipart, a->np)
	if (!UseSelect || IS_SELECT(ipart))
      {

      /*  Count number of particles written  */
      NumWrite++;

      /*  Print tag type and atom number  */
      fprintf (OutputFile, "HETATM%6i", NumWrite);

      /*  Print type name  */
      NameStr = a->TypeNames[a->type[ipart]];
      if (NameStr == NULL)
         fprintf (OutputFile, "%-3i", a->type[ipart]+1);
      else
         fprintf (OutputFile, "%-3s", NameStr);

      /*  Print stuff  */
#if 0
      fprintf (OutputFile, "  AUT     1    ");
#endif
      fprintf (OutputFile, "          1    ");

      /*  Print coordinates in angstroms  */
      fprintf
         (
         OutputFile,
#if 0
         "%8.3lf%8.3lf%8.3lf  0.00  0.00\n",
#endif
         "%8.3lf%8.3lf%8.3lf\n",
         CM * a->cur[NDIR*ipart+X],
         CM * a->cur[NDIR*ipart+Y],
         CM * a->cur[NDIR*ipart+Z]
         );

      }

   /*  Print end tag  */
   fprintf (OutputFile, "END\n");

   }


/*  Calculate average, minimum and maximum values of an array  */
void CalcAvgMinMax
(
double *Array,
int    NumValue,
int    NumDim,
double *Avg,
double *Min,
double *Max,
BOOLEAN UseSelect
)
	{
	int idim;
	int ival;
	int NumAvg;
	double  Value;
	BOOLEAN IsMinMaxInitialized = FALSE;


	/*  Empty array, return  */
	if (NumValue==0 || NumDim==0 || Array==NULL)
		return;
	
	ASSERT(Avg!=NULL)
	ASSERT(Min!=NULL)
	ASSERT(Max!=NULL)

	/*  Initialize output to zero  */
	LOOP (idim, NumDim)
		{
		Avg[idim] = 0.0;
		}


	/*  Loop through array  */
	NumAvg=0;
	LOOP (ival, NumValue)
	if   (!UseSelect || IS_SELECT(ival))
		{

		/*  Treat each dimension of current array group  */
		LOOP (idim, NumDim  )
			{
	
			/*  Get value of this array element */
			Value = Array[ival*NumDim + idim];

			/*  Add to average sum  */
			Avg[idim] += Value;

			/*  If Min,Max uninitialized, then intialize them to current value */
			if (!IsMinMaxInitialized)
				{
				Min[idim] = Value;
				Max[idim] = Value;
				}

			/*  Compare current value to Min,Max  */
			else
				{
				if (Value < Min[idim])
					{
					Min[idim] = Value;
					}
				if (Value > Max[idim])
					{
					Max[idim] = Value;
					}
				}
			}

		/*  Increment number of average sums  */
		NumAvg++;

		/*  Min and Max have been initialized  */
		IsMinMaxInitialized = TRUE;
		}

	/*  Done if no particles selected  */
	if (NumAvg==0)
		return;

	/*  Take average  */
	LOOP (idim, NumDim)
		{
		Avg[idim] /= NumAvg;
		}
	}


void CalcRdf
(
WORD ipart,
WORD *NeighList,
double *NeighRadSqr,
WORD NumNeigh
)
	{
	int ineigh;
	int jpart;
	int Index;
	BOOLEAN TypesOK;
	double  Radius;

	ASSERT (NeighRadSqr!=NULL)
	ASSERT (NeighList!=NULL)

	TypesOK = TRUE;

	LOOP (ineigh, NumNeigh)
		{

		/*
		Using only pairs of specific type
		*/
		if (UseTypes_m)
			{
			jpart = NeighList[ineigh];
			TypesOK =
				(Type_m[ipart]==Type1_m && Type_m[jpart]==Type2_m) ||
				(Type_m[jpart]==Type1_m && Type_m[ipart]==Type2_m);
		
			}

		if (TypesOK)
			{
			Radius = sqrt(NeighRadSqr[ineigh]);
			Index  = RdfTableFactor_m*(Radius-MinRadius_m);
			if (Index>=0 && Index<NumTable_m)
				RdfTable_m[Index]++;
			}
		}
	}



/*  
Convert symbol (t,m,x,y,z,vx,vy,vz,ep)  
to corresponding default print format
*/
char *get_pfmt (char *symbol) {

	if (!strcmp("i",symbol)) return "%d";
	if (!strcmp("g",symbol)) return "%d";
	if (!strcmp("t",symbol)) return "%d";
	if (!strcmp("m",symbol)) return "%10.4lf";
	if (!strcmp("x",symbol)) return "%10.4lf";
	if (!strcmp("y",symbol)) return "%10.4lf";
	if (!strcmp("z",symbol)) return "%10.4lf";
	if (!strcmp("vx",symbol)) return "%12.4e";
	if (!strcmp("vy",symbol)) return "%12.4e";
	if (!strcmp("vz",symbol)) return "%12.4e";
	if (!strcmp("ep",symbol)) return "%13.6e";
	/*  Default format  */
	return "%13.6e";
}


/*  
Convert symbol (t,m,x,y,z,vx,vy,vz,ep), index and format into
a pointer static formatted string.
*/
char *get_pptr (Particle_t *a, char *symbol, int index, char *fmt) {
	static char buffer[NSTR];

	/*  Index  */
	if (!strcmp("i",symbol)) return  sprintf(buffer,fmt,index+1),buffer;

	/*  Type  */
	if (!strcmp("t",symbol)) return  PFMT(a->type, a->type[index]+1);

	/*  Tag   */
	if (!strcmp("g",symbol)) return  PFMT(a->tag, a->tag[index]+1);

	/*  Mass  */
	if (!strcmp("m",symbol)) return  PFMT(a->mass, AVAG*a->mass[index]);

	/*  Position  */
	if (!strcmp("x",symbol)) return  PFMT(a->cur, 1e8*a->cur[NDIR*index+X]);
	if (!strcmp("y",symbol)) return  PFMT(a->cur, 1e8*a->cur[NDIR*index+Y]);
	if (!strcmp("z",symbol)) return  PFMT(a->cur, 1e8*a->cur[NDIR*index+Z]);

	/*  Velocity */
	if (!strcmp("vx",symbol)) return PFMT(a->v, a->v[NDIR*index+X]/s->dtime);
	if (!strcmp("vy",symbol)) return PFMT(a->v, a->v[NDIR*index+Y]/s->dtime);
	if (!strcmp("vz",symbol)) return PFMT(a->v, a->v[NDIR*index+Z]/s->dtime);

	/*  Eatom  */
	if (!strcmp("ep",symbol)) return PFMT(a->eatom, s->eunit*a->eatom[index]);

	/*  Ekin  */
	if (!strcmp("ek",symbol))  {
		double vx,vy,vz;
		if (a->v==NULL) return "0";
		if (a->mass==NULL) return "0";
		vx=a->v[NDIR*index+X];
		vy=a->v[NDIR*index+Y];
		vz=a->v[NDIR*index+Z];
		sprintf (buffer, "%e", 
			s->eunit*0.5*a->mass[index]*(vx*vx+vy*vy+vz*vz)/(s->dtime*s->dtime));
		return buffer;
	}

	/*  Temperature  */
	if (!strcmp("tm",symbol)) {
		double vx,vy,vz,df;
		if (a->v==NULL || IS_FIX(index) || a->mass==NULL)  return "0";
		vx=a->v[NDIR*index+X];
		vy=a->v[NDIR*index+Y];
		vz=a->v[NDIR*index+Z];
		df = (a->cons && a->cons[index]) ? 
				((a->cons[index]->type==CONSTR_LINE) ? 1 : 2) : 3;
		sprintf (buffer, "%e", 
			a->mass[index]*(vx*vx+vy*vy+vz*vz)/(s->dtime*s->dtime)/(df*BK));
		return buffer;
	}

	/*  Degree of freedom  */
	if (!strcmp("df",symbol))  {
		sprintf (buffer, "%d",
			(a->cons && a->cons[index]) ?  
				((a->cons[index]->type==0) ? 1 : 2) : 3);
		return buffer;
	}

	return "0";
}