ptmc_data.c

/*! \file ptmc_data.c
 * \brief This program processes the data produced by ptmc.c
 * 
 * Input	:	The program requieres 1 or 2 files depending on wether 
 * 				the radial distribution function will be calculated or
 * 				not. These files should be passes at the time of 
 * 				execution. It is assumed that this program will be run 
 * 				in the same directory in which the output of files of
 * 				ptmc.c are located. \n\n\n 
 * 
 * 				The first file that should be passed is the same file
 * 				that was passed as first argument to the main program
 * 				ptmc.c . See the documentation of ptmc.c and ptmc.h
 * 
 * 				The second file should only be passed if the macro
 * 				CALCULATE_RDF was set to a nonzero value in the file
 * 				ptmc.h when the ptmc.c was compiled, which implies
 * 				that ptmc.c also calculated values of the radial
 * 				distribution function. In such case the file that should
 * 				be passed is the 3 file that was passed as an argument 
 * 				to ptmc.c, i.e., the configuration file for the
 * 				radial distribution functions.
 * 
 * Output	:	Assuming the ptmc_data is run in the same directory 
 * 				in which the output of the ptmc is the program will 
 * 				process the data files generated by ptmc.
 * 
 * 				The most important file it will generate is 
 * 				data.dat this file will contain the following 
 * 				information in columns: \n
 * 				1. Replica number, i. \n
 * 				2. Temperature of the replica, ti.\n
 * 				3. Internal energy of the replica ui. \n
 * 				4. Constant volume heat capacity, C_vi. \n
 * 				5. Number of times the replica i tried to move out of 
 * 				of the simulation box. \n
 * 				6. % of rejected moves. \n
 * 				7. Number exchanges attempted between replica i and i+1.
 * 				\n
 * 				8. Number of exchanges between replica i and i+1.\n
 * 				9. Ratio of exchanges and attempted exchanges between 
 * 				replica i and replica i+1.\n\n
 * 				
 * 				If the macro NUMBER_OF_DOPANTS was set to a nonzero value 
 * 				then the following extra column will appear:\n
 * 				10. Ratio of accepted exchanges between the two 
 * 				different types of atoms.\n\n
 * 
 * 				Then the program will also generate the following files
 * 				for each replica i:\n\n
 * 				
 * 				energy[i].dat: the first column of this file will contain
 * 				the energy and the second the probability having that 
 * 				energy for the ith replica.\n\n
 * 				
 * 				If CALCULATE_RDF was set to a nonzero value then also
 * 				the following files will be generated:\n\n
 * 
 * 				radialA[i].dat: the first column of this file will contain
 * 				the distance from the center of mass and the second the
 * 				probability of the particles of the matrix being at that
 * 				distance from the center of mass.\n\n
 * 
 * 				If also NUMBER_OF_DOPANTS is set to a non zero then
 * 				also the following files will be generated:
 * 				radialB[i].dat: the first column of this file will contain
 * 				the distance from the center of mass and the second the
 * 				probability of the particles of the dopant being at that
 * 				distance from the center of mass.\n\n
 * 
 * 				If the macro CALCULATE_RDF is set to a nonzero value
 * 				then this program will also generate a file called
 * 				rdf_data.dat that will contain the following:\n
 * 				1. Replica number, i. \n
 * 				2. Temperature of the replica, ti.\n
 * 				3. Average position of the matrix atoms measured from the 
 * 				geometric center of the distribution.\n
 * 				4. Standard deviation of the main position of the matrix	
 * 				atoms from the geometric center of the distribution.\n
 * 
 * 				If also NUMBER_OF_DOPANTS is non zero the the two extra 
 * 				colums will be produced:\n
 * 
 * 				5. Average position of the dopant atoms measured from the 
 * 				geometric center of the distribution.\n
 * 				6. Standard deviation of the main position of the dopant 	
 * 				atoms from the geometric center of the distribution.\n  
 */ 









#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#include"ptmc.h"


