/*

  Usage:   thisprogram   Cgra.Dos.val Cgra.Dos.vec 

  output: Cgra.DOS.Tetrahedron
          Cgra.PDOS.Tetrahedron

*/

 
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "Inputtools.h"

#define BINARYWRITE 

typedef struct {
  double r,i;
}  dcomplex; 

#define eV2Hartree 27.211386
#define PI 3.1415926535897932384626

#define YOUSO10  100

void *malloc_multidimarray(char *type, int N, int *size);
void OrderE0(double *e, int n);
void OrderE(double *e,double *a, int n);
void ATM_Dos(double *et, double *e, double *dos);
void ATM_Spectrum(double *et,double *at, double *e, double *spectrum);

double ***EigenE;
double ****PDOS;
int **Msize1;

/*
   input s0
   output r0
   if input='../test.Dos.val' output='test'
*/

char *Delete_path_extension(char *s0, char *r0)
{
  char *c;
  /* char s[YOUSO10]; */

  /* find '/' */
  c=rindex(s0,'/');
  if (c)  {
    strcpy(r0,c+1);
  }
  else {
    strcpy(r0,s0);
  }
  printf("<%s>\n",r0);

  if (strlen(r0)==0 ) { return  NULL; }
  
  c =index(r0,'.');
  if (c) {
    *c='\0';
  }
  printf("<%s>\n",r0);
  return r0;

}


void Dos_Histgram( char *basename, int Dos_N,
		   int knum_i,int knum_j, int knum_k,
		   int SpinP_switch, double *****EIGEN, int neg,
		   double Dos_emin, double Dos_emax, 
		   int iemin, int iemax )
{
  int ie,spin,i,j,k,ieg,iecenter;
  double *DosE, **Dos;
  double eg,x,factor;
  int N_Dos , i_Dos[10];

  char file_Dos[YOUSO10];
  FILE *fp_Dos;

  printf("<Dos_Histgram> start\n");

  if (strlen(basename)==0) {
    strcpy(file_Dos,"DOS.Histgram");
  }
  else {
    sprintf(file_Dos,"%s.DOS.Histgram",basename);
  }
  if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
    printf("can not open a file %s\n",file_Dos);
    return;
  }

  printf("<Dos_Histgram> make %s\n",file_Dos);

  /* allocation */
  DosE=(double*)malloc(sizeof(double)*Dos_N);
  N_Dos=2; i_Dos[0]=SpinP_switch+1; i_Dos[1]=Dos_N;
  Dos=(double**)malloc_multidimarray("double",N_Dos,i_Dos);


  /* initialize */
  for (ie=0;ie<Dos_N;ie++) {
    DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
  }
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ie=0;ie<Dos_N;ie++) {
      Dos[spin][ie]=0.0;
    }
  }

  for (spin=0;spin<=SpinP_switch;spin++) {
    for (i=0; i<=(knum_i-1); i++){
      for (j=0; j<=(knum_j-1); j++){
	for (k=0; k<=(knum_k-1); k++){
	  for (ieg=0;ieg<neg; ieg++) {

	    eg = EIGEN[spin][i][j][k][ieg] ;
	    x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
	    iecenter = (int)x ;

	    if (iecenter>=0 && iecenter <Dos_N )
	      Dos[spin][iecenter] = Dos[spin][iecenter] + 1.0;
	  } /* ieg */
	} /* k */
      } /* j */
    } /* i */
  } /* spin */

  /* normalize */
  factor =  eV2Hartree * knum_i * knum_j * knum_k * (Dos_emax-Dos_emin)/(Dos_N-1); 
  factor = 1.0/factor;
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ie=0;ie<Dos_N;ie++) {
      Dos[spin][ie] = Dos[spin][ie] * factor;
    }
  }

  if (SpinP_switch==1) {
    for (ie=0;ie<Dos_N;ie++) {
      fprintf(fp_Dos,"%lf %lf %lf\n",(DosE[ie])*eV2Hartree, Dos[0][ie], -Dos[1][ie]);
    }
  }
  else {
    for (ie=0;ie<Dos_N;ie++) {
      fprintf(fp_Dos,"%lf %lf\n",(DosE[ie])*eV2Hartree, Dos[0][ie]*2.0);
    }
  }

  if (fp_Dos) fclose(fp_Dos);


  free(DosE);
  free_multidimarray(Dos,N_Dos,i_Dos);


}


void Dos_Gaussian( char *basename, int Dos_N, double Dos_Gaussian, 
		   int knum_i,int knum_j, int knum_k,
		   int SpinP_switch, double *****EIGEN, double *******ev,
                   int neg,
		   double Dos_emin, double Dos_emax, 
		   int iemin, int iemax,
                   int *WhatSpecies, int *Spe_Total_CNO, int atomnum )
{

  /*****************************************************************
                            Method = Gaussian 
  *****************************************************************/
  /*     exp( -(x/a)^2 ) 
	 a= Dos_Gaussian 
	 x=-3a : 3a is counted */
  double pi2;
  int iewidth,ie;

  int spin,i,j,k,ieg,iecenter;
  int wanA,GA_AN,tnoA,i1;
  double eg,x,xa,factor,wval;

  double *DosE, **Dos;
  int N_Dos,i_Dos[10];

  char file_Dos[YOUSO10];
  FILE *fp_Dos;

  printf("<Dos_Gaussian> start\n");

  if (strlen(basename)==0) {
    strcpy(file_Dos,"DOS.Gaussian");
  }
  else {
    sprintf(file_Dos,"%s.DOS.Gaussian",basename);
  }
  if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
    printf("can not open a file %s\n",file_Dos);
    return;
  }
  
  printf("<Dos_Gaussian> make %s\n",file_Dos);

  /* allocation */
  DosE=(double*)malloc(sizeof(double)*Dos_N);
  N_Dos=2; i_Dos[0]=SpinP_switch+1; i_Dos[1]=Dos_N;
  Dos=(double**)malloc_multidimarray("double",N_Dos,i_Dos);

  /* initialize */
  for (ie=0;ie<Dos_N;ie++) {
    DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
  }
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ie=0;ie<Dos_N;ie++) {
      Dos[spin][ie]=0.0;
    }
  }

  pi2 = sqrt(PI);
  iewidth = Dos_Gaussian*3.0/(Dos_emax-Dos_emin)*(Dos_N-1)+3;
  printf("Gaussian width=%d\n",iewidth); 

  for (spin=0;spin<=SpinP_switch;spin++) {
    for (i=0; i<=(knum_i-1); i++){
      for (j=0; j<=(knum_j-1); j++){
	for (k=0; k<=(knum_k-1); k++){
	  for (ieg=0;ieg<neg; ieg++) {

            /* generalization for spin non-collinear */
            wval = 0.0;
   	    for (GA_AN=0; GA_AN<atomnum; GA_AN++){
  	      wanA = WhatSpecies[GA_AN];
	      tnoA = Spe_Total_CNO[wanA];
  	      for (i1=0; i1<tnoA; i1++){
                wval += ev[spin][i][j][k][GA_AN][i1][ieg];
	      }
	    }

	    eg = EIGEN[spin][i][j][k][ieg] ;
	    x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
	    iecenter = (int)x ; 

	    iemin = iecenter - iewidth;
	    iemax = iecenter + iewidth ;

	    if (iemin<0) iemin=0;
	    if (iemax>=Dos_N) iemax=Dos_N-1;
	    if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax <Dos_N ) {

	      for (ie=iemin;ie<=iemax;ie++) {
		xa = (eg-DosE[ie])/Dos_Gaussian; 

		/*
            printf("Dos_N=%3d iemin=%3d iemax=%3d Dos_Gaussian=%15.12f %15.12f\n",
                    Dos_N,iemin,iemax,Dos_Gaussian,wval*exp( -xa*xa)/(Dos_Gaussian*pi2));
		*/


		Dos[spin][ie] += wval*exp(- xa*xa)/(Dos_Gaussian*pi2);
	      }
	    } 
	  } /* ieg */
	} /* k */
      } /* j */
    } /* i */
  } /* spin */
 
  /* normalize */
  factor = 1.0/(double)( eV2Hartree * knum_i * knum_j * knum_k );
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ie=0;ie<Dos_N;ie++) {
      Dos[spin][ie] = Dos[spin][ie] * factor;
    }
  }

  if (SpinP_switch==1) {
    for (ie=0;ie<Dos_N;ie++) {
      fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie], -Dos[1][ie]);
    }
  }
  else {
    for (ie=0;ie<Dos_N;ie++) {
      fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie]*2.0);
    }
  }

  if (fp_Dos) fclose(fp_Dos);

  free(DosE);
  free_multidimarray(Dos,N_Dos,i_Dos);
}



