cdwrite.c

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

						(OutputFile,
						"%10.4f %10.4f %10.4f\n",
						1e8*a->disp[X+3*i],
						1e8*a->disp[Y+3*i],
						1e8*a->disp[Z+3*i]
					);
				}
			}
      }


   else if (!strcmpi("EPOT", tokstr) || !strcmpi("ENERGY", tokstr) )
      {
		double Epot = 0.0;
		double NumPart = 0;

		/*  calculate latest energy  */
      energy_all (a, s, UseSelect);

		/*  test for mass and velocity present  */
		if (a->eatom!=NULL && a->np!=0)
			{
			LOOP (i, a->np)
				{
				if (!UseSelect || IS_SELECT(i))
					{
					Epot += a->eatom[i];
					NumPart++;
					}
				}
			}

		/*  prepare energy  */
		if (NumPart==0)
			Epot = 0.0;
		else
   		Epot = s->eunit*Epot/NumPart;

  		fprintf (OutputFile, "EPOT %17.9e\n", Epot);
      }

   else if (!strcmpi("EATOM", tokstr))
      {
		double ekin;
		double (*Vel)[NDIR] = (double (*)[NDIR]) a->v;

      /*  Check for s->dtime zero  */
      if (s->dtime==0.0)
         {
         printf
			("WARNING:  Time step is zero, cannot calculate kinetic energy.\n");
         IncrementNumberOfWarnings();
			}

		/*  Get latest energy  */
      energy_all (a, s, UseSelect);

		/*  Print results  */
      fprintf (OutputFile, "EATOM %i\n", NumPartUsed);

		/*  Print verbose comments  */
		if (s->UseVerboseOutputFile)
			{
			fprintf (OutputFile, 
   	    "# type          x          y          z epot (%6s) ekin (%6s)\n",
				s->eulabel, s->eulabel);
			fprintf (OutputFile, "#----- ---------- ----------");
			fprintf (OutputFile, " ---------- ------------- -------------\n");
			}

		/*  Print energies  */
		LOOP(i, a->np)
		if (!UseSelect || IS_SELECT(i))
         {
         fprintf
				(
				OutputFile, "%6i %10.4lf %10.4lf %10.4lf %+13.6e",
            a->type[i]+1,
            1e8*a->cur[X+3*i],
            1e8*a->cur[Y+3*i],
            1e8*a->cur[Z+3*i],
				a->eatom[i]*s->eunit
				);

			/*  Calculate kinetic energy if appropriate  */
			if (s->dtime==0 || Vel==NULL || IS_FIX(i))
				fprintf (OutputFile, "\n");
			else
				{
				ekin = 0.5 * a->mass[i] *
					(SQR(Vel[i][X]) + SQR(Vel[i][Y]) + SQR(Vel[i][Z]))
					/ SQR(s->dtime);
				fprintf (OutputFile, " %+13.6e\n", ekin*s->eunit);
				}
         }
      }

   else if (!strcmpi("ETOT", tokstr))
      {
      fprintf (OutputFile, "ETOT %14.6e\n", (a->epot+a->ekin)*s->eunit/a->np);
      }

	/*  Temperature of particles which are both SELECTED and FREE  */
   else if (!strcmpi("TEMP", tokstr) )
		{
 	   fprintf 
			( 
			OutputFile, 
			"TEMP %14.6e\n", 
			GetTemperature(a,s,UseSelect)
			);
		}

	else if (!strcmpi("EKIN", tokstr) )
      {
		
      /*  KE of selected particle per atom  */
		if (UseSelect) {
			Value = s->eunit*GetKE(a, s, UseSelect);
		/*  Kinetic Energy of PARTICLES and BOX divided by ALL Part */
		} else {
			CalcKineticEnergy(a,s);
			Value = s->eunit*a->ekin/a->np;
		}

  	   fprintf (OutputFile, "EKIN %14.6e\n", Value);
      }

   else if (!strcmpi("DTIME", tokstr))
      {
      fprintf (OutputFile, "DTIME %14.6e\n", s->dtime);
      }

   else if (!strcmpi("CSTEP", tokstr))
      {
      fprintf (OutputFile, "CSTEP %14.6e\n", s->cstep);
      }

   else if (!strcmpi("COR", tokstr))
      {
      FILE *fout;

      tokstr = strhed (&instr);

		/*  Check for TAG option */
		if (!strcmpi("TAG",tokstr)) {
			tokstr = strhed (&instr);
			UseTag = TRUE;
			}

      if (*tokstr=='+')
         fout = fopen (tokstr+1, "ab");
      else
         fout = fopen (tokstr, "wb");
      if (fout==NULL)
         {
         printf ("ERROR:  Cannot open file %s\n", tokstr);
         return;
         }
      WriteCorFile (fout, a, UseSelect, UseTag);
      fclose (fout);
      }


   else if (!strcmpi("STRESS", tokstr))
      {
      BOOLEAN  DeallocForce;
		BOOLEAN  OldStressState;

      /*  Determine if force should be deallocated  */
      DeallocForce = (a->f == NULL);

		/*  Save previous stress output state  */
		OldStressState = s->StressStepOutput.Save;
		s->StressStepOutput.Save = TRUE;

		/*  Set stress select flag  */
		s->UseSelectStress = UseSelect;


      /*  Allocate space if needed  */
		if (a->f==NULL)
			{
			ALLOCATE (a->f, double, NDIR*a->NumPartAlloc)
			}

      /*  Calculate force  */
      ( (PotForce_f *) s->forcepnt ) (a, s);

      /*  Print Stress */
		fprintf
			(
			OutputFile,
			"STRESS %10.3le %10.3le %10.3le %10.3le %10.3le %10.3le\n",
			a->Stress[XX]/a->np,
			a->Stress[YY]/a->np,
			a->Stress[ZZ]/a->np,
			a->Stress[YZ]/a->np,
			a->Stress[ZX]/a->np,
			a->Stress[XY]/a->np
			);

      /*  Deallocate force if needed  */
      if (DeallocForce)
         FREE(a->f)

		/*  Reset old stress state  */
		s->StressStepOutput.Save = OldStressState;

		/*  Reset stress select flag  */
		s->UseSelectStress = FALSE;
      }


   else if (!strcmpi("FORCE", tokstr))
      {
      BOOLEAN  DeallocForce;

      /*  Determine if force should be deallocated  */
      DeallocForce = (a->f == NULL);

      /*  Allocate space if needed  */
		if (a->f==NULL)
			{
			ALLOCATE (a->f, double, NDIR*a->NumPartAlloc)
			}

      /*  Calculate force  */
      ( (PotForce_f *) s->forcepnt ) (a, s);

		/*  Print position statistics  */
		if (UseAvg || UseMin || UseMax)
			{
			CalcAvgMinMax (a->f, a->np, NDIR, Avg, Min, Max, UseSelect);

			if (UseAvg)
				fprintf
					(
					OutputFile,
					"AVG_FORCE %12.4le %12.4le %12.4le\n",
					Avg[X],
					Avg[Y],
					Avg[Z]
					);

			if (UseMin)
				fprintf
					(
					OutputFile,
					"MIN_FORCE %12.4le %12.4le %12.4le\n",
					Min[X],
					Min[Y],
					Min[Z]
					);

			if (UseMax)
				fprintf
					(
					OutputFile,
					"MAX_FORCE %12.4le %12.4le %12.4le\n",
					Max[X],
					Max[Y],
					Max[Z]
					);
			}


		/*  Print individual velocities  */
		else
			{
	      /*  Print force  */
	      fprintf (OutputFile, "FORCE %i\n", NumPartUsed);
			LOOP (i, a->np)
				if (!UseSelect || IS_SELECT(i))
	            fprintf (OutputFile, "%4i %12.4le %12.4le %12.4le\n",
	                    a->type[i]+1,
	                    a->f[X + 3*i],
	                    a->f[Y + 3*i],
	                    a->f[Z + 3*i]);
			}

      /*  Deallocate force if needed  */
      if (DeallocForce)
         FREE(a->f)
      }

   else if (!strcmpi("EXTFORCE", tokstr))
      {
      /*  Test for existance of external force  */
      if (a->ForcePtrList== NULL)
         {
         printf ("WARNING: There is no applied external force.\n");
         IncrementNumberOfWarnings();
         }
      else
         {
			fprintf (OutputFile, "#  Force units are dynes\n");
         fprintf (OutputFile, "EXTFORCE %i\n", NumPartUsed);
         PrintExternalVectorList
            (
            OutputFile,
            a->ForcePtrList,
            a->type,
            a->Selection,
            a->np,
            UseSelect
            );
         }
      }


   else if (!strcmpi("EXTSPRING", tokstr))
      {
      /*  Test for existance of external spring  */
      if (a->SpringPtrList== NULL)
         {
         printf ("WARNING: There is no applied external spring.\n");
         IncrementNumberOfWarnings();
         }
      else
         {
			fprintf (OutputFile, "#  Spring constant units are erg/cm^2\n");
         fprintf (OutputFile, "EXTSPRING %i\n", NumPartUsed);
         PrintExternalVectorList
            (
            OutputFile,
            a->SpringPtrList,
            a->type,
            a->Selection,
            a->np,
            UseSelect
            );
         }
      }

   /*  Write data in PDB (Protein Data Bank) format  */
   else if (!strcmpi("PDB", 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  */
      WritePDB (fout, a, UseSelect);
      fclose (fout);
      }


   else if (!strcmpi("STATE", tokstr))
		{
		/*  Get first token for state file name  */
		tokstr = strhed (&instr);
		
		/*  Write state file  */
      writestate (a, s, tokstr);

		/*  
		Update Neighbor List so that subsequent trajectory
		will agree with trajectory that follows from this 
		state file
		*/
		UpdateNeighborList(a);
		}

   else if (!strcmpi("STEP", tokstr))
      fprintf (OutputFile, "STEP %li\n", a->step);

   else if (!strcmpi("RCV", 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;
         }
      writercv (fout, a, UseSelect);
      fclose (fout);
      }

   else if (!strcmpi("RDF", tokstr))
      {
		double   fac;
		double   r;
      BOOLEAN rep[NDIR];
      long   CumulativeCount;

      NumTable_m  = intstrf (&instr);
      MinRadius_m = 1e-8 * dblstrf (&instr);
      MaxRadius_m = 1e-8 * dblstrf (&instr);

      if (IsBlank(instr))
         UseTypes_m = FALSE;
      else
         {
         UseTypes_m = TRUE;
         Type1_m = intstrf (&instr)-1;
         Type2_m = intstrf (&instr)-1;
         }

      /*  Test for type in range  */
      if (UseTypes_m && (Type1_m < 0  ||  Type2_m < 0))
         {
         printf ("ERROR:  One or more types are out of range (<0).\n");
         CleanAfterError();
         }

      /*  ALLOCATE TABLE  */
      if (NumTable_m<=0)
         return;
		ALLOCATE (RdfTable_m, long, NumTable_m)
      rep[X] = !a->surf[X];
      rep[Y] = !a->surf[Y];
      rep[Z] = !a->surf[Z];

#if 0
/*
Old code for use with module bondsub.c
*/
      /*  Call RDF without types if type==0  */
      if (UseTypes_m)
         CalcTypeRDF
            (
            a->cur, a->type, a->np, rmin, rmax,
            type1, type2, a->bcur, rep, table, ntable
            );
      else
         CalcRDF
            (
            a->cur, a->np, rmin, rmax, a->bcur, rep, table, ntable
            );
#endif

/*
New code (dec 12, 1997) for use with CalcNeighborPerAtom()
*/

		/*
		Store values for recall via CallBack Function
		*/
		Type_m      = a->type;
		RdfTableFactor_m = NumTable_m / (MaxRadius_m-MinRadius_m);

		/*  Call neighbor routine (which calls callback function)  */
		CalcNeighborPerAtom
			(
			(double (*)[NDIR]) a->cur,
			a->Selection,
			a->np,
			UseSelect,
			MinRadius_m,
			MaxRadius_m,
			a->bcur,
			rep,
			TRUE,
			CalcRdf
			);
			

      fprintf
			(
			OutputFile,
			"RDF  %i  %lf  %lf\n",
			NumTable_m,
			1e8*MinRadius_m,
			1e8*MaxRadius_m
			);
		fac = (MaxRadius_m - MinRadius_m) / NumTable_m;
      CumulativeCount = 0;
		LOOP (i, NumTable_m)
         {
         r = 1e8 * (MinRadius_m + (i+0.5)*fac);
         CumulativeCount += RdfTable_m[i];
         fprintf
				(
				OutputFile,
				"%8.2lf %8li %8li\n",
				r,
				RdfTable_m[i],
				CumulativeCount
				);
         }

      /*  FREE STORAGE  */
		FREE (RdfTable_m)
      }

   else if (!strcmpi("RUN", tokstr))
      fprintf (OutputFile, "RUN %li\n", a->run);

   else if (!strcmpi("ILIST", 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;
            }

         }