1 #include <math.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include "Inputtools.h"
6
7 #define asize10 100
8 #define BohrR 0.529177249
9 #define print_std 0
10
11 static void Input_NEGF(char *argv[]);
12 static void SpeciesString2int(int p);
13 static int Species2int(char Species[asize10]);
14 static int SEQ(char str1[asize10], char str2[asize10]);
15 static void fnjoint2(char name1[asize10], char name2[asize10],
16 char name3[asize10], char name4[asize10]);
17 static int Find_Size_LeftLeads();
18 static int Find_Size_RightLeads();
19 static void Make_new_inputfile(char *argv[], int NL, int NR);
20 static void FreeArrays();
21
22 static int atomnum,Catomnum;
23 static int Latomnum,Ratomnum;
24 static int SpeciesNum;
25 static int coordinates_unit;
26 static int unitvector_unit;
27 static int Solver,Lbuffer,Rbuffer;
28
29 static int *C_WhatSpecies;
30 static int *L_WhatSpecies;
31 static int *R_WhatSpecies;
32 static int **Spe_Num_Basis;
33 static int **Spe_Num_CBasis;
34 static int *Spe_MaxL_Basis;
35
36 static double tv[4][4];
37 static double Left_tv[4][4];
38 static double Right_tv[4][4];
39 static double **Cxyz;
40 static double **Lxyz;
41 static double **Rxyz;
42 static double **Cdens;
43 static double **Ldens;
44 static double **Rdens;
45 static double *Spe_Atom_Cut1;
46
47 static char **SpeName;
48 static char **SpeBasis;
49 static char **SpeBasisName;
50 static char **SpeVPS;
51
main(int argc,char * argv[])52 int main(int argc, char *argv[])
53 {
54 static int NL,NR;
55
56 Input_NEGF(argv);
57 NL = Find_Size_LeftLeads();
58 NR = Find_Size_RightLeads();
59 Make_new_inputfile(argv,NL,NR);
60 FreeArrays();
61 }
62
Make_new_inputfile(char * argv[],int NL,int NR)63 void Make_new_inputfile(char *argv[], int NL, int NR)
64 {
65
66 FILE *fp,*fp1,*fp2;
67 char *sp0,*sp1,*sp2;
68 char *sp3,*sp4,*sp5;
69 char fname2[asize10];
70 char command0[3*asize10];
71 static int po,po1,i,num,cell,spe;
72 static double x,y,z;
73
74 sprintf(fname2,"%s_NEGF1",argv[1]);
75
76 printf("\n");
77 printf("The new input file, %s, was generated.\n",fname2);
78 printf("\n");
79
80 if ((fp2 = fopen(fname2,"w")) != NULL){
81 if ((fp = fopen(argv[1],"r")) != NULL){
82
83 po = 0;
84 do {
85 if (fgets(command0,asize10,fp)!=NULL){
86 command0[strlen(command0)-1] = '\0';
87
88 sp0 = strstr(command0,"LeftLeadAtoms.Number");
89 sp1 = strstr(command0,"<LeftLeadAtoms.SpeciesAndCoordinates");
90 sp2 = strstr(command0,"<LeftLeadAtoms.UnitVectors");
91 sp3 = strstr(command0,"RightLeadAtoms.Number");
92 sp4 = strstr(command0,"<RightLeadAtoms.SpeciesAndCoordinates");
93 sp5 = strstr(command0,"<RightLeadAtoms.UnitVectors");
94
95 /****************************************************
96 LeftLeadAtoms.Number
97 ****************************************************/
98
99 if (sp0!=NULL){
100 fprintf(fp2,"LeftLeadAtoms.Number %i\n",NL*Latomnum);
101 }
102
103 /****************************************************
104 LeftLeadAtoms.SpeciesAndCoordinates
105 ****************************************************/
106
107 else if (sp1!=NULL){
108 fprintf(fp2,"%s\n",command0);
109
110 for (i=0; i<=Latomnum; i++){
111 fgets(command0,asize10,fp);
112 command0[strlen(command0)-1] = '\0';
113 }
114
115 num = 0;
116 for (cell=0; cell<NL; cell++){
117 for (i=1; i<=Latomnum; i++){
118
119 spe = L_WhatSpecies[i];
120
121 x = Lxyz[i][1] + (double)cell*Left_tv[3][1];
122 y = Lxyz[i][2] + (double)cell*Left_tv[3][2];
123 z = Lxyz[i][3] + (double)cell*Left_tv[3][3];
124 num++;
125
126 if (coordinates_unit==0){
127 fprintf(fp2," %4d %3s %15.12f %15.12f %15.12f %6.3f %6.3f\n",
128 num,SpeName[spe],x*BohrR,y*BohrR,z*BohrR,
129 Ldens[i][0],Ldens[i][1]);
130 }
131 else{
132 fprintf(fp2," %4d %3s %15.12f %15.12f %15.12f %6.3f %6.3f\n",
133 num,SpeName[spe],x,y,z,Ldens[i][0],Ldens[i][1]);
134 }
135
136 }
137 }
138 fprintf(fp2,"LeftLeadAtoms.SpeciesAndCoordinates>\n");
139
140 }
141
142 /****************************************************
143 LeftLeadAtoms.UnitVectors
144 ****************************************************/
145
146 else if (sp2!=NULL){
147 fprintf(fp2,"%s\n",command0);
148
149 for (i=1; i<=2; i++){
150 fgets(command0,asize10,fp);
151 command0[strlen(command0)-1] = '\0';
152 fprintf(fp2,"%s\n",command0);
153 }
154 for (i=1; i<=2; i++){
155 fgets(command0,asize10,fp);
156 command0[strlen(command0)-1] = '\0';
157 }
158
159 if (coordinates_unit==0){
160 fprintf(fp2," %15.12f %15.12f %15.12f\n",
161 (double)NL*Left_tv[3][1]*BohrR,
162 (double)NL*Left_tv[3][2]*BohrR,
163 (double)NL*Left_tv[3][3]*BohrR);
164 }
165 else {
166 fprintf(fp2," %15.12f %15.12f %15.12f\n",
167 (double)NL*Left_tv[3][1],
168 (double)NL*Left_tv[3][2],
169 (double)NL*Left_tv[3][3]);
170 }
171
172 fprintf(fp2,"LeftLeadAtoms.UnitVectors>\n");
173 }
174
175 /****************************************************
176 RightLeadAtoms.Number
177 ****************************************************/
178
179 else if (sp3!=NULL){
180 fprintf(fp2,"RightLeadAtoms.Number %i\n",NR*Ratomnum);
181 }
182
183 /****************************************************
184 RightLeadAtoms.SpeciesAndCoordinates
185 ****************************************************/
186
187 else if (sp4!=NULL){
188 fprintf(fp2,"%s\n",command0);
189
190 for (i=0; i<=Ratomnum; i++){
191 fgets(command0,asize10,fp);
192 command0[strlen(command0)-1] = '\0';
193 }
194
195 num = 0;
196 for (cell=0; cell<NR; cell++){
197 for (i=1; i<=Ratomnum; i++){
198
199 spe = R_WhatSpecies[i];
200
201 x = Rxyz[i][1] + (double)cell*Right_tv[3][1];
202 y = Rxyz[i][2] + (double)cell*Right_tv[3][2];
203 z = Rxyz[i][3] + (double)cell*Right_tv[3][3];
204 num++;
205
206 if (coordinates_unit==0){
207 fprintf(fp2," %4d %3s %15.12f %15.12f %15.12f %6.3f %6.3f\n",
208 num,SpeName[spe],x*BohrR,y*BohrR,z*BohrR,Rdens[i][0],Rdens[i][1]);
209 }
210 else{
211 fprintf(fp2," %4d %3s %15.12f %15.12f %15.12f %6.3f %6.3f\n",
212 num,SpeName[spe],x,y,z,Rdens[i][0],Rdens[i][1]);
213 }
214
215 }
216 }
217 fprintf(fp2,"RightLeadAtoms.SpeciesAndCoordinates>\n");
218
219 }
220
221 /****************************************************
222 RightLeadAtoms.UnitVectors
223 ****************************************************/
224
225
226 else if (sp5!=NULL){
227
228 fprintf(fp2,"%s\n",command0);
229
230 for (i=1; i<=2; i++){
231 fgets(command0,asize10,fp);
232 command0[strlen(command0)-1] = '\0';
233 fprintf(fp2,"%s\n",command0);
234 }
235 for (i=1; i<=2; i++){
236 fgets(command0,asize10,fp);
237 command0[strlen(command0)-1] = '\0';
238 }
239
240 if (coordinates_unit==0){
241 fprintf(fp2," %15.12f %15.12f %15.12f\n",
242 (double)NR*Right_tv[3][1]*BohrR,
243 (double)NR*Right_tv[3][2]*BohrR,
244 (double)NR*Right_tv[3][3]*BohrR);
245 }
246 else {
247 fprintf(fp2," %15.12f %15.12f %15.12f\n",
248 (double)NR*Right_tv[3][1],
249 (double)NR*Right_tv[3][2],
250 (double)NR*Right_tv[3][3]);
251 }
252
253 fprintf(fp2,"RightLeadAtoms.UnitVectors>\n");
254 }
255
256 /****************************************************
257 Others
258 ****************************************************/
259
260 else{
261 fprintf(fp2,"%s\n",command0);
262 }
263
264 }
265 else{
266 po = 1;
267 }
268 }while(po==0);
269
270 fclose(fp);
271 }
272 else{
273 printf("Could not find %s\n",argv[1]);
274 exit(0);
275 }
276
277 fclose(fp2);
278 }
279 else{
280 printf("Could not open %s\n",fname2);
281 exit(0);
282 }
283
284 }
285
Find_Size_LeftLeads()286 int Find_Size_LeftLeads()
287 {
288
289 static int cell,i,j,po,spe1,spe2;
290 static int Num_LCcell,Num_LLcell,Num_Lcell;
291 static int TFNAN_p,TFNAN_c;
292 static double dx,dy,dz;
293 static double R,r1,r2;
294
295 /****************************************************
296 interations between central and left regions
297 ****************************************************/
298
299 po = 0;
300 Num_LCcell = 0;
301 TFNAN_c = 0;
302
303 do{
304
305 Num_LCcell++;
306
307 TFNAN_p = TFNAN_c;
308 TFNAN_c = 0;
309
310 for (i=1; i<=Catomnum; i++){
311 spe1 = C_WhatSpecies[i];
312 r1 = Spe_Atom_Cut1[spe1];
313
314 for (cell=0; cell<Num_LCcell; cell++){
315 for (j=1; j<=Latomnum; j++){
316
317 spe2 = L_WhatSpecies[j];
318 r2 = Spe_Atom_Cut1[spe2];
319
320 dx = Cxyz[i][1] - Lxyz[j][1] - (double)cell*Left_tv[3][1];
321 dy = Cxyz[i][2] - Lxyz[j][2] - (double)cell*Left_tv[3][2];
322 dz = Cxyz[i][3] - Lxyz[j][3] - (double)cell*Left_tv[3][3];
323 R = sqrt(dx*dx + dy*dy + dz*dz);
324 if ( R<(r1+r2) ) TFNAN_c++;
325 }
326 }
327 }
328
329 if ((TFNAN_c-TFNAN_p)==0) po = 1;
330
331 }while(po==0);
332
333 Num_LCcell--;
334
335 /****************************************************
336 interations between left0 and left regions
337 ****************************************************/
338
339 po = 0;
340 Num_LLcell = 0;
341 TFNAN_c = 0;
342
343 do{
344
345 Num_LLcell++;
346
347 TFNAN_p = TFNAN_c;
348 TFNAN_c = 0;
349
350 for (i=1; i<=Latomnum; i++){
351 spe1 = L_WhatSpecies[i];
352 r1 = Spe_Atom_Cut1[spe1];
353
354 for (cell=0; cell<Num_LLcell; cell++){
355 for (j=1; j<=Latomnum; j++){
356 spe2 = L_WhatSpecies[j];
357 r2 = Spe_Atom_Cut1[spe2];
358
359 dx = Lxyz[i][1] - Lxyz[j][1] - (double)cell*Left_tv[3][1];
360 dy = Lxyz[i][2] - Lxyz[j][2] - (double)cell*Left_tv[3][2];
361 dz = Lxyz[i][3] - Lxyz[j][3] - (double)cell*Left_tv[3][3];
362 R = sqrt(dx*dx + dy*dy + dz*dz);
363 if ( R<(r1+r2) ) TFNAN_c++;
364 }
365 }
366 }
367
368 if ((TFNAN_c-TFNAN_p)==0) po = 1;
369
370 }while(po==0);
371
372 Num_LLcell = Num_LLcell - 2;
373
374 if (Num_LLcell<Num_LCcell) Num_Lcell = Num_LCcell;
375 else Num_Lcell = Num_LLcell;
376
377 printf("\n");
378 printf("************ Left lead ************\n");
379 printf(" Minimum cell size of the left lead\n");
380 printf(" for the central region = %i\n",Num_LCcell);
381 printf(" in the left lead = %i\n",Num_LLcell);
382 printf(" Requested buffer size = %i\n",Lbuffer);
383 printf(" Size of an adjusted left lead = %i\n",Num_Lcell+Lbuffer);
384
385 Num_Lcell = Num_Lcell + Lbuffer;
386
387 return Num_Lcell;
388 }
389
390
Find_Size_RightLeads()391 int Find_Size_RightLeads()
392 {
393
394 static int cell,i,j,po,spe1,spe2;
395 static int Num_RCcell,Num_RRcell,Num_Rcell;
396 static int TFNAN_p,TFNAN_c;
397 static double dx,dy,dz;
398 static double R,r1,r2;
399
400 /****************************************************
401 interations between central and right regions
402 ****************************************************/
403
404 po = 0;
405 Num_RCcell = 0;
406 TFNAN_c = 0;
407
408 do{
409
410 Num_RCcell++;
411
412 TFNAN_p = TFNAN_c;
413 TFNAN_c = 0;
414
415 for (i=1; i<=Catomnum; i++){
416 spe1 = C_WhatSpecies[i];
417 r1 = Spe_Atom_Cut1[spe1];
418
419 for (cell=0; cell<Num_RCcell; cell++){
420 for (j=1; j<=Ratomnum; j++){
421
422 spe2 = R_WhatSpecies[j];
423 r2 = Spe_Atom_Cut1[spe2];
424
425 dx = Cxyz[i][1] - Rxyz[j][1] - (double)cell*Right_tv[3][1];
426 dy = Cxyz[i][2] - Rxyz[j][2] - (double)cell*Right_tv[3][2];
427 dz = Cxyz[i][3] - Rxyz[j][3] - (double)cell*Right_tv[3][3];
428 R = sqrt(dx*dx + dy*dy + dz*dz);
429 if ( R<(r1+r2) ) TFNAN_c++;
430 }
431 }
432 }
433
434 if ((TFNAN_c-TFNAN_p)==0) po = 1;
435
436 }while(po==0);
437
438 Num_RCcell--;
439
440 /****************************************************
441 interations between right0 and right regions
442 ****************************************************/
443
444 po = 0;
445 Num_RRcell = 0;
446 TFNAN_c = 0;
447
448 do{
449
450 Num_RRcell++;
451
452 TFNAN_p = TFNAN_c;
453 TFNAN_c = 0;
454
455 for (i=1; i<=Ratomnum; i++){
456 spe1 = R_WhatSpecies[i];
457 r1 = Spe_Atom_Cut1[spe1];
458
459 for (cell=0; cell<Num_RRcell; cell++){
460 for (j=1; j<=Ratomnum; j++){
461 spe2 = R_WhatSpecies[j];
462 r2 = Spe_Atom_Cut1[spe2];
463
464 dx = Rxyz[i][1] - Rxyz[j][1] - (double)cell*Right_tv[3][1];
465 dy = Rxyz[i][2] - Rxyz[j][2] - (double)cell*Right_tv[3][2];
466 dz = Rxyz[i][3] - Rxyz[j][3] - (double)cell*Right_tv[3][3];
467 R = sqrt(dx*dx + dy*dy + dz*dz);
468 if ( R<(r1+r2) ) TFNAN_c++;
469 }
470 }
471 }
472
473 if ((TFNAN_c-TFNAN_p)==0) po = 1;
474
475 }while(po==0);
476
477 Num_RRcell = Num_RRcell - 2;
478
479 if (Num_RRcell<Num_RCcell) Num_Rcell = Num_RCcell;
480 else Num_Rcell = Num_RRcell;
481
482 printf("\n");
483 printf("************ Right lead ************\n");
484 printf(" Minimum cell size of the right lead\n");
485 printf(" for the central region = %i\n",Num_RCcell);
486 printf(" in the right lead = %i\n",Num_RRcell);
487 printf(" Requested buffer size = %i\n",Rbuffer);
488 printf(" Size of an adjusted left lead = %i\n",Num_Rcell+Rbuffer);
489
490 Num_Rcell = Num_Rcell + Rbuffer;
491
492 return Num_Rcell;
493 }
494
495
496
497
Input_NEGF(char * argv[])498 void Input_NEGF(char *argv[])
499 {
500 FILE *fp;
501 int i,j,ct_AN,spe;
502 int po=0; /* error count */
503 double r_vec[20],r_vec2[20];
504 int i_vec[20],i_vec2[20];
505 char *s_vec[20];
506 char Species[20];
507 char c;
508 double tmp0,tmp1;
509 char DirPAO[asize10];
510 char ExtPAO[asize10] = ".pao";
511 char FN_PAO[asize10];
512
513 printf("\n*******************************************\n");
514 printf(" This program adjusts the size of leads. \n");
515 printf("*******************************************\n\n");
516
517 /****************************************************
518 open a file
519 ****************************************************/
520
521 input_open(argv[1]);
522
523 /****************************************************
524 set DirPAO and DirVPS
525 ****************************************************/
526
527 sprintf(DirPAO,"%s/PAO/","/usr/local/share/openmx/DFT_DATA13");
528
529 /****************************************************
530 Number of Atomic Species
531 ****************************************************/
532
533 input_int("Species.Number",&SpeciesNum,0);
534
535 if (SpeciesNum<=0){
536 printf("Species.Number may be wrong.\n");
537 po++;
538 }
539
540 /****************************************************
541 allocation of arrays:
542
543 char SpeName[SpeciesNum][asize10]
544 char SpeBasis[SpeciesNum][asize10]
545 char SpeBasisName[SpeciesNum][asize10]
546 char SpeVPS[SpeciesNum][asize10]
547 int Spe_Num_Basis[SpeciesNum][6];
548 int Spe_Num_CBasis[SpeciesNum][6];
549 int Spe_MaxL_Basis[SpeciesNum];
550 ****************************************************/
551
552 SpeName = (char**)malloc(sizeof(char*)*SpeciesNum);
553 for (i=0; i<SpeciesNum; i++){
554 SpeName[i] = (char*)malloc(sizeof(long)*asize10);
555 }
556
557 SpeBasis = (char**)malloc(sizeof(char*)*SpeciesNum);
558 for (i=0; i<SpeciesNum; i++){
559 SpeBasis[i] = (char*)malloc(sizeof(char)*asize10);
560 }
561
562 SpeBasisName = (char**)malloc(sizeof(char*)*SpeciesNum);
563 for (i=0; i<SpeciesNum; i++){
564 SpeBasisName[i] = (char*)malloc(sizeof(char)*asize10);
565 }
566
567 SpeVPS = (char**)malloc(sizeof(char*)*SpeciesNum);
568 for (i=0; i<SpeciesNum; i++){
569 SpeVPS[i] = (char*)malloc(sizeof(char)*asize10);
570 }
571
572 Spe_Num_Basis = (int**)malloc(sizeof(int*)*SpeciesNum);
573 for (i=0; i<SpeciesNum; i++){
574 Spe_Num_Basis[i] = (int*)malloc(sizeof(int)*6);
575 }
576
577 Spe_Num_CBasis = (int**)malloc(sizeof(int*)*SpeciesNum);
578 for (i=0; i<SpeciesNum; i++){
579 Spe_Num_CBasis[i] = (int*)malloc(sizeof(int)*6);
580 }
581
582 Spe_MaxL_Basis = (int*)malloc(sizeof(int)*SpeciesNum);
583 Spe_Atom_Cut1 = (double*)malloc(sizeof(double)*SpeciesNum);
584
585 /****************************************************
586 definition of Atomic Species
587 ****************************************************/
588
589 if (fp=input_find("<Definition.of.Atomic.Species")) {
590 for (i=0; i<SpeciesNum; i++){
591 fscanf(fp,"%s %s %s",SpeName[i],SpeBasis[i],SpeVPS[i]);
592 SpeciesString2int(i);
593 }
594 if (! input_last("Definition.of.Atomic.Species>")) {
595 /* format error */
596 printf("Format error for Definition.of.Atomic.Species\n");
597 po++;
598 }
599 }
600
601 if (print_std==1){
602 for (i=0; i<SpeciesNum; i++){
603 printf("<Input_std> %i Name %s\n",i,SpeName[i]);
604 printf("<Input_std> %i Basis %s\n",i,SpeBasis[i]);
605 printf("<Input_std> %i VPS %s\n",i,SpeVPS[i]);
606 }
607 }
608
609 /****************************************************
610 eigensolver
611 ****************************************************/
612
613 s_vec[0]="Recursion"; s_vec[1]="Cluster"; s_vec[2]="Band";
614 s_vec[3]="NEGF"; s_vec[4]="DC"; s_vec[5]="GDC";
615 s_vec[6]="Cluster-DIIS"; s_vec[7]="Krylov";
616
617 i_vec[0]=1; i_vec[1]=2; i_vec[2]=3;
618 i_vec[3]=4; i_vec[4]=5; i_vec[5]=6;
619 i_vec[6]=7; i_vec[7]=8;
620
621 input_string2int("scf.EigenvalueSolver", &Solver, 8, s_vec,i_vec);
622
623 /****************************************************
624 Atoms
625 ****************************************************/
626
627 /* except for NEGF */
628
629 if (Solver!=4){
630 printf("Plese set NEGF for scf.EigenvalueSolve\n");
631 exit(0);
632 } /* if (Solver!=4){ */
633
634 /* NEGF */
635
636 else{
637
638 /* center */
639 input_int("Atoms.Number",&Catomnum,0);
640 if (Catomnum<=0){
641 printf("Atoms.Number may be wrong.\n");
642 po++;
643 }
644
645 /* left */
646 input_int("LeftLeadAtoms.Number",&Latomnum,0);
647 if (Latomnum<=0){
648 printf("LeftLeadAtoms.Number may be wrong.\n");
649 po++;
650 }
651
652 /* right */
653 input_int("RightLeadAtoms.Number",&Ratomnum,0);
654 if (Ratomnum<=0){
655 printf("RightLeadAtoms.Number may be wrong.\n");
656 po++;
657 }
658
659 atomnum = Catomnum + Latomnum + Ratomnum;
660
661 /****************************************************
662 allocation of arrays:
663
664 int C_WhatSpecies[Catomnum+1];
665 int L_WhatSpecies[Latomnum+1];
666 int R_WhatSpecies[Ratomnum+1];
667 double Cxyz[Catomnum+1][4];
668 double Lxyz[Latomnum+1][4];
669 double Rxyz[Ratomnum+1][4];
670 double Cdens[Catomnum+1][2];
671 double Ldens[Latomnum+1][2];
672 double Rdens[Ratomnum+1][2];
673 ****************************************************/
674
675 C_WhatSpecies = (int*)malloc(sizeof(int)*(Catomnum+1));
676 L_WhatSpecies = (int*)malloc(sizeof(int)*(Latomnum+1));
677 R_WhatSpecies = (int*)malloc(sizeof(int)*(Ratomnum+1));
678
679 Cxyz = (double**)malloc(sizeof(double*)*(Catomnum+1));
680 for (ct_AN=0; ct_AN<=Catomnum; ct_AN++){
681 Cxyz[ct_AN] = (double*)malloc(sizeof(double)*4);
682 }
683
684 Lxyz = (double**)malloc(sizeof(double*)*(Latomnum+1));
685 for (ct_AN=0; ct_AN<=Latomnum; ct_AN++){
686 Lxyz[ct_AN] = (double*)malloc(sizeof(double)*4);
687 }
688
689 Rxyz = (double**)malloc(sizeof(double*)*(Ratomnum+1));
690 for (ct_AN=0; ct_AN<=Ratomnum; ct_AN++){
691 Rxyz[ct_AN] = (double*)malloc(sizeof(double)*4);
692 }
693
694 Cdens = (double**)malloc(sizeof(double*)*(Catomnum+1));
695 for (ct_AN=0; ct_AN<=Catomnum; ct_AN++){
696 Cdens[ct_AN] = (double*)malloc(sizeof(double)*2);
697 }
698
699 Ldens = (double**)malloc(sizeof(double*)*(Latomnum+1));
700 for (ct_AN=0; ct_AN<=Latomnum; ct_AN++){
701 Ldens[ct_AN] = (double*)malloc(sizeof(double)*2);
702 }
703
704 Rdens = (double**)malloc(sizeof(double*)*(Ratomnum+1));
705 for (ct_AN=0; ct_AN<=Ratomnum; ct_AN++){
706 Rdens[ct_AN] = (double*)malloc(sizeof(double)*2);
707 }
708
709 /****************************************************
710 coordinates unit
711 ****************************************************/
712
713 s_vec[0]="Ang"; s_vec[1]="AU";
714 i_vec[0]= 0; i_vec[1]=1;
715 input_string2int("Atoms.SpeciesAndCoordinates.Unit",
716 &coordinates_unit,2,s_vec,i_vec);
717
718 /****************************************************
719 coordinates of center
720 ****************************************************/
721
722 if (fp=input_find("<Atoms.SpeciesAndCoordinates") ) {
723 for (i=1; i<=Catomnum; i++){
724 fscanf(fp,"%i %s %lf %lf %lf %lf %lf",&j,&Species,
725 &Cxyz[i][1],
726 &Cxyz[i][2],
727 &Cxyz[i][3],
728 &Cdens[i][0],
729 &Cdens[i][1]);
730
731 C_WhatSpecies[i] = Species2int(Species);
732
733 if (i!=j){
734 printf("Error of sequential number %i in <Atoms.SpeciesAndCoordinates\n",j);
735 po++;
736 }
737
738 if (print_std==1){
739 printf("<Input_std> ct_AN=%2d C_WhatSpecies=%2d\n",i,C_WhatSpecies[i]);
740 }
741 }
742
743 if (!input_last("Atoms.SpeciesAndCoordinates>")) {
744 /* format error */
745 printf("Format error for Atoms.SpeciesAndCoordinates\n");
746 po++;
747 }
748 }
749
750 /****************************************************
751 coordinates of left
752 ****************************************************/
753
754 if (fp=input_find("<LeftLeadAtoms.SpeciesAndCoordinates") ) {
755 for (i=1; i<=Latomnum; i++){
756 fscanf(fp,"%i %s %lf %lf %lf %lf %lf",&j,&Species,
757 &Lxyz[i][1],
758 &Lxyz[i][2],
759 &Lxyz[i][3],
760 &Ldens[i][0],
761 &Ldens[i][1]);
762
763 L_WhatSpecies[i] = Species2int(Species);
764
765 if (i!=j){
766 printf("Error of sequential number %i in <LeftLeadAtoms.SpeciesAndCoordinates\n",j);
767 po++;
768 }
769
770 if (print_std==1){
771 printf("<Input_std> L_AN=%2d L_WhatSpecies=%2d\n",
772 i,L_WhatSpecies[i]);
773 }
774
775 }
776
777 if (!input_last("LeftLeadAtoms.SpeciesAndCoordinates>")) {
778 /* format error */
779 printf("Format error for LeftLeadAtoms.SpeciesAndCoordinates\n");
780 po++;
781 }
782 }
783
784 /****************************************************
785 coordinates of right
786 ****************************************************/
787
788 if (fp=input_find("<RightLeadAtoms.SpeciesAndCoordinates") ) {
789 for (i=1; i<=Ratomnum; i++){
790 fscanf(fp,"%i %s %lf %lf %lf %lf %lf",&j,&Species,
791 &Rxyz[i][1],
792 &Rxyz[i][2],
793 &Rxyz[i][3],
794 &Rdens[i][0],
795 &Rdens[i][1]);
796
797 R_WhatSpecies[i] = Species2int(Species);
798
799 if (i!=j){
800 printf("Error of sequential number %i in <RightLeadAtoms.SpeciesAndCoordinates\n",j);
801 po++;
802 }
803
804 if (print_std==1){
805 printf("<Input_std> R_AN=%2d R_WhatSpecies=%2d\n",
806 i,R_WhatSpecies[i]);
807 }
808
809 }
810
811 if (!input_last("RightLeadAtoms.SpeciesAndCoordinates>")) {
812 /* format error */
813 printf("Format error for RightLeadAtoms.SpeciesAndCoordinates\n");
814 po++;
815 }
816 }
817
818 if (coordinates_unit==0){
819 for (i=1; i<=Catomnum; i++){
820 Cxyz[i][1] = Cxyz[i][1]/BohrR;
821 Cxyz[i][2] = Cxyz[i][2]/BohrR;
822 Cxyz[i][3] = Cxyz[i][3]/BohrR;
823 }
824
825 for (i=1; i<=Latomnum; i++){
826 Lxyz[i][1] = Lxyz[i][1]/BohrR;
827 Lxyz[i][2] = Lxyz[i][2]/BohrR;
828 Lxyz[i][3] = Lxyz[i][3]/BohrR;
829 }
830
831 for (i=1; i<=Latomnum; i++){
832 Rxyz[i][1] = Rxyz[i][1]/BohrR;
833 Rxyz[i][2] = Rxyz[i][2]/BohrR;
834 Rxyz[i][3] = Rxyz[i][3]/BohrR;
835 }
836 }
837
838 /****************************************************
839 Unit cell
840 ****************************************************/
841
842 s_vec[0]="Ang"; s_vec[1]="AU";
843 i_vec[1]=0; i_vec[1]=1;
844 input_string2int("Atoms.UnitVectors.Unit",&unitvector_unit,2,s_vec,i_vec);
845
846 /* center */
847 if (fp=input_find("<Atoms.Unitvectors")) {
848 for (i=1; i<=3; i++){
849 fscanf(fp,"%lf %lf %lf",&tv[i][1],&tv[i][2],&tv[i][3]);
850 }
851 if ( ! input_last("Atoms.Unitvectors>") ) {
852 /* format error */
853 printf("Format error for Atoms.Unitvectors\n");
854 po++;
855 }
856 }
857
858 /* left */
859 if (fp=input_find("<LeftLeadAtoms.Unitvectors")) {
860 for (i=1; i<=3; i++){
861 fscanf(fp,"%lf %lf %lf",&Left_tv[i][1],&Left_tv[i][2],&Left_tv[i][3]);
862 }
863 if ( ! input_last("LeftLeadAtoms.Unitvectors>") ) {
864 /* format error */
865 printf("Format error for LeftLeadAtoms.Unitvectors\n");
866 po++;
867 }
868 }
869
870 /* right */
871 if (fp=input_find("<RightLeadAtoms.Unitvectors")) {
872 for (i=1; i<=3; i++){
873 fscanf(fp,"%lf %lf %lf",&Right_tv[i][1],&Right_tv[i][2],&Right_tv[i][3]);
874 }
875 if ( ! input_last("RightLeadAtoms.Unitvectors>") ) {
876 /* format error */
877 printf("Format error for RightLeadAtoms.Unitvectors\n");
878 po++;
879 }
880 }
881
882 if (unitvector_unit==0){
883 for (i=1; i<=3; i++){
884 tv[i][1] = tv[i][1]/BohrR;
885 tv[i][2] = tv[i][2]/BohrR;
886 tv[i][3] = tv[i][3]/BohrR;
887 }
888
889 for (i=1; i<=3; i++){
890 Left_tv[i][1] = Left_tv[i][1]/BohrR;
891 Left_tv[i][2] = Left_tv[i][2]/BohrR;
892 Left_tv[i][3] = Left_tv[i][3]/BohrR;
893 }
894
895 for (i=1; i<=3; i++){
896 Right_tv[i][1] = Right_tv[i][1]/BohrR;
897 Right_tv[i][2] = Right_tv[i][2]/BohrR;
898 Right_tv[i][3] = Right_tv[i][3]/BohrR;
899 }
900 }
901
902 /****************************************************
903 buffer size
904 ****************************************************/
905
906 input_int("Buffer.Size.LeftLead", &Lbuffer, 1);
907 input_int("Buffer.Size.RightLead", &Rbuffer, 1);
908
909 } /* else{ */
910
911
912 /****************************************************
913 input_close
914 ****************************************************/
915
916 input_close();
917
918 if (po>0 || input_errorCount()>0) {
919 printf("errors in the inputfile\n");
920 exit(1);
921 }
922
923 /****************************************************
924 print out to std
925 ****************************************************/
926
927 printf("<Input_NEGF> Your input file was normally read.\n");
928 printf("<Input_NEGF> The system includes %i species and %i atoms.\n",
929 SpeciesNum,atomnum);
930
931 /****************************************************
932 find rcut of pao
933 ****************************************************/
934
935 for (spe=0; spe<SpeciesNum; spe++){
936 fnjoint2(DirPAO,SpeBasisName[spe],ExtPAO,FN_PAO);
937 if ((fp = fopen(FN_PAO,"r")) != NULL){
938 input_open(FN_PAO);
939 input_double("radial.cutoff.pao",&Spe_Atom_Cut1[spe],(double)0.0);
940 input_close();
941 fclose(fp);
942 printf("<Input_NEGF> PAOs of species %s were normally found.\n",SpeName[spe]);
943 printf("<Input_NEGF> rcut of PAOs of species %s = %5.3f (Bohr)\n",
944 SpeName[spe],Spe_Atom_Cut1[spe]);
945 }
946 else{
947 printf("Could not find %s\n",FN_PAO);
948 exit(1);
949 }
950 }
951
952 }
953
SpeciesString2int(int p)954 void SpeciesString2int(int p)
955 {
956 static int i,l,n,po;
957 char c,cstr[asize10];
958
959 /* Get basis name */
960
961 i = 0;
962 po = 0;
963 while((c=SpeBasis[p][i])!='\0'){
964 if (c=='-') po = 1;
965 if (po==0) SpeBasisName[p][i] = SpeBasis[p][i];
966 i++;
967 }
968
969 if (print_std==1){
970 printf("<Input_std> SpeBasisName=%s\n",SpeBasisName[p]);
971 }
972
973 /* Get basis type */
974
975 for (l=0; l<5; l++){
976 Spe_Num_Basis[p][l] = 0;
977 }
978
979 i = 0;
980 po = 0;
981
982 while((c=SpeBasis[p][i])!='\0'){
983 if (po==1){
984 if (c=='s'){ l=0; n=0; }
985 else if (c=='p'){ l=1; n=0; }
986 else if (c=='d'){ l=2; n=0; }
987 else if (c=='f'){ l=3; n=0; }
988 else{
989
990 if (n==0){
991 cstr[0] = c;
992 cstr[1] = '\0';
993 Spe_Num_Basis[p][l] = atoi(cstr);
994 Spe_Num_CBasis[p][l] = atoi(cstr);
995 n++;
996 }
997 else if (n==1){
998 cstr[0] = c;
999 cstr[1] = '\0';
1000 Spe_Num_CBasis[p][l] = atoi(cstr);
1001 if (Spe_Num_Basis[p][l]<Spe_Num_CBasis[p][l]){
1002 printf("# of contracted orbitals are larger than # of primitive oribitals\n");
1003 exit(1);
1004 }
1005
1006 n++;
1007 }
1008 else {
1009 printf("Format error in Definition of Atomic Species\n");
1010 exit(1);
1011 }
1012 }
1013 }
1014
1015 if (SpeBasis[p][i]=='-') po = 1;
1016 i++;
1017 }
1018
1019 for (l=0; l<5; l++){
1020 if (Spe_Num_Basis[p][l]!=0) Spe_MaxL_Basis[p] = l;
1021
1022 if (print_std==1){
1023 printf("<Input_std> p=%2d l=%2d %2d %2d\n",
1024 p,l,Spe_Num_Basis[p][l],Spe_Num_CBasis[p][l]);
1025 }
1026 }
1027 }
1028
Species2int(char Species[asize10])1029 int Species2int(char Species[asize10])
1030 {
1031 static int i,po;
1032
1033 i = 0;
1034 po = 0;
1035 while (i<SpeciesNum && po==0){
1036 if (SEQ(Species,SpeName[i])==1){
1037 po = 1;
1038 }
1039 if (po==0) i++;
1040 };
1041
1042 if (po==1) return i;
1043 else {
1044 printf("%s is an invalid species name in Atoms.SpeciesAndCoordinates\n",
1045 Species);
1046 printf("Please check your input file\n");
1047 exit(1);
1048 }
1049 }
1050
1051
SEQ(char str1[asize10],char str2[asize10])1052 static int SEQ(char str1[asize10],
1053 char str2[asize10])
1054 {
1055
1056 int i,result,l1,l2;
1057
1058 l1 = strlen(str1);
1059 l2 = strlen(str2);
1060
1061 result = 1;
1062 if (l1 == l2){
1063 for (i=0; i<=l1-1; i++){
1064 if (str1[i]!=str2[i]) result = 0;
1065 }
1066 }
1067 else
1068 result = 0;
1069
1070 return result;
1071 }
1072
1073
fnjoint2(char name1[asize10],char name2[asize10],char name3[asize10],char name4[asize10])1074 void fnjoint2(char name1[asize10], char name2[asize10],
1075 char name3[asize10], char name4[asize10])
1076 {
1077 char *f1 = name1,
1078 *f2 = name2,
1079 *f3 = name3,
1080 *f4 = name4;
1081
1082 while(*f1)
1083 {
1084 *f4 = *f1;
1085 *f1++;
1086 *f4++;
1087 }
1088 while(*f2)
1089 {
1090 *f4 = *f2;
1091 *f2++;
1092 *f4++;
1093 }
1094 while(*f3)
1095 {
1096 *f4 = *f3;
1097 *f3++;
1098 *f4++;
1099 }
1100 *f4 = *f3;
1101 }
1102
FreeArrays()1103 void FreeArrays()
1104 {
1105 static int i,ct_AN,spe;
1106
1107 free(C_WhatSpecies);
1108 free(L_WhatSpecies);
1109 free(R_WhatSpecies);
1110
1111 for (i=0; i<SpeciesNum; i++){
1112 free(Spe_Num_Basis[i]);
1113 }
1114 free(Spe_Num_Basis);
1115
1116 for (i=0; i<SpeciesNum; i++){
1117 free(Spe_Num_CBasis[i]);
1118 }
1119 free(Spe_Num_CBasis);
1120
1121 free(Spe_MaxL_Basis);
1122
1123 for (ct_AN=0; ct_AN<=Catomnum; ct_AN++){
1124 free(Cxyz[ct_AN]);
1125 }
1126 free(Cxyz);
1127
1128 for (ct_AN=0; ct_AN<=Latomnum; ct_AN++){
1129 free(Lxyz[ct_AN]);
1130 }
1131 free(Lxyz);
1132
1133 for (ct_AN=0; ct_AN<=Ratomnum; ct_AN++){
1134 free(Rxyz[ct_AN]);
1135 }
1136 free(Rxyz);
1137
1138 for (ct_AN=0; ct_AN<=Catomnum; ct_AN++){
1139 free(Cdens[ct_AN]);
1140 }
1141 free(Cdens);
1142
1143 for (ct_AN=0; ct_AN<=Latomnum; ct_AN++){
1144 free(Ldens[ct_AN]);
1145 }
1146 free(Ldens);
1147
1148 for (ct_AN=0; ct_AN<=Ratomnum; ct_AN++){
1149 free(Rdens[ct_AN]);
1150 }
1151 free(Rdens);
1152
1153 free(Spe_Atom_Cut1);
1154
1155 for (i=0; i<SpeciesNum; i++){
1156 free(SpeName[i]);
1157 }
1158 free(SpeName);
1159
1160 for (i=0; i<SpeciesNum; i++){
1161 free(SpeBasis[i]);
1162 }
1163 free(SpeBasis);
1164
1165 for (i=0; i<SpeciesNum; i++){
1166 free(SpeBasisName[i]);
1167 }
1168 free(SpeBasisName);
1169
1170 for (i=0; i<SpeciesNum; i++){
1171 free(SpeVPS[i]);
1172 }
1173 free(SpeVPS);
1174
1175 }
1176
1177
1178
1179
1180