void DosDC_Gaussian( char *basename, int atomnum, 
                     int Dos_N, double Dos_Gaussian, 
		     int SpinP_switch,
  		     double Dos_emin, double Dos_emax, 
		     int iemin, int iemax,
                     int *WhatSpecies, int *Spe_Total_CNO )
{

  /*****************************************************************
                          Method = Gaussian 
  *****************************************************************/
  /********************************
          exp( -(x/a)^2 ) 
	  a= Dos_Gaussian 
 	  x=-3a : 3a is counted 
  ********************************/

  double pi2;
  int iewidth,ie;
  int GA_AN,wanA,tnoA;

  int spin,i,j,k,ieg,iecenter;
  double eg,x,xa,factor,sum;

  double *DosE, **Dos;
  int N_Dos,i_Dos[10];

  char file_Dos[YOUSO10];
  FILE *fp_Dos;

  printf("<Dos_Gaussian> start\n");

  if (strlen(basename)==0) {
    strcpy(file_Dos,"DOS.Gaussian");
  }
  else {
    sprintf(file_Dos,"%s.DOS.Gaussian",basename);
  }
  if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
    printf("can not open a file %s\n",file_Dos);
    return;
  }
  
  printf("<Dos_Gaussian> make %s\n",file_Dos);

  /* allocation */
  DosE=(double*)malloc(sizeof(double)*Dos_N);
  N_Dos=2; i_Dos[0]=SpinP_switch+1; i_Dos[1]=Dos_N;
  Dos=(double**)malloc_multidimarray("double",N_Dos,i_Dos);


  /* initialize */
  for (ie=0;ie<Dos_N;ie++) {
    DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
  }
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ie=0;ie<Dos_N;ie++) {
      Dos[spin][ie]=0.0;
    }
  }

  pi2=sqrt(PI);
  iewidth = Dos_Gaussian*3.0/(Dos_emax-Dos_emin)*(Dos_N-1)+3;
  printf("Gaussian width=%d\n",iewidth); 

  for (spin=0;spin<=SpinP_switch;spin++) {

    for (GA_AN=0; GA_AN<atomnum; GA_AN++){

      wanA = WhatSpecies[GA_AN];
      tnoA = Spe_Total_CNO[wanA];

      for (ieg=0; ieg<Msize1[spin][GA_AN]; ieg++) {

        eg = EigenE[spin][GA_AN][ieg];
        x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
        iecenter = (int)x ; 
        iemin = iecenter - iewidth;
        iemax = iecenter + iewidth ;

        if (iemin<0) iemin=0;
        if (iemax>=Dos_N) iemax=Dos_N-1;

        if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax <Dos_N ) {
          sum = 0.0;
          for (i=0; i<tnoA; i++) sum += PDOS[spin][GA_AN][ieg][i];
          for (ie=iemin; ie<=iemax; ie++) {
            xa = (eg-DosE[ie])/Dos_Gaussian; 
            Dos[spin][ie] += sum*exp(- xa*xa)/(Dos_Gaussian*pi2);
	  }
        } 

      }   /* ieg */
    }     /* GA_AN */
  }       /* spin */
 
  /* normalize */
  factor = 1.0/eV2Hartree;
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ie=0;ie<Dos_N;ie++) {
      Dos[spin][ie] = Dos[spin][ie] * factor;
    }
  }

  if (SpinP_switch==1) {
    for (ie=0;ie<Dos_N;ie++) {
      fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie], -Dos[1][ie]);
    }
  }
  else {
    for (ie=0;ie<Dos_N;ie++) {
      fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie]*2.0);
    }
  }

  if (fp_Dos) fclose(fp_Dos);

  free(DosE);
  free_multidimarray(Dos,N_Dos,i_Dos);

}


void Dos_Tetrahedron( char *basename, int Dos_N, 
		      int knum_i,int knum_j, int knum_k,
		      int SpinP_switch, double *****EIGEN, double *******ev,
                      int neg, 
		      double Dos_emin, double Dos_emax, 
		      int iemin, int iemax,
                      int *WhatSpecies, int *Spe_Total_CNO, int atomnum )
{
  int spin,ieg,i,j,k,ie,i1,tnoA,wanA,GA_AN;
  int i_in,j_in,k_in,ic,itetra; 
  double x,result,factor,wval; 
  double cell_e[8], tetra_e[4];
  static int tetra_id[6][4]= { {0,1,2,5}, {1,2,3,5}, {2,3,5,7},
			       {0,2,4,5}, {2,4,5,6}, {2,5,6,7} };

  double *DosE, **Dos;
  int N_Dos,i_Dos[10];

  char file_Dos[YOUSO10];
  FILE *fp_Dos;

  printf("<Dos_Tetrahedron> start\n");

  if (strlen(basename)==0) {
    strcpy(file_Dos,"DOS.Tetrahedron");
  }
  else {
    sprintf(file_Dos,"%s.DOS.Tetrahedron",basename);
  }
  if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
    printf("can not open a file %s\n",file_Dos);
    return;
  }

  printf("<Dos_Tetrahedron> make %s\n",file_Dos);

  /* allocation */
  DosE=(double*)malloc(sizeof(double)*Dos_N);
  N_Dos=2; i_Dos[0]=SpinP_switch+1; i_Dos[1]=Dos_N;
  Dos=(double**)malloc_multidimarray("double",N_Dos,i_Dos);
   
  /* initialize */
  for (ie=0;ie<Dos_N;ie++) {
    DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
  }
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ie=0;ie<Dos_N;ie++) {
      Dos[spin][ie]=0.0;
    }
  }

  /********************************************************************
                           Method = tetrahedron 
  *******************************************************************/

  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ieg=0;ieg<neg; ieg++) {
      for (i=0; i<=(knum_i-1); i++){
	for (j=0; j<=(knum_j-1); j++){
	  for (k=0; k<=(knum_k-1); k++){

            /* generalization for spin non-collinear */
            wval = 0.0;
   	    for (GA_AN=0; GA_AN<atomnum; GA_AN++){
  	      wanA = WhatSpecies[GA_AN];
	      tnoA = Spe_Total_CNO[wanA];
  	      for (i1=0; i1<tnoA; i1++){
                wval += ev[spin][i][j][k][GA_AN][i1][ieg];
	      }
	    }

	    for (i_in=0;i_in<2;i_in++) {
	      for (j_in=0;j_in<2;j_in++) {
		for (k_in=0;k_in<2;k_in++) {
		  cell_e[i_in*4+j_in*2+k_in] = 
		    EIGEN[spin][ (i+i_in)%knum_i ][ (j+j_in)%knum_j ][ (k+k_in)%knum_k ][ieg] ;
		}
	      }
	    }
	    for (itetra=0;itetra<6;itetra++) {
	      for (ic=0;ic<4;ic++) {
		tetra_e[ic]=cell_e[ tetra_id[itetra][ic] ];
	      }
	      OrderE0(tetra_e,4);
	      x = (tetra_e[0]-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1)-1.0;
	      iemin=(int)x;
	      x = (tetra_e[3]-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1)+1.0;
	      iemax=(int)x;
	      if (iemin<0) { iemin=0; }
	      if (iemax>=Dos_N) {iemax=Dos_N-1; }
	      if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax<Dos_N ) {
		for (ie=iemin;ie<=iemax;ie++) {
		  ATM_Dos( tetra_e, &DosE[ie], &result);
		  Dos[spin][ie] += wval*result;
		}
	      }

	    } /* itetra */
	  } /* k */
	} /* j */
      } /* i */
    } /* ieg */
  } /* spin */

  /* normalize */
  factor = 1.0/(double)( eV2Hartree * knum_i * knum_j * knum_k * 6 );
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ie=0;ie<Dos_N;ie++) {
      Dos[spin][ie] = Dos[spin][ie] * factor;
    }
  }

  if (SpinP_switch==1) {
    for (ie=0;ie<Dos_N;ie++) {
      fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie], -Dos[1][ie]);
    }
  } 
  else {
    for (ie=0;ie<Dos_N;ie++) {
      fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie]*2.0);
    }
  }

  if (fp_Dos) fclose(fp_Dos);

  free(DosE);
  free_multidimarray(Dos,N_Dos,i_Dos);


}
   