int
main (int argc, char *argv[])
{

  char charU[MAX_LENGTH_CHAR], average_save[MAX_LENGTH_CHAR],
    charE[MAX_LENGTH_CHAR];

  int gg, i, j, nU;
  int nprocs;
  int n, swap_freq, mc_steps, save_freq, eq_steps;
  double t0, tf, U0, Uf, dU, box_radius;
  double tried, accepted;
  double counter, countrej, outmoves, normU, U, U2;
  double t;
  FILE *np;
  FILE *exc;
  FILE *in;
  FILE *out;
  FILE *dat;
#if CALCULATE_RDF
  FILE *rdf_dat;
#endif


#if CALCULATE_RDF
  char rpf1char[MAX_LENGTH_CHAR], ordchar[MAX_LENGTH_CHAR],
    rdf1char[MAX_LENGTH_CHAR], pos[MAX_LENGTH_CHAR];
  int rdf_freq, nradial;
  double radialmin, radialmax, dradial, r;
  double rpf1a, rpf1a2, normrpf1;
#if NUMBER_OF_DOPANTS
  char rpf2char[MAX_LENGTH_CHAR], rdf2char[MAX_LENGTH_CHAR];
  double rpf2a, rpf2a2, normrpf2;
#endif
#endif

#if NUMBER_OF_DOPANTS
  double swap, swap_rejected;
#endif

  dat = fopen ("data.dat", "w");
#if CALCULATE_RDF  
  rdf_dat=fopen ("rdf_data.dat", "w");
#endif
  np = fopen ("nprocs.dat", "r");
  exc = fopen ("exchange.dat", "r");
  gg = fscanf (np, "%d", &nprocs);
  fprintf (stdout, "%d\n", nprocs);

  in = fopen (argv[1], "r");
  if (in == NULL)
    {
      fprintf (stdout, "The input file %s was not found\n", argv[1]);
      exit (0);
    }
  else
    {
      gg = fscanf (in, "%d %d %d %d %d %lf %lf %lf %lf %lf %lf",
		   &n, &swap_freq, &mc_steps, &save_freq, &eq_steps, &t0, &tf,
		   &U0, &Uf, &dU, &box_radius);
      fclose (in);
    }

#if CALCULATE_RDF
  in = fopen (argv[2], "r");
  if (in == NULL)
    {
      fprintf (stdout, "The input file %s was not found\n", argv[2]);
      exit (0);
    }
  else
    {
      gg =
	fscanf (in, "%d %lf %lf %lf", &rdf_freq, &radialmin, &radialmax,
		&dradial);
      fclose (in);
    }
  nradial = (int) ((radialmax - radialmin) / dradial);
#endif

  nU = (int) ((Uf - U0) / dU);

  for (i = 1; i < nprocs; i++)
    {
    /****************************************************************
     * Write the names of the input files in their respective strings
     ***************************************************************/

      sprintf (charU, "e%d.dat", i);
      sprintf (charE, "energy%d.dat", i);
      sprintf (average_save, "avs%d.dat", i);
#if CALCULATE_RDF
      sprintf (rpf1char, "rpfA%d.dat", i);
      sprintf (rdf1char, "radialA%d.dat", i);

#if NUMBER_OF_DOPANTS
      sprintf (rpf2char, "rpfB%d.dat", i);
      sprintf (rdf2char, "radialB%d.dat", i);
#endif
      sprintf (ordchar, "ord%d.dat", i);
      sprintf (pos, "pos%d.dat", i);
#endif

    /*******************************************************************
     * Generates energy Historams with their proper origin
     ******************************************************************/

      in = fopen (charU, "r");
      out = fopen (charE, "w");
      for (j = 0; j < nU; j++)
	{
	  gg = fscanf (in, "%lf", &U);
	  fprintf (out, "%lf %lf\n", U0 + j * dU + dU / 2, U / (nU));
	}
      fclose (in);
      fclose (out);

    /*******************************************************************
     * Calculates the main observables such as internal energy and
     * C_v for each replica running at a different temperature
     ******************************************************************/

      gg = fscanf (np, "%lf", &t);
      gg = fscanf (exc, "%lf %lf", &tried, &accepted);
      in = fopen (average_save, "r");
      if (in == NULL)
	{
	  fprintf (stdout, "The file %s was not found\n", average_save);
	}
      else
	{
	  while (!feof (in))
	    {
#if NUMBER_OF_DOPANTS
	      gg = fscanf (in, "%lf %lf %lf %lf %lf %lf %lf %lf",
			   &counter, &countrej, &outmoves, &normU, &U, &U2,
			   &swap, &swap_rejected);
#else
	      gg = fscanf (in, "%lf %lf %lf %lf %lf %lf",
			   &counter, &countrej, &outmoves, &normU, &U, &U2);
#endif
	    }

	  fclose (in);
	  fprintf (dat, "%d %lf %lf %lf %lf %lf %lf %lf %lf ", i, t, U / normU,
		   ((U2 / normU) - (U / normU) * (U / normU)) / (t * t),
		   outmoves, countrej / counter, tried, accepted,
		   accepted / tried);
#if NUMBER_OF_DOPANTS
	  fprintf (dat, "%lf ", swap_rejected / swap);
#endif
	  fprintf (dat, "\n");
	}

	/*******************************************************************
     * Calculates the observables related to the RDFS
     * for each replica running at a different temperature
     ******************************************************************/


#if CALCULATE_RDF
      in = fopen (rpf1char, "r");
      out = fopen (rdf1char, "w");
      for (j = 0; j < nradial; j++)
	{
	  gg = fscanf (in, "%lf", &r);
	  fprintf (out, "%lf %lf\n", radialmin + j * dradial + dradial / 2,
		   r / (nradial));
	}
      fclose (in);
      fclose (out);
#if NUMBER_OF_DOPANTS
      in = fopen (rpf2char, "r");
      out = fopen (rdf2char, "w");
      for (j = 0; j < nradial; j++)
	{
	  gg = fscanf (in, "%lf", &r);
	  fprintf (out, "%lf %lf\n", radialmin + j * dradial + dradial / 2,
		   r / (nradial));
	}
      fclose (in);
      fclose (out);

#endif
#endif


#if CALCULATE_RDF
      in = fopen (ordchar, "r");
      if (in == NULL)
	{
	  fprintf (stdout, "The file %s was not found\n", ordchar);
	}
      else
	{
	  while (!feof (in))
	    {
#if NUMBER_OF_DOPANTS
	      gg =
		fscanf (in, "%lf %lf %lf %lf %lf %lf", &normrpf1, &rpf1a,
			&rpf1a2, &normrpf2, &rpf2a, &rpf2a2);
#else
	      gg = fscanf (in, "%lf %lf %lf", &normrpf1, &rpf1a, &rpf1a2);
#endif
	    }
#endif
	  fclose (in);


	}
#if CALCULATE_RDF
      fprintf (rdf_dat, "%d %lf %lf %lf ", i, t, rpf1a / normrpf1,
	       sqrt (rpf1a2 / normrpf1 -
		     (rpf1a / normrpf1) * (rpf1a / normrpf1)));
#if NUMBER_OF_DOPANTS
      fprintf (rdf_dat, "%lf %lf", rpf2a / normrpf2,
	       sqrt (rpf2a2 / normrpf2 -
		     (rpf2a / normrpf2) * (rpf2a / normrpf2)));
#endif
      fprintf (rdf_dat, "\n");
#endif


    }
  fclose (np);
  fclose (dat);
#if CALCULATE_RDF  
  fclose (rdf_dat);
#endif
  return 0;
}