1 /*
2
3 Usage: thisprogram Cgra.Dos.val Cgra.Dos.vec
4
5 output: Cgra.DOS.Tetrahedron
6 Cgra.PDOS.Tetrahedron
7
8 */
9
10
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <string.h>
14 #include <math.h>
15 #include "Inputtools.h"
16
17 #define BINARYWRITE
18
19 #define eV2Hartree 27.2113845
20 #define PI 3.1415926535897932384626
21 #define fp_bsize 1048576 /* buffer size for setvbuf */
22
23 #define YOUSO10 100
24
25 typedef struct DCOMPLEX{double r,i;} dcomplex;
26
27 double ***EigenE;
28 double ****PDOS;
29 int **Msize1;
30
31 void *malloc_multidimarray(char *type, int N, int *size);
32 void OrderE0(double *e, int n);
33 void OrderE(double *e,double *a, int n);
34 void ATM_Dos(double *et, double *e, double *dos);
35 void ATM_Spectrum(double *et,double *at, double *e, double *spectrum);
36 char *Delete_path_extension(char *s0, char *r0);
37 void Dos_Gaussian( char *basename, int Dos_N, double Dos_Gaussian,
38 int knum_i,int knum_j, int knum_k,
39 int SpinP_switch, double *****EIGEN, double *******ev,
40 int neg,
41 double Dos_emin, double Dos_emax,
42 int iemin, int iemax,
43 int *WhatSpecies, int *Spe_Total_CNO, int atomnum );
44 void DosDC_Gaussian( char *basename, int atomnum,
45 int Dos_N, double Dos_Gaussian,
46 int SpinP_switch,
47 double Dos_emin, double Dos_emax,
48 int iemin, int iemax,
49 int *WhatSpecies, int *Spe_Total_CNO );
50 void Dos_NEGF( char *basename,
51 char *extension,
52 int Dos_N,
53 int SpinP_switch, double *****EIGEN, double *******ev,
54 int neg,
55 double Dos_emin, double Dos_emax,
56 int iemin, int iemax,
57 int *WhatSpecies, int *Spe_Total_CNO, int atomnum );
58 void Dos_Tetrahedron( char *basename, int Dos_N,
59 int knum_i,int knum_j, int knum_k,
60 int SpinP_switch, double *****EIGEN, double *******ev,
61 int neg,
62 double Dos_emin, double Dos_emax,
63 int iemin, int iemax,
64 int *WhatSpecies, int *Spe_Total_CNO, int atomnum );
65 void Spectra_Gaussian( int pdos_n, int *pdos_atoms, char *basename, int Dos_N,
66 double Dos_Gaussian,
67 int knum_i, int knum_j, int knum_k,
68 int SpinP_switch, double *****EIGEN, double *******ev, int neg,
69 double Dos_emin, double Dos_emax,
70 int iemin, int iemax,
71 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
72 int MaxL, int **Spe_Num_CBasis, int **Spe_Num_Relation);
73 void SpectraDC_Gaussian( int pdos_n, int *pdos_atoms, char *basename, int Dos_N,
74 double Dos_Gaussian,
75 int knum_i,int knum_j, int knum_k,
76 int SpinP_switch, double *****EIGEN, double *******ev, int neg,
77 double Dos_emin, double Dos_emax,
78 int iemin, int iemax,
79 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
80 int MaxL,int **Spe_Num_CBasis, int **Spe_Num_Relation);
81 void Spectra_Tetrahedron( int pdos_n, int *pdos_atoms, char *basename, int Dos_N,
82 int knum_i,int knum_j, int knum_k,
83 int SpinP_switch, double *****EIGEN, double *******ev, int neg,
84 double Dos_emin, double Dos_emax,
85 int iemin, int iemax,
86 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
87 int MaxL,int **Spe_Num_CBasis, int **Spe_Num_Relation);
88 void Spectra_NEGF( int pdos_n, int *pdos_atoms, char *basename, char *extension, int Dos_N,
89 int SpinP_switch, double *****EIGEN, double *******ev, int neg,
90 double Dos_emin, double Dos_emax,
91 int iemin, int iemax,
92 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
93 int MaxL, int **Spe_Num_CBasis, int **Spe_Num_Relation);
94 void Spe_Num_CBasis2Relation(
95 int SpeciesNum,int MaxL,int *Spe_Total_CNO,int **Spe_Num_CBasis,
96 int **Spe_Num_Relation);
97 void input_file_eg(char *file_eg,
98 int *mode, int *NonCol,
99 int *n, double *emin, double *emax, int *iemin,int *iemax,
100 int *SpinP_switch, int Kgrid[3], int *atomnum, int **WhatSpecies,
101 int *SpeciesNum, int **Spe_Total_CNO, int * Max_tnoA, int *MaxL,
102 int ***Spe_Num_CBasis, int ***Spe_Num_Relation, double *ChemP,
103 double **Angle0_Spin, double **Angle1_Spin,
104 double ******ko);
105 void input_file_ev( char *file_ev, int mode, int NonCol, int Max_tnoA,
106 int Kgrid[3], int SpinP_switch, int iemin, int iemax,
107 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
108 double ********ev, double ChemP);
109 void input_main( int mode, int Kgrid[3], int atomnum,
110 int *method, int *todo,
111 double *gaussian,
112 int *pdos_n, int **pdos_atoms);
113
114
115
main(int argc,char ** argv)116 int main(int argc, char **argv)
117 {
118
119 char *file_eg, *file_ev;
120
121 int mode,NonCol;
122 int n;
123 double emax,emin;
124 int iemax,iemin;
125 int SpinP_switch;
126 int Kgrid[3];
127 int atomnum,SpeciesNum;
128 double ChemP;
129 int *WhatSpecies, *Spe_Total_CNO,MaxL;
130
131 double *Angle0_Spin,*Angle1_Spin;
132 double *****ko;
133 double *******ev;
134 int **Spe_Num_CBasis;
135 int **Spe_Num_Relation;
136
137 double gaussian;
138 int Dos_N;
139 int Max_tnoA;
140
141 char basename[YOUSO10];
142
143 /**************************************************
144 read data from file_eg
145 ***************************************************/
146
147 file_eg=argv[1];
148 file_ev=argv[2];
149
150 input_file_eg(file_eg,
151 &mode, &NonCol, &n, &emin, &emax, &iemin,&iemax,
152 &SpinP_switch, Kgrid, &atomnum, &WhatSpecies,
153 &SpeciesNum, &Spe_Total_CNO, &Max_tnoA, &MaxL,
154 &Spe_Num_CBasis, &Spe_Num_Relation, &ChemP,
155 &Angle0_Spin, &Angle1_Spin,
156 &ko);
157
158
159 input_file_ev( file_ev, mode, NonCol, Max_tnoA,
160 Kgrid, SpinP_switch, iemin, iemax,
161 atomnum, WhatSpecies, Spe_Total_CNO,
162 &ev, ChemP);
163
164
165 Delete_path_extension(file_eg,basename);
166
167 {
168 int todo,method;
169 int pdos_n, *pdos_atoms;
170
171 input_main( mode, Kgrid, atomnum,
172 &method, &todo,
173 &gaussian,
174 &pdos_n, &pdos_atoms);
175
176 /* set a suitable Dos_N */
177
178 if (method==1){
179 Dos_N = 500;
180 }
181 else if (method==4){
182 Dos_N = iemax - iemin + 1;
183 }
184 else{
185 Dos_N = (emax-emin)/gaussian*5;
186 }
187
188
189 /* Total DOS */
190 if (todo==1) {
191
192 switch (method) {
193 case 1:
194 Dos_Tetrahedron( basename,Dos_N,
195 Kgrid[0],Kgrid[1],Kgrid[2],
196 SpinP_switch, ko, ev, iemax-iemin+1,
197 emin, emax,
198 iemin, iemax, WhatSpecies, Spe_Total_CNO, atomnum );
199 break;
200 case 2:
201
202 /* DC */
203 if (mode==5){
204
205 DosDC_Gaussian( basename, atomnum, Dos_N, gaussian,
206 SpinP_switch, emin, emax,
207 iemin, iemax, WhatSpecies, Spe_Total_CNO );
208 }
209 else{
210
211 Dos_Gaussian( basename, Dos_N, gaussian,
212 Kgrid[0],Kgrid[1],Kgrid[2],
213 SpinP_switch, ko, ev, iemax-iemin+1,
214 emin, emax,
215 iemin, iemax, WhatSpecies, Spe_Total_CNO, atomnum );
216
217 }
218
219 break;
220 #if 0
221 case 3:
222 Dos_Histgram( basename,Dos_N,
223 Kgrid[0],Kgrid[1],Kgrid[2],
224 SpinP_switch, ko, iemax-iemin+1,
225 emin, emax,
226 iemin, iemax );
227 break;
228 #endif
229
230 case 4:
231
232 if (mode==6){
233
234 Dos_NEGF( basename, "NEGF", Dos_N,
235 SpinP_switch, ko, ev, iemax-iemin+1,
236 emin, emax,
237 iemin, iemax, WhatSpecies, Spe_Total_CNO, atomnum );
238 }
239
240 else if (mode==7){
241
242 Dos_NEGF( basename, "Gaussian", Dos_N,
243 SpinP_switch, ko, ev, iemax-iemin+1,
244 emin, emax,
245 iemin, iemax, WhatSpecies, Spe_Total_CNO, atomnum );
246 }
247
248 break;
249 }
250 }
251
252 /* Projected DOS */
253
254 else if (todo==2) {
255
256 switch(method) {
257 case 1:
258 Spectra_Tetrahedron( pdos_n, pdos_atoms, basename,Dos_N,
259 Kgrid[0],Kgrid[1],Kgrid[2],
260 SpinP_switch, ko,ev, iemax-iemin+1,
261 emin, emax,
262 iemin, iemax,
263 atomnum, WhatSpecies, Spe_Total_CNO,
264 MaxL,Spe_Num_CBasis, Spe_Num_Relation );
265 break;
266 case 2:
267
268 /* DC */
269 if (mode==5){
270
271 SpectraDC_Gaussian( pdos_n, pdos_atoms, basename,Dos_N, gaussian,
272 Kgrid[0],Kgrid[1],Kgrid[2],
273 SpinP_switch, ko,ev, iemax-iemin+1,
274 emin, emax,
275 iemin, iemax,
276 atomnum, WhatSpecies, Spe_Total_CNO,
277 MaxL,Spe_Num_CBasis, Spe_Num_Relation );
278 }
279 else {
280
281 Spectra_Gaussian( pdos_n, pdos_atoms, basename,Dos_N, gaussian,
282 Kgrid[0],Kgrid[1],Kgrid[2],
283 SpinP_switch, ko,ev, iemax-iemin+1,
284 emin, emax,
285 iemin, iemax,
286 atomnum, WhatSpecies, Spe_Total_CNO,
287 MaxL,Spe_Num_CBasis, Spe_Num_Relation );
288 }
289
290 break;
291
292 case 4:
293
294 if (mode==6){
295
296 Spectra_NEGF( pdos_n, pdos_atoms, basename, "NEGF", Dos_N,
297 SpinP_switch, ko,ev, iemax-iemin+1,
298 emin, emax,
299 iemin, iemax,
300 atomnum, WhatSpecies, Spe_Total_CNO,
301 MaxL,Spe_Num_CBasis, Spe_Num_Relation );
302 }
303
304 else if (mode==7){
305
306 Spectra_NEGF( pdos_n, pdos_atoms, basename, "Gaussian", Dos_N,
307 SpinP_switch, ko,ev, iemax-iemin+1,
308 emin, emax,
309 iemin, iemax,
310 atomnum, WhatSpecies, Spe_Total_CNO,
311 MaxL,Spe_Num_CBasis, Spe_Num_Relation );
312 }
313
314 break;
315
316 }
317 }
318 }
319
320 exit(0);
321 }
322
323
324
325
326
327 /*
328 input s0
329 output r0
330 if input='../test.Dos.val' output='test'
331 */
332
Delete_path_extension(char * s0,char * r0)333 char *Delete_path_extension(char *s0, char *r0)
334 {
335 char *c;
336 /* char s[YOUSO10]; */
337
338 /* find '/' */
339 c=rindex(s0,'/');
340 if (c) {
341 strcpy(r0,c+1);
342 }
343 else {
344 strcpy(r0,s0);
345 }
346 printf("<%s>\n",r0);
347
348 if (strlen(r0)==0 ) { return NULL; }
349
350 c =index(r0,'.');
351 if (c) {
352 *c='\0';
353 }
354 printf("<%s>\n",r0);
355 return r0;
356
357 }
358
359
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)360 void Dos_Histgram( char *basename, int Dos_N,
361 int knum_i,int knum_j, int knum_k,
362 int SpinP_switch, double *****EIGEN, int neg,
363 double Dos_emin, double Dos_emax,
364 int iemin, int iemax )
365 {
366 int ie,spin,i,j,k,ieg,iecenter;
367 double *DosE, **Dos;
368 double eg,x,factor;
369 int N_Dos , i_Dos[10];
370 char file_Dos[YOUSO10];
371 FILE *fp_Dos;
372 char fp_buf[fp_bsize]; /* setvbuf */
373
374 printf("<Dos_Histgram> start\n");
375
376 if (strlen(basename)==0) {
377 strcpy(file_Dos,"DOS.Histgram");
378 }
379 else {
380 sprintf(file_Dos,"%s.DOS.Histgram",basename);
381 }
382 if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
383
384 #ifdef xt3
385 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
386 #endif
387
388 printf("can not open a file %s\n",file_Dos);
389 return;
390 }
391
392 printf("<Dos_Histgram> make %s\n",file_Dos);
393
394 /* allocation */
395 DosE=(double*)malloc(sizeof(double)*Dos_N);
396 N_Dos=2; i_Dos[0]=SpinP_switch+1; i_Dos[1]=Dos_N;
397 Dos=(double**)malloc_multidimarray("double",N_Dos,i_Dos);
398
399 /* initialize */
400 for (ie=0;ie<Dos_N;ie++) {
401 DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
402 }
403 for (spin=0;spin<=SpinP_switch;spin++) {
404 for (ie=0;ie<Dos_N;ie++) {
405 Dos[spin][ie]=0.0;
406 }
407 }
408
409 for (spin=0;spin<=SpinP_switch;spin++) {
410 for (i=0; i<=(knum_i-1); i++){
411 for (j=0; j<=(knum_j-1); j++){
412 for (k=0; k<=(knum_k-1); k++){
413 for (ieg=0;ieg<neg; ieg++) {
414
415 eg = EIGEN[spin][i][j][k][ieg] ;
416 x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
417 iecenter = (int)x ;
418
419 if (iecenter>=0 && iecenter <Dos_N )
420 Dos[spin][iecenter] = Dos[spin][iecenter] + 1.0;
421 } /* ieg */
422 } /* k */
423 } /* j */
424 } /* i */
425 } /* spin */
426
427 /* normalize */
428 factor = eV2Hartree * knum_i * knum_j * knum_k * (Dos_emax-Dos_emin)/(Dos_N-1);
429 factor = 1.0/factor;
430 for (spin=0;spin<=SpinP_switch;spin++) {
431 for (ie=0;ie<Dos_N;ie++) {
432 Dos[spin][ie] = Dos[spin][ie] * factor;
433 }
434 }
435
436 if (SpinP_switch==1) {
437 for (ie=0;ie<Dos_N;ie++) {
438 fprintf(fp_Dos,"%lf %lf %lf\n",(DosE[ie])*eV2Hartree, Dos[0][ie], -Dos[1][ie]);
439 }
440 }
441 else {
442 for (ie=0;ie<Dos_N;ie++) {
443 fprintf(fp_Dos,"%lf %lf\n",(DosE[ie])*eV2Hartree, Dos[0][ie]*2.0);
444 }
445 }
446
447 if (fp_Dos) fclose(fp_Dos);
448
449
450 free(DosE);
451
452 for (i=0; i<i_Dos[0]; i++){
453 free(Dos[i]);
454 }
455 free(Dos);
456
457 }
458
459
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)460 void Dos_Gaussian( char *basename, int Dos_N, double Dos_Gaussian,
461 int knum_i,int knum_j, int knum_k,
462 int SpinP_switch, double *****EIGEN, double *******ev,
463 int neg,
464 double Dos_emin, double Dos_emax,
465 int iemin, int iemax,
466 int *WhatSpecies, int *Spe_Total_CNO, int atomnum )
467 {
468
469 /*****************************************************************
470 Method = Gaussian
471 *****************************************************************/
472 /* exp( -(x/a)^2 )
473 a= Dos_Gaussian
474 x=-3a : 3a is counted */
475 double pi2;
476 int iewidth,ie;
477
478 int spin,i,j,k,ieg,iecenter;
479 int wanA,GA_AN,tnoA,i1;
480 double eg,x,xa,factor,wval;
481
482 double *DosE, **Dos;
483 int N_Dos,i_Dos[10];
484
485 char file_Dos[YOUSO10];
486 FILE *fp_Dos;
487 char fp_buf[fp_bsize]; /* setvbuf */
488
489 /* sawada */
490
491 int q, sw;
492 double h,s1,s2;
493 double ssum[2][Dos_N];
494
495 printf("<Dos_Gaussian> start\n");
496
497 if (strlen(basename)==0) {
498 strcpy(file_Dos,"DOS.Gaussian");
499 }
500 else {
501 sprintf(file_Dos,"%s.DOS.Gaussian",basename);
502 }
503 if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
504 printf("can not open a file %s\n",file_Dos);
505 return;
506 }
507
508 printf("<Dos_Gaussian> make %s\n",file_Dos);
509
510 /* allocation */
511 DosE=(double*)malloc(sizeof(double)*Dos_N);
512 N_Dos=2; i_Dos[0]=SpinP_switch+1; i_Dos[1]=Dos_N;
513 Dos=(double**)malloc_multidimarray("double",N_Dos,i_Dos);
514
515 /* initialize */
516 for (ie=0;ie<Dos_N;ie++) {
517 DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
518 }
519 for (spin=0;spin<=SpinP_switch;spin++) {
520 for (ie=0;ie<Dos_N;ie++) {
521 Dos[spin][ie]=0.0;
522 }
523 }
524
525 pi2 = sqrt(PI);
526 iewidth = Dos_Gaussian*3.0/(Dos_emax-Dos_emin)*(Dos_N-1)+3;
527 printf("Gaussian width=%d\n",iewidth);
528
529 for (spin=0;spin<=SpinP_switch;spin++) {
530 for (i=0; i<=(knum_i-1); i++){
531 for (j=0; j<=(knum_j-1); j++){
532 for (k=0; k<=(knum_k-1); k++){
533 for (ieg=0;ieg<neg; ieg++) {
534
535 /* generalization for spin non-collinear */
536 wval = 0.0;
537 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
538 wanA = WhatSpecies[GA_AN];
539 tnoA = Spe_Total_CNO[wanA];
540 for (i1=0; i1<tnoA; i1++){
541 wval += ev[spin][i][j][k][GA_AN][i1][ieg];
542 }
543 }
544
545 eg = EIGEN[spin][i][j][k][ieg] ;
546 x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
547 iecenter = (int)x ;
548
549 iemin = iecenter - iewidth;
550 iemax = iecenter + iewidth;
551
552 if (iemin<0) iemin=0;
553 if (iemax>=Dos_N) iemax=Dos_N-1;
554 if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax <Dos_N ) {
555
556 for (ie=iemin;ie<=iemax;ie++) {
557 xa = (eg-DosE[ie])/Dos_Gaussian;
558
559 /*
560 printf("Dos_N=%3d iemin=%3d iemax=%3d Dos_Gaussian=%15.12f %15.12f\n",
561 Dos_N,iemin,iemax,Dos_Gaussian,wval*exp( -xa*xa)/(Dos_Gaussian*pi2));
562 */
563
564 Dos[spin][ie] += wval*exp(- xa*xa)/(Dos_Gaussian*pi2);
565 }
566 }
567 } /* ieg */
568 } /* k */
569 } /* j */
570 } /* i */
571 } /* spin */
572
573 /* normalize */
574 factor = 1.0/(double)( eV2Hartree * knum_i * knum_j * knum_k );
575 for (spin=0;spin<=SpinP_switch;spin++) {
576 for (ie=0;ie<Dos_N;ie++) {
577 Dos[spin][ie] = Dos[spin][ie] * factor;
578 }
579
580 /* sawada */
581
582 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
583 sw = 2;
584
585 if (sw == 1) {
586
587 /* Trapezoidal rule */
588
589 for (q=0;q<Dos_N;q++) {
590 s1 = 0.0;
591 for (ie=1;ie<q;ie++) {
592 s1 += Dos[spin][ie];
593 }
594 ssum[spin][q] = (0.5*(Dos[spin][0]+Dos[spin][q])+s1)*h;
595 }
596 } else {
597
598 /* Simpson's rule */
599
600 for (q=0;q<Dos_N;q++) {
601 s1 = 0.0;
602 s2 = 0.0;
603 for (ie=1;ie<q;ie+=2) {
604 s1 += Dos[spin][ie];
605 }
606 for (ie=2;ie<q;ie+=2) {
607 s2 += Dos[spin][ie];
608 }
609 ssum[spin][q] = (Dos[spin][0]+4.0*s1+2.0*s2+Dos[spin][q])*h/3.0;
610 }
611 }
612
613 }
614
615 if (SpinP_switch==1) {
616 for (ie=0;ie<Dos_N;ie++) {
617 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie], -Dos[1][ie], ssum[0][ie], ssum[1][ie]);
618 }
619 }
620
621 else {
622 for (ie=0;ie<Dos_N;ie++) {
623 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie]*2.0, ssum[0][ie]*2.0);
624 }
625 }
626
627 if (fp_Dos) fclose(fp_Dos);
628
629 free(DosE);
630
631 for (i=0; i<i_Dos[0]; i++){
632 free(Dos[i]);
633 }
634 free(Dos);
635 }
636
637
638
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)639 void DosDC_Gaussian( char *basename, int atomnum,
640 int Dos_N, double Dos_Gaussian,
641 int SpinP_switch,
642 double Dos_emin, double Dos_emax,
643 int iemin, int iemax,
644 int *WhatSpecies, int *Spe_Total_CNO )
645 {
646
647 /*****************************************************************
648 Method = Gaussian
649 *****************************************************************/
650 /********************************
651 exp( -(x/a)^2 )
652 a= Dos_Gaussian
653 x=-3a : 3a is counted
654 ********************************/
655
656 double pi2;
657 int iewidth,ie;
658 int GA_AN,wanA,tnoA;
659
660 int spin,i,j,k,ieg,iecenter;
661 double eg,x,xa,factor,sum;
662
663 double *DosE, **Dos;
664 int N_Dos,i_Dos[10];
665 char file_Dos[YOUSO10];
666 FILE *fp_Dos;
667 char fp_buf[fp_bsize]; /* setvbuf */
668
669 /* sawada */
670
671 int q, sw;
672 double h,s1,s2;
673 double ssum[2][Dos_N];
674
675 printf("<Dos_Gaussian> start\n");
676
677 if (strlen(basename)==0) {
678 strcpy(file_Dos,"DOS.Gaussian");
679 }
680 else {
681 sprintf(file_Dos,"%s.DOS.Gaussian",basename);
682 }
683 if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
684
685 #ifdef xt3
686 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
687 #endif
688
689 printf("can not open a file %s\n",file_Dos);
690 return;
691 }
692
693 printf("<Dos_Gaussian> make %s\n",file_Dos);
694
695 /* allocation */
696 DosE=(double*)malloc(sizeof(double)*Dos_N);
697 N_Dos=2; i_Dos[0]=SpinP_switch+1; i_Dos[1]=Dos_N;
698 Dos=(double**)malloc_multidimarray("double",N_Dos,i_Dos);
699
700
701 /* initialize */
702 for (ie=0;ie<Dos_N;ie++) {
703 DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
704 }
705 for (spin=0;spin<=SpinP_switch;spin++) {
706 for (ie=0;ie<Dos_N;ie++) {
707 Dos[spin][ie]=0.0;
708 }
709 }
710
711 pi2=sqrt(PI);
712 iewidth = Dos_Gaussian*3.0/(Dos_emax-Dos_emin)*(Dos_N-1)+3;
713 printf("Gaussian width=%d\n",iewidth);
714
715 for (spin=0;spin<=SpinP_switch;spin++) {
716
717 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
718
719 wanA = WhatSpecies[GA_AN];
720 tnoA = Spe_Total_CNO[wanA];
721
722 for (ieg=0; ieg<Msize1[spin][GA_AN]; ieg++) {
723
724 eg = EigenE[spin][GA_AN][ieg];
725 x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
726 iecenter = (int)x ;
727 iemin = iecenter - iewidth;
728 iemax = iecenter + iewidth ;
729
730 if (iemin<0) iemin=0;
731 if (iemax>=Dos_N) iemax=Dos_N-1;
732
733 if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax <Dos_N ) {
734 sum = 0.0;
735 for (i=0; i<tnoA; i++) sum += PDOS[spin][GA_AN][ieg][i];
736 for (ie=iemin; ie<=iemax; ie++) {
737 xa = (eg-DosE[ie])/Dos_Gaussian;
738 Dos[spin][ie] += sum*exp(- xa*xa)/(Dos_Gaussian*pi2);
739 }
740 }
741
742 } /* ieg */
743 } /* GA_AN */
744 } /* spin */
745
746 /* normalize */
747 factor = 1.0/eV2Hartree;
748 for (spin=0;spin<=SpinP_switch;spin++) {
749 for (ie=0;ie<Dos_N;ie++) {
750 Dos[spin][ie] = Dos[spin][ie] * factor;
751 }
752
753 /* sawada */
754
755 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
756 sw = 2;
757
758 if (sw == 1) {
759
760 /* Trapezoidal rule */
761
762 for (q=0;q<Dos_N;q++) {
763 s1 = 0.0;
764 for (ie=1;ie<q;ie++) {
765 s1 += Dos[spin][ie];
766 }
767 ssum[spin][q] = (0.5*(Dos[spin][0]+Dos[spin][q])+s1)*h;
768 }
769 } else {
770
771 /* Simpson's rule */
772
773 for (q=0;q<Dos_N;q++) {
774 s1 = 0.0;
775 s2 = 0.0;
776 for (ie=1;ie<q;ie+=2) {
777 s1 += Dos[spin][ie];
778 }
779 for (ie=2;ie<q;ie+=2) {
780 s2 += Dos[spin][ie];
781 }
782 ssum[spin][q] = (Dos[spin][0]+4.0*s1+2.0*s2+Dos[spin][q])*h/3.0;
783 }
784 }
785
786 }
787
788 if (SpinP_switch==1) {
789 for (ie=0;ie<Dos_N;ie++) {
790 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie], -Dos[1][ie], ssum[0][ie], ssum[1][ie]);
791 }
792 }
793 else {
794 for (ie=0;ie<Dos_N;ie++) {
795 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie]*2.0, ssum[0][ie]*2.0);
796 }
797 }
798
799 if (fp_Dos) fclose(fp_Dos);
800
801 free(DosE);
802
803 for (i=0; i<i_Dos[0]; i++){
804 free(Dos[i]);
805 }
806 free(Dos);
807 }
808
809
810
811
812
813
Dos_NEGF(char * basename,char * extension,int Dos_N,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)814 void Dos_NEGF( char *basename,
815 char *extension,
816 int Dos_N,
817 int SpinP_switch, double *****EIGEN, double *******ev,
818 int neg,
819 double Dos_emin, double Dos_emax,
820 int iemin, int iemax,
821 int *WhatSpecies, int *Spe_Total_CNO, int atomnum )
822 {
823 /*****************************************************************
824 Method = NEGF
825 *****************************************************************/
826
827 double pi2;
828 int iewidth,ie;
829
830 int spin,i,j,k,ieg,iecenter;
831 int wanA,GA_AN,tnoA,i1;
832 double eg,x,xa,factor,wval;
833
834 double *DosE, **Dos;
835 int N_Dos,i_Dos[10];
836
837 char file_Dos[YOUSO10];
838 FILE *fp_Dos;
839 char fp_buf[fp_bsize]; /* setvbuf */
840
841 /* sawada */
842
843 int q,sw;
844 double h,s1,s2;
845 double ssum[2][Dos_N];
846
847 printf("<Dos_%s> start\n",extension);
848
849 if (strlen(basename)==0) {
850 sprintf(file_Dos,"DOS.%s",extension);
851 }
852 else {
853 sprintf(file_Dos,"%s.DOS.%s",basename,extension);
854 }
855 if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
856 printf("can not open a file %s\n",file_Dos);
857 return;
858 }
859
860 printf("<Dos_%s> make %s\n",extension,file_Dos);
861
862 /* allocation */
863 DosE=(double*)malloc(sizeof(double)*Dos_N);
864 N_Dos=2; i_Dos[0]=SpinP_switch+1; i_Dos[1]=Dos_N;
865 Dos=(double**)malloc_multidimarray("double",N_Dos,i_Dos);
866
867 /* note: Dos_N == neg */
868
869 /* initialize */
870
871 for (ie=0; ie<Dos_N; ie++) {
872 DosE[ie] = EIGEN[0][0][0][0][ie];
873 }
874
875 for (spin=0; spin<=SpinP_switch; spin++){
876 for (ie=0; ie<Dos_N; ie++){
877 Dos[spin][ie] = 0.0;
878 }
879 }
880
881 for (spin=0; spin<=SpinP_switch; spin++) {
882 for (ieg=0; ieg<neg; ieg++) {
883
884 i = 0;
885 j = 0;
886 k = 0;
887
888 wval = 0.0;
889
890 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
891 wanA = WhatSpecies[GA_AN];
892 tnoA = Spe_Total_CNO[wanA];
893 for (i1=0; i1<tnoA; i1++){
894 wval += ev[spin][i][j][k][GA_AN][i1][ieg];
895 }
896 }
897
898 Dos[spin][ieg] = wval/eV2Hartree;
899
900 } /* ieg */
901
902 /* sawada */
903
904 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
905 sw = 2;
906
907 if (sw == 1) {
908
909 /* Trapezoidal rule */
910
911 for (q=0;q<Dos_N;q++) {
912 s1 = 0.0;
913 for (ie=1;ie<q;ie++) {
914 s1 += Dos[spin][ie];
915 }
916 ssum[spin][q] = (0.5*(Dos[spin][0]+Dos[spin][q])+s1)*h;
917 }
918 } else {
919
920 /* Simpson's rule */
921
922 for (q=0;q<Dos_N;q++) {
923 s1 = 0.0;
924 s2 = 0.0;
925 for (ie=1;ie<q;ie+=2) {
926 s1 += Dos[spin][ie];
927 }
928 for (ie=2;ie<q;ie+=2) {
929 s2 += Dos[spin][ie];
930 }
931 ssum[spin][q] = (Dos[spin][0]+4.0*s1+2.0*s2+Dos[spin][q])*h/3.0;
932 }
933 }
934
935 } /* spin */
936
937 /* save */
938
939 if (SpinP_switch==1) {
940 for (ie=0; ie<Dos_N; ie++) {
941 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie], -Dos[1][ie], ssum[0][ie], ssum[1][ie]);
942 }
943 }
944 else {
945 for (ie=0; ie<Dos_N; ie++) {
946 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie]*2.0, ssum[0][ie]*2.0);
947 }
948 }
949
950 if (fp_Dos) fclose(fp_Dos);
951
952 /* free arrays */
953
954 free(DosE);
955
956 for (i=0; i<i_Dos[0]; i++){
957 free(Dos[i]);
958 }
959 free(Dos);
960 }
961
962
963
964
965
966
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)967 void Dos_Tetrahedron( char *basename, int Dos_N,
968 int knum_i,int knum_j, int knum_k,
969 int SpinP_switch, double *****EIGEN, double *******ev,
970 int neg,
971 double Dos_emin, double Dos_emax,
972 int iemin, int iemax,
973 int *WhatSpecies, int *Spe_Total_CNO, int atomnum )
974 {
975 int spin,ieg,i,j,k,ie,i1,tnoA,wanA,GA_AN;
976 int i_in,j_in,k_in,ic,itetra;
977 double x,result,factor,wval;
978 double cell_e[8], tetra_e[4];
979 static int tetra_id[6][4]= { {0,1,2,5}, {1,2,3,5}, {2,3,5,7},
980 {0,2,4,5}, {2,4,5,6}, {2,5,6,7} };
981
982 double *DosE, **Dos;
983 int N_Dos,i_Dos[10];
984 char file_Dos[YOUSO10];
985 FILE *fp_Dos;
986 char fp_buf[fp_bsize]; /* setvbuf */
987
988 /* sawada */
989
990 int q, sw;
991 double h,s1,s2;
992 double ssum[2][Dos_N];
993
994 printf("<Dos_Tetrahedron> start\n");
995
996 if (strlen(basename)==0) {
997 strcpy(file_Dos,"DOS.Tetrahedron");
998 }
999 else {
1000 sprintf(file_Dos,"%s.DOS.Tetrahedron",basename);
1001 }
1002 if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
1003
1004 #ifdef xt3
1005 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
1006 #endif
1007
1008 printf("can not open a file %s\n",file_Dos);
1009 return;
1010 }
1011
1012 printf("<Dos_Tetrahedron> make %s\n",file_Dos);
1013
1014 /* allocation */
1015 DosE=(double*)malloc(sizeof(double)*Dos_N);
1016 N_Dos=2; i_Dos[0]=SpinP_switch+1; i_Dos[1]=Dos_N;
1017 Dos=(double**)malloc_multidimarray("double",N_Dos,i_Dos);
1018
1019 /* initialize */
1020 for (ie=0;ie<Dos_N;ie++) {
1021 DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
1022 }
1023 for (spin=0;spin<=SpinP_switch;spin++) {
1024 for (ie=0;ie<Dos_N;ie++) {
1025 Dos[spin][ie]=0.0;
1026 }
1027 }
1028
1029 /********************************************************************
1030 Method = tetrahedron
1031 *******************************************************************/
1032
1033 for (spin=0;spin<=SpinP_switch;spin++) {
1034 for (ieg=0;ieg<neg; ieg++) {
1035 for (i=0; i<=(knum_i-1); i++){
1036 for (j=0; j<=(knum_j-1); j++){
1037 for (k=0; k<=(knum_k-1); k++){
1038
1039 /* generalization for spin non-collinear */
1040 wval = 0.0;
1041 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
1042 wanA = WhatSpecies[GA_AN];
1043 tnoA = Spe_Total_CNO[wanA];
1044 for (i1=0; i1<tnoA; i1++){
1045 wval += ev[spin][i][j][k][GA_AN][i1][ieg];
1046 }
1047 }
1048
1049 for (i_in=0;i_in<2;i_in++) {
1050 for (j_in=0;j_in<2;j_in++) {
1051 for (k_in=0;k_in<2;k_in++) {
1052 cell_e[i_in*4+j_in*2+k_in] =
1053 EIGEN[spin][ (i+i_in)%knum_i ][ (j+j_in)%knum_j ][ (k+k_in)%knum_k ][ieg] ;
1054 }
1055 }
1056 }
1057 for (itetra=0;itetra<6;itetra++) {
1058 for (ic=0;ic<4;ic++) {
1059 tetra_e[ic]=cell_e[ tetra_id[itetra][ic] ];
1060 }
1061 OrderE0(tetra_e,4);
1062 x = (tetra_e[0]-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1)-1.0;
1063 iemin=(int)x;
1064 x = (tetra_e[3]-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1)+1.0;
1065 iemax=(int)x;
1066 if (iemin<0) { iemin=0; }
1067 if (iemax>=Dos_N) {iemax=Dos_N-1; }
1068 if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax<Dos_N ) {
1069 for (ie=iemin;ie<=iemax;ie++) {
1070 ATM_Dos( tetra_e, &DosE[ie], &result);
1071 Dos[spin][ie] += wval*result;
1072 }
1073 }
1074
1075 } /* itetra */
1076 } /* k */
1077 } /* j */
1078 } /* i */
1079 } /* ieg */
1080 } /* spin */
1081
1082 /* normalize */
1083 factor = 1.0/(double)( eV2Hartree * knum_i * knum_j * knum_k * 6 );
1084 for (spin=0;spin<=SpinP_switch;spin++) {
1085 for (ie=0;ie<Dos_N;ie++) {
1086 Dos[spin][ie] = Dos[spin][ie] * factor;
1087 }
1088
1089 /* sawada */
1090
1091 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
1092 sw = 2;
1093
1094 if (sw == 1) {
1095
1096 /* Trapezoidal rule */
1097
1098 for (q=0;q<Dos_N;q++) {
1099 s1 = 0.0;
1100 for (ie=1;ie<q;ie++) {
1101 s1 += Dos[spin][ie];
1102 }
1103 ssum[spin][q] = (0.5*(Dos[spin][0]+Dos[spin][q])+s1)*h;
1104 }
1105 } else {
1106
1107 /* Simpson's rule */
1108
1109 for (q=0;q<Dos_N;q++) {
1110 s1 = 0.0;
1111 s2 = 0.0;
1112 for (ie=1;ie<q;ie+=2) {
1113 s1 += Dos[spin][ie];
1114 }
1115 for (ie=2;ie<q;ie+=2) {
1116 s2 += Dos[spin][ie];
1117 }
1118 ssum[spin][q] = (Dos[spin][0]+4.0*s1+2.0*s2+Dos[spin][q])*h/3.0;
1119 }
1120 }
1121
1122 }
1123
1124 if (SpinP_switch==1) {
1125 for (ie=0;ie<Dos_N;ie++) {
1126 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie], -Dos[1][ie], ssum[0][ie], ssum[1][ie]);
1127 }
1128 }
1129 else {
1130 for (ie=0;ie<Dos_N;ie++) {
1131 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, Dos[0][ie]*2.0, ssum[0][ie]*2.0);
1132 }
1133 }
1134
1135 if (fp_Dos) fclose(fp_Dos);
1136
1137 free(DosE);
1138
1139 for (i=0; i<i_Dos[0]; i++){
1140 free(Dos[i]);
1141 }
1142 free(Dos);
1143 }
1144
1145
1146
1147
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)1148 void Spectra_Gaussian( int pdos_n, int *pdos_atoms, char *basename, int Dos_N,
1149 double Dos_Gaussian,
1150 int knum_i, int knum_j, int knum_k,
1151 int SpinP_switch, double *****EIGEN, double *******ev, int neg,
1152 double Dos_emin, double Dos_emax,
1153 int iemin, int iemax,
1154 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
1155 int MaxL, int **Spe_Num_CBasis, int **Spe_Num_Relation)
1156 {
1157
1158 /*****************************************************************
1159 Method = Gaussian
1160 *****************************************************************/
1161 /* exp( -(x/a)^2 )
1162 a= Dos_Gaussian
1163 x=-3a : 3a is counted */
1164
1165 int Max_tnoA;
1166 int i,spin,j,k;
1167 int wanA, GA_AN,tnoA,i1,ieg;
1168 double eg,x,rval,xa;
1169 double *DosE, ****Dos;
1170 int N_Dos, i_Dos[10];
1171 double pi2,factor;
1172 int iewidth,ie,iecenter;
1173 int iatom,L,M,LM;
1174 static char *Lname[5]={"s","p","d","f","g"};
1175 double dossum[2],MulP[5],dE;
1176
1177 char file_Dos[YOUSO10];
1178 FILE *fp_Dos;
1179 char fp_buf[fp_bsize]; /* setvbuf */
1180 static char *extension="PDOS.Gaussian";
1181
1182 /* sawada */
1183
1184 int q, sw;
1185 double h,s1,s2;
1186 double ssum[2][Dos_N], Doss[2][Dos_N];
1187
1188 printf("<Spectra_Gaussian> start\n");
1189
1190 pi2=sqrt(PI);
1191 iewidth = Dos_Gaussian*3.0/(Dos_emax-Dos_emin)*(Dos_N-1)+3;
1192
1193
1194 Max_tnoA=0;
1195 for (i=0;i<atomnum;i++) {
1196 wanA=WhatSpecies[i];
1197 if (Max_tnoA<Spe_Total_CNO[wanA]) Max_tnoA=Spe_Total_CNO[wanA];
1198 }
1199 /* allocation */
1200 DosE=(double*)malloc(sizeof(double)*Dos_N);
1201
1202 N_Dos=4; i_Dos[0]=SpinP_switch+1;
1203 i_Dos[1] = atomnum; i_Dos[2]= Max_tnoA;
1204 i_Dos[3]=Dos_N;
1205 Dos=(double****)malloc_multidimarray("double",N_Dos,i_Dos);
1206
1207
1208 /* initialize */
1209 for (ie=0;ie<Dos_N;ie++) {
1210 DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
1211 }
1212 for (spin=0;spin<=SpinP_switch;spin++) {
1213 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
1214 wanA = WhatSpecies[GA_AN];
1215 tnoA = Spe_Total_CNO[wanA];
1216 for (i1=0; i1<=(tnoA-1); i1++){
1217
1218 for (ie=0;ie<Dos_N;ie++) {
1219 Dos[spin][GA_AN][i1][ie]=0.0;
1220 }
1221 }
1222 }
1223 }
1224
1225 /********************************************************************
1226 Method = gaussian
1227 *******************************************************************/
1228
1229 factor = 1.0/(double)(knum_i * knum_j * knum_k); /* normalization factor */
1230
1231 for (spin=0; spin<=SpinP_switch; spin++) {
1232 for (i=0; i<=(knum_i-1); i++){
1233 for (j=0; j<=(knum_j-1); j++){
1234 for (k=0; k<=(knum_k-1); k++){
1235 for (ieg=0;ieg<neg; ieg++) {
1236
1237 eg = EIGEN[spin][i][j][k][ieg];
1238
1239 x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
1240 iecenter = (int)x ;
1241 iemin = iecenter - iewidth;
1242 iemax = iecenter + iewidth ;
1243 if (iemin<0) iemin=0;
1244 if (iemax>=Dos_N) iemax=Dos_N-1;
1245
1246 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
1247 wanA = WhatSpecies[GA_AN];
1248 tnoA = Spe_Total_CNO[wanA];
1249 for (i1=0; i1<=(tnoA-1); i1++){
1250
1251 rval = ev[spin][i][j][k][GA_AN][i1][ieg];
1252
1253 if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax <Dos_N ) {
1254 for (ie=iemin; ie<=iemax; ie++) {
1255 xa = (eg-DosE[ie])/Dos_Gaussian;
1256
1257 Dos[spin][GA_AN][i1][ie] += factor*rval* exp(-xa*xa)/(Dos_Gaussian*pi2*eV2Hartree);
1258 }
1259 } /* if */
1260 } /* i1 */
1261 } /* GA_AN */
1262 } /* ieg */
1263 } /* k */
1264 } /* j */
1265 } /* i */
1266 } /* spin */
1267
1268 /* orbital projection */
1269 for (iatom=0;iatom<pdos_n;iatom++) {
1270 GA_AN = pdos_atoms[iatom]-1;
1271 wanA = WhatSpecies[GA_AN];
1272 tnoA = Spe_Total_CNO[wanA];
1273 for (L=0;L<=MaxL; L++) {
1274 if (Spe_Num_CBasis[wanA][L]>0) {
1275 for (M=0;M<2*L+1;M++) {
1276 LM=100*L+M+1;
1277 sprintf(file_Dos,"%s.%s.atom%d.%s%d",basename,extension,pdos_atoms[iatom],Lname[L],M+1);
1278 if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL ) {
1279
1280 #ifdef xt3
1281 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
1282 #endif
1283
1284 printf("can not open %s\n",file_Dos);
1285 }
1286 else {
1287 printf("make %s\n",file_Dos);
1288
1289 /* sawada */
1290
1291 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
1292 sw = 2;
1293
1294 for (spin=0;spin<=SpinP_switch;spin++) {
1295 for (ie=1;ie<Dos_N;ie++) {
1296 Doss[spin][ie] = 0.0;
1297 }
1298 }
1299
1300 for (spin=0;spin<=SpinP_switch;spin++) {
1301 if (sw == 1) {
1302 /* Trapezoidal rule */
1303 for (q=0;q<Dos_N;q++) {
1304 s1 = 0.0;
1305 for (ie=1;ie<q;ie++) {
1306 Doss[spin][ie] = 0.0;
1307 for (i1=0;i1< tnoA; i1++) {
1308 if (LM== Spe_Num_Relation[wanA][i1]) {
1309 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1310 }
1311 }
1312 s1 += Doss[spin][ie];
1313 }
1314 ssum[spin][q] = (0.5*(Doss[spin][0]+Doss[spin][q])+s1)*h;
1315 }
1316 } else {
1317 /* Simpson's rule */
1318 for (q=0;q<Dos_N;q++) {
1319 s1 = 0.0;
1320 s2 = 0.0;
1321 for (ie=1;ie<q;ie+=2) {
1322 Doss[spin][ie] = 0.0;
1323 for (i1=0;i1< tnoA; i1++) {
1324 if (LM== Spe_Num_Relation[wanA][i1]) {
1325 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1326 }
1327 }
1328 s1 += Doss[spin][ie];
1329 }
1330 for (ie=2;ie<q;ie+=2) {
1331 Doss[spin][ie] = 0.0;
1332 for (i1=0;i1< tnoA; i1++) {
1333 if (LM== Spe_Num_Relation[wanA][i1]) {
1334 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1335 }
1336 }
1337 s2 += Doss[spin][ie];
1338 }
1339 ssum[spin][q] = (Doss[spin][0]+4.0*s1+2.0*s2+Doss[spin][q])*h/3.0;
1340 }
1341 }
1342 }
1343
1344 for (ie=0;ie<Dos_N;ie++) {
1345 for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
1346 for (i1=0;i1< tnoA; i1++) {
1347 if (LM== Spe_Num_Relation[wanA][i1]) {
1348 for (spin=0;spin<=SpinP_switch;spin++) {
1349 dossum[spin]+= Dos[spin][GA_AN][i1][ie];
1350 }
1351 }
1352 }
1353
1354 if (SpinP_switch==1) {
1355 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree,
1356 dossum[0],-dossum[1], ssum[0][ie], ssum[1][ie]);
1357 }
1358 else {
1359 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0, ssum[0][ie]*2.0);
1360 }
1361
1362 }
1363 }
1364 if (fp_Dos) fclose(fp_Dos);
1365 } /* M */
1366 }
1367 } /* L */
1368 } /* iatom */
1369
1370 /* atom projection */
1371 for (iatom=0; iatom<pdos_n; iatom++) {
1372
1373 GA_AN = pdos_atoms[iatom]-1;
1374 wanA = WhatSpecies[GA_AN];
1375 tnoA = Spe_Total_CNO[wanA];
1376
1377 sprintf(file_Dos,"%s.%s.atom%d",basename,extension,pdos_atoms[iatom]);
1378
1379 if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL ) {
1380
1381 #ifdef xt3
1382 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
1383 #endif
1384
1385 printf("can not open %s\n",file_Dos);
1386 }
1387 else {
1388 printf("make %s\n",file_Dos);
1389
1390 /*
1391 for (spin=0; spin<=SpinP_switch; spin++){
1392 MulP[spin] = 0.0;
1393 }
1394 */
1395
1396 /* sawada */
1397
1398 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
1399 sw = 2;
1400
1401 for (spin=0;spin<=SpinP_switch;spin++) {
1402 for (ie=1;ie<Dos_N;ie++) {
1403 Doss[spin][ie] = 0.0;
1404 }
1405 }
1406
1407 for (spin=0;spin<=SpinP_switch;spin++) {
1408 if (sw == 1) {
1409 /* Trapezoidal rule */
1410 for (q=0;q<Dos_N;q++) {
1411 s1 = 0.0;
1412 for (ie=1;ie<q;ie++) {
1413 Doss[spin][ie] = 0.0;
1414 for (i1=0;i1< tnoA; i1++) {
1415 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1416 }
1417 s1 += Doss[spin][ie];
1418 }
1419 ssum[spin][q] = (0.5*(Doss[spin][0]+Doss[spin][q])+s1)*h;
1420 }
1421 } else {
1422 /* Simpson's rule */
1423 for (q=0;q<Dos_N;q++) {
1424 s1 = 0.0;
1425 s2 = 0.0;
1426 for (ie=1;ie<q;ie+=2) {
1427 Doss[spin][ie] = 0.0;
1428 for (i1=0;i1< tnoA; i1++) {
1429 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1430 }
1431 s1 += Doss[spin][ie];
1432 }
1433 for (ie=2;ie<q;ie+=2) {
1434 Doss[spin][ie] = 0.0;
1435 for (i1=0;i1< tnoA; i1++) {
1436 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1437 }
1438 s2 += Doss[spin][ie];
1439 }
1440 ssum[spin][q] = (Doss[spin][0]+4.0*s1+2.0*s2+Doss[spin][q])*h/3.0;
1441 }
1442 }
1443 }
1444
1445
1446 for (ie=0;ie<Dos_N;ie++) {
1447 for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
1448 for (i1=0;i1< tnoA; i1++) {
1449 for (spin=0;spin<=SpinP_switch;spin++)
1450 dossum[spin]+= Dos[spin][GA_AN][i1][ie];
1451 }
1452
1453 /* integration */
1454 /*
1455 if (DosE[ie]<=0.00){
1456 dE = DosE[1] - DosE[0];
1457 for (spin=0; spin<=SpinP_switch; spin++){
1458 MulP[spin] += dossum[spin]*dE;
1459 }
1460 }
1461 */
1462
1463 if (SpinP_switch==1) {
1464 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree,
1465 dossum[0],-dossum[1], ssum[0][ie], ssum[1][ie]);
1466 }
1467 else {
1468 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0, ssum[0][ie]*2.0);
1469 }
1470
1471 } /* ie */
1472 } /* fp_Dos */
1473 if (fp_Dos) fclose(fp_Dos);
1474
1475
1476 /*
1477 for (spin=0; spin<=SpinP_switch; spin++){
1478 printf("spin=%2d MulP=%15.12f\n",spin,MulP[spin]*eV2Hartree);
1479 }
1480 */
1481
1482
1483 } /* iatom */
1484
1485
1486 /*********************************************************
1487 close and free
1488 *********************************************************/
1489
1490 free(DosE);
1491
1492 for (i=0; i<i_Dos[0]; i++){
1493 for (j=0; j<i_Dos[1]; j++){
1494 for (k=0; k<i_Dos[2]; k++){
1495 free(Dos[i][j][k]);
1496 }
1497 free(Dos[i][j]);
1498 }
1499 free(Dos[i]);
1500 }
1501 free(Dos);
1502 }
1503
1504
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)1505 void SpectraDC_Gaussian( int pdos_n, int *pdos_atoms, char *basename, int Dos_N,
1506 double Dos_Gaussian,
1507 int knum_i,int knum_j, int knum_k,
1508 int SpinP_switch, double *****EIGEN, double *******ev, int neg,
1509 double Dos_emin, double Dos_emax,
1510 int iemin, int iemax,
1511 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
1512 int MaxL,int **Spe_Num_CBasis, int **Spe_Num_Relation)
1513 {
1514
1515 /*****************************************************************
1516 Method = Gaussian
1517 *****************************************************************/
1518 /* exp( -(x/a)^2 )
1519 a= Dos_Gaussian
1520 x=-3a : 3a is counted */
1521
1522
1523 int Max_tnoA;
1524 int i,spin,j,k;
1525 int wanA,GA_AN,tnoA,i1,ieg;
1526 double eg,x,rval,xa,tmp1;
1527 double *DosE, ****Dos;
1528 int N_Dos, i_Dos[10];
1529 double pi2;
1530 int iewidth,ie,iecenter;
1531 int iatom,L,M,LM;
1532 static char *Lname[5]={"s","p","d","f","g"};
1533 double dossum[2];
1534
1535 char file_Dos[YOUSO10];
1536 FILE *fp_Dos;
1537 char fp_buf[fp_bsize]; /* setvbuf */
1538 static char *extension="PDOS.Gaussian";
1539
1540 /* sawada */
1541
1542 int q, sw;
1543 double h,s1,s2;
1544 double ssum[2][Dos_N], Doss[2][Dos_N];
1545
1546 printf("<Spectra_Gaussian> start\n");
1547
1548 pi2=sqrt(PI);
1549 iewidth = Dos_Gaussian*3.0/(Dos_emax-Dos_emin)*(Dos_N-1)+3;
1550
1551 Max_tnoA=0;
1552 for (i=0;i<atomnum;i++) {
1553 wanA=WhatSpecies[i];
1554 if (Max_tnoA<Spe_Total_CNO[wanA]) Max_tnoA=Spe_Total_CNO[wanA];
1555 }
1556 /* allocation */
1557 DosE=(double*)malloc(sizeof(double)*Dos_N);
1558
1559 Dos = (double****)malloc(sizeof(double***)*(SpinP_switch+1));
1560 for (spin=0; spin<=SpinP_switch; spin++) {
1561 Dos[spin] =(double***)malloc(sizeof(double**)*pdos_n);
1562 for (iatom=0;iatom<pdos_n;iatom++) {
1563 GA_AN = pdos_atoms[iatom]-1;
1564 wanA = WhatSpecies[GA_AN];
1565 tnoA = Spe_Total_CNO[wanA];
1566 Dos[spin][iatom] =(double**)malloc(sizeof(double*)*tnoA);
1567 for (i=0; i<tnoA; i++){
1568 Dos[spin][iatom][i] =(double*)malloc(sizeof(double)*Dos_N);
1569 for (ie=0; ie<Dos_N; ie++) {
1570 Dos[spin][iatom][i][ie] =0.0;
1571 }
1572 }
1573 }
1574 }
1575
1576 /* initialize */
1577 for (ie=0;ie<Dos_N;ie++) {
1578 DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
1579 }
1580
1581 /********************************************************************
1582 Method = gaussian
1583 *******************************************************************/
1584
1585 for (spin=0;spin<=SpinP_switch;spin++) {
1586 for (iatom=0;iatom<pdos_n;iatom++) {
1587
1588 GA_AN = pdos_atoms[iatom]-1;
1589 wanA = WhatSpecies[GA_AN];
1590 tnoA = Spe_Total_CNO[wanA];
1591
1592 for (ieg=0; ieg<Msize1[spin][GA_AN]; ieg++) {
1593
1594 eg = EigenE[spin][GA_AN][ieg];
1595 x = (eg-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1);
1596 iecenter = (int)x ;
1597 iemin = iecenter - iewidth;
1598 iemax = iecenter + iewidth ;
1599 if (iemin<0) iemin=0;
1600 if (iemax>=Dos_N) iemax=Dos_N-1;
1601
1602 if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax <Dos_N ) {
1603 for (ie=iemin; ie<=iemax; ie++) {
1604 xa = (eg-DosE[ie])/Dos_Gaussian;
1605
1606 tmp1 = exp( -xa*xa)/(Dos_Gaussian*pi2*eV2Hartree);
1607
1608 for (i=0; i<tnoA; i++){
1609 Dos[spin][iatom][i][ie] += PDOS[spin][GA_AN][ieg][i]*tmp1;
1610 }
1611 }
1612 }
1613
1614 } /* ieg */
1615 } /* iatom */
1616 } /* spin */
1617
1618 /* orbital projection */
1619 for (iatom=0;iatom<pdos_n;iatom++) {
1620 GA_AN = pdos_atoms[iatom]-1;
1621 wanA = WhatSpecies[GA_AN];
1622 tnoA = Spe_Total_CNO[wanA];
1623 for (L=0;L<=MaxL; L++) {
1624 if (Spe_Num_CBasis[wanA][L]>0) {
1625 for (M=0;M<2*L+1;M++) {
1626 LM=100*L+M+1;
1627 sprintf(file_Dos,"%s.%s.atom%d.%s%d",basename,extension,pdos_atoms[iatom],
1628 Lname[L],M+1);
1629 if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL ) {
1630
1631 #ifdef xt3
1632 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
1633 #endif
1634
1635 printf("can not open %s\n",file_Dos);
1636 }
1637 else {
1638 printf("make %s\n",file_Dos);
1639
1640 /* sawada */
1641
1642 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
1643 sw = 2;
1644
1645 for (spin=0;spin<=SpinP_switch;spin++) {
1646 for (ie=1;ie<Dos_N;ie++) {
1647 Doss[spin][ie] = 0.0;
1648 }
1649 }
1650
1651 for (spin=0;spin<=SpinP_switch;spin++) {
1652 if (sw == 1) {
1653 /* Trapezoidal rule */
1654 for (q=0;q<Dos_N;q++) {
1655 s1 = 0.0;
1656 for (ie=1;ie<q;ie++) {
1657 Doss[spin][ie] = 0.0;
1658 for (i1=0;i1< tnoA; i1++) {
1659 if (LM== Spe_Num_Relation[wanA][i1]) {
1660 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1661 }
1662 }
1663 s1 += Doss[spin][ie];
1664 }
1665 ssum[spin][q] = (0.5*(Doss[spin][0]+Doss[spin][q])+s1)*h;
1666 }
1667 } else {
1668 /* Simpson's rule */
1669 for (q=0;q<Dos_N;q++) {
1670 s1 = 0.0;
1671 s2 = 0.0;
1672 for (ie=1;ie<q;ie+=2) {
1673 Doss[spin][ie] = 0.0;
1674 for (i1=0;i1< tnoA; i1++) {
1675 if (LM== Spe_Num_Relation[wanA][i1]) {
1676 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1677 }
1678 }
1679 s1 += Doss[spin][ie];
1680 }
1681 for (ie=2;ie<q;ie+=2) {
1682 Doss[spin][ie] = 0.0;
1683 for (i1=0;i1< tnoA; i1++) {
1684 if (LM== Spe_Num_Relation[wanA][i1]) {
1685 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1686 }
1687 }
1688 s2 += Doss[spin][ie];
1689 }
1690 ssum[spin][q] = (Doss[spin][0]+4.0*s1+2.0*s2+Doss[spin][q])*h/3.0;
1691 }
1692 }
1693 }
1694
1695 for (ie=0;ie<Dos_N;ie++) {
1696 for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
1697 for (i1=0;i1< tnoA; i1++) {
1698 if (LM== Spe_Num_Relation[wanA][i1]) {
1699 for (spin=0;spin<=SpinP_switch;spin++)
1700 dossum[spin]+= Dos[spin][iatom][i1][ie];
1701 }
1702 }
1703 if (SpinP_switch==1) {
1704 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree,
1705 dossum[0],-dossum[1], ssum[0][ie], ssum[1][ie]);
1706 }
1707 else {
1708 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0, ssum[0][ie]*2.0);
1709 }
1710
1711 }
1712 }
1713 if (fp_Dos) fclose(fp_Dos);
1714 } /* M */
1715 }
1716 } /* L */
1717 } /* iatom */
1718
1719 /* atom projection */
1720 for (iatom=0;iatom<pdos_n;iatom++) {
1721 GA_AN = pdos_atoms[iatom]-1;
1722 wanA = WhatSpecies[GA_AN];
1723 tnoA = Spe_Total_CNO[wanA];
1724
1725 sprintf(file_Dos,"%s.%s.atom%d",basename,extension,pdos_atoms[iatom]);
1726 if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL ) {
1727
1728 #ifdef xt3
1729 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
1730 #endif
1731
1732 printf("can not open %s\n",file_Dos);
1733 }
1734 else {
1735 printf("make %s\n",file_Dos);
1736
1737 /* sawada */
1738
1739 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
1740 sw = 2;
1741
1742 for (spin=0;spin<=SpinP_switch;spin++) {
1743 for (ie=1;ie<Dos_N;ie++) {
1744 Doss[spin][ie] = 0.0;
1745 }
1746 }
1747
1748 for (spin=0;spin<=SpinP_switch;spin++) {
1749 if (sw == 1) {
1750 /* Trapezoidal rule */
1751 for (q=0;q<Dos_N;q++) {
1752 s1 = 0.0;
1753 for (ie=1;ie<q;ie++) {
1754 Doss[spin][ie] = 0.0;
1755 for (i1=0;i1< tnoA; i1++) {
1756 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1757 }
1758 s1 += Doss[spin][ie];
1759 }
1760 ssum[spin][q] = (0.5*(Doss[spin][0]+Doss[spin][q])+s1)*h;
1761 }
1762 } else {
1763 /* Simpson's rule */
1764 for (q=0;q<Dos_N;q++) {
1765 s1 = 0.0;
1766 s2 = 0.0;
1767 for (ie=1;ie<q;ie+=2) {
1768 Doss[spin][ie] = 0.0;
1769 for (i1=0;i1< tnoA; i1++) {
1770 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1771 }
1772 s1 += Doss[spin][ie];
1773 }
1774 for (ie=2;ie<q;ie+=2) {
1775 Doss[spin][ie] = 0.0;
1776 for (i1=0;i1< tnoA; i1++) {
1777 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
1778 }
1779 s2 += Doss[spin][ie];
1780 }
1781 ssum[spin][q] = (Doss[spin][0]+4.0*s1+2.0*s2+Doss[spin][q])*h/3.0;
1782 }
1783 }
1784 }
1785
1786 for (ie=0;ie<Dos_N;ie++) {
1787 for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
1788 for (i1=0;i1< tnoA; i1++) {
1789 for (spin=0;spin<=SpinP_switch;spin++)
1790 dossum[spin]+= Dos[spin][iatom][i1][ie];
1791 }
1792 if (SpinP_switch==1) {
1793 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree,
1794 dossum[0],-dossum[1], ssum[0][ie], ssum[1][ie]);
1795 }
1796 else {
1797 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0, ssum[0][ie]*2.0);
1798 }
1799
1800 } /* ie */
1801 } /* fp_Dos */
1802 if (fp_Dos) fclose(fp_Dos);
1803 } /* iatom */
1804
1805
1806 /*********************************************************
1807 close and free
1808 *********************************************************/
1809
1810 free(DosE);
1811
1812 for (spin=0; spin<=SpinP_switch; spin++) {
1813 for (iatom=0;iatom<pdos_n;iatom++) {
1814 GA_AN = pdos_atoms[iatom]-1;
1815 wanA = WhatSpecies[GA_AN];
1816 tnoA = Spe_Total_CNO[wanA];
1817 for (i=0; i<tnoA; i++){
1818 free(Dos[spin][iatom][i]);
1819 }
1820 free(Dos[spin][iatom]);
1821 }
1822 free(Dos[spin]);
1823 }
1824 free(Dos);
1825
1826 }
1827
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)1828 void Spectra_Tetrahedron( int pdos_n, int *pdos_atoms, char *basename, int Dos_N,
1829 int knum_i,int knum_j, int knum_k,
1830 int SpinP_switch, double *****EIGEN, double *******ev, int neg,
1831 double Dos_emin, double Dos_emax,
1832 int iemin, int iemax,
1833 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
1834 int MaxL,int **Spe_Num_CBasis, int **Spe_Num_Relation)
1835 {
1836 int spin,ieg,i,j,k,ie;
1837 int i_in,j_in,k_in,ic,itetra;
1838 double x,result,factor;
1839 double cell_e[8], tetra_e[4],tetra_a[4];
1840 static int tetra_id[6][4]= { {0,1,2,5}, {1,2,3,5}, {2,3,5,7},
1841 {0,2,4,5}, {2,4,5,6}, {2,5,6,7} };
1842 int wanA,Max_tnoA,GA_AN,tnoA,i1,iatom;
1843 int L,M,LM;
1844 static char *Lname[5]={"s","p","d","f","g"};
1845 double dossum[2];
1846
1847 double *DosE, ****Dos, rval;
1848 int N_Dos,i_Dos[10];
1849 double ***cell_a;
1850 int N_cell_a, i_cell_a[10];
1851
1852 char file_Dos[YOUSO10];
1853 FILE *fp_Dos;
1854 char fp_buf[fp_bsize]; /* setvbuf */
1855
1856 /* sawada */
1857
1858 int q, sw;
1859 double h,s1,s2;
1860 double ssum[2][Dos_N], Doss[2][Dos_N];
1861
1862 static char *extension="PDOS.Tetrahedron";
1863
1864 printf("<Spectra_Tetrahedron> start\n");
1865
1866 #if 0
1867 if (strlen(basename)==0) {
1868 strcpy(file_Dos,extension);
1869 }
1870 else {
1871 sprintf(file_Dos,"%s.%s",basename,extension);
1872 }
1873
1874 if ( (fp_Dos=fopen(file_Dos,"w"))==NULL ) {
1875
1876 #ifdef xt3
1877 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
1878 #endif
1879
1880 printf("can not open a file %s\n",file_Dos);
1881 return;
1882 }
1883
1884 printf("<Spectra_Tetrahedron> make %s\n",file_Dos);
1885 #endif
1886
1887 Max_tnoA=Spe_Total_CNO[0];
1888 for (i=0;i<atomnum;i++) {
1889 wanA=WhatSpecies[i];
1890 if (Max_tnoA<Spe_Total_CNO[wanA]) Max_tnoA=Spe_Total_CNO[wanA];
1891 }
1892
1893
1894
1895 /* allocation */
1896 DosE=(double*)malloc(sizeof(double)*Dos_N);
1897
1898 N_Dos=4; i_Dos[0]=SpinP_switch+1;
1899 i_Dos[1] = atomnum; i_Dos[2]= Max_tnoA;
1900 i_Dos[3]=Dos_N;
1901 Dos=(double****)malloc_multidimarray("double",N_Dos,i_Dos);
1902
1903
1904 N_cell_a=3; i_cell_a[0]= atomnum; i_cell_a[1]=Max_tnoA; i_cell_a[2]= 8;
1905 cell_a = (double***) malloc_multidimarray("double",N_cell_a,i_cell_a);
1906
1907 /* initialize */
1908 for (ie=0;ie<Dos_N;ie++) {
1909 DosE[ie] = Dos_emin+(Dos_emax-Dos_emin)*(double)ie/(double)(Dos_N-1);
1910 }
1911 for (spin=0;spin<=SpinP_switch;spin++) {
1912 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
1913 wanA = WhatSpecies[GA_AN];
1914 tnoA = Spe_Total_CNO[wanA];
1915 for (i1=0; i1<=(tnoA-1); i1++){
1916
1917 for (ie=0;ie<Dos_N;ie++) {
1918 Dos[spin][GA_AN][i1][ie]=0.0;
1919 }
1920 }
1921 }
1922 }
1923
1924 /********************************************************************
1925 Method = tetrahedron
1926 *******************************************************************/
1927
1928 for (spin=0;spin<=SpinP_switch;spin++) {
1929 for (ieg=0;ieg<neg; ieg++) {
1930 for (i=0; i<=(knum_i-1); i++){
1931 for (j=0; j<=(knum_j-1); j++){
1932 for (k=0; k<=(knum_k-1); k++){
1933
1934 for (i_in=0;i_in<2;i_in++) {
1935 for (j_in=0;j_in<2;j_in++) {
1936 for (k_in=0;k_in<2;k_in++) {
1937 cell_e[i_in*4+j_in*2+k_in] =
1938 EIGEN[spin][ (i+i_in)%knum_i ][ (j+j_in)%knum_j ][ (k+k_in)%knum_k ][ieg] ;
1939 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
1940 wanA = WhatSpecies[GA_AN];
1941 tnoA = Spe_Total_CNO[wanA];
1942 for (i1=0; i1<=(tnoA-1); i1++){
1943 rval= ev[spin][ (i+i_in)%knum_i ][ (j+j_in)%knum_j ][ (k+k_in)%knum_k ][GA_AN][i1][ieg];
1944 cell_a[GA_AN][i1][i_in*4+j_in*2+k_in] = rval;
1945 }
1946 }
1947 }
1948 }
1949 }
1950
1951 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
1952 wanA = WhatSpecies[GA_AN];
1953 tnoA = Spe_Total_CNO[wanA];
1954 for (i1=0; i1<=(tnoA-1); i1++){
1955
1956 for (itetra=0;itetra<6;itetra++) {
1957 for (ic=0;ic<4;ic++) {
1958 tetra_e[ic]=cell_e[ tetra_id[itetra][ic] ];
1959 tetra_a[ic]=cell_a[ GA_AN ][i1][ tetra_id[itetra][ic] ];
1960 }
1961 OrderE(tetra_e,tetra_a,4);
1962
1963 for (ie=0;ie<4;ie++)
1964 #if 0
1965 printf("%lf(%lf) ",tetra_e[ie],tetra_a[ie]);
1966 printf("\n");
1967 #endif
1968
1969 x = (tetra_e[0]-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1)-1.0;
1970 iemin=(int)x;
1971 x = (tetra_e[3]-Dos_emin)/(Dos_emax-Dos_emin)*(Dos_N-1)+1.0;
1972 iemax=(int)x;
1973 if (iemin<0) { iemin=0; }
1974 if (iemax>=Dos_N) {iemax=Dos_N-1; }
1975 #if 0
1976 printf("%lf %lf %lf %lf\n", DosE[iemin],tetra_e[0],tetra_e[3], DosE[iemax] );
1977 #endif
1978 if ( 0<=iemin && iemin<Dos_N && 0<=iemax && iemax<Dos_N ) {
1979 for (ie=iemin;ie<=iemax;ie++) {
1980 ATM_Spectrum( tetra_e,tetra_a, &DosE[ie], &result);
1981 Dos[spin][GA_AN][i1][ie] += result;
1982 #if 0
1983 printf("%lf-> %lf\n",DosE[ie],result);
1984 #endif
1985 }
1986 }
1987 } /* itetra */
1988 }
1989 }
1990
1991 } /* k */
1992 } /* j */
1993 } /* i */
1994 } /* ieg */
1995 } /* spin */
1996
1997 /* normalize */
1998 factor = 1.0/(double)( eV2Hartree * knum_i * knum_j * knum_k * 6 );
1999 for (spin=0;spin<=SpinP_switch;spin++) {
2000 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
2001 wanA = WhatSpecies[GA_AN];
2002 tnoA = Spe_Total_CNO[wanA];
2003 for (i1=0; i1<=(tnoA-1); i1++){
2004
2005 for (ie=0;ie<Dos_N;ie++) {
2006 Dos[spin][GA_AN][i1][ie] = Dos[spin][GA_AN][i1][ie] * factor;
2007 }
2008 }
2009 }
2010 }
2011
2012 /***************************************************
2013 output
2014 ***************************************************/
2015
2016 #if 0
2017 if (SpinP_switch==1) {
2018 for (ie=0;ie<Dos_N;ie++) {
2019 for (sum=0.0,sum2=0.0,GA_AN=0; GA_AN<atomnum; GA_AN++){
2020 wanA = WhatSpecies[GA_AN];
2021 tnoA = Spe_Total_CNO[wanA];
2022 for (i1=0; i1<=(tnoA-1); i1++){
2023 sum += Dos[0][GA_AN][i1][ie];
2024 sum2 += Dos[1][GA_AN][i1][ie];
2025 }
2026 }
2027
2028 fprintf(fp_Dos,"%lf %lf %lf", (DosE[ie])*eV2Hartree,
2029 sum, -sum2 );
2030 }
2031 }
2032 else {
2033 for (ie=0;ie<Dos_N;ie++) {
2034 for (sum=0.0,sum2=0.0,GA_AN=0; GA_AN<atomnum; GA_AN++){
2035 wanA = WhatSpecies[GA_AN];
2036 tnoA = Spe_Total_CNO[wanA];
2037 for (i1=0; i1<=(tnoA-1); i1++){
2038 sum += Dos[0][GA_AN][i1][ie];
2039 }
2040 }
2041
2042 fprintf(fp_Dos,"%lf %lf\n", (DosE[ie])*eV2Hartree, sum*2.0);
2043 }
2044 }
2045 #else
2046
2047 /*
2048 for (iatom=0;iatom<pdos_n;iatom++) {
2049 GA_AN = pdos_atoms[iatom]-1;
2050 wanA = WhatSpecies[GA_AN];
2051 tnoA = Spe_Total_CNO[wanA];
2052 for (i1=0;i1< tnoA; i1++) {
2053 printf("Spe_Num_Relation %d %d %d\n",wanA,i1,Spe_Num_Relation[wanA][i1]);
2054 }
2055 }
2056 */
2057
2058 /* orbital projection */
2059 for (iatom=0;iatom<pdos_n;iatom++) {
2060 GA_AN = pdos_atoms[iatom]-1;
2061 wanA = WhatSpecies[GA_AN];
2062 tnoA = Spe_Total_CNO[wanA];
2063 for (L=0;L<=MaxL; L++) {
2064 if (Spe_Num_CBasis[wanA][L]>0) {
2065 for (M=0;M<2*L+1;M++) {
2066 LM=100*L+M+1;
2067 sprintf(file_Dos,"%s.%s.atom%d.%s%d",basename,extension,pdos_atoms[iatom],
2068 Lname[L],M+1);
2069 if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL ) {
2070
2071 #ifdef xt3
2072 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
2073 #endif
2074
2075 printf("can not open %s\n",file_Dos);
2076 }
2077 else {
2078 printf("make %s\n",file_Dos);
2079
2080 /* sawada */
2081
2082 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
2083 sw = 2;
2084
2085 for (spin=0;spin<=SpinP_switch;spin++) {
2086 for (ie=1;ie<Dos_N;ie++) {
2087 Doss[spin][ie] = 0.0;
2088 }
2089 }
2090
2091 for (spin=0;spin<=SpinP_switch;spin++) {
2092 if (sw == 1) {
2093 /* Trapezoidal rule */
2094 for (q=0;q<Dos_N;q++) {
2095 s1 = 0.0;
2096 for (ie=1;ie<q;ie++) {
2097 Doss[spin][ie] = 0.0;
2098 for (i1=0;i1< tnoA; i1++) {
2099 if (LM== Spe_Num_Relation[wanA][i1]) {
2100 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2101 }
2102 }
2103 s1 += Doss[spin][ie];
2104 }
2105 ssum[spin][q] = (0.5*(Doss[spin][0]+Doss[spin][q])+s1)*h;
2106 }
2107 } else {
2108 /* Simpson's rule */
2109 for (q=0;q<Dos_N;q++) {
2110 s1 = 0.0;
2111 s2 = 0.0;
2112 for (ie=1;ie<q;ie+=2) {
2113 Doss[spin][ie] = 0.0;
2114 for (i1=0;i1< tnoA; i1++) {
2115 if (LM== Spe_Num_Relation[wanA][i1]) {
2116 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2117 }
2118 }
2119 s1 += Doss[spin][ie];
2120 }
2121 for (ie=2;ie<q;ie+=2) {
2122 Doss[spin][ie] = 0.0;
2123 for (i1=0;i1< tnoA; i1++) {
2124 if (LM== Spe_Num_Relation[wanA][i1]) {
2125 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2126 }
2127 }
2128 s2 += Doss[spin][ie];
2129 }
2130 ssum[spin][q] = (Doss[spin][0]+4.0*s1+2.0*s2+Doss[spin][q])*h/3.0;
2131 }
2132 }
2133 }
2134
2135 for (ie=0;ie<Dos_N;ie++) {
2136 for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
2137 for (i1=0;i1< tnoA; i1++) {
2138 /* if (i1==4) {
2139 printf("%d %d\n",LM,Spe_Num_Relation[wanA][i1]);
2140 }
2141 */
2142 if (LM== Spe_Num_Relation[wanA][i1]) {
2143 for (spin=0;spin<=SpinP_switch;spin++)
2144 dossum[spin]+= Dos[spin][GA_AN][i1][ie];
2145 }
2146 }
2147 if (SpinP_switch==1) {
2148 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n",
2149 (DosE[ie])*eV2Hartree,
2150 dossum[0],-dossum[1],ssum[0][ie],ssum[1][ie]);
2151 }
2152 else {
2153 fprintf(fp_Dos,"%lf %lf %lf\n",
2154 (DosE[ie])*eV2Hartree,
2155 dossum[0]*2.0,ssum[0][ie]*2.0);
2156 }
2157
2158 }
2159 }
2160 } /* M */
2161 }
2162 } /* L */
2163 } /* iatom */
2164
2165
2166 /* atom projection */
2167 for (iatom=0;iatom<pdos_n;iatom++) {
2168 GA_AN = pdos_atoms[iatom]-1;
2169 wanA = WhatSpecies[GA_AN];
2170 tnoA = Spe_Total_CNO[wanA];
2171
2172 sprintf(file_Dos,"%s.%s.atom%d",basename,extension,pdos_atoms[iatom]);
2173 if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL ) {
2174
2175 #ifdef xt3
2176 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
2177 #endif
2178
2179 printf("can not open %s\n",file_Dos);
2180 }
2181 else {
2182 printf("make %s\n",file_Dos);
2183
2184 /* sawada */
2185
2186 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
2187 sw = 2;
2188
2189 for (spin=0;spin<=SpinP_switch;spin++) {
2190 for (ie=1;ie<Dos_N;ie++) {
2191 Doss[spin][ie] = 0.0;
2192 }
2193 }
2194
2195 for (spin=0;spin<=SpinP_switch;spin++) {
2196 if (sw == 1) {
2197 /* Trapezoidal rule */
2198 for (q=0;q<Dos_N;q++) {
2199 s1 = 0.0;
2200 for (ie=1;ie<q;ie++) {
2201 Doss[spin][ie] = 0.0;
2202 for (i1=0;i1< tnoA; i1++) {
2203 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2204 }
2205 s1 += Doss[spin][ie];
2206 }
2207 ssum[spin][q] = (0.5*(Doss[spin][0]+Doss[spin][q])+s1)*h;
2208 }
2209 } else {
2210 /* Simpson's rule */
2211 for (q=0;q<Dos_N;q++) {
2212 s1 = 0.0;
2213 s2 = 0.0;
2214 for (ie=1;ie<q;ie+=2) {
2215 Doss[spin][ie] = 0.0;
2216 for (i1=0;i1< tnoA; i1++) {
2217 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2218 }
2219 s1 += Doss[spin][ie];
2220 }
2221 for (ie=2;ie<q;ie+=2) {
2222 Doss[spin][ie] = 0.0;
2223 for (i1=0;i1< tnoA; i1++) {
2224 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2225 }
2226 s2 += Doss[spin][ie];
2227 }
2228 ssum[spin][q] = (Doss[spin][0]+4.0*s1+2.0*s2+Doss[spin][q])*h/3.0;
2229 }
2230 }
2231 }
2232
2233 for (ie=0;ie<Dos_N;ie++) {
2234 for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
2235 for (i1=0;i1< tnoA; i1++) {
2236 for (spin=0;spin<=SpinP_switch;spin++)
2237 dossum[spin]+= Dos[spin][GA_AN][i1][ie];
2238 }
2239 if (SpinP_switch==1) {
2240 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree,
2241 dossum[0],-dossum[1],ssum[0][ie],ssum[1][ie]);
2242 }
2243 else {
2244 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree,
2245 dossum[0]*2.0, ssum[0][ie]*2.0);
2246 }
2247
2248 } /* ie */
2249 } /* fp_Dos */
2250 } /* iatom */
2251
2252
2253
2254 #endif
2255
2256
2257 /*********************************************************
2258 close and free
2259 *********************************************************/
2260
2261 #if 0
2262 if (fp_Dos) fclose(fp_Dos);
2263 #endif
2264
2265 free(DosE);
2266
2267 for (spin=0; spin<i_Dos[0]; spin++) {
2268 for (iatom=0; iatom<i_Dos[1]; iatom++) {
2269 for (i=0; i<i_Dos[2]; i++){
2270 free(Dos[spin][iatom][i]);
2271 }
2272 free(Dos[spin][iatom]);
2273 }
2274 free(Dos[spin]);
2275 }
2276 free(Dos);
2277
2278 for (i=0; i<i_cell_a[0]; i++){
2279 for (j=0; j<i_cell_a[1]; j++){
2280 free(cell_a[i][j]);
2281 }
2282 free(cell_a[i]);
2283 }
2284 free(cell_a);
2285
2286 }
2287
2288
2289
2290
2291
Spectra_NEGF(int pdos_n,int * pdos_atoms,char * basename,char * extension,int Dos_N,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)2292 void Spectra_NEGF( int pdos_n, int *pdos_atoms, char *basename, char *extension, int Dos_N,
2293 int SpinP_switch, double *****EIGEN, double *******ev, int neg,
2294 double Dos_emin, double Dos_emax,
2295 int iemin, int iemax,
2296 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
2297 int MaxL, int **Spe_Num_CBasis, int **Spe_Num_Relation)
2298 {
2299 /*****************************************************************
2300 Method = NEGF
2301 *****************************************************************/
2302
2303 int Max_tnoA;
2304 int i,spin,j,k;
2305 int wanA, GA_AN,tnoA,i1,ieg;
2306 double eg,x,rval,xa;
2307 double *DosE, ****Dos;
2308 int N_Dos, i_Dos[10];
2309 double pi2,factor;
2310 int iewidth,ie,iecenter;
2311 int iatom,L,M,LM;
2312 static char *Lname[5]={"s","p","d","f","g"};
2313 double dossum[2],MulP[5],dE;
2314
2315 char file_Dos[YOUSO10];
2316 FILE *fp_Dos;
2317 char fp_buf[fp_bsize]; /* setvbuf */
2318
2319 /* sawada */
2320
2321 int q, sw;
2322 double h,s1,s2;
2323 double ssum[2][Dos_N], Doss[2][Dos_N];
2324
2325 printf("<Spectra_%s> start\n",extension);
2326
2327 Max_tnoA=0;
2328 for (i=0; i<atomnum; i++) {
2329 wanA=WhatSpecies[i];
2330 if (Max_tnoA<Spe_Total_CNO[wanA]) Max_tnoA=Spe_Total_CNO[wanA];
2331 }
2332 /* allocation */
2333 DosE=(double*)malloc(sizeof(double)*Dos_N);
2334
2335 N_Dos=4; i_Dos[0]=SpinP_switch+1;
2336 i_Dos[1] = atomnum; i_Dos[2]= Max_tnoA;
2337 i_Dos[3]=Dos_N;
2338 Dos=(double****)malloc_multidimarray("double",N_Dos,i_Dos);
2339
2340 /* note: Dos_N == neg */
2341
2342 /* initialize */
2343
2344 for (ie=0; ie<Dos_N; ie++) {
2345 DosE[ie] = EIGEN[0][0][0][0][ie];
2346 }
2347
2348 /********************************************************************
2349 Method = NEGF
2350 *******************************************************************/
2351
2352 for (spin=0;spin<=SpinP_switch;spin++) {
2353 for (ieg=0; ieg<neg; ieg++) {
2354 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
2355 wanA = WhatSpecies[GA_AN];
2356 tnoA = Spe_Total_CNO[wanA];
2357 for (i1=0; i1<tnoA; i1++){
2358 rval = ev[spin][0][0][0][GA_AN][i1][ieg];
2359 Dos[spin][GA_AN][i1][ieg] = rval/eV2Hartree;
2360 } /* i1 */
2361 } /* GA_AN */
2362 } /* ieg */
2363 } /* spin */
2364
2365 /* orbital projection */
2366
2367 for (iatom=0; iatom<pdos_n; iatom++) {
2368
2369 GA_AN = pdos_atoms[iatom]-1;
2370 wanA = WhatSpecies[GA_AN];
2371 tnoA = Spe_Total_CNO[wanA];
2372
2373 for (L=0;L<=MaxL; L++) {
2374 if (Spe_Num_CBasis[wanA][L]>0) {
2375 for (M=0;M<2*L+1;M++) {
2376 LM=100*L+M+1;
2377
2378 sprintf(file_Dos,"%s.PDOS.%s.atom%d.%s%d",
2379 basename,extension,pdos_atoms[iatom],
2380 Lname[L],M+1);
2381
2382 if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL ) {
2383
2384 #ifdef xt3
2385 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
2386 #endif
2387 printf("can not open %s\n",file_Dos);
2388 }
2389
2390 else {
2391 printf("make %s\n",file_Dos);
2392
2393 /* sawada */
2394
2395 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
2396 sw = 2;
2397
2398 for (spin=0;spin<=SpinP_switch;spin++) {
2399 for (ie=1;ie<Dos_N;ie++) {
2400 Doss[spin][ie] = 0.0;
2401 }
2402 }
2403
2404 for (spin=0;spin<=SpinP_switch;spin++) {
2405 if (sw == 1) {
2406 /* Trapezoidal rule */
2407 for (q=0;q<Dos_N;q++) {
2408 s1 = 0.0;
2409 for (ie=1;ie<q;ie++) {
2410 Doss[spin][ie] = 0.0;
2411 for (i1=0;i1< tnoA; i1++) {
2412 if (LM== Spe_Num_Relation[wanA][i1]) {
2413 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2414 }
2415 }
2416 s1 += Doss[spin][ie];
2417 }
2418 ssum[spin][q] = (0.5*(Doss[spin][0]+Doss[spin][q])+s1)*h;
2419 }
2420 } else {
2421 /* Simpson's rule */
2422 for (q=0;q<Dos_N;q++) {
2423 s1 = 0.0;
2424 s2 = 0.0;
2425 for (ie=1;ie<q;ie+=2) {
2426 Doss[spin][ie] = 0.0;
2427 for (i1=0;i1< tnoA; i1++) {
2428 if (LM== Spe_Num_Relation[wanA][i1]) {
2429 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2430 }
2431 }
2432 s1 += Doss[spin][ie];
2433 }
2434 for (ie=2;ie<q;ie+=2) {
2435 Doss[spin][ie] = 0.0;
2436 for (i1=0;i1< tnoA; i1++) {
2437 if (LM== Spe_Num_Relation[wanA][i1]) {
2438 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2439 }
2440 }
2441 s2 += Doss[spin][ie];
2442 }
2443 ssum[spin][q] = (Doss[spin][0]+4.0*s1+2.0*s2+Doss[spin][q])*h/3.0;
2444 }
2445 }
2446 }
2447
2448 for (ie=0;ie<Dos_N;ie++) {
2449 for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
2450 for (i1=0;i1< tnoA; i1++) {
2451 if (LM== Spe_Num_Relation[wanA][i1]) {
2452 for (spin=0;spin<=SpinP_switch;spin++) {
2453 dossum[spin]+= Dos[spin][GA_AN][i1][ie];
2454 }
2455 }
2456 }
2457 if (SpinP_switch==1) {
2458 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree,
2459 dossum[0],-dossum[1],ssum[0][ie],ssum[1][ie]);
2460 }
2461 else {
2462 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0, ssum[0][ie]*2.0);
2463 }
2464
2465 }
2466 }
2467 if (fp_Dos) fclose(fp_Dos);
2468 } /* M */
2469 }
2470 } /* L */
2471 } /* iatom */
2472
2473 /* atom projection */
2474
2475 for (iatom=0;iatom<pdos_n;iatom++) {
2476 GA_AN = pdos_atoms[iatom]-1;
2477 wanA = WhatSpecies[GA_AN];
2478 tnoA = Spe_Total_CNO[wanA];
2479
2480 sprintf(file_Dos,"%s.PDOS.%s.atom%d",basename,extension,pdos_atoms[iatom]);
2481 if ( (fp_Dos=fopen(file_Dos,"w")) ==NULL ) {
2482
2483 #ifdef xt3
2484 setvbuf(fp_Dos,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
2485 #endif
2486
2487 printf("can not open %s\n",file_Dos);
2488 }
2489 else {
2490
2491 printf("make %s\n",file_Dos);
2492
2493 /*
2494 for (spin=0; spin<=SpinP_switch; spin++){
2495 MulP[spin] = 0.0;
2496 }
2497 */
2498
2499 /* sawada */
2500
2501 h = (DosE[Dos_N-1] - DosE[0])/(Dos_N-1)*eV2Hartree;
2502 sw = 2;
2503
2504 for (spin=0;spin<=SpinP_switch;spin++) {
2505 for (ie=1;ie<Dos_N;ie++) {
2506 Doss[spin][ie] = 0.0;
2507 }
2508 }
2509
2510 for (spin=0;spin<=SpinP_switch;spin++) {
2511 if (sw == 1) {
2512 /* Trapezoidal rule */
2513 for (q=0;q<Dos_N;q++) {
2514 s1 = 0.0;
2515 for (ie=1;ie<q;ie++) {
2516 Doss[spin][ie] = 0.0;
2517 for (i1=0;i1< tnoA; i1++) {
2518 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2519 }
2520 s1 += Doss[spin][ie];
2521 }
2522 ssum[spin][q] = (0.5*(Doss[spin][0]+Doss[spin][q])+s1)*h;
2523 }
2524 } else {
2525 /* Simpson's rule */
2526 for (q=0;q<Dos_N;q++) {
2527 s1 = 0.0;
2528 s2 = 0.0;
2529 for (ie=1;ie<q;ie+=2) {
2530 Doss[spin][ie] = 0.0;
2531 for (i1=0;i1< tnoA; i1++) {
2532 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2533 }
2534 s1 += Doss[spin][ie];
2535 }
2536 for (ie=2;ie<q;ie+=2) {
2537 Doss[spin][ie] = 0.0;
2538 for (i1=0;i1< tnoA; i1++) {
2539 Doss[spin][ie]+= Dos[spin][GA_AN][i1][ie];
2540 }
2541 s2 += Doss[spin][ie];
2542 }
2543 ssum[spin][q] = (Doss[spin][0]+4.0*s1+2.0*s2+Doss[spin][q])*h/3.0;
2544 }
2545 }
2546 }
2547
2548 for (ie=0;ie<Dos_N;ie++) {
2549 for (spin=0;spin<=SpinP_switch;spin++) dossum[spin]=0.0;
2550 for (i1=0;i1< tnoA; i1++) {
2551 for (spin=0;spin<=SpinP_switch;spin++)
2552 dossum[spin]+= Dos[spin][GA_AN][i1][ie];
2553 }
2554
2555 /* integration */
2556 /*
2557 if (DosE[ie]<=0.00){
2558 dE = DosE[1] - DosE[0];
2559 for (spin=0; spin<=SpinP_switch; spin++){
2560 MulP[spin] += dossum[spin]*dE;
2561 }
2562 }
2563 */
2564
2565 if (SpinP_switch==1) {
2566 fprintf(fp_Dos,"%lf %lf %lf %lf %lf\n", (DosE[ie])*eV2Hartree,
2567 dossum[0],-dossum[1],ssum[0][ie],ssum[1][ie]);
2568 }
2569 else {
2570 fprintf(fp_Dos,"%lf %lf %lf\n", (DosE[ie])*eV2Hartree, dossum[0]*2.0, ssum[0][ie]*2.0);
2571 }
2572
2573 } /* ie */
2574 } /* fp_Dos */
2575 if (fp_Dos) fclose(fp_Dos);
2576
2577 /*
2578 for (spin=0; spin<=SpinP_switch; spin++){
2579 printf("spin=%2d MulP=%15.12f\n",spin,MulP[spin]*eV2Hartree);
2580 }
2581 */
2582
2583 } /* iatom */
2584
2585
2586 /*********************************************************
2587 close and free
2588 *********************************************************/
2589
2590 free(DosE);
2591
2592 for (i=0; i<i_Dos[0]; i++){
2593 for (j=0; j<i_Dos[1]; j++){
2594 for (k=0; k<i_Dos[2]; k++){
2595 free(Dos[i][j][k]);
2596 }
2597 free(Dos[i][j]);
2598 }
2599 free(Dos[i]);
2600 }
2601 free(Dos);
2602 }
2603
2604
2605
2606
2607
2608 /* input int SpeciesNum,int MaxL,int *Spe_Total_CNO,int **Spe_Num_CBasis,
2609 output Spe_Num_Relation
2610 */
Spe_Num_CBasis2Relation(int SpeciesNum,int MaxL,int * Spe_Total_CNO,int ** Spe_Num_CBasis,int ** Spe_Num_Relation)2611 void Spe_Num_CBasis2Relation(
2612 int SpeciesNum,int MaxL,int *Spe_Total_CNO,int **Spe_Num_CBasis,
2613 int **Spe_Num_Relation)
2614 {
2615 int i,j,k,l,id;
2616
2617 for (i=0;i<SpeciesNum;i++) {
2618 id=0;
2619 for (j=0;j<=MaxL;j++) {
2620 if (Spe_Num_CBasis[i][j]>0) {
2621 for (k=0;k<Spe_Num_CBasis[i][j];k++) {
2622 for (l=0;l<2*j+1;l++) {
2623 Spe_Num_Relation[i][id++]=100*j+l+1;
2624 }
2625 }
2626 }
2627 }
2628 for (j=0;j<Spe_Total_CNO[i];j++) {
2629 printf("%d ",Spe_Num_Relation[i][j]);
2630 }
2631 printf("\n");
2632 }
2633
2634
2635 }
2636
2637
2638
2639
2640
input_file_eg(char * file_eg,int * mode,int * NonCol,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 ** Angle0_Spin,double ** Angle1_Spin,double ****** ko)2641 void input_file_eg(char *file_eg,
2642 int *mode, int *NonCol,
2643 int *n, double *emin, double *emax, int *iemin,int *iemax,
2644 int *SpinP_switch, int Kgrid[3], int *atomnum, int **WhatSpecies,
2645 int *SpeciesNum, int **Spe_Total_CNO, int * Max_tnoA, int *MaxL,
2646 int ***Spe_Num_CBasis, int ***Spe_Num_Relation, double *ChemP,
2647 double **Angle0_Spin, double **Angle1_Spin,
2648 double ******ko)
2649 {
2650 double r_vec[10], r_vec2[10];
2651 int i_vec[10], i_vec2[10];
2652 FILE *fp;
2653 int N_Spe_Num_CBasis, i_Spe_Num_CBasis[10];
2654 int j,k;
2655 int N_Spe_Num_Relation, i_Spe_Num_Relation[10];
2656 int N_ko, i_ko[10];
2657 int i,i1,j1,k1,ie,spin;
2658
2659 input_open(file_eg);
2660
2661 /*****************************************
2662 mode
2663
2664 1: normal for cluster and band
2665 2:
2666 3:
2667 4:
2668 5: DC, GDC, and Krylov
2669 6: NEGF
2670 7: Gaussian DOS for collinear case
2671 *****************************************/
2672
2673 input_int("mode",mode,1);
2674
2675 /*****************************************
2676 NonCol
2677
2678 0: collinear calc.
2679 1: non-collinear calc.
2680 *****************************************/
2681
2682 input_int("NonCol",NonCol,0);
2683
2684 /*****************************************
2685 In case of divide-conquer
2686 *****************************************/
2687
2688 if (*mode==5){
2689
2690 input_int("Nspin",SpinP_switch,0);
2691
2692 r_vec2[0]=0.0; r_vec2[1]=0.0;
2693 input_doublev("Erange",2,r_vec,r_vec2);
2694 *emin=r_vec[0]; *emax=r_vec[1];
2695
2696 input_int("atomnum",atomnum,0);
2697
2698 *WhatSpecies = (int*)malloc(sizeof(int)*(*atomnum));
2699 if (fp=input_find("<WhatSpecies") ) {
2700 for (i=0;i<*atomnum;i++) {
2701 fscanf( fp, "%d", &(*WhatSpecies)[i]);
2702 }
2703 if (!input_last("WhatSpecies>")) {
2704 printf("format error near WhatSpecies>\n");
2705 exit(10);
2706 }
2707 }
2708 else {
2709 printf("<WhatSpecies is not found\n");
2710 exit(10);
2711 }
2712
2713 input_int("SpeciesNum",SpeciesNum,0);
2714 *Spe_Total_CNO= (int*)malloc(sizeof(int)*(*SpeciesNum));
2715 if (fp=input_find("<Spe_Total_CNO") ) {
2716 for (i=0;i<*SpeciesNum;i++) {
2717 fscanf( fp, "%d", &(*Spe_Total_CNO)[i]);
2718 }
2719 if (!input_last("Spe_Total_CNO>")) {
2720 printf("format error near Spe_Total_CNO>\n");
2721 exit(10);
2722 }
2723 }
2724 else {
2725 printf("<Spe_Total_CNO is not found\n");
2726 exit(10);
2727 }
2728
2729 *Max_tnoA = (*Spe_Total_CNO)[0];
2730 for (i=1;i<*SpeciesNum;i++) {
2731 if ( *Max_tnoA < (*Spe_Total_CNO)[i] ) *Max_tnoA=(*Spe_Total_CNO)[i];
2732 }
2733 printf("Max of Spe_Total_CNO = %d\n",*Max_tnoA);
2734
2735 input_int("MaxL",MaxL,4);
2736 N_Spe_Num_CBasis=2; i_Spe_Num_CBasis[0]=*SpeciesNum; i_Spe_Num_CBasis[1]=*MaxL+1;
2737 *Spe_Num_CBasis=(int**)malloc_multidimarray("int",N_Spe_Num_CBasis, i_Spe_Num_CBasis);
2738 if (fp=input_find("<Spe_Num_CBasis") ) {
2739 for (i=0;i<*SpeciesNum;i++) {
2740 for (j=0;j<=*MaxL;j++) {
2741 fscanf(fp,"%d",&(*Spe_Num_CBasis)[i][j]);
2742 }
2743 }
2744 if (!input_last("Spe_Num_CBasis>")) {
2745 printf("format error near Spe_Num_CBasis>\n");
2746 exit(10);
2747 }
2748 }
2749 else {
2750 printf("<Spe_Num_CBasis is not found\n");
2751 exit(10);
2752 }
2753
2754 /* relation, index to orbital, 0->s, 1->p ... and so on. */
2755 N_Spe_Num_Relation=2; i_Spe_Num_Relation[0]=*SpeciesNum; i_Spe_Num_Relation[1]=*Max_tnoA;
2756 *Spe_Num_Relation=(int**)malloc_multidimarray("int",N_Spe_Num_Relation, i_Spe_Num_Relation);
2757 Spe_Num_CBasis2Relation(*SpeciesNum,*MaxL,*Spe_Total_CNO,*Spe_Num_CBasis, *Spe_Num_Relation);
2758
2759 input_double("ChemP",ChemP,0.0);
2760
2761 if (*NonCol==1){
2762 *Angle0_Spin = (double*)malloc(sizeof(double)*(*atomnum));
2763 *Angle1_Spin = (double*)malloc(sizeof(double)*(*atomnum));
2764 if (fp=input_find("<SpinAngle") ) {
2765 for (i=0; i<*atomnum; i++) {
2766 fscanf( fp, "%lf %lf", &(*Angle0_Spin)[i],&(*Angle0_Spin)[i]);
2767 }
2768 if (!input_last("SpinAngle>")) {
2769 printf("format error near SpinAngle>\n");
2770 exit(10);
2771 }
2772 }
2773 else {
2774 printf("<SpinAngle is not found\n");
2775 exit(10);
2776 }
2777 }
2778 }
2779
2780 /*****************************************
2781 in case of cluster and band
2782 *****************************************/
2783
2784 else {
2785
2786 input_int("N",n,0);
2787 r_vec2[0]=0.0; r_vec2[1]=0.0;
2788 input_doublev("Erange",2,r_vec,r_vec2);
2789
2790 *emin=r_vec[0]; *emax=r_vec[1];
2791 i_vec2[0]=0; i_vec2[1]=0; i_vec2[2]=0;
2792
2793 input_intv("irange" , 2,i_vec,i_vec2);
2794 *iemin=i_vec[0]; *iemax= i_vec[1];
2795
2796 input_int("Nspin",SpinP_switch,0);
2797
2798 i_vec2[0]=i_vec2[1]=i_vec2[2]=0;
2799 input_intv("Kgrid",3, Kgrid,i_vec2);
2800
2801 input_int("atomnum",atomnum,0);
2802 *WhatSpecies = (int*)malloc(sizeof(int)* (*atomnum));
2803 if (fp=input_find("<WhatSpecies") ) {
2804 for (i=0;i<*atomnum;i++) {
2805 fscanf( fp, "%d", &(*WhatSpecies)[i]);
2806 }
2807 if (!input_last("WhatSpecies>")) {
2808 printf("format error near WhatSpecies>\n");
2809 exit(10);
2810 }
2811 }
2812 else {
2813 printf("<WhatSpecies is not found\n");
2814 exit(10);
2815 }
2816
2817 input_int("SpeciesNum",SpeciesNum,0);
2818 *Spe_Total_CNO= (int*)malloc(sizeof(int)*(*SpeciesNum));
2819 if (fp=input_find("<Spe_Total_CNO") ) {
2820 for (i=0;i<*SpeciesNum;i++) {
2821 fscanf( fp, "%d", &(*Spe_Total_CNO)[i]);
2822 }
2823 if (!input_last("Spe_Total_CNO>")) {
2824 printf("format error near Spe_Total_CNO>\n");
2825 exit(10);
2826 }
2827 }
2828 else {
2829 printf("<Spe_Total_CNO is not found\n");
2830 exit(10);
2831 }
2832
2833 *Max_tnoA = (*Spe_Total_CNO)[0];
2834 for (i=1;i<*SpeciesNum;i++) {
2835 if ( *Max_tnoA < (*Spe_Total_CNO)[i] ) *Max_tnoA=(*Spe_Total_CNO)[i];
2836 }
2837 printf("Max of Spe_Total_CNO = %d\n",*Max_tnoA);
2838
2839 input_int("MaxL",MaxL,4);
2840 N_Spe_Num_CBasis=2; i_Spe_Num_CBasis[0]=*SpeciesNum; i_Spe_Num_CBasis[1]=*MaxL+1;
2841 *Spe_Num_CBasis=(int**)malloc_multidimarray("int",N_Spe_Num_CBasis, i_Spe_Num_CBasis);
2842 if (fp=input_find("<Spe_Num_CBasis") ) {
2843 for (i=0;i<*SpeciesNum;i++) {
2844 for (j=0;j<=*MaxL;j++) {
2845 fscanf(fp,"%d",&(*Spe_Num_CBasis)[i][j]);
2846 }
2847 }
2848 if (!input_last("Spe_Num_CBasis>")) {
2849 printf("format error near Spe_Num_CBasis>\n");
2850 exit(10);
2851 }
2852 }
2853 else {
2854 printf("<Spe_Num_CBasis is not found\n");
2855 exit(10);
2856 }
2857
2858 /* relation, index to orbital, 0->s, 1->p ... and so on. */
2859 N_Spe_Num_Relation=2; i_Spe_Num_Relation[0]=*SpeciesNum; i_Spe_Num_Relation[1]=*Max_tnoA;
2860 *Spe_Num_Relation=(int**)malloc_multidimarray("int",N_Spe_Num_Relation, i_Spe_Num_Relation);
2861 Spe_Num_CBasis2Relation(*SpeciesNum,*MaxL,*Spe_Total_CNO,*Spe_Num_CBasis, *Spe_Num_Relation);
2862
2863 input_double("ChemP",ChemP,0.0);
2864
2865 if (*NonCol==1){
2866 *Angle0_Spin = (double*)malloc(sizeof(double)*(*atomnum));
2867 *Angle1_Spin = (double*)malloc(sizeof(double)*(*atomnum));
2868 if (fp=input_find("<SpinAngle") ) {
2869 for (i=0; i<*atomnum; i++) {
2870 fscanf( fp, "%lf %lf", &(*Angle0_Spin)[i],&(*Angle0_Spin)[i]);
2871 }
2872 if (!input_last("SpinAngle>")) {
2873 printf("format error near SpinAngle>\n");
2874 exit(10);
2875 }
2876 }
2877 else {
2878 printf("<SpinAngle is not found\n");
2879 exit(10);
2880 }
2881 }
2882
2883 if (fp=input_find("<Eigenvalues") ) {
2884 N_ko=5;
2885 i=0; i_ko[i++]=*SpinP_switch+1;
2886 i_ko[i++]=Kgrid[0]; i_ko[i++]=Kgrid[1]; i_ko[i++]=Kgrid[2]; i_ko[i++]= *iemax-*iemin+1;
2887 *ko=(double*****)malloc_multidimarray("double", N_ko,i_ko);
2888
2889 for (i=0;i<Kgrid[0];i++) {
2890 for (j=0;j<Kgrid[1];j++) {
2891 for (k=0;k<Kgrid[2];k++) {
2892 for (spin=0;spin<=*SpinP_switch;spin++) {
2893 fscanf(fp,"%d%d%d",&i1,&j1,&k1);
2894
2895 /*
2896 This part was generalized for the parallel calculation.
2897 if (i1!=i || j1!=j || k1!=k) {
2898 printf("format error in Eigenvalues\n");
2899 exit(10);
2900 }
2901 */
2902
2903 for (ie=*iemin;ie<=*iemax;ie++) {
2904 fscanf(fp,"%lf",&(*ko)[spin][i1][j1][k1][ie-*iemin]);
2905 /* printf("%lf ",ko[spin][i1][j1][k1][ie-iemin]); */
2906 (*ko)[spin][i1][j1][k1][ie-*iemin]=(*ko)[spin][i1][j1][k1][ie-*iemin]-*ChemP;
2907
2908 }
2909 } /* spin */
2910 } /* k*/
2911 } /* j */
2912 } /* i */
2913 if (!input_last("Eigenvalues>")) {
2914 printf("format error near Eigenvalues>\n");
2915 exit(10);
2916 }
2917 } /* if */
2918 else {
2919 printf("can not find <Eigenvalues\n");
2920 exit(10);
2921 }
2922 }
2923
2924 input_close();
2925 }
2926
2927
2928
2929
2930
2931
input_file_ev(char * file_ev,int mode,int NonCol,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)2932 void input_file_ev( char *file_ev, int mode, int NonCol, int Max_tnoA,
2933 int Kgrid[3], int SpinP_switch, int iemin, int iemax,
2934 int atomnum, int *WhatSpecies, int *Spe_Total_CNO,
2935 double ********ev, double ChemP)
2936
2937 {
2938 int N_ev, i_ev[10];
2939 int i,j,k,spin,ie,i1,j1,k1,n;
2940 int ii,jj,kk;
2941 int GA_AN, wanA, tnoA;
2942 int i_vec[10];
2943 float *fSD;
2944 double sum;
2945 char keyword[YOUSO10];
2946 char fp_buf[fp_bsize]; /* setvbuf */
2947 FILE *fp;
2948
2949 /*****************************************
2950 In case of divide-conquer
2951 *****************************************/
2952
2953 if (mode==5){
2954
2955 input_open(file_ev);
2956
2957 Msize1 = (int**)malloc(sizeof(int*)*(SpinP_switch+1));
2958 EigenE = (double***)malloc(sizeof(double**)*(SpinP_switch+1));
2959 PDOS = (double****)malloc(sizeof(double***)*(SpinP_switch+1));
2960
2961 if (NonCol==0){
2962
2963 for (spin=0; spin<=SpinP_switch; spin++){
2964
2965 Msize1[spin] = (int*)malloc(sizeof(int)*atomnum);
2966 EigenE[spin] = (double**)malloc(sizeof(double*)*atomnum);
2967 PDOS[spin] = (double***)malloc(sizeof(double**)*atomnum);
2968
2969 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
2970
2971 printf("read GA_AN=%4d\n",GA_AN);
2972
2973 sprintf(keyword,"<AN%dAN%d",GA_AN+1,spin);
2974
2975 if (fp=input_find(keyword) ) {
2976 wanA = WhatSpecies[GA_AN];
2977 tnoA = Spe_Total_CNO[wanA];
2978
2979 fscanf(fp,"%d",&Msize1[spin][GA_AN]);
2980
2981 EigenE[spin][GA_AN] = (double*)malloc(sizeof(double)*Msize1[spin][GA_AN]);
2982 PDOS[spin][GA_AN] = (double**)malloc(sizeof(double*)*Msize1[spin][GA_AN]);
2983
2984 for (n=0; n<Msize1[spin][GA_AN]; n++){
2985
2986 PDOS[spin][GA_AN][n] = (double*)malloc(sizeof(double)*tnoA);
2987
2988 fscanf(fp,"%d", &j);
2989 fscanf(fp,"%lf",&EigenE[spin][GA_AN][n]);
2990 EigenE[spin][GA_AN][n] = EigenE[spin][GA_AN][n] - ChemP;
2991
2992 for (i=0; i<tnoA; i++) {
2993 fscanf(fp,"%lf",&PDOS[spin][GA_AN][n][i]);
2994 }
2995 }
2996
2997 sprintf(keyword,"AN%dAN%d>",GA_AN+1,spin);
2998 if (!input_last(keyword)) {
2999 printf("format error near AN%dAN>\n",GA_AN+1);
3000 exit(10);
3001 }
3002
3003 }
3004 else {
3005 printf("<AN%dAN is not found\n",GA_AN+1);
3006 exit(10);
3007 }
3008 }
3009 }
3010 }
3011
3012 else{
3013
3014 Msize1[0] = (int*)malloc(sizeof(int)*atomnum);
3015 Msize1[1] = (int*)malloc(sizeof(int)*atomnum);
3016 EigenE[0] = (double**)malloc(sizeof(double*)*atomnum);
3017 EigenE[1] = (double**)malloc(sizeof(double*)*atomnum);
3018 PDOS[0] = (double***)malloc(sizeof(double**)*atomnum);
3019 PDOS[1] = (double***)malloc(sizeof(double**)*atomnum);
3020
3021
3022 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
3023
3024 printf("read GA_AN=%4d\n",GA_AN);
3025
3026 sprintf(keyword,"<AN%d",GA_AN+1);
3027
3028 if (fp=input_find(keyword) ) {
3029
3030 wanA = WhatSpecies[GA_AN];
3031 tnoA = Spe_Total_CNO[wanA];
3032
3033 fscanf(fp,"%d %d",&Msize1[0][GA_AN],&Msize1[1][GA_AN]);
3034
3035 EigenE[0][GA_AN] = (double*)malloc(sizeof(double)*Msize1[0][GA_AN]);
3036 EigenE[1][GA_AN] = (double*)malloc(sizeof(double)*Msize1[1][GA_AN]);
3037 PDOS[0][GA_AN] = (double**)malloc(sizeof(double*)*Msize1[0][GA_AN]);
3038 PDOS[1][GA_AN] = (double**)malloc(sizeof(double*)*Msize1[1][GA_AN]);
3039
3040 for (n=0; n<Msize1[0][GA_AN]; n++){
3041
3042 PDOS[0][GA_AN][n] = (double*)malloc(sizeof(double)*tnoA);
3043 PDOS[1][GA_AN][n] = (double*)malloc(sizeof(double)*tnoA);
3044
3045 fscanf(fp,"%d", &j);
3046 fscanf(fp,"%lf %lf",&EigenE[0][GA_AN][n],&EigenE[1][GA_AN][n]);
3047 EigenE[0][GA_AN][n] = EigenE[0][GA_AN][n] - ChemP;
3048 EigenE[1][GA_AN][n] = EigenE[1][GA_AN][n] - ChemP;
3049
3050 for (i=0; i<tnoA; i++) {
3051 fscanf(fp,"%lf %lf",&PDOS[0][GA_AN][n][i],&PDOS[1][GA_AN][n][i]);
3052 }
3053 }
3054
3055 sprintf(keyword,"AN%d>",GA_AN+1);
3056 if (!input_last(keyword)) {
3057 printf("format error near AN%dAN>\n",GA_AN+1);
3058 exit(10);
3059 }
3060
3061 }
3062 else {
3063 printf("<AN%dAN is not found\n",GA_AN+1);
3064 exit(10);
3065 }
3066 }
3067
3068 }
3069
3070 input_close();
3071
3072 }
3073
3074 /*****************************************
3075 In case of cluster and band
3076 *****************************************/
3077
3078 else {
3079
3080 if ( (fp=fopen(file_ev,"r"))==NULL ) {
3081
3082 #ifdef xt3
3083 setvbuf(fp,fp_buf,_IOFBF,fp_bsize); /* setvbuf */
3084 #endif
3085
3086 printf("cannot open a file %s\n", file_ev);
3087 exit(10);
3088 }
3089
3090 N_ev=7;
3091 i=0; i_ev[i++]=SpinP_switch+1;
3092 i_ev[i++]=Kgrid[0]; i_ev[i++]=Kgrid[1]; i_ev[i++]=Kgrid[2];
3093 i_ev[i++] = atomnum ; i_ev[i++] = Max_tnoA;
3094 i_ev[i++]= iemax-iemin+1;
3095 (*ev)=(double*******)malloc_multidimarray("double", N_ev,i_ev);
3096
3097 fSD=(float*)malloc(sizeof(float)*Max_tnoA*2);
3098
3099 if (NonCol==0){
3100
3101 for (i=0;i<Kgrid[0];i++) {
3102 for (j=0;j<Kgrid[1];j++) {
3103 for (k=0;k<Kgrid[2];k++) {
3104 for (spin=0;spin<=SpinP_switch;spin++) {
3105 for (ie=iemin;ie<=iemax;ie++) {
3106
3107 fread(i_vec,sizeof(int),3,fp);
3108 ii=i_vec[0]; jj=i_vec[1]; kk=i_vec[2];
3109
3110 /*
3111 This part was generalized for the parallel calculation.
3112 if (ii!=i || jj!=j || kk!=k) {
3113 printf("parameter error in reading Eigenvectors\n");
3114 exit(10);
3115 }
3116 */
3117
3118 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
3119 wanA = WhatSpecies[GA_AN];
3120 tnoA = Spe_Total_CNO[wanA];
3121
3122 fread(fSD,sizeof(float),tnoA,fp);
3123 for (i1=0; i1<=(tnoA-1); i1++){
3124 (*ev)[spin][ii][jj][kk][GA_AN][i1][ie-iemin] = fSD[i1];
3125 }
3126 }
3127
3128 /* printf("\n"); */
3129 /* check normalization */
3130 sum=0.0;
3131 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
3132 wanA = WhatSpecies[GA_AN];
3133 tnoA = Spe_Total_CNO[wanA];
3134 for (i1=0; i1<=(tnoA-1); i1++){
3135 sum+= (*ev)[spin][ii][jj][kk][GA_AN][i1][ie-iemin];
3136 }
3137 }
3138
3139 /*
3140 if (fabs(sum-1.0) > 1.0e-3) {
3141 printf("%d %d %d %d %d, sum=%lf\n",spin,ii,jj,kk,ie,sum);
3142 }
3143 */
3144
3145 }
3146 }
3147 }
3148 }
3149 }
3150 }
3151
3152 else{
3153
3154 for (i=0;i<Kgrid[0];i++) {
3155 for (j=0;j<Kgrid[1];j++) {
3156 for (k=0;k<Kgrid[2];k++) {
3157
3158 for (ie=iemin; ie<=iemax; ie++) {
3159
3160 fread(i_vec,sizeof(int),3,fp);
3161 ii=i_vec[0]; jj=i_vec[1]; kk=i_vec[2];
3162
3163 /*
3164 This part was generalized for the parallel calculation.
3165 if (ii!=i || jj!=j || kk!=k) {
3166 printf("parameter error in reading Eigenvectors\n");
3167 exit(10);
3168 }
3169 */
3170
3171 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
3172 wanA = WhatSpecies[GA_AN];
3173 tnoA = Spe_Total_CNO[wanA];
3174
3175 fread(fSD,sizeof(float),2*tnoA,fp);
3176
3177 for (i1=0; i1<tnoA; i1++){
3178 (*ev)[0][ii][jj][kk][GA_AN][i1][ie-iemin] = fSD[i1];
3179 }
3180
3181 for (i1=0; i1<tnoA; i1++){
3182 (*ev)[1][ii][jj][kk][GA_AN][i1][ie-iemin] = fSD[tnoA+i1];
3183 }
3184 }
3185
3186 /* printf("\n"); */
3187 /* check normalization */
3188 /*
3189 sum=0.0;
3190 for (GA_AN=0; GA_AN<atomnum; GA_AN++){
3191 wanA = WhatSpecies[GA_AN];
3192 tnoA = Spe_Total_CNO[wanA];
3193 for (i1=0; i1<=(tnoA-1); i1++){
3194 sum+= (*ev)[0][ii][jj][kk][GA_AN][i1][ie-iemin];
3195 }
3196 }
3197 */
3198
3199 /*
3200 if (fabs(sum-1.0) > 1.0e-3) {
3201 printf("%d %d %d %d %d, sum=%lf\n",spin,ii,jj,kk,ie,sum);
3202 }
3203 */
3204
3205 }
3206 }
3207 }
3208 }
3209 }
3210
3211 free(fSD);
3212 fclose(fp);
3213 }
3214
3215
3216 }
3217
3218
3219
input_main(int mode,int Kgrid[3],int atomnum,int * method,int * todo,double * gaussian,int * pdos_n,int ** pdos_atoms)3220 void input_main( int mode, int Kgrid[3], int atomnum,
3221 int *method, int *todo,
3222 double *gaussian,
3223 int *pdos_n, int **pdos_atoms)
3224 {
3225 char buf[1000],buf2[1000],*c;
3226 int i;
3227
3228 if (mode==5){
3229 printf("The input data caluculated by the divide-conquer method\n");
3230 printf("Gaussian Broadening is employed\n");
3231 *method=2;
3232 }
3233
3234 else if (mode==6){
3235 printf("The input data caluculated by the NEGF method\n");
3236 *method=4;
3237 }
3238
3239 else if (mode==7){
3240 printf("The input data caluculated by the Gaussian broadening method\n");
3241 *method=4;
3242 }
3243
3244 else{
3245 if (Kgrid[0]==1 && Kgrid[1]==1 && Kgrid[2]==1 ) {
3246 printf("Kgrid= 1 1 1 => Gaussian Broadening is employed\n");
3247 *method=2;
3248 }
3249 else {
3250 printf("Which method do you use?, Tetrahedron(1), Gaussian Broadening(2)\n");
3251 fgets(buf,1000,stdin); sscanf(buf,"%d",method);
3252 }
3253 }
3254
3255 if (*method==2) {
3256 printf("Please input a value of gaussian (double) (eV)\n");
3257 fgets(buf,1000,stdin); sscanf(buf,"%lf",gaussian); *gaussian = *gaussian/eV2Hartree;
3258 }
3259
3260 printf("Do you want Dos(1) or PDos(2)?\n");
3261 fgets(buf,1000,stdin);sscanf(buf,"%d",todo);
3262
3263 if (*todo==2) {
3264 printf("\nNumber of atoms=%d\n",atomnum);
3265 if (atomnum==1) {
3266 (*pdos_atoms)=(int*)malloc(sizeof(int)*1);
3267 *pdos_n=1;
3268 (*pdos_atoms)[0]=1;
3269 }
3270 else {
3271 printf("Which atoms for PDOS : (1,...,%d), ex 1 2\n",atomnum);
3272 fgets(buf,1000,stdin);
3273 strcpy(buf2,buf);
3274 /* printf("<%s>\n",buf); */
3275 *pdos_n=0;
3276 c=strtok(buf," ") ;
3277 while ( c ) {
3278 /* printf("<%s> ",c); */
3279 if (sscanf(c,"%d",&i)==1) (*pdos_n)++;
3280 c=strtok(NULL," ");
3281 }
3282 /* printf("pdos_n=%d\n",*pdos_n); */
3283 *pdos_atoms=(int*)malloc(sizeof(int)*(*pdos_n));
3284 for (c=strtok(buf2," ") , i=0;i<*pdos_n;c=strtok(NULL," "),i++) {
3285 /* printf("<%s> ",c); */
3286 sscanf(c,"%d",&(*pdos_atoms)[i]);
3287 }
3288 }
3289 printf("pdos_n=%d\n",*pdos_n);
3290 for (i=0;i<*pdos_n;i++) {
3291 printf("%d ",(*pdos_atoms)[i]);
3292 }
3293 printf("\n");
3294 }
3295
3296 }
3297
3298