void Spectra_Gaussian( int pdos_n, int *pdos_atoms, char *basename, int Dos_N, 
		       double Dos_Gaussian, 
		       int knum_i,int knum_j, int knum_k,
		       int SpinP_switch, double *****EIGEN, double *******ev, int neg,
		       double Dos_emin, double Dos_emax,
		       int iemin, int iemax,
		       int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
		       int MaxL,int **Spe_Num_CBasis, int **Spe_Num_Relation)
{

  /*****************************************************************
                           Method = Gaussian 
  *****************************************************************/
  /*     exp( -(x/a)^2 ) 
         a= Dos_Gaussian 
         x=-3a : 3a is counted */


  int Max_tnoA;
  int i,spin,j,k;
  int wanA, GA_AN,tnoA,i1,ieg;
  double eg,x,rval,xa;
  double *DosE, ****Dos;
  int N_Dos, i_Dos[10];
  double pi2;
  int  iewidth,ie,iecenter;
  int iatom,L,M,LM;
  static char *Lname[4]={"s","p", "d","f"};
  double dossum[2],MulP[5],dE;

  char file_Dos[YOUSO10];
  FILE *fp_Dos;
  static char *extension="PDOS.Gaussian";

  printf("<Spectra_Gaussian> start\n");

  pi2=sqrt(PI);
  iewidth = Dos_Gaussian*3.0/(Dos_emax-Dos_emin)*(Dos_N-1)+3;



  Max_tnoA=0;
  for (i=0;i<atomnum;i++) {
    wanA=WhatSpecies[i];
    if (Max_tnoA<Spe_Total_CNO[wanA]) Max_tnoA=Spe_Total_CNO[wanA];
  }
  /* allocation */
  DosE=(double*)malloc(sizeof(double)*Dos_N);

  N_Dos=4; i_Dos[0]=SpinP_switch+1; 
  i_Dos[1] = atomnum; i_Dos[2]= Max_tnoA;
  i_Dos[3]=Dos_N;
  Dos=(double****)malloc_multidimarray("double",N_Dos,i_Dos);


  /* initialize */
  for (ie=0;ie<Dos_N;ie++) {
    DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
  }
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (GA_AN=0; GA_AN<atomnum; GA_AN++){
      wanA = WhatSpecies[GA_AN];
      tnoA = Spe_Total_CNO[wanA];
      for (i1=0; i1<=(tnoA-1); i1++){

        for (ie=0;ie<Dos_N;ie++) {
          Dos[spin][GA_AN][i1][ie]=0.0;
        }
      }
    }
  }

  /********************************************************************
                           Method = gaussian
  *******************************************************************/

  i=j=k=0;
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ieg=0;ieg<neg; ieg++) {
      eg=EIGEN[spin][i][j][k][ieg];

      x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
      iecenter = (int)x ;
      iemin = iecenter - iewidth;
      iemax = iecenter + iewidth ;
      if (iemin<0) iemin=0;
      if (iemax>=Dos_N) iemax=Dos_N-1;


         
      for (GA_AN=0; GA_AN<atomnum; GA_AN++){
	wanA = WhatSpecies[GA_AN];
	tnoA = Spe_Total_CNO[wanA];
	for (i1=0; i1<=(tnoA-1); i1++){
 
	  rval = ev[spin][i][j][k][GA_AN][i1][ieg];

	  if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax <Dos_N ) {
	    for (ie=iemin;ie<=iemax;ie++) {
	      xa = (eg-DosE[ie])/Dos_Gaussian; 

	      Dos[spin][GA_AN][i1][ie] += rval* exp(- xa*xa)/(Dos_Gaussian*pi2*eV2Hartree);
	    }
	  } /* if */
	} /* i1 */
      } /* GA_AN */
    } /* ieg */
  } /* spin */


  /* orbital projection */
  for (iatom=0;iatom<pdos_n;iatom++) {
    GA_AN = pdos_atoms[iatom]-1;
    wanA = WhatSpecies[GA_AN];
    tnoA = Spe_Total_CNO[wanA];
    for (L=0;L<=MaxL; L++) {
      if (Spe_Num_CBasis[wanA][L]>0) {
        for (M=0;M<2*L+1;M++) {
          LM=100*L+M+1;
	  sprintf(file_Dos,"%s.%s.atom%d.%s%d",basename,extension,pdos_atoms[iatom],
		  Lname[L],M+1);
	  if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL )  {
	    printf("can not open %s\n",file_Dos);
	  }
	  else {
	    printf("make %s\n",file_Dos);
	    for (ie=0;ie<Dos_N;ie++) {
	      for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
	      for (i1=0;i1< tnoA; i1++) {
		if (LM== Spe_Num_Relation[wanA][i1]) {
		  for (spin=0;spin<=SpinP_switch;spin++) 
		    dossum[spin]+=  Dos[spin][GA_AN][i1][ie];
		}
	      }
	      if (SpinP_switch==1) {
		fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree,
                        dossum[0],-dossum[1]);
	      }
	      else {
		fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0);
	      }

	    }
	  }
	  if (fp_Dos) fclose(fp_Dos);
	} /* M */
      }
    } /* L */
  } /* iatom */

  /* atom projection */
  for (iatom=0;iatom<pdos_n;iatom++) {
    GA_AN = pdos_atoms[iatom]-1;
    wanA = WhatSpecies[GA_AN];
    tnoA = Spe_Total_CNO[wanA];

    sprintf(file_Dos,"%s.%s.atom%d",basename,extension,pdos_atoms[iatom],
	    Lname[L]);
    if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL )  {
      printf("can not open %s\n",file_Dos);
    }
    else {
      printf("make %s\n",file_Dos);

      /*
      for (spin=0; spin<=SpinP_switch; spin++){
        MulP[spin] = 0.0;
      }
      */

      for (ie=0;ie<Dos_N;ie++) {
	for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
	for (i1=0;i1< tnoA; i1++) {
	  for (spin=0;spin<=SpinP_switch;spin++)
	    dossum[spin]+=  Dos[spin][GA_AN][i1][ie];
	}

        /* integration */
	/*
        if (DosE[ie]<=0.00){
          dE = DosE[1] - DosE[0];
          for (spin=0; spin<=SpinP_switch; spin++){
            MulP[spin] += dossum[spin]*dE;
          }
	}
	*/

	if (SpinP_switch==1) {
	  fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree,
		  dossum[0],-dossum[1]);
	}
	else {
	  fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0);
	}

      } /* ie */
    } /* fp_Dos */
    if (fp_Dos) fclose(fp_Dos);


    /*
    for (spin=0; spin<=SpinP_switch; spin++){
      printf("spin=%2d MulP=%15.12f\n",spin,MulP[spin]*eV2Hartree);
    }
    */


  } /* iatom */
  

  /*********************************************************
                        close and free
  *********************************************************/

  free(DosE);
  free_multidimarray(Dos,N_Dos,i_Dos);


}


