/*************************************************************************
 *                                                                       *
 *               Script for calculating density of states                *
 *                        from CPMD 3.4.1 output                         *
 *                        Veronika Brazdova 2002                         *
 *                        vb@chemie.hu-berlin.de                         *
 *                                                                       *
 *               Based on dosscript and dos98.f written by               *
 *                  Schwarz and Wittenberg (FHI Berlin)                  *
 *                                                                       *
 *                   compilation: cc dos.c -lm -o dos                    *
 *                                                                       *
 *************************************************************************/

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

#define MAX_NUM_EIG 3000 
#define NUM_M 1000 /* max. number of k points */
#define PI 3.141592654
#define BUFFER 150

double broad(double, double, double, int, double);

/// Return codes for function main:
//  0 successful run
//  1 no command line arguments or file not found
//  2 incorrect command line arguments

int main(int argc, char *argv[])
{
  int    i, j, k, type, found;
  int    no_kpt, no_states, no_eig;
  double width, mu, min_eig, max_eig, delta_energy;
  double e_sum, t_dos, energy;
  double weight[NUM_M], dos[2][MAX_NUM_EIG];
  double eig[MAX_NUM_EIG], eig_weight[MAX_NUM_EIG];
  char   line[BUFFER];
  FILE   *fp;
  double difference;

  found = 0;

  // argv[0] name-of-the-executable
  // argv[1] input file
  // argv[2] broadening function type: l=Lorentz, g=Gauss
  // argv[3] width of the broadening function

  /// Command line arguments
  // 
  if (argc < 3)
  {
    fprintf(stderr,
      "dos:  constructs density of states from CPMD 3.4.1 output\n");
    fprintf(stderr,"usage:    dos cpmd_output_file function width\n\n");
    fprintf(stderr,"function: the broadening function; possible values\n");
    fprintf(stderr,"          l for Lorentz function\n");
    fprintf(stderr,"          g for Gauss function\n");
    fprintf(stderr,"width:    width of the broadening function (in eV)\n"); 
    fprintf(stderr,"default width: 0.05 eV for l, 0.1 eV for g\n");
    return 1;
  }
  
  if ((fp = fopen(argv[1],"r")) == NULL )
  {
    fprintf(stderr,"%s : cannot open\n",argv[1]);
    return 1;
  } 
  

  /// Check ig argv[2] is of an allowed value, because function 
  /// 'broad' doesn't do it.
  if (!strcmp(argv[2],"l"))
    type = 1; 
  else
  {
    if (!strcmp(argv[2],"g"))
      type = 0;
    else 
    {
      fprintf(stderr,
        "%s is not an allowed value for broadening function type\n",
        argv[2]);
      return 2;
    }
  }

  /// Default width: 0.1 for Gaussian broadening, 0.05 for Lorentz 
  if (argc < 4)
    if(type == 0)
      width = 0.1; 
    else
      width = 0.05;
  else
    width = atof(argv[3]);

  /// Read in number of states

  while (fgets(line,BUFFER,fp))
    if (strstr(line,"NUMBER OF STATES") != NULL)
    {
      sscanf(line," NUMBER OF STATES: %i",&no_states);
      break;
    }

  /// Read in number of k-points and their weight
  while (fgets(line,BUFFER,fp))
    if (strstr(line,"KPOINTS") != NULL && 
       (strstr(line,"CARTESIAN")) != NULL)
    {
      found = 1;
      sscanf(line," KPOINTS (IN CARTESIAN COORDINATES AS INPUT): %i",
        &no_kpt);
      fgets(line,BUFFER,fp);
      for (i = 0; i < no_kpt; i++)
      {
        fgets(line,BUFFER,fp);
        sscanf(line,"    %*i %*lf %*lf %*lf %lf %*i",&weight[i]);
      }
      break;
    }

  if (!found)
  {
    fprintf(stderr,"%s\n","Number of k points not specified!");
    return 2;
  }

  no_eig = no_kpt * no_states;

  /// Read in the final eigenvalues for each k-point 
  /// and then the chemical potential

  while (fgets(line,BUFFER,fp))
    if (strstr(line,"FINAL RESULTS") != NULL)
      break;

  i = 0; k = 0;

  while (fgets(line,BUFFER,fp))
  {
    if (strstr(line,"K POINT") != NULL)
    {
      for(j = i; j < (i+no_states); j += 2)
      {
        fgets(line,BUFFER,fp);
        sscanf(line," %*i %lf %*lf %*i %lf %*lf",&eig[j],&eig[j+1]);
        eig_weight[j] = weight[k];
        eig_weight[j+1] = weight[k];
      }
      i = i + no_states;
      k = k + 1;
    }
    if (strstr(line,"CHEMICAL POTENTIAL") != NULL)
    {
      sscanf(line," CHEMICAL POTENTIAL = %lf %*s",&mu);
      break;
    }
  }

  /// Find the min and max eigenvalue
  min_eig = eig[0];
  max_eig = eig[0];
  
  for (i = 1; i < no_eig; i++)
  {
    if (eig[i] < min_eig) min_eig = eig[i];
    if (eig[i] > max_eig) max_eig = eig[i];
  }

  delta_energy = max_eig - min_eig;
  
  /// Construct DOS

  e_sum = 0;
  for(i = 0; i < NUM_M; i++)
  {
    t_dos = 0;
    energy = min_eig + (1.2 * i/NUM_M - 0.1) * delta_energy;      
    for (j = 0; j < no_eig; j++)
      t_dos = t_dos + broad(energy,eig[j],width,type,eig_weight[j]);
    dos[0][i] = energy - mu;
    dos[1][i] = t_dos;
    e_sum = e_sum + t_dos;
    printf("%lf %14.10lf\n",dos[0][i],dos[1][i]);
  } 
  
  /// Check the integral of DOS

  e_sum = e_sum * 1.2 * delta_energy/NUM_M;
  difference = (e_sum-2 * no_states)/(2 * no_states);
  printf("# TOTAL DOS = %lf\n",e_sum);   
  printf("# difference : %lf\n",difference);
  if (difference < 0.0)
    difference = -difference;
  if (difference > 0.001)
    printf("# The integral of DOS does not have a correct value.\n");
      
}

/// Function broad broadens _energy_ by a Lorentz or Gauss function
//  of width _width_, centered at _eigenvalue_, normed by _e_weight_
//  of the _eigenvalue_ (this is in fact the weight of a given k-point
//  type 0: Gauss function
//  type 1: Lorentz function

double broad(double energy, double eigenvalue, double width, int type,\
             double e_weight)
{
  double value, tmp;
         
  if (type) 
  { 
    tmp = pow((energy-eigenvalue),2);
    value = (2 * e_weight/PI) * width/(tmp + pow(width,2));
    return value;   
  }
  else
  {
    tmp = pow(((energy - eigenvalue)/width),2);
    value=(2 * e_weight/sqrt(PI)/width) * exp(-tmp); 
    return value;
  }
}