/* Usage: thisprogram Cgra.Dos.val Cgra.Dos.vec output: Cgra.DOS.Tetrahedron Cgra.PDOS.Tetrahedron */ #include #include #include #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(" 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(" 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=0 && iecenter 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(" 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) iemax=Dos_N-1; if ( 0<=iemin && iemin 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(" 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) iemax=Dos_N-1; if ( 0<=iemin && iemin 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(" 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) {iemax=Dos_N-1; } if ( 0<=iemin && iemin 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=Dos_N) iemax=Dos_N-1; for (GA_AN=0; GA_AN0) { 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 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=Dos_N) iemax=Dos_N-1; if ( 0<=iemin && iemin0) { 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 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(" make %s\n",file_Dos); #endif Max_tnoA=Spe_Total_CNO[0]; for (i=0;i=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 %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_AN0) { 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;ie0) { for (k=0;k")) { printf("format error near WhatSpecies>\n"); exit(10); } } else { printf("")) { printf("format error near Spe_Total_CNO>\n"); exit(10); } } else { printf("")) { printf("format error near Spe_Num_CBasis>\n"); exit(10); } } else { printf("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("")) { printf("format error near WhatSpecies>\n"); exit(10); } } else { printf("")) { printf("format error near Spe_Total_CNO>\n"); exit(10); } } else { printf("")) { printf("format error near Spe_Num_CBasis>\n"); exit(10); } } else { printf("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("")) { printf("format error near Eigenvalues>\n"); exit(10); } } /* if */ else { printf("can not find ",GA_AN+1,spin); if (!input_last(keyword)) { printf("format error near AN%dAN>\n",GA_AN+1); exit(10); } } else { printf(" 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); }