void SpectraDC_Gaussian( int pdos_n, int *pdos_atoms, char *basename, int Dos_N, 
  		         double Dos_Gaussian, 
		         int knum_i,int knum_j, int knum_k,
		         int SpinP_switch, double *****EIGEN, double *******ev, int neg,
		         double Dos_emin, double Dos_emax,
		         int iemin, int iemax,
		         int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
		         int MaxL,int **Spe_Num_CBasis, int **Spe_Num_Relation)
{

  /*****************************************************************
                           Method = Gaussian 
  *****************************************************************/
  /*     exp( -(x/a)^2 ) 
         a= Dos_Gaussian 
         x=-3a : 3a is counted */


  int Max_tnoA;
  int i,spin,j,k;
  int wanA,GA_AN,tnoA,i1,ieg;
  double eg,x,rval,xa,tmp1;
  double *DosE, ****Dos;
  int N_Dos, i_Dos[10];
  double pi2;
  int  iewidth,ie,iecenter;
  int iatom,L,M,LM;
  static char *Lname[4]={"s","p", "d","f"};
  double dossum[2];

  char file_Dos[YOUSO10];
  FILE *fp_Dos;
  static char *extension="PDOS.Gaussian";

  printf("<Spectra_Gaussian> start\n");

  pi2=sqrt(PI);
  iewidth = Dos_Gaussian*3.0/(Dos_emax-Dos_emin)*(Dos_N-1)+3;

  Max_tnoA=0;
  for (i=0;i<atomnum;i++) {
    wanA=WhatSpecies[i];
    if (Max_tnoA<Spe_Total_CNO[wanA]) Max_tnoA=Spe_Total_CNO[wanA];
  }
  /* allocation */
  DosE=(double*)malloc(sizeof(double)*Dos_N);

  Dos = (double****)malloc(sizeof(double***)*(SpinP_switch+1));
  for (spin=0; spin<=SpinP_switch; spin++) {
    Dos[spin] =(double***)malloc(sizeof(double**)*pdos_n);
    for (iatom=0;iatom<pdos_n;iatom++) {
      GA_AN = pdos_atoms[iatom]-1;
      wanA = WhatSpecies[GA_AN];
      tnoA = Spe_Total_CNO[wanA];
      Dos[spin][iatom] =(double**)malloc(sizeof(double*)*tnoA);
      for (i=0; i<tnoA; i++){
        Dos[spin][iatom][i] =(double*)malloc(sizeof(double)*Dos_N);
        for (ie=0; ie<Dos_N; ie++) {
          Dos[spin][iatom][i][ie] =0.0;
	}
      }
    }
  }

  /* initialize */
  for (ie=0;ie<Dos_N;ie++) {
    DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
  }

  /********************************************************************
                           Method = gaussian
  *******************************************************************/

  for (spin=0;spin<=SpinP_switch;spin++) {
    for (iatom=0;iatom<pdos_n;iatom++) {

      GA_AN = pdos_atoms[iatom]-1;
      wanA = WhatSpecies[GA_AN];
      tnoA = Spe_Total_CNO[wanA];

      for (ieg=0; ieg<Msize1[spin][GA_AN]; ieg++) {

        eg = EigenE[spin][GA_AN][ieg];
        x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
        iecenter = (int)x ; 
        iemin = iecenter - iewidth;
        iemax = iecenter + iewidth ;
        if (iemin<0) iemin=0;
        if (iemax>=Dos_N) iemax=Dos_N-1;

        if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax <Dos_N ) {
          for (ie=iemin; ie<=iemax; ie++) {
            xa = (eg-DosE[ie])/Dos_Gaussian; 

            tmp1 = exp( -xa*xa)/(Dos_Gaussian*pi2*eV2Hartree);

            for (i=0; i<tnoA; i++){
              Dos[spin][iatom][i][ie] += PDOS[spin][GA_AN][ieg][i]*tmp1;
	    }
	  }
	}

      } /* ieg */
    }   /* iatom */
   }    /* spin */

  /* orbital projection */
  for (iatom=0;iatom<pdos_n;iatom++) {
    GA_AN = pdos_atoms[iatom]-1;
    wanA = WhatSpecies[GA_AN];
    tnoA = Spe_Total_CNO[wanA];
    for (L=0;L<=MaxL; L++) {
      if (Spe_Num_CBasis[wanA][L]>0) {
        for (M=0;M<2*L+1;M++) {
          LM=100*L+M+1;
	  sprintf(file_Dos,"%s.%s.atom%d.%s%d",basename,extension,pdos_atoms[iatom],
		  Lname[L],M+1);
	  if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL )  {
	    printf("can not open %s\n",file_Dos);
	  }
	  else {
	    printf("make %s\n",file_Dos);
	    for (ie=0;ie<Dos_N;ie++) {
	      for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
	      for (i1=0;i1< tnoA; i1++) {
		if (LM== Spe_Num_Relation[wanA][i1]) {
		  for (spin=0;spin<=SpinP_switch;spin++) 
		    dossum[spin]+=  Dos[spin][iatom][i1][ie];
		}
	      }
	      if (SpinP_switch==1) {
		fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree,
                        dossum[0],-dossum[1]);
	      }
	      else {
		fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0);
	      }

	    }
	  }
	  if (fp_Dos) fclose(fp_Dos);
	} /* M */
      }
    }     /* L */
  }       /* iatom */

  /* atom projection */
  for (iatom=0;iatom<pdos_n;iatom++) {
    GA_AN = pdos_atoms[iatom]-1;
    wanA = WhatSpecies[GA_AN];
    tnoA = Spe_Total_CNO[wanA];

    sprintf(file_Dos,"%s.%s.atom%d",basename,extension,pdos_atoms[iatom],
	    Lname[L]);
    if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL )  {
      printf("can not open %s\n",file_Dos);
    }
    else {
      printf("make %s\n",file_Dos);
      for (ie=0;ie<Dos_N;ie++) {
	for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
	for (i1=0;i1< tnoA; i1++) {
	  for (spin=0;spin<=SpinP_switch;spin++)
	    dossum[spin]+=  Dos[spin][iatom][i1][ie];
	}
	if (SpinP_switch==1) {
	  fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree,
		  dossum[0],-dossum[1]);
	}
	else {
	  fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0);
	}

      } /* ie */
    } /* fp_Dos */
    if (fp_Dos) fclose(fp_Dos);
  } /* iatom */
  

  /*********************************************************
                        close and free
  *********************************************************/

  free(DosE);

  for (spin=0; spin<=SpinP_switch; spin++) {
    for (iatom=0;iatom<pdos_n;iatom++) {
      GA_AN = pdos_atoms[iatom]-1;
      wanA = WhatSpecies[GA_AN];
      tnoA = Spe_Total_CNO[wanA];
      for (i=0; i<tnoA; i++){
        free(Dos[spin][iatom][i]);
      }
      free(Dos[spin][iatom]);
    }
    free(Dos[spin]);
  }
  free(Dos);

}

void Spectra_Tetrahedron( int pdos_n, int *pdos_atoms, char *basename, int Dos_N, 
			  int knum_i,int knum_j, int knum_k,
			  int SpinP_switch, double *****EIGEN, double *******ev, int neg, 
			  double Dos_emin, double Dos_emax, 
			  int iemin, int iemax,
			  int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
                          int MaxL,int **Spe_Num_CBasis, int **Spe_Num_Relation)
{
  int spin,ieg,i,j,k,ie;
  int i_in,j_in,k_in,ic,itetra; 
  double x,result,factor; 
  double cell_e[8], tetra_e[4],tetra_a[4];
  static int tetra_id[6][4]= { {0,1,2,5}, {1,2,3,5}, {2,3,5,7},
			       {0,2,4,5}, {2,4,5,6}, {2,5,6,7} };
  int wanA,Max_tnoA,GA_AN,tnoA,i1,iatom;
  int L,M,LM;
  static char *Lname[4]={"s","p", "d","f"};
  double dossum[2];

  double *DosE, ****Dos, rval;
  int N_Dos,i_Dos[10];
  double ***cell_a;
  int N_cell_a, i_cell_a[10];

  char file_Dos[YOUSO10];
  FILE *fp_Dos;

  static char *extension="PDOS.Tetrahedron";

  printf("<Spectra_Tetrahedron> start\n");

#if 0
  if (strlen(basename)==0) {
    strcpy(file_Dos,extension);
  }
  else {
    sprintf(file_Dos,"%s.%s",basename,extension);
  }

  if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
    printf("can not open a file %s\n",file_Dos);
    return;
  }

  printf("<Spectra_Tetrahedron> make %s\n",file_Dos);
#endif

  Max_tnoA=Spe_Total_CNO[0];
  for (i=0;i<atomnum;i++) {
    wanA=WhatSpecies[i];
    if (Max_tnoA<Spe_Total_CNO[wanA]) Max_tnoA=Spe_Total_CNO[wanA];
  }



  /* allocation */
  DosE=(double*)malloc(sizeof(double)*Dos_N);

  N_Dos=4; i_Dos[0]=SpinP_switch+1; 
  i_Dos[1] = atomnum; i_Dos[2]= Max_tnoA;
  i_Dos[3]=Dos_N;
  Dos=(double****)malloc_multidimarray("double",N_Dos,i_Dos);


  N_cell_a=3; i_cell_a[0]= atomnum; i_cell_a[1]=Max_tnoA; i_cell_a[2]= 8;
  cell_a = (double***) malloc_multidimarray("double",N_cell_a,i_cell_a);
   
  /* initialize */
  for (ie=0;ie<Dos_N;ie++) {
    DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
  }
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (GA_AN=0; GA_AN<atomnum; GA_AN++){
      wanA = WhatSpecies[GA_AN];
      tnoA = Spe_Total_CNO[wanA];
      for (i1=0; i1<=(tnoA-1); i1++){

	for (ie=0;ie<Dos_N;ie++) {
	  Dos[spin][GA_AN][i1][ie]=0.0;
	}
      }
    }
  }

  /********************************************************************
                          Method = tetrahedron 
  *******************************************************************/
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (ieg=0;ieg<neg; ieg++) {
      for (i=0; i<=(knum_i-1); i++){
	for (j=0; j<=(knum_j-1); j++){
	  for (k=0; k<=(knum_k-1); k++){

	    for (i_in=0;i_in<2;i_in++) {
	      for (j_in=0;j_in<2;j_in++) {
		for (k_in=0;k_in<2;k_in++) {
		  cell_e[i_in*4+j_in*2+k_in] = 
		    EIGEN[spin][ (i+i_in)%knum_i ][ (j+j_in)%knum_j ][ (k+k_in)%knum_k ][ieg] ;
                  for (GA_AN=0; GA_AN<atomnum; GA_AN++){
                    wanA = WhatSpecies[GA_AN];
                    tnoA = Spe_Total_CNO[wanA];
                    for (i1=0; i1<=(tnoA-1); i1++){
                      rval= ev[spin][ (i+i_in)%knum_i ][ (j+j_in)%knum_j ][ (k+k_in)%knum_k ][GA_AN][i1][ieg];
                      cell_a[GA_AN][i1][i_in*4+j_in*2+k_in] = rval;
                    }
                  }
		}
	      }
	    }

	    for (GA_AN=0; GA_AN<atomnum; GA_AN++){
	      wanA = WhatSpecies[GA_AN];
	      tnoA = Spe_Total_CNO[wanA];
	      for (i1=0; i1<=(tnoA-1); i1++){

		for (itetra=0;itetra<6;itetra++) {
		  for (ic=0;ic<4;ic++) {
		    tetra_e[ic]=cell_e[ tetra_id[itetra][ic] ];
		    tetra_a[ic]=cell_a[ GA_AN ][i1][ tetra_id[itetra][ic] ];
		  }
		  OrderE(tetra_e,tetra_a,4);

		  for (ie=0;ie<4;ie++) 
#if 0
 	          printf("%lf(%lf) ",tetra_e[ie],tetra_a[ie]);
		  printf("\n");
#endif

		  x = (tetra_e[0]-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1)-1.0;
		  iemin=(int)x;
		  x = (tetra_e[3]-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1)+1.0;
		  iemax=(int)x;
		  if (iemin<0) { iemin=0; }
		  if (iemax>=Dos_N) {iemax=Dos_N-1; }
#if 0
		  printf("%lf %lf %lf %lf\n", DosE[iemin],tetra_e[0],tetra_e[3], DosE[iemax] );
#endif
		  if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax<Dos_N ) {
		    for (ie=iemin;ie<=iemax;ie++) {
		      ATM_Spectrum( tetra_e,tetra_a, &DosE[ie], &result);
		      Dos[spin][GA_AN][i1][ie] += result;
#if 0
		      printf("%lf-> %lf\n",DosE[ie],result);
#endif
		    }
		  }
		} /* itetra */
	      } 
	    }

	  } /* k */
	} /* j */
      } /* i */
    } /* ieg */
  } /* spin */

  /* normalize */
  factor = 1.0/(double)( eV2Hartree * knum_i * knum_j * knum_k * 6 );
  for (spin=0;spin<=SpinP_switch;spin++) {
    for (GA_AN=0; GA_AN<atomnum; GA_AN++){
      wanA = WhatSpecies[GA_AN];
      tnoA = Spe_Total_CNO[wanA];
      for (i1=0; i1<=(tnoA-1); i1++){

	for (ie=0;ie<Dos_N;ie++) {
	  Dos[spin][GA_AN][i1][ie] = Dos[spin][GA_AN][i1][ie] * factor;
	}
      }
    }
  }

  /***************************************************
                         output 
  ***************************************************/

#if 0
  if (SpinP_switch==1) {
    for (ie=0;ie<Dos_N;ie++) {
      for (sum=0.0,sum2=0.0,GA_AN=0; GA_AN<atomnum; GA_AN++){
	wanA = WhatSpecies[GA_AN];
	tnoA = Spe_Total_CNO[wanA];
	for (i1=0; i1<=(tnoA-1); i1++){
	  sum += Dos[0][GA_AN][i1][ie];
	  sum2 += Dos[1][GA_AN][i1][ie];
	}
      }

      fprintf(fp_Dos,"%lf %lf %lf", (DosE[ie])*eV2Hartree,
	      sum, -sum2 );
    }
  } 
  else {
    for (ie=0;ie<Dos_N;ie++) {
      for (sum=0.0,sum2=0.0,GA_AN=0; GA_AN<atomnum; GA_AN++){
	wanA = WhatSpecies[GA_AN];
	tnoA = Spe_Total_CNO[wanA];
	for (i1=0; i1<=(tnoA-1); i1++){
	  sum += Dos[0][GA_AN][i1][ie];
	}
      }

      fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree, sum*2.0);
    }
  }
#else

  for (iatom=0;iatom<pdos_n;iatom++) {
    GA_AN = pdos_atoms[iatom]-1;
    wanA = WhatSpecies[GA_AN];
    tnoA = Spe_Total_CNO[wanA];
    for (i1=0;i1< tnoA; i1++) {
      printf("Spe_Num_Relation %d %d %d\n",wanA,i1,Spe_Num_Relation[wanA][i1]);
    }
  }

  /* orbital projection */
  for (iatom=0;iatom<pdos_n;iatom++) {
    GA_AN = pdos_atoms[iatom]-1;
    wanA = WhatSpecies[GA_AN];
    tnoA = Spe_Total_CNO[wanA];
    for (L=0;L<=MaxL; L++) {
      if (Spe_Num_CBasis[wanA][L]>0) {
        for (M=0;M<2*L+1;M++) {
          LM=100*L+M+1;
	  sprintf(file_Dos,"%s.%s.atom%d.%s%d",basename,extension,pdos_atoms[iatom],
		  Lname[L],M+1);
	  if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL )  {
	    printf("can not open %s\n",file_Dos);
	  }
	  else {
	    printf("make %s\n",file_Dos);
	    for (ie=0;ie<Dos_N;ie++) {
	      for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
	      for (i1=0;i1< tnoA; i1++) {
		/*          if (i1==4) {
			    printf("%d %d\n",LM,Spe_Num_Relation[wanA][i1]);
			    } 
		*/
		if (LM== Spe_Num_Relation[wanA][i1]) {
		  for (spin=0;spin<=SpinP_switch;spin++) 
		    dossum[spin]+=  Dos[spin][GA_AN][i1][ie];
		}
	      }
	      if (SpinP_switch==1) {
		fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, 
			dossum[0],dossum[1]);
	      }
	      else {
		fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree,
			dossum[0]*2.0);
	      }

	    }
	  }
	} /* M */
      }
    } /* L */
  } /* iatom */
 

  /* atom projection */
  for (iatom=0;iatom<pdos_n;iatom++) {
    GA_AN = pdos_atoms[iatom]-1;
    wanA = WhatSpecies[GA_AN];
    tnoA = Spe_Total_CNO[wanA];

    sprintf(file_Dos,"%s.%s.atom%d",basename,extension,pdos_atoms[iatom],
	    Lname[L]);
    if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL )  {
      printf("can not open %s\n",file_Dos);
    }
    else {
      printf("make %s\n",file_Dos);
      for (ie=0;ie<Dos_N;ie++) {
	for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
	for (i1=0;i1< tnoA; i1++) {
	  for (spin=0;spin<=SpinP_switch;spin++)
	    dossum[spin]+=  Dos[spin][GA_AN][i1][ie];
	}
	if (SpinP_switch==1) {
	  fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree,
		  dossum[0],dossum[1]);
	}
	else {
	  fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree,
		  dossum[0]*2.0);
	}

      } /* ie */
    } /* fp_Dos */
  } /* iatom */



#endif


  /*********************************************************
                       close and free
  *********************************************************/

#if 0
  if (fp_Dos) fclose(fp_Dos);
#endif

  free(DosE);
  free_multidimarray(Dos,N_Dos,i_Dos);
  free_multidimarray(cell_a,N_cell_a,i_cell_a);



}
   

/* input    int SpeciesNum,int MaxL,int *Spe_Total_CNO,int **Spe_Num_CBasis,
   output Spe_Num_Relation
*/
void  Spe_Num_CBasis2Relation(
			      int SpeciesNum,int MaxL,int *Spe_Total_CNO,int **Spe_Num_CBasis,
			      int **Spe_Num_Relation)
{
  int i,j,k,l,id;

  for (i=0;i<SpeciesNum;i++) {
    id=0;
    for (j=0;j<=MaxL;j++) {
      if (Spe_Num_CBasis[i][j]>0) {
	for (k=0;k<Spe_Num_CBasis[i][j];k++) {
	  for (l=0;l<2*j+1;l++) {
	    Spe_Num_Relation[i][id++]=100*j+l+1;
	  }
	}
      } 
    } 
    for (j=0;j<Spe_Total_CNO[i];j++) {
      printf("%d ",Spe_Num_Relation[i][j]);
    }
    printf("\n");
  }


}


void input_file_eg(char *file_eg,
		   int *mode,  int *n, double *emin, double *emax, int *iemin,int *iemax, 
		   int *SpinP_switch, int Kgrid[3], int *atomnum, int **WhatSpecies,
		   int *SpeciesNum, int **Spe_Total_CNO, int * Max_tnoA, int *MaxL,
		   int ***Spe_Num_CBasis, int ***Spe_Num_Relation, double *ChemP, 
		   double ******ko)
{
  double r_vec[10], r_vec2[10];
  int i_vec[10], i_vec2[10];
  FILE *fp;
  int N_Spe_Num_CBasis, i_Spe_Num_CBasis[10];
  int j,k;
  int N_Spe_Num_Relation, i_Spe_Num_Relation[10];
  int N_ko, i_ko[10];
  int i,i1,j1,k1,ie,spin;

  input_open(file_eg);

  input_int("mode",mode,1);

  /*****************************************
          In case of divide-conquer
  *****************************************/

  if (*mode==5){

    input_int("Nspin",SpinP_switch,0);

    r_vec2[0]=0.0; r_vec2[1]=0.0;
    input_doublev("Erange",2,r_vec,r_vec2);
    *emin=r_vec[0]; *emax=r_vec[1];

    input_int("atomnum",atomnum,0);

    *WhatSpecies = (int*)malloc(sizeof(int)* (*atomnum));
    if (fp=input_find("<WhatSpecies") ) {
      for (i=0;i<*atomnum;i++) {
	fscanf( fp, "%d", &(*WhatSpecies)[i]);
      }
      if (!input_last("WhatSpecies>")) {
	printf("format error near WhatSpecies>\n");
	exit(10);
      }
    }
    else {
      printf("<WhatSpecies is not found\n");
      exit(10);
    }

    input_int("SpeciesNum",SpeciesNum,0);
    *Spe_Total_CNO= (int*)malloc(sizeof(int)*(*SpeciesNum));
    if (fp=input_find("<Spe_Total_CNO") ) {
      for (i=0;i<*SpeciesNum;i++) {
	fscanf( fp, "%d", &(*Spe_Total_CNO)[i]);
      }
      if (!input_last("Spe_Total_CNO>")) {
	printf("format error near Spe_Total_CNO>\n");
	exit(10);
      }
    }
    else {
      printf("<Spe_Total_CNO is not found\n");
      exit(10);
    }

    *Max_tnoA = (*Spe_Total_CNO)[0];
    for (i=1;i<*SpeciesNum;i++) {
      if ( *Max_tnoA  < (*Spe_Total_CNO)[i] ) *Max_tnoA=(*Spe_Total_CNO)[i];
    }
    printf("Max of Spe_Total_CNO = %d\n",*Max_tnoA);

    input_int("MaxL",MaxL,3);
    N_Spe_Num_CBasis=2; i_Spe_Num_CBasis[0]=*SpeciesNum; i_Spe_Num_CBasis[1]=*MaxL+1;
    *Spe_Num_CBasis=(int**)malloc_multidimarray("int",N_Spe_Num_CBasis, i_Spe_Num_CBasis);
    if (fp=input_find("<Spe_Num_CBasis") ) {
      for (i=0;i<*SpeciesNum;i++) {
	for (j=0;j<=*MaxL;j++) {
	  fscanf(fp,"%d",&(*Spe_Num_CBasis)[i][j]);
	}
      }
      if (!input_last("Spe_Num_CBasis>")) {
	printf("format error near Spe_Num_CBasis>\n");
	exit(10);
      }
    }
    else {
      printf("<Spe_Num_CBasis is not found\n");
      exit(10);
    }

    /* relation, index to orbital, 0->s, 1->p ... and so on. */
    N_Spe_Num_Relation=2; i_Spe_Num_Relation[0]=*SpeciesNum; i_Spe_Num_Relation[1]=*Max_tnoA;
    *Spe_Num_Relation=(int**)malloc_multidimarray("int",N_Spe_Num_Relation, i_Spe_Num_Relation);
    Spe_Num_CBasis2Relation(*SpeciesNum,*MaxL,*Spe_Total_CNO,*Spe_Num_CBasis, *Spe_Num_Relation);


    input_double("ChemP",ChemP,0.0);

  }

  /*****************************************
         In case of cluster and band
  *****************************************/

  else {

    input_int("N",n,0);
    r_vec2[0]=0.0; r_vec2[1]=0.0;
    input_doublev("Erange",2,r_vec,r_vec2);

    *emin=r_vec[0]; *emax=r_vec[1];
    i_vec2[0]=0; i_vec2[1]=0; i_vec2[2]=0;

    input_intv("irange" , 2,i_vec,i_vec2);
    *iemin=i_vec[0];  *iemax= i_vec[1];

    input_int("Nspin",SpinP_switch,0);

    i_vec2[0]=i_vec2[1]=i_vec2[2]=0;
    input_intv("Kgrid",3, Kgrid,i_vec2);
  
    input_int("atomnum",atomnum,0);
    *WhatSpecies = (int*)malloc(sizeof(int)* (*atomnum));
    if (fp=input_find("<WhatSpecies") ) {
      for (i=0;i<*atomnum;i++) {
	fscanf( fp, "%d", &(*WhatSpecies)[i]);
      }
      if (!input_last("WhatSpecies>")) {
	printf("format error near WhatSpecies>\n");
	exit(10);
      }
    }
    else {
      printf("<WhatSpecies is not found\n");
      exit(10);
    }
  
    input_int("SpeciesNum",SpeciesNum,0);
    *Spe_Total_CNO= (int*)malloc(sizeof(int)*(*SpeciesNum));
    if (fp=input_find("<Spe_Total_CNO") ) {
      for (i=0;i<*SpeciesNum;i++) {
	fscanf( fp, "%d", &(*Spe_Total_CNO)[i]);
      }
      if (!input_last("Spe_Total_CNO>")) {
	printf("format error near Spe_Total_CNO>\n");
	exit(10);
      }
    }
    else {
      printf("<Spe_Total_CNO is not found\n");
      exit(10);
    }

    *Max_tnoA = (*Spe_Total_CNO)[0];
    for (i=1;i<*SpeciesNum;i++) {
      if ( *Max_tnoA  < (*Spe_Total_CNO)[i] ) *Max_tnoA=(*Spe_Total_CNO)[i];
    }
    printf("Max of Spe_Total_CNO = %d\n",*Max_tnoA);

    input_int("MaxL",MaxL,3);
    N_Spe_Num_CBasis=2; i_Spe_Num_CBasis[0]=*SpeciesNum; i_Spe_Num_CBasis[1]=*MaxL+1;
    *Spe_Num_CBasis=(int**)malloc_multidimarray("int",N_Spe_Num_CBasis, i_Spe_Num_CBasis);
    if (fp=input_find("<Spe_Num_CBasis") ) {
      for (i=0;i<*SpeciesNum;i++) {
	for (j=0;j<=*MaxL;j++) {
	  fscanf(fp,"%d",&(*Spe_Num_CBasis)[i][j]);
	}
      }
      if (!input_last("Spe_Num_CBasis>")) {
	printf("format error near Spe_Num_CBasis>\n");
	exit(10);
      }
    }
    else {
      printf("<Spe_Num_CBasis is not found\n");
      exit(10);
    }

    /* relation, index to orbital, 0->s, 1->p ... and so on. */
    N_Spe_Num_Relation=2; i_Spe_Num_Relation[0]=*SpeciesNum; i_Spe_Num_Relation[1]=*Max_tnoA;
    *Spe_Num_Relation=(int**)malloc_multidimarray("int",N_Spe_Num_Relation, i_Spe_Num_Relation);
    Spe_Num_CBasis2Relation(*SpeciesNum,*MaxL,*Spe_Total_CNO,*Spe_Num_CBasis, *Spe_Num_Relation);
  
    input_double("ChemP",ChemP,0.0);

    if (fp=input_find("<Eigenvalues") ) {
      N_ko=5; 
      i=0; i_ko[i++]=*SpinP_switch+1; 
      i_ko[i++]=Kgrid[0]; i_ko[i++]=Kgrid[1]; i_ko[i++]=Kgrid[2]; i_ko[i++]= *iemax-*iemin+1; 
      *ko=(double*****)malloc_multidimarray("double", N_ko,i_ko);

      for (i=0;i<Kgrid[0];i++) {
	for (j=0;j<Kgrid[1];j++) {
	  for (k=0;k<Kgrid[2];k++) {
	    for (spin=0;spin<=*SpinP_switch;spin++) {
	      fscanf(fp,"%d%d%d",&i1,&j1,&k1); 
	      if (i1!=i || j1!=j || k1!=k) {
		printf("format error in Eigenvalues\n");
		exit(10);
	      }
	      for (ie=*iemin;ie<=*iemax;ie++) {
		fscanf(fp,"%lf",&(*ko)[spin][i][j][k][ie-*iemin]);
		/*   printf("%lf ",ko[spin][i][j][k][ie-iemin]); */
		(*ko)[spin][i][j][k][ie-*iemin]=(*ko)[spin][i][j][k][ie-*iemin]-*ChemP;

	      }
	    } /* spin */
	  } /* k*/
	} /* j */
      } /* i */
      if (!input_last("Eigenvalues>")) {
	printf("format error near Eigenvalues>\n");
	exit(10);
      }
    } /* if */
    else {
      printf("can not find <Eigenvalues\n");
      exit(10);
    }

  }

  input_close();
}


void input_file_ev( char *file_ev, int mode,  int Max_tnoA, 
		    int Kgrid[3], int SpinP_switch, int iemin, int iemax, 
                    int atomnum,  int *WhatSpecies, int *Spe_Total_CNO, 
		    double ********ev, double ChemP)

{

  int N_ev, i_ev[10];
  int i,j,k,spin,ie,i1,j1,k1,n;
  int GA_AN, wanA, tnoA;
  int i_vec[10];
  float *fSD;
  double sum;
  char keyword[YOUSO10];

  FILE *fp; 
  
  /*****************************************
          In case of divide-conquer
  *****************************************/
    
  if (mode==5){

    input_open(file_ev);

    Msize1 = (int**)malloc(sizeof(int*)*(SpinP_switch+1));
    EigenE = (double***)malloc(sizeof(double**)*(SpinP_switch+1));
    PDOS = (double****)malloc(sizeof(double***)*(SpinP_switch+1));

    for (spin=0; spin<=SpinP_switch; spin++){

      Msize1[spin] = (int*)malloc(sizeof(int)*atomnum);
      EigenE[spin] = (double**)malloc(sizeof(double*)*atomnum);
      PDOS[spin] = (double***)malloc(sizeof(double**)*atomnum);

      for (GA_AN=0; GA_AN<atomnum; GA_AN++){

	sprintf(keyword,"<AN%dAN%d",GA_AN+1,spin);

	if (fp=input_find(keyword) ) {
	  wanA = WhatSpecies[GA_AN];
	  tnoA = Spe_Total_CNO[wanA];

	  fscanf(fp,"%d",&Msize1[spin][GA_AN]);

	  EigenE[spin][GA_AN] = (double*)malloc(sizeof(double)*Msize1[spin][GA_AN]);
	  PDOS[spin][GA_AN] = (double**)malloc(sizeof(double*)*Msize1[spin][GA_AN]);

	  for (n=0; n<Msize1[spin][GA_AN]; n++){

	    PDOS[spin][GA_AN][n] = (double*)malloc(sizeof(double)*tnoA);

	    fscanf(fp,"%d", &j);
	    fscanf(fp,"%lf",&EigenE[spin][GA_AN][n]);
            EigenE[spin][GA_AN][n] = EigenE[spin][GA_AN][n] - ChemP;

	    for (i=0; i<tnoA; i++) {
	      fscanf(fp,"%lf",&PDOS[spin][GA_AN][n][i]);
	    }
	  }

	  sprintf(keyword,"AN%dAN%d>",GA_AN+1,spin);
	  if (!input_last(keyword)) {
	    printf("format error near AN%dAN>\n",GA_AN+1);
	    exit(10);
	  }

	}
	else {
	  printf("<AN%dAN is not found\n",GA_AN+1);
	  exit(10);
	}
      }
    }

    input_close();
  
  }  
  
  /*****************************************
         In case of cluster and band
  *****************************************/

  else {

    if ( (fp=fopen(file_ev,"r"))==NULL ) {
      printf("cannot open a file %s\n", file_ev);
      exit(10);
    }

    N_ev=7;
    i=0; i_ev[i++]=SpinP_switch+1;
    i_ev[i++]=Kgrid[0]; i_ev[i++]=Kgrid[1]; i_ev[i++]=Kgrid[2]; 
    i_ev[i++] = atomnum ;  i_ev[i++] = Max_tnoA;
    i_ev[i++]= iemax-iemin+1;
    (*ev)=(double*******)malloc_multidimarray("double", N_ev,i_ev);

    fSD=(float*)malloc(sizeof(float)*Max_tnoA);

    for (i=0;i<Kgrid[0];i++) {
      for (j=0;j<Kgrid[1];j++) {
	for (k=0;k<Kgrid[2];k++) {
	  for (spin=0;spin<=SpinP_switch;spin++) {
	    for (ie=iemin;ie<=iemax;ie++) {

#ifdef BINARYWRITE
	      fread(i_vec,sizeof(int),3,fp);
	      i1=i_vec[0]; j1=i_vec[1]; k1=i_vec[2];
#else
	      fscanf(fp,"%d%d%d",&i1,&j1,&k1);
#endif
	      if (i1!=i || j1!=j || k1!=k) {
		printf("parameter error in reading Eigenvectors\n");
		exit(10);
	      }


	      for (GA_AN=0; GA_AN<atomnum; GA_AN++){
		wanA = WhatSpecies[GA_AN];
		tnoA = Spe_Total_CNO[wanA];
#ifdef BINARYWRITE
		fread(fSD,sizeof(float),tnoA,fp);
		for (i1=0; i1<=(tnoA-1); i1++){
		  (*ev)[spin][i][j][k][GA_AN][i1][ie-iemin] = fSD[i1];
		}
#else
		for (i1=0; i1<=(tnoA-1); i1++){
		  fscanf(fp,"%lf",&(*ev)[spin][i][j][k][GA_AN][i1][ie-iemin]);
		  /* printf("%lf ",ev[spin][i][j][k][GA_AN][i1][ie-iemin]); */
		}
#endif
	      }
	      /*   printf("\n"); */
	      /* check normalization */
	      sum=0.0;
	      for (GA_AN=0; GA_AN<atomnum; GA_AN++){
		wanA = WhatSpecies[GA_AN];
		tnoA = Spe_Total_CNO[wanA];
		for (i1=0; i1<=(tnoA-1); i1++){
		  sum+= (*ev)[spin][i][j][k][GA_AN][i1][ie-iemin];
		}
	      }

	      /*
	      if (fabs(sum-1.0) > 1.0e-3) {
                printf("%d %d %d %d %d, sum=%lf\n",spin,i,j,k,ie,sum);
              }
	      */

	    }
	  }
	}
      }
    } 
    free(fSD);
    fclose(fp);
  }


}



void input_main( int mode, int Kgrid[3], int atomnum, 
		 int *method, int *todo, 
		 double *gaussian, 
		 int *pdos_n, int **pdos_atoms)
{
  char buf[1000],buf2[1000],*c;
  int i;

  if (mode==5){
    printf("The input data caluculated by the divide-conquer method\n"); 
    printf("Gaussian Broadeninig is employed\n"); 
    *method=2;
  }
  else{
    if (Kgrid[0]==1 && Kgrid[1]==1 && Kgrid[2]==1 ) {
      printf("Kgrid= 1 1 1 => Gaussian Boardening is employed\n");
      *method=2;
    }
    else {
      printf("Which method do you use?, Tetrahedron(1), Gaussian Broadeninig(2)\n"); 
      fgets(buf,1000,stdin); sscanf(buf,"%d",method);
    }
  }

  if (*method==2) {
    printf("Please input a value of gaussian (double) (eV)\n"); 
    fgets(buf,1000,stdin); sscanf(buf,"%lf",gaussian); *gaussian = *gaussian/eV2Hartree; 
  }
    
  printf("Do you want Dos(1) or PDos(2)?\n");
  fgets(buf,1000,stdin);sscanf(buf,"%d",todo);
  
  if (*todo==2) {
    printf("\nNumber of atoms=%d\n",atomnum);
    if (atomnum==1) {
      (*pdos_atoms)=(int*)malloc(sizeof(int)*1);
      *pdos_n=1;
      (*pdos_atoms)[0]=1;
    }
    else {
      printf("Which atoms for PDOS : (1,...,%d), ex 1 2\n",atomnum);
      fgets(buf,1000,stdin);
      strcpy(buf2,buf);
      /*	printf("<%s>\n",buf); */
      *pdos_n=0;
      c=strtok(buf," ")  ;
      while ( c  ) {
	/*  printf("<%s> ",c); */
	if (sscanf(c,"%d",&i)==1)    (*pdos_n)++;
	c=strtok(NULL," ");
      }
      /*  printf("pdos_n=%d\n",*pdos_n); */
      *pdos_atoms=(int*)malloc(sizeof(int)*(*pdos_n));
      for (c=strtok(buf2," ") , i=0;i<*pdos_n;c=strtok(NULL," "),i++)  {
	/*  printf("<%s> ",c); */
	sscanf(c,"%d",&(*pdos_atoms)[i]);
      }
    }
    printf("pdos_n=%d\n",*pdos_n);
    for (i=0;i<*pdos_n;i++) {
      printf("%d ",(*pdos_atoms)[i]);
    }
    printf("\n");
  }

}

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

  char *file_eg, *file_ev;
 
  int mode;
  int n;
  double emax,emin;
  int iemax,iemin;
  int SpinP_switch;
  int Kgrid[3];
  int atomnum,SpeciesNum;
  double ChemP;
  int *WhatSpecies, *Spe_Total_CNO,MaxL;

  double *****ko;
  double *******ev; 
  int   **Spe_Num_CBasis; 
  int **Spe_Num_Relation; 

  double gaussian; 
  int Dos_N;
  int Max_tnoA;

  char basename[YOUSO10];

  /**************************************************
	     read data  from file_eg
  ***************************************************/

  file_eg=argv[1];
  file_ev=argv[2];

  input_file_eg(file_eg,
		&mode,  &n, &emin, &emax, &iemin,&iemax,
		&SpinP_switch,  Kgrid, &atomnum, &WhatSpecies,
		&SpeciesNum, &Spe_Total_CNO, &Max_tnoA, &MaxL,
		&Spe_Num_CBasis, &Spe_Num_Relation, &ChemP,
		&ko); 


  input_file_ev( file_ev, mode, Max_tnoA,
		 Kgrid,  SpinP_switch,  iemin,  iemax,
		 atomnum,  WhatSpecies,  Spe_Total_CNO,
		 &ev, ChemP);

 
    

  Delete_path_extension(file_eg,basename);

  {
    int todo,method;
    int pdos_n, *pdos_atoms;
 
    input_main(  mode, Kgrid,  atomnum,
		 &method, &todo,
		 &gaussian,
		 &pdos_n, &pdos_atoms); 

    /* set a suitable Dos_N */

    if (method==1){
      Dos_N = 300;
    }
    else{
      Dos_N = (emax-emin)/gaussian*5;
    }

    /* Total DOS */
    if (todo==1) {

      switch (method) {
      case 1:
	Dos_Tetrahedron(  basename,Dos_N,
			  Kgrid[0],Kgrid[1],Kgrid[2], 
			  SpinP_switch, ko, ev, iemax-iemin+1, 
			  emin, emax, 
			  iemin,  iemax, WhatSpecies, Spe_Total_CNO, atomnum );
	break;
      case 2: 

        /* DC */
        if (mode==5){

	DosDC_Gaussian(  basename, atomnum, Dos_N, gaussian,
		         SpinP_switch, emin, emax,
		         iemin, iemax, WhatSpecies, Spe_Total_CNO );
	}
	else{

	Dos_Gaussian(  basename,Dos_N,  gaussian,  
		       Kgrid[0],Kgrid[1],Kgrid[2],
		       SpinP_switch, ko, ev, iemax-iemin+1,
		       emin, emax,
		       iemin,  iemax, WhatSpecies, Spe_Total_CNO, atomnum );

	}

	break;
#if 0
      case 3:
	Dos_Histgram(  basename,Dos_N,
		       Kgrid[0],Kgrid[1],Kgrid[2],
		       SpinP_switch, ko, iemax-iemin+1,
		       emin, emax,
		       iemin,  iemax );
	break;
#endif
      }
    }

    /* Projected DOS */
    else if (todo==2) {

      switch(method) {
      case 1:
	Spectra_Tetrahedron( pdos_n, pdos_atoms, basename,Dos_N,
			     Kgrid[0],Kgrid[1],Kgrid[2],
			     SpinP_switch, ko,ev, iemax-iemin+1,
			     emin, emax,
			     iemin,  iemax,
			     atomnum, WhatSpecies, Spe_Total_CNO, 
			     MaxL,Spe_Num_CBasis, Spe_Num_Relation );
	break;
      case 2:

        /* DC */
        if (mode==5){

	SpectraDC_Gaussian( pdos_n, pdos_atoms, basename,Dos_N, gaussian, 
			    Kgrid[0],Kgrid[1],Kgrid[2],
			    SpinP_switch, ko,ev, iemax-iemin+1,
			    emin, emax,
			    iemin,  iemax,
			    atomnum, WhatSpecies, Spe_Total_CNO,
			    MaxL,Spe_Num_CBasis, Spe_Num_Relation );
	}
        else {

	Spectra_Gaussian( pdos_n, pdos_atoms, basename,Dos_N, gaussian, 
			  Kgrid[0],Kgrid[1],Kgrid[2],
			  SpinP_switch, ko,ev, iemax-iemin+1,
			  emin, emax,
			  iemin,  iemax,
			  atomnum, WhatSpecies, Spe_Total_CNO,
			  MaxL,Spe_Num_CBasis, Spe_Num_Relation );
        }

	break;
      }

    }

  }

  exit(0);
}