1 /**********************************************************************
2 TRAN_Main_Analysis_NC.c:
3
4 TRAN_Main_Analysis_NC.c is a subroutine to analyze transport properties
5 such as electronic transmission, current, eigen channel, and current
6 distribution in real space based on the non-colliear density functional
7 theories and the non-equilibrium Green's function method.
8
9 Log of TRAN_Main_Analysis_NC.c:
10
11 15/Dec./2011 released by Y. Xiao
12 2/June/2015 integrated in OpenMX by T.Ozaki
13
14 ***********************************************************************/
15
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <math.h>
19 #include <time.h>
20 #include <string.h>
21 #include "Inputtools.h"
22 #include "mpi.h"
23 #include "tran_prototypes.h"
24
25 #define eV2Hartree 27.2113845
26 #define kB 0.00008617251324000000 /* eV/K */
27 #define PI 3.1415926535897932384626
28
29 #define Host_ID 0
30 #define PrintLevel 0
31
32 /*
33 #define SCC_ref(i,j) ( ((j)-1)*NUM_c + (i)-1 )
34 #define SCL_ref(i,j) ( ((j)-1)*NUM_c + (i)-1 )
35 #define SCR_ref(i,j) ( ((j)-1)*NUM_c + (i)-1 )
36 #define S00l_ref(i,j) ( ((j)-1)*NUM_e[0]+(i)-1 )
37 #define S00r_ref(i,j) ( ((j)-1)*NUM_e[1]+(i)-1 )
38 */
39
40 static int SpinP_switch,SpinP_switch2;
41 static int NUM_c, NUM_e[2];
42 static double E_Temp;
43
44 /* the center region */
45
46 static double ChemP;
47 static double *****OLP;
48 static double *****H;
49 static double *****iHNL;
50 static double *****H0;
51 static int atomnum;
52 static int SpeciesNum;
53 static int *WhatSpecies;
54 static int *Spe_Total_CNO;
55 static int *FNAN;
56 static int **natn;
57 static int **ncn;
58 static int **atv_ijk;
59 static int Max_FSNAN;
60 static double ScaleSize;
61 static int TCpyCell;
62 static int *TRAN_region;
63 static int *TRAN_Original_Id;
64
65 /* the leads region */
66
67 static double ChemP_e[2];
68 static double *****OLP_e[2];
69 static double *****H_e[2];
70 static double *****iHNL_e[2];
71 static int atomnum_e[2];
72 static int SpeciesNum_e[2];
73 static int *WhatSpecies_e[2];
74 static int *Spe_Total_CNO_e[2];
75 static int *FNAN_e[2];
76 static int **natn_e[2];
77 static int **ncn_e[2];
78 static int **atv_ijk_e[2];
79 static int Max_FSNAN_e[2];
80 static double ScaleSize_e[2];
81 static int TCpyCell_e[2];
82
83 /* k-dependent matrices */
84
85 static dcomplex *S00_e[2];
86 static dcomplex *S01_e[2];
87 static dcomplex *H00_e[2];
88 static dcomplex *H01_e[2];
89
90 static dcomplex *SCC;
91 static dcomplex *SCL;
92 static dcomplex *SCR;
93 static dcomplex *HCC;
94 static dcomplex *HCL;
95 static dcomplex *HCR;
96
97 static dcomplex ****tran_transmission;
98 static dcomplex **tran_transmission_iv;
99
100 static int tran_surfgreen_iteration_max;
101 static double tran_surfgreen_eps;
102
103 static int Tran_current_num_step;
104 static double Tran_current_energy_step,Tran_current_cutoff;
105 static double Tran_current_lower_bound,Tran_current_im_energy;
106
107 static int SCF_tran_bias_apply;
108 static int tran_transmission_on;
109 static double tran_transmission_energyrange[3];
110 static int tran_transmission_energydiv;
111 static int tran_interpolate;
112 static char interpolate_filename1[YOUSO10];
113 static char interpolate_filename2[YOUSO10];
114 static double interpolate_c1,interpolate_c2;
115 static int TRAN_TKspace_grid2,TRAN_TKspace_grid3;
116 static double ***current;
117 static int Order_Lead_Side[2];
118
119 static char filepath[100];
120 static char filename[100];
121
122 /* S MitsuakiKAWAMURA */
123
124 /* Eigenchannel analysis */
125 static int TRAN_Channel;
126 static int TRAN_CurrentDensity;
127 static int TRAN_OffDiagonalCurrent;
128 static int TRAN_Channel_Num;
129 static int TRAN_Channel_Nkpoint;
130 static int TRAN_Channel_Nenergy;
131 static double **TRAN_Channel_kpoint;
132 static double *TRAN_Channel_energy;
133
134 /* Current Density */
135 static dcomplex *VCC;
136 static double ****JLocSym, ***JLocASym, ***RhoNL, ****Jmat;
137
138 #define BUFSIZE 200
139 /* E MitsuakiKAWAMURA */
140 void Make_Comm_Worlds(
141 MPI_Comm MPI_Curret_Comm_WD,
142 int myid0,
143 int numprocs0,
144 int Num_Comm_World,
145 int *myworld1,
146 MPI_Comm *MPI_CommWD, /* size: Num_Comm_World */
147 int *NPROCS1_ID, /* size: numprocs0 */
148 int *Comm_World1, /* size: numprocs0 */
149 int *NPROCS1_WD, /* size: Num_Comm_World */
150 int *Comm_World_StartID /* size: Num_Comm_World */
151 );
152
153
154 static void MTRAN_Read_Tran_HS(MPI_Comm comm1, char *filepath, char *filename, char *ext ) ;
155
156 static void MTRAN_Transmission(
157 MPI_Comm comm1,
158 int numprocs,
159 int myid,
160 int SpinP_switch,
161 double ChemP_e[2],
162 int NUM_c,
163 int NUM_e[2],
164 dcomplex *H00_e[2],
165 dcomplex *S00_e[2],
166 dcomplex *H01_e[2],
167 dcomplex *S01_e[2],
168 dcomplex *HCC,
169 dcomplex *HCL,
170 dcomplex *HCR,
171 dcomplex *SCC,
172 dcomplex *SCL,
173 dcomplex *SCR,
174 double tran_surfgreen_iteration_max,
175 double tran_surfgreen_eps,
176 double tran_transmission_energyrange[3],
177 int tran_transmission_energydiv,
178 dcomplex **tran_transmission);
179
180 /* S MitsuakiKAWAMURA*/
181
182 static void MTRAN_Current(
183 MPI_Comm comm1,
184 int numprocs,
185 int myid,
186 int SpinP_switch,
187 double ChemP_e[2],
188 double E_Temp,
189 int NUM_c,
190 int NUM_e[2],
191 dcomplex *H00_e[2],
192 dcomplex *S00_e[2],
193 dcomplex *H01_e[2],
194 dcomplex *S01_e[2],
195 dcomplex *HCC,
196 dcomplex *HCL,
197 dcomplex *HCR,
198 dcomplex *SCC,
199 dcomplex *SCL,
200 dcomplex *SCR,
201 double tran_surfgreen_iteration_max,
202 double tran_surfgreen_eps,
203 double *current,
204 double k2,
205 double k3);
206
207 void TRAN_Calc_Sinv(
208 int NUM_c,
209 dcomplex *SCC,
210 dcomplex *Sinv);
211
212 void TRAN_Calc_CurrentDensity_NC(
213 int NUM_c,
214 dcomplex *GC,
215 dcomplex *SigmaL,
216 dcomplex *SigmaR,
217 dcomplex *VCC,
218 dcomplex *Sinv,
219 double *kvec,
220 double fL,
221 double fR,
222 double Tran_current_energy_step,
223 double ****JLocSym,
224 double ***JLocAsym,
225 double ***RhoNLNL,
226 double ****Jmat
227 );
228
229 void TRAN_CDen_Main(
230 int NUM_c,
231 int *MP,
232 double ****JLocSym,
233 double ***JLocASym,
234 double ***Rho,
235 double ****Jmat,
236 dcomplex *SCC,
237 int TRAN_OffDiagonalCurrent);
238
239 void MTRAN_EigenChannel_NC(
240 MPI_Comm comm1,
241 int numprocs,
242 int myid,
243 int myid0,
244 int SpinP_switch,
245 double ChemP_e[2],
246 int NUM_c,
247 int NUM_e[2],
248 dcomplex *H00_e[2],
249 dcomplex *S00_e[2],
250 dcomplex *H01_e[2],
251 dcomplex *S01_e[2],
252 dcomplex *HCC,
253 dcomplex *HCL,
254 dcomplex *HCR,
255 dcomplex *SCC,
256 dcomplex *SCL,
257 dcomplex *SCR,
258 double tran_surfgreen_iteration_max,
259 double tran_surfgreen_eps,
260 double tran_transmission_energyrange[3],
261 int TRAN_Channel_Nenergy,
262 double *TRAN_Channel_energy,
263 int TRAN_Channel_Num,
264 int kloop,
265 double *TRAN_Channel_kpoint,
266 dcomplex ***EChannel,
267 double **eigentrans,
268 double **eigentrans_sum
269 ); /* void MTRAN_EigenChannel_NC */
270
271 void TRAN_Output_eigentrans_sum(
272 int TRAN_Channel_Nkpoint,
273 int TRAN_Channel_Nenergy,
274 double ***eigentrans_sum);
275
276 void TRAN_Output_ChannelCube(
277 int kloop,
278 int iw,
279 int ispin,
280 int orbit,
281 int NUM_c,
282 double *TRAN_Channel_kpoint,
283 dcomplex *EChannel,
284 int *MP,
285 double eigentrans,
286 double TRAN_Channel_energy
287 ); /* void TRAN_Output_ChannelCube */
288
289 static void MTRAN_Free_All_NC();
290 /* E MitsuakiKAWAMURA */
291
292 static void MTRAN_Output_Transmission(
293 MPI_Comm comm1,
294 int ID,
295 char *fname,
296 double k2,
297 double k3,
298 int SpinP_switch,
299 int tran_transmission_energydiv,
300 double tran_transmission_energyrange[3],
301 dcomplex **tran_transmission
302 );
303
304 static void MTRAN_Output_Current(
305 MPI_Comm comm1,
306 char *fname,
307 int TRAN_TKspace_grid2,
308 int TRAN_TKspace_grid3,
309 int SpinP_switch,
310 double ***current
311 );
312
313 static void MTRAN_Output_Conductance(
314 MPI_Comm comm1,
315 char *fname,
316 int *T_k_ID,
317 int T_knum,
318 int TRAN_TKspace_grid2,
319 int TRAN_TKspace_grid3,
320 int *T_IGrids2,
321 int *T_IGrids3,
322 double *T_KGrids2,
323 double *T_KGrids3,
324 int SpinP_switch,
325 int tran_transmission_energydiv,
326 double tran_transmission_energyrange[3],
327 dcomplex ****tran_transmission
328 );
329
330 static void MTRAN_Input(
331 MPI_Comm comm1,
332 int argc,
333 char *fname,
334 double ChemP_e[2],
335 int *TRAN_TKspace_grid2,
336 int *TRAN_TKspace_grid3,
337 int *SpinP_switch,
338 double *E_Temp,
339 int *tran_surfgreen_iteration_max,
340 double *tran_surfgreen_eps,
341 int *tran_transmission_on,
342 double tran_transmission_energyrange[3],
343 int *tran_transmission_energydiv,
344 dcomplex ****(*tran_transmission)
345 );
346
347
348 static void MTRAN_Input_Sys(
349 int argc,char *file, char *filepath, char *filename,
350 int *tran_interpolate,
351 char *interpolate_filename1,
352 char *interpolate_filename2);
353
354 static void MTRAN_Set_MP(
355 int job,
356 int anum, int *WhatSpecies, int *Spe_Total_CNO,
357 int *NUM, /* output */
358 int *MP /* output */
359 );
360
361 static void MTRAN_Set_SurfOverlap(
362 char *position,
363 double k2,
364 double k3,
365 int SpinP_switch,
366 int atomnum_e[2],
367 double *****OLP_e[2],
368 double *****H_e[2],
369 double *****iHNL_e[2],
370 int SpeciesNum_e[2],
371 int *WhatSpecies_e[2],
372 int *Spe_Total_CNO_e[2],
373 int *FNAN_e[2],
374 int **natn_e[2],
375 int **ncn_e[2],
376 int **atv_ijk_e[2],
377 dcomplex *S00_e[2],
378 dcomplex *S01_e[2],
379 dcomplex *H00_e[2],
380 dcomplex *H01_e[2]
381 );
382
383 static void MTRAN_Set_CentOverlap(
384 int job,
385 int SpinP_switch,
386 double k2,
387 double k3,
388 int NUM_c,
389 int NUM_e[2],
390 double *****H,
391 double *****iHNL,
392 double *****OLP,
393 int atomnum,
394 int atomnum_e[2],
395 int *WhatSpecies,
396 int *WhatSpecies_e[2],
397 int *Spe_Total_CNO,
398 int *Spe_Total_CNO_e[2],
399 int *FNAN,
400 int **natn,
401 int **ncn,
402 int **atv_ijk,
403 int *TRAN_region,
404 int *TRAN_Original_Id
405 );
406
407 static void MTRAN_Allocate_HS(
408 int NUM_c,
409 int NUM_e[2],
410 int SpinP_switch);
411
412
413
414
415
TRAN_Main_Analysis_NC(MPI_Comm comm1,int argc,char * argv[],int Matomnum,int * M2G,int * GridN_Atom,int ** GridListAtom,int ** CellListAtom,Type_Orbs_Grid *** Orbs_Grid,int TNumGrid)416 void TRAN_Main_Analysis_NC( MPI_Comm comm1,
417 int argc, char *argv[],
418 int Matomnum, int *M2G,
419 int *GridN_Atom,
420 int **GridListAtom,
421 int **CellListAtom,
422 Type_Orbs_Grid ***Orbs_Grid,
423 int TNumGrid )
424 {
425 int i,j,i2,i3,ii2,ii3,k,iw;
426 int Gc_AN,h_AN,Gh_AN;
427 int iside,tno0,tno1,Cwan,Hwan;
428 int kloop0,kloop,S_knum,E_knum,k_op;
429 /* S MitsuakiKAWAMURA */
430 int myworld1, myworld2, parallel_mode, num_kloop0;
431 int myid0, numprocs0, myid1, myid2, numprocs1, numprocs2;
432 /* E MitsuakiKAWAMURA */
433 int ID, myid_tmp, numprocs_tmp, T_knum;
434 int **op_flag,*T_op_flag,*T_k_ID;
435 double *T_KGrids2,*T_KGrids3;
436 int *T_IGrids2,*T_IGrids3;
437 /* S MitsuakiKAWAMURA */
438 int Num_Comm_World1, Num_Comm_World2;
439 int *NPROCS_ID1, *NPROCS_ID2;
440 int *Comm_World1, *Comm_World2;
441 int *NPROCS_WD1, *NPROCS_WD2;
442 int *Comm_World_StartID1, *Comm_World_StartID2;
443 MPI_Comm *MPI_CommWD1, *MPI_CommWD2;
444 double ***eigentrans;
445 double ***eigentrans_sum;
446 dcomplex ****EChannel;
447 int *MP;
448 int NUM_cs;
449 /* E MitsuakiKAWAMURA */
450 double k2,k3,tmp;
451 MPI_Comm comm_tmp;
452 char fnameout[100];
453
454 MPI_Comm_size(comm1,&numprocs0);
455 MPI_Comm_rank(comm1,&myid0);
456
457 /**********************************************
458 show something
459 **********************************************/
460
461 if (myid0==Host_ID){
462 printf("\n*******************************************************\n");
463 printf("*******************************************************\n");
464 printf(" Welcome to TranMain_NC \n");
465 printf(" This is a post-processing code of OpenMX to calculate \n");
466 printf(" electronic transmission and current. \n");
467 printf(" Copyright (C), 2002-2013, H.Kino and T.Ozaki \n");
468 printf(" TranMain_NC comes with ABSOLUTELY NO WARRANTY. \n");
469 printf(" This is free software, and you are welcome to \n");
470 printf(" redistribute it under the constitution of the GNU-GPL.\n");
471 printf("*******************************************************\n");
472 printf("*******************************************************\n\n");
473 }
474
475 /**********************************************
476 read system
477 **********************************************/
478
479 MTRAN_Input_Sys(argc,argv[1],filepath,filename,&tran_interpolate,interpolate_filename1,interpolate_filename2);
480
481 /**********************************************
482 read tranb file
483 **********************************************/
484
485 MTRAN_Read_Tran_HS(comm1, filepath, filename, "tranb" );
486
487 /**********************************************
488 read input file
489 **********************************************/
490
491 MTRAN_Input(comm1,
492 argc,
493 argv[1],
494 ChemP_e,
495 /* output */
496 &TRAN_TKspace_grid2,
497 &TRAN_TKspace_grid3,
498 &SpinP_switch2,
499 &E_Temp,
500 &tran_surfgreen_iteration_max,
501 &tran_surfgreen_eps,
502 &tran_transmission_on,
503 tran_transmission_energyrange,
504 &tran_transmission_energydiv,
505 &tran_transmission );
506
507 if (SpinP_switch!=SpinP_switch2) {
508 printf("SpinP_switch conflicts\n");fflush(stdout);
509 printf("SpinP_switch=%d SpinP_switch2=%d\n", SpinP_switch,SpinP_switch2);fflush(stdout);
510 exit(0);
511 }
512
513 MTRAN_Allocate_HS(NUM_c,NUM_e,SpinP_switch);
514
515 /**********************************************
516 calculate transmission
517 **********************************************/
518
519 if (tran_transmission_on) {
520
521 /* allocation of arrays */
522
523 current = (double***)malloc(sizeof(double**)*TRAN_TKspace_grid2);
524 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
525 current[i2] = (double**)malloc(sizeof(double*)*TRAN_TKspace_grid3);
526 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
527 current[i2][i3] = (double*)malloc(sizeof(double)*3);
528 }
529 }
530
531 op_flag = (int**)malloc(sizeof(int*)*TRAN_TKspace_grid2);
532 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
533 op_flag[i2] = (int*)malloc(sizeof(int)*TRAN_TKspace_grid3);
534 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
535 op_flag[i2][i3] = 1;
536 }
537 }
538
539 /***********************************
540 one-dimentionalize for MPI
541 ************************************/
542
543 T_knum = 0;
544 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
545 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
546 if (0<op_flag[i2][i3]) T_knum++;
547 }
548 }
549
550 T_KGrids2 = (double*)malloc(sizeof(double)*T_knum);
551 T_KGrids3 = (double*)malloc(sizeof(double)*T_knum);
552 T_IGrids2 = (int*)malloc(sizeof(int)*T_knum);
553 T_IGrids3 = (int*)malloc(sizeof(int)*T_knum);
554 T_op_flag = (int*)malloc(sizeof(int)*T_knum);
555 T_k_ID = (int*)malloc(sizeof(int)*T_knum);
556
557 T_knum = 0;
558
559 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
560
561 k2 = -0.5 + (2.0*(double)i2+1.0)/(2.0*(double)TRAN_TKspace_grid2) + Shift_K_Point;
562
563 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
564
565 k3 = -0.5 + (2.0*(double)i3+1.0)/(2.0*(double)TRAN_TKspace_grid3) - Shift_K_Point;
566
567 if (0<op_flag[i2][i3]){
568
569 T_KGrids2[T_knum] = k2;
570 T_KGrids3[T_knum] = k3;
571 T_IGrids2[T_knum] = i2;
572 T_IGrids3[T_knum] = i3;
573 T_op_flag[T_knum] = op_flag[i2][i3];
574
575 T_knum++;
576 }
577 }
578 }
579
580 /***************************************************
581 allocate calculations of k-points into processors
582 ***************************************************/
583
584 if (numprocs0<T_knum){
585
586 /* set parallel_mode */
587 parallel_mode = 0;
588
589 /* allocation of kloop to ID */
590
591 for (ID=0; ID<numprocs0; ID++){
592
593 tmp = (double)T_knum/(double)numprocs0;
594 S_knum = (int)((double)ID*(tmp+1.0e-12));
595 E_knum = (int)((double)(ID+1)*(tmp+1.0e-12)) - 1;
596 if (ID==(numprocs0-1)) E_knum = T_knum - 1;
597 if (E_knum<0) E_knum = 0;
598
599 for (k=S_knum; k<=E_knum; k++){
600 /* ID in the first level world */
601 T_k_ID[k] = ID;
602 }
603 }
604
605 /* find own informations */
606
607 tmp = (double)T_knum/(double)numprocs0;
608 S_knum = (int)((double)myid0*(tmp+1.0e-12));
609 E_knum = (int)((double)(myid0+1)*(tmp+1.0e-12)) - 1;
610 if (myid0==(numprocs0-1)) E_knum = T_knum - 1;
611 if (E_knum<0) E_knum = 0;
612
613 num_kloop0 = E_knum - S_knum + 1;
614
615 }
616
617 else {
618
619 /* set parallel_mode */
620 parallel_mode = 1;
621 num_kloop0 = 1;
622
623 Num_Comm_World1 = T_knum;
624
625 NPROCS_ID1 = (int*)malloc(sizeof(int)*numprocs0);
626 Comm_World1 = (int*)malloc(sizeof(int)*numprocs0);
627 NPROCS_WD1 = (int*)malloc(sizeof(int)*Num_Comm_World1);
628 Comm_World_StartID1 = (int*)malloc(sizeof(int)*Num_Comm_World1);
629 MPI_CommWD1 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World1);
630
631 Make_Comm_Worlds(comm1, myid0, numprocs0, Num_Comm_World1, &myworld1, MPI_CommWD1,
632 NPROCS_ID1, Comm_World1, NPROCS_WD1, Comm_World_StartID1);
633
634 MPI_Comm_size(MPI_CommWD1[myworld1],&numprocs1);
635 MPI_Comm_rank(MPI_CommWD1[myworld1],&myid1);
636
637 S_knum = myworld1;
638
639 /* allocate k-points into processors */
640
641 for (k=0; k<T_knum; k++){
642 /* ID in the first level world */
643 T_k_ID[k] = Comm_World_StartID1[k];
644 }
645
646 }
647
648 /***********************************************************
649 start "kloop0"
650 ***********************************************************/
651
652 /* S MitsuakiKAWAMURA 2*/
653 NUM_cs = NUM_c / 2;
654
655 if (fabs(ChemP_e[0] - ChemP_e[1]) < 0.000001) TRAN_CurrentDensity = 0;
656
657 JLocSym = (double****)malloc(sizeof(double***) * (SpinP_switch + 1));
658 Jmat = (double****)malloc(sizeof(double***) * (SpinP_switch + 1));
659 JLocASym = (double***)malloc(sizeof(double**) * (SpinP_switch + 1));
660 RhoNL = (double***)malloc(sizeof(double**) * (SpinP_switch + 1));
661 for (k = 0; k < SpinP_switch + 1; k++) {
662 JLocSym[k] = (double***)malloc(sizeof(double**) * 3);
663 for (iside = 0; iside < 3; iside++) {
664 JLocSym[k][iside] = (double**)malloc(sizeof(double*) * NUM_cs);
665 for (i = 0; i < NUM_cs; i++) {
666 JLocSym[k][iside][i] = (double*)malloc(sizeof(double) * NUM_cs);
667 for (j = 0; j < NUM_cs; j++) JLocSym[k][iside][i][j] = 0.0;
668 } /* for (i = 0; i < NUM_cs; i++) */
669 } /* for (iside = 0; iside < 3; iside++) */
670
671 Jmat[k] = (double***)malloc(sizeof(double**) * 2);
672 for (iside = 0; iside < 2; iside++) {
673 Jmat[k][iside] = (double**)malloc(sizeof(double*) * NUM_cs);
674 for (i = 0; i < NUM_cs; i++) {
675 Jmat[k][iside][i] = (double*)malloc(sizeof(double) * NUM_cs);
676 for (j = 0; j < NUM_cs; j++) Jmat[k][iside][i][j] = 0.0;
677 } /* for (i = 0; i < NUM_cs; i++) */
678 } /* for (iside = 0; iside < 2; iside++) */
679
680 JLocASym[k] = (double**)malloc(sizeof(double*) * NUM_cs);
681 RhoNL[k] = (double**)malloc(sizeof(double*) * NUM_cs);
682 for (i = 0; i < NUM_cs; i++) {
683 JLocASym[k][i] = (double*)malloc(sizeof(double) * NUM_cs);
684 RhoNL[k][i] = (double*)malloc(sizeof(double) * NUM_cs);
685 for (j = 0; j < NUM_cs; j++) {
686 JLocASym[k][i][j] = 0.0;
687 RhoNL[k][i][j] = 0.0;
688 } /* for (j = 0; j < NUM_cs; j++) */
689 } /* for (i = 0; i < NUM_cs; i++) */
690 } /* for (k = 0; k < SpinP_switch + 1; k++) */
691 /* E MitsuakiKAWAMURA 2*/
692
693 if (myid0==Host_ID) printf("\n calculating...\n\n"); fflush(stdout);
694 MPI_Barrier(comm1);
695
696 for (kloop0=0; kloop0<num_kloop0; kloop0++){
697
698 kloop = S_knum + kloop0;
699
700 k2 = T_KGrids2[kloop];
701 k3 = T_KGrids3[kloop];
702 i2 = T_IGrids2[kloop];
703 i3 = T_IGrids3[kloop];
704 k_op = T_op_flag[kloop];
705
706 printf(" myid0=%2d i2=%2d i3=%2d k2=%8.4f k3=%8.4f\n",myid0,i2,i3,k2,k3); fflush(stdout);
707
708 if (parallel_mode){
709 comm_tmp = MPI_CommWD1[myworld1];
710 numprocs_tmp = numprocs1;
711 myid_tmp = myid1;
712 }
713 else{
714 comm_tmp = comm1;
715 numprocs_tmp = 1;
716 myid_tmp = 0;
717 }
718
719 /* set Hamiltonian and overlap matrices of left and right leads */
720
721 MTRAN_Set_SurfOverlap( "left", k2,k3,SpinP_switch,atomnum_e,
722 OLP_e,H_e,iHNL_e,SpeciesNum_e,WhatSpecies_e,Spe_Total_CNO_e,
723 FNAN_e,natn_e,ncn_e,atv_ijk_e,S00_e,S01_e,H00_e,H01_e );
724
725 MTRAN_Set_SurfOverlap( "right", k2,k3,SpinP_switch,atomnum_e,
726 OLP_e,H_e,iHNL_e,SpeciesNum_e,WhatSpecies_e,Spe_Total_CNO_e,
727 FNAN_e,natn_e,ncn_e,atv_ijk_e,S00_e,S01_e,H00_e,H01_e );
728
729 /* set CC, CL and CR */
730
731 MTRAN_Set_CentOverlap(3,
732 SpinP_switch,
733 k2,
734 k3,
735 NUM_c,
736 NUM_e,
737 H,
738 iHNL,
739 OLP,
740 atomnum,
741 atomnum_e,
742 WhatSpecies,
743 WhatSpecies_e,
744 Spe_Total_CNO,
745 Spe_Total_CNO_e,
746 FNAN,
747 natn,
748 ncn,
749 atv_ijk,
750 TRAN_region,
751 TRAN_Original_Id
752 );
753
754 /* calculate transmission */
755
756 MTRAN_Transmission(comm_tmp,
757 numprocs_tmp,
758 myid_tmp,
759 SpinP_switch,
760 ChemP_e,
761 NUM_c,
762 NUM_e,
763 H00_e,
764 S00_e,
765 H01_e,
766 S01_e,
767 HCC,
768 HCL,
769 HCR,
770 SCC,
771 SCL,
772 SCR,
773 tran_surfgreen_iteration_max,
774 tran_surfgreen_eps,
775 tran_transmission_energyrange,
776 tran_transmission_energydiv,
777 /* output */
778 tran_transmission[i2][i3]
779 );
780
781 /* calculate current */
782
783 MTRAN_Current(comm_tmp,
784 numprocs_tmp,
785 myid_tmp,
786 SpinP_switch,
787 ChemP_e,
788 E_Temp,
789 NUM_c,
790 NUM_e,
791 H00_e,
792 S00_e,
793 H01_e,
794 S01_e,
795 HCC,
796 HCL,
797 HCR,
798 SCC,
799 SCL,
800 SCR,
801 tran_surfgreen_iteration_max,
802 tran_surfgreen_eps,
803 current[i2][i3],
804 k2, k3
805 );
806
807
808 } /* kloop0 */
809 } /* if (tran_transmission_on) */
810
811 /**********************************************
812 MPI: current
813 **********************************************/
814
815 for (k=0; k<T_knum; k++){
816
817 ID = T_k_ID[k];
818
819 i2 = T_IGrids2[k];
820 i3 = T_IGrids3[k];
821
822 MPI_Bcast(current[i2][i3], 2, MPI_DOUBLE, ID, comm1);
823 MPI_Barrier(comm1);
824
825 }
826
827 /**********************************************
828 output transmission
829 **********************************************/
830
831 if (tran_transmission_on) {
832
833 MPI_Barrier(comm1);
834 if (myid0==Host_ID) printf("\nTransmission: files\n\n");fflush(stdout);
835 MPI_Barrier(comm1);
836
837 for (k=0; k<T_knum; k++){
838
839 ID = T_k_ID[k];
840
841 k2 = T_KGrids2[k];
842 k3 = T_KGrids3[k];
843 i2 = T_IGrids2[k];
844 i3 = T_IGrids3[k];
845
846 sprintf(fnameout,"%s%s.tran%i_%i",filepath,filename,i2,i3);
847
848 MTRAN_Output_Transmission(
849 comm1,
850 ID,
851 fnameout,
852 k2,
853 k3,
854 SpinP_switch,
855 tran_transmission_energydiv,
856 tran_transmission_energyrange,
857 tran_transmission[i2][i3]
858 );
859
860 }
861
862 }
863
864 /**********************************************
865 output current
866 **********************************************/
867
868 if (tran_transmission_on) {
869
870 MPI_Barrier(comm1);
871 if (myid0==Host_ID) printf("\nCurrent: file\n\n");fflush(stdout);
872 MPI_Barrier(comm1);
873
874 sprintf(fnameout,"%s%s.current",filepath,filename);
875
876 MTRAN_Output_Current(
877 comm1,
878 fnameout,
879 TRAN_TKspace_grid2,
880 TRAN_TKspace_grid3,
881 SpinP_switch,
882 current
883 );
884
885 if (myid0==Host_ID) printf("\n");fflush(stdout);
886 }
887
888 /**********************************************
889 output conductance
890 **********************************************/
891
892 if (tran_transmission_on) {
893
894 MPI_Barrier(comm1);
895 if (myid0==Host_ID) printf("\nConductance: file\n\n");fflush(stdout);
896 MPI_Barrier(comm1);
897
898 sprintf(fnameout,"%s%s.conductance",filepath,filename);
899
900 MTRAN_Output_Conductance(
901 comm1,
902 fnameout,
903 T_k_ID,
904 T_knum,
905 TRAN_TKspace_grid2,
906 TRAN_TKspace_grid3,
907 T_IGrids2,
908 T_IGrids3,
909 T_KGrids2,
910 T_KGrids3,
911 SpinP_switch,
912 tran_transmission_energydiv,
913 tran_transmission_energyrange,
914 tran_transmission);
915
916 if (myid0==Host_ID) printf("\n");fflush(stdout);
917 }
918
919 /*S MitsuakiKAWAMURA2*/
920 /********************************************
921 Compute & Output Current density
922 *********************************************/
923
924 if (TRAN_CurrentDensity == 1) {
925
926 if (myid0 == Host_ID) printf("\nCurrentdensity: file\n\n"); fflush(stdout);
927
928 MP = (int*)malloc(sizeof(int)*(atomnum + 1));
929 MTRAN_Set_MP(1, atomnum, WhatSpecies, Spe_Total_CNO, &i, MP);
930
931 for (k = 0; k < SpinP_switch + 1; k++) {
932 for (iside = 0; iside < 3; iside++) {
933 for (i = 0; i < NUM_cs; i++) {
934 for (j = 0; j < NUM_cs; j++)
935 JLocSym[k][iside][i][j] = JLocSym[k][iside][i][j]
936 / (double)(TRAN_TKspace_grid2 * TRAN_TKspace_grid3);
937 MPI_Allreduce(MPI_IN_PLACE, JLocSym[k][iside][i], NUM_cs,
938 MPI_DOUBLE_PRECISION, MPI_SUM, comm1);
939 } /* for (i = 0; i < NUM_cs; i++) */
940 } /* for (iside = 0; iside < 3; iside++) */
941
942 for (iside = 0; iside < 2; iside++) {
943 for (i = 0; i < NUM_cs; i++) {
944 for (j = 0; j < NUM_cs; j++)
945 Jmat[k][iside][i][j] = Jmat[k][iside][i][j]
946 / (double)(TRAN_TKspace_grid2 * TRAN_TKspace_grid3);
947 MPI_Allreduce(MPI_IN_PLACE, Jmat[k][iside][i], NUM_cs,
948 MPI_DOUBLE_PRECISION, MPI_SUM, comm1);
949 } /* for (i = 0; i < NUM_cs; i++) */
950 } /* for (iside = 0; iside < 2; iside++) */
951
952 for (i = 0; i < NUM_cs; i++) {
953 for (j = 0; j < NUM_cs; j++) {
954 JLocASym[k][i][j] = JLocASym[k][i][j]
955 / (double)(TRAN_TKspace_grid2 * TRAN_TKspace_grid3);
956 RhoNL[k][i][j] = RhoNL[k][i][j]
957 / (double)(TRAN_TKspace_grid2 * TRAN_TKspace_grid3);
958 } /* for (j = 0; j < NUM_cs; j++) */
959 MPI_Allreduce(MPI_IN_PLACE, JLocASym[k][i], NUM_cs,
960 MPI_DOUBLE_PRECISION, MPI_SUM, comm1);
961 MPI_Allreduce(MPI_IN_PLACE, RhoNL[k][i], NUM_cs,
962 MPI_DOUBLE_PRECISION, MPI_SUM, comm1);
963 } /* for (i = 0; i < NUM_cs; i++) */
964 } /* for (k = 0; k < SpinP_switch + 1; k++) */
965
966 TRAN_CDen_Main(NUM_cs, MP, JLocSym, JLocASym, RhoNL, Jmat, SCC, TRAN_OffDiagonalCurrent);
967 free(MP);
968 } /*if (TRAN_CurrentDensity == 1)*/
969
970 for (k = 0; k < SpinP_switch + 1; k++) {
971 for (iside = 0; iside < 3; iside++) {
972 for (i = 0; i < NUM_cs; i++) {
973 free(JLocSym[k][iside][i]);
974 } /* for (i = 0; i < NUM_cs; i++) */
975 free(JLocSym[k][iside]);
976 } /* for (iside = 0; iside < 3; iside++) */
977 free(JLocSym[k]);
978
979 for (iside = 0; iside < 2; iside++) {
980 for (i = 0; i < NUM_cs; i++) {
981 free(Jmat[k][iside][i]);
982 } /* for (i = 0; i < NUM_cs; i++) */
983 free(Jmat[k][iside]);
984 } /* for (iside = 0; iside < 2; iside++) */
985 free(Jmat[k]);
986
987 for (i = 0; i < NUM_cs; i++) {
988 free(JLocASym[k][i]);
989 free(RhoNL[k][i]);
990 } /* for (i = 0; i < NUM_cs; i++) */
991 free(JLocASym[k]);
992 free(RhoNL[k]);
993 } /* for (k = 0; k < SpinP_switch + 1; k++) */
994 free(JLocSym);
995 free(Jmat);
996 free(JLocASym);
997 free(RhoNL);
998 /*E MitsuakiKAWAMURA2*/
999
1000 /* S MitsuakiKAWAMURA*/
1001 /********************************************
1002 Start EigenChannel Analysis
1003 *********************************************/
1004
1005 /***********************************************************
1006 start "kloop0"
1007 ***********************************************************/
1008 if (TRAN_Channel_Nenergy * TRAN_Channel_Nkpoint <= 0) TRAN_Channel == 0;
1009
1010 if (TRAN_Channel == 1){
1011
1012 if (myid0 == Host_ID){
1013 printf("\n**************************************************\n");
1014 printf(" Calculation of transmission eigenchannels starts\n");
1015 printf("**************************************************\n\n");
1016 printf(" File index : *.traneval#k_#E_#spin *.tranevec#k_#E_#spin \n\n");
1017 }
1018
1019 fflush(stdout);
1020 MPI_Barrier(comm1);
1021
1022 eigentrans = (double***)malloc(sizeof(double**) * TRAN_Channel_Nkpoint);
1023 eigentrans_sum = (double***)malloc(sizeof(double**) * TRAN_Channel_Nkpoint);
1024 EChannel = (dcomplex****)malloc(sizeof(dcomplex***) * TRAN_Channel_Nkpoint);
1025 for (kloop = 0; kloop < TRAN_Channel_Nkpoint; kloop++){
1026 eigentrans[kloop] = (double**)malloc(sizeof(double*) * TRAN_Channel_Nenergy);
1027 eigentrans_sum[kloop] = (double**)malloc(sizeof(double*) * TRAN_Channel_Nenergy);
1028 EChannel[kloop] = (dcomplex***)malloc(sizeof(dcomplex**) * TRAN_Channel_Nenergy);
1029 for (iw = 0; iw < TRAN_Channel_Nenergy; iw++){
1030 eigentrans[kloop][iw] = (double*)malloc(sizeof(double) * (TRAN_Channel_Num + 1));
1031 eigentrans_sum[kloop][iw] = (double*)malloc(sizeof(double) * 1);
1032 eigentrans_sum[kloop][iw][0] = 0.0;
1033 EChannel[kloop][iw] = (dcomplex**)malloc(sizeof(dcomplex*) * (TRAN_Channel_Num + 1));
1034 for (i = 0; i < TRAN_Channel_Num + 1; i++){
1035 eigentrans[kloop][iw][i] = 0.0;
1036 EChannel[kloop][iw][i] = (dcomplex*)malloc(sizeof(dcomplex) * NUM_c);
1037 for (j = 0; j < NUM_c; j++){
1038 EChannel[kloop][iw][i][j].r = 0.0;
1039 EChannel[kloop][iw][i][j].i = 0.0;
1040 }
1041 }
1042 }
1043 }
1044
1045 /* Set k parallelization mode*/
1046
1047 free(T_k_ID);
1048 T_k_ID = (int*)malloc(sizeof(int)*TRAN_Channel_Nkpoint);
1049
1050 if (numprocs0<TRAN_Channel_Nkpoint){
1051
1052 /* set parallel_mode */
1053 parallel_mode = 0;
1054
1055 /* allocation of kloop to ID */
1056
1057 for (ID = 0; ID<numprocs0; ID++){
1058
1059 tmp = (double)TRAN_Channel_Nkpoint / (double)numprocs0;
1060 S_knum = (int)((double)ID*(tmp + 1.0e-12));
1061 E_knum = (int)((double)(ID + 1)*(tmp + 1.0e-12)) - 1;
1062 if (ID == (numprocs0 - 1)) E_knum = TRAN_Channel_Nkpoint - 1;
1063 if (E_knum<0) E_knum = 0;
1064
1065 for (k = S_knum; k <= E_knum; k++){
1066 /* ID in the first level world */
1067 T_k_ID[k] = ID;
1068 }
1069 }
1070
1071 /* find own informations */
1072
1073 tmp = (double)TRAN_Channel_Nkpoint / (double)numprocs0;
1074 S_knum = (int)((double)myid0*(tmp + 1.0e-12));
1075 E_knum = (int)((double)(myid0 + 1)*(tmp + 1.0e-12)) - 1;
1076 if (myid0 == (numprocs0 - 1)) E_knum = TRAN_Channel_Nkpoint - 1;
1077 if (E_knum<0) E_knum = 0;
1078
1079 num_kloop0 = E_knum - S_knum + 1;
1080
1081 } /* if (numprocs0<TRAN_Channel_Nkpoint) */
1082 else {
1083
1084 /* set parallel_mode */
1085 parallel_mode = 1;
1086 num_kloop0 = 1;
1087
1088 Num_Comm_World2 = TRAN_Channel_Nkpoint;
1089
1090 NPROCS_ID2 = (int*)malloc(sizeof(int)*numprocs0);
1091 Comm_World2 = (int*)malloc(sizeof(int)*numprocs0);
1092 NPROCS_WD2 = (int*)malloc(sizeof(int)*Num_Comm_World2);
1093 Comm_World_StartID2 = (int*)malloc(sizeof(int)*Num_Comm_World2);
1094 MPI_CommWD2 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World2);
1095
1096 Make_Comm_Worlds(comm1, myid0, numprocs0, Num_Comm_World2, &myworld2, MPI_CommWD2,
1097 NPROCS_ID2, Comm_World2, NPROCS_WD2, Comm_World_StartID2);
1098
1099 MPI_Comm_size(MPI_CommWD2[myworld2], &numprocs2);
1100 MPI_Comm_rank(MPI_CommWD2[myworld2], &myid2);
1101
1102 S_knum = myworld2;
1103
1104 /* allocate k-points into processors */
1105
1106 for (k = 0; k<TRAN_Channel_Nkpoint; k++){
1107 /* ID in the first level world */
1108 T_k_ID[k] = Comm_World_StartID2[k];
1109 }
1110
1111 } /* else if (numprocs0>=TRAN_Channel_Nkpoint) */
1112
1113 for (kloop0 = 0; kloop0 < num_kloop0; kloop0++){
1114
1115 kloop = kloop0 + S_knum;
1116
1117 k2 = TRAN_Channel_kpoint[kloop][0];
1118 k3 = TRAN_Channel_kpoint[kloop][1];
1119
1120 if (parallel_mode){
1121 comm_tmp = MPI_CommWD2[myworld2];
1122 numprocs_tmp = numprocs2;
1123 myid_tmp = myid2;
1124 }
1125 else{
1126 comm_tmp = comm1;
1127 numprocs_tmp = 1;
1128 myid_tmp = 0;
1129 }
1130
1131 /* set Hamiltonian and overlap matrices of left and right leads */
1132
1133 MTRAN_Set_SurfOverlap("left", k2, k3, SpinP_switch, atomnum_e,
1134 OLP_e, H_e, iHNL_e, SpeciesNum_e, WhatSpecies_e, Spe_Total_CNO_e,
1135 FNAN_e, natn_e, ncn_e, atv_ijk_e, S00_e, S01_e, H00_e, H01_e);
1136
1137 MTRAN_Set_SurfOverlap("right", k2, k3, SpinP_switch, atomnum_e,
1138 OLP_e, H_e, iHNL_e, SpeciesNum_e, WhatSpecies_e, Spe_Total_CNO_e,
1139 FNAN_e, natn_e, ncn_e, atv_ijk_e, S00_e, S01_e, H00_e, H01_e);
1140
1141 /* set CC, CL and CR */
1142
1143 MTRAN_Set_CentOverlap(3,
1144 SpinP_switch,
1145 k2,
1146 k3,
1147 NUM_c,
1148 NUM_e,
1149 H,
1150 iHNL,
1151 OLP,
1152 atomnum,
1153 atomnum_e,
1154 WhatSpecies,
1155 WhatSpecies_e,
1156 Spe_Total_CNO,
1157 Spe_Total_CNO_e,
1158 FNAN,
1159 natn,
1160 ncn,
1161 atv_ijk,
1162 TRAN_region,
1163 TRAN_Original_Id
1164 );
1165
1166 /* calculate transmission */
1167
1168 MTRAN_EigenChannel_NC(
1169 comm_tmp,
1170 numprocs_tmp,
1171 myid_tmp,
1172 myid0,
1173 SpinP_switch,
1174 ChemP_e,
1175 NUM_c,
1176 NUM_e,
1177 H00_e,
1178 S00_e,
1179 H01_e,
1180 S01_e,
1181 HCC,
1182 HCL,
1183 HCR,
1184 SCC,
1185 SCL,
1186 SCR,
1187 tran_surfgreen_iteration_max,
1188 tran_surfgreen_eps,
1189 tran_transmission_energyrange,
1190 TRAN_Channel_Nenergy,
1191 TRAN_Channel_energy,
1192 TRAN_Channel_Num,
1193 kloop,
1194 TRAN_Channel_kpoint[kloop],
1195 EChannel[kloop],
1196 eigentrans[kloop],
1197 eigentrans_sum[kloop]
1198 );
1199
1200 } /* for (kloop0 = 0; kloop0 < num_kloop0; kloop0++) */
1201
1202 fflush(stdout);
1203
1204 MP = (int*)malloc(sizeof(int)*(atomnum + 1));
1205 MTRAN_Set_MP(1, atomnum, WhatSpecies, Spe_Total_CNO, &i, MP);
1206
1207 for (kloop = 0; kloop < TRAN_Channel_Nkpoint; kloop++){
1208 for (iw = 0; iw < TRAN_Channel_Nenergy; iw++){
1209 MPI_Allreduce(MPI_IN_PLACE, eigentrans[kloop][iw], TRAN_Channel_Num,
1210 MPI_DOUBLE_PRECISION, MPI_SUM, comm1);
1211 MPI_Allreduce(MPI_IN_PLACE, eigentrans_sum[kloop][iw], 1,
1212 MPI_DOUBLE_PRECISION, MPI_SUM, comm1);
1213 for (i = 0; i < TRAN_Channel_Num; i++){
1214 MPI_Allreduce(MPI_IN_PLACE, EChannel[kloop][iw][i], NUM_c,
1215 MPI_DOUBLE_COMPLEX, MPI_SUM, comm1);
1216 }
1217 }
1218 }
1219
1220 TRAN_Output_eigentrans_sum(TRAN_Channel_Nkpoint, TRAN_Channel_Nenergy, eigentrans_sum);
1221
1222 if (myid0 == Host_ID) {
1223 printf("\n Eigenchannel calculation finished \n");
1224 printf(" They are written in plottable files. \n");
1225 printf(" File index : *.tranec#k_#E_#spin_#branch_r.xsf \n");
1226 printf(" *.tranec#k_#E_#spin_#branch_i.xsf \n\n");
1227 }
1228
1229 for (kloop = 0; kloop < TRAN_Channel_Nkpoint; kloop++){
1230 for (iw = 0; iw < TRAN_Channel_Nenergy; iw++){
1231 for (i = 0; i < TRAN_Channel_Num; i++){
1232 for (k = 0; k < 2; k++){
1233 TRAN_Output_ChannelCube(
1234 kloop,
1235 iw,
1236 k,
1237 i,
1238 NUM_c / 2,
1239 TRAN_Channel_kpoint[kloop],
1240 &EChannel[kloop][iw][i][k * NUM_c / 2],
1241 MP,
1242 eigentrans[kloop][iw][i],
1243 TRAN_Channel_energy[iw]);
1244 }
1245 }
1246 }
1247
1248 }
1249
1250 free(MP);
1251 for (kloop = 0; kloop < TRAN_Channel_Nkpoint; kloop++){
1252 for (iw = 0; iw < TRAN_Channel_Nenergy; iw++){
1253 for (i = 0; i < TRAN_Channel_Num + 1; i++){
1254 free(EChannel[kloop][iw][i]);
1255 }
1256 free(eigentrans[kloop][iw]);
1257 free(EChannel[kloop][iw]);
1258 }
1259 free(eigentrans[kloop]);
1260 free(EChannel[kloop]);
1261 }
1262 free(eigentrans);
1263 free(EChannel);
1264
1265 free(NPROCS_ID2);
1266 free(Comm_World2);
1267 free(NPROCS_WD2);
1268 free(Comm_World_StartID2);
1269 free(MPI_CommWD2);
1270 } /* if (TRAN_Channel == 1) */
1271
1272 /********************************************
1273 Finish EigenChannel Analysis
1274 *********************************************/
1275 /* E MitsuakiKAWAMURA*/
1276
1277 /**********************************************
1278 freeing of arrays
1279 **********************************************/
1280
1281 if (tran_transmission_on) {
1282
1283 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
1284 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
1285 free(current[i2][i3]);
1286 }
1287 free(current[i2]);
1288 }
1289 free(current);
1290
1291 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
1292 free(op_flag[i2]);
1293 }
1294 free(op_flag);
1295
1296 free(T_KGrids2);
1297 free(T_KGrids3);
1298 free(T_IGrids2);
1299 free(T_IGrids3);
1300 free(T_op_flag);
1301 free(T_k_ID);
1302
1303 if (T_knum<=numprocs0){
1304
1305 if (Num_Comm_World1<=numprocs0){
1306 MPI_Comm_free(&MPI_CommWD1[myworld1]);
1307 }
1308
1309 free(NPROCS_ID1);
1310 free(Comm_World1);
1311 free(NPROCS_WD1);
1312 free(Comm_World_StartID1);
1313 free(MPI_CommWD1);
1314 }
1315 }
1316
1317 /* MPI_Finalize() and exit(0) should not be done here.
1318 MitsuakiKawamura
1319
1320 MPI_Finalize();
1321 exit(0);
1322 */
1323 }
1324
1325
1326
1327
1328
1329
1330
MTRAN_Read_Tran_HS(MPI_Comm comm1,char * filepath,char * filename,char * ext)1331 void MTRAN_Read_Tran_HS(MPI_Comm comm1, char *filepath, char *filename, char *ext )
1332 {
1333 FILE *fp;
1334 int iv[20];
1335 double v[20];
1336 double vtmp[50];
1337 int i,j,k,id;
1338 int Gc_AN,Cwan,tno0,tno1;
1339 int h_AN,Gh_AN,Hwan;
1340 int iside,spin;
1341 int size1;
1342 int *ia_vec;
1343 char fname[300];
1344 int inconsistency;
1345 int numprocs0,myid0;
1346
1347 MPI_Comm_size(comm1,&numprocs0);
1348 MPI_Comm_rank(comm1,&myid0);
1349
1350 if (tran_interpolate==0){
1351 sprintf(fname,"%s%s.%s",filepath,filename,ext);
1352 }
1353 else if (tran_interpolate==1){
1354 sprintf(fname,"%s%s",filepath,interpolate_filename1);
1355 }
1356
1357 /**********************************************
1358 read the first Hamiltonian regardless
1359 of performing the interpolation
1360 **********************************************/
1361
1362 if ( (fp=fopen(fname,"r"))==NULL ) {
1363 printf("cannot open %s\n",fname);
1364 printf("in MTRAN_Read_Tran_HS\n");
1365 exit(0);
1366 }
1367
1368 /* SpinP_switch, NUM_c, and NUM_e */
1369
1370 i=0;
1371 fread(iv, sizeof(int),4,fp);
1372 SpinP_switch = iv[i++];
1373 NUM_c = iv[i++];
1374 NUM_e[0] = iv[i++];
1375 NUM_e[1] = iv[i++];
1376
1377 if (PrintLevel){
1378 printf("spin=%d NUM_c=%d NUM_e=%d %d\n", SpinP_switch, NUM_c, NUM_e[0], NUM_e[1]);
1379 }
1380
1381 /* chemical potential */
1382
1383 i=0;
1384 fread(v,sizeof(double),3,fp);
1385 ChemP = v[i++];
1386 ChemP_e[0]= v[i++];
1387 ChemP_e[1]= v[i++];
1388
1389 if (myid0==0 && tran_interpolate==0){
1390 printf("\n");fflush(stdout);
1391 printf("Chemical potentials used in the SCF calculation\n"); fflush(stdout);
1392 printf(" Left lead: %15.12f (eV)\n",ChemP_e[0]*eV2Hartree);fflush(stdout);
1393 printf(" Right lead: %15.12f (eV)\n",ChemP_e[1]*eV2Hartree);fflush(stdout);
1394 }
1395
1396 /* SCF_tran_bias_apply */
1397
1398 fread(iv, sizeof(int),1,fp);
1399 SCF_tran_bias_apply = iv[0];
1400
1401 /* the number of atoms */
1402
1403 i=0;
1404 fread(iv, sizeof(int),3,fp);
1405 atomnum = iv[i++];
1406 atomnum_e[0] = iv[i++];
1407 atomnum_e[1] = iv[i++];
1408
1409 if (PrintLevel){
1410 printf("atomnum=%d atomnum_e=%d %d\n", atomnum, atomnum_e[0], atomnum_e[1]);
1411 }
1412
1413 /* the number of species */
1414
1415 i=0;
1416 fread(iv, sizeof(int),3,fp);
1417 SpeciesNum = iv[i++];
1418 SpeciesNum_e[0] = iv[i++];
1419 SpeciesNum_e[1] = iv[i++];
1420
1421 if (PrintLevel){
1422 printf("SpeciesNum=%d SpeciesNum_e=%d %d\n",SpeciesNum, SpeciesNum_e[0], SpeciesNum_e[1]);
1423 }
1424
1425 /* TCpyCell */
1426
1427 i=0;
1428 fread(iv, sizeof(int),3,fp);
1429 TCpyCell = iv[i++];
1430 TCpyCell_e[0] = iv[i++];
1431 TCpyCell_e[1] = iv[i++];
1432
1433 if (PrintLevel){
1434 printf("TCpyCell=%d TCpyCell_e=%d %d\n",TCpyCell, TCpyCell_e[0], TCpyCell_e[1]);
1435 }
1436
1437 /* TRAN_region */
1438
1439 TRAN_region = (int*)malloc(sizeof(int)*(atomnum+1));
1440 fread(TRAN_region,sizeof(int),atomnum+1,fp);
1441
1442 /* TRAN_Original_Id */
1443
1444 TRAN_Original_Id = (int*)malloc(sizeof(int)*(atomnum+1));
1445 fread(TRAN_Original_Id,sizeof(int),atomnum+1,fp);
1446
1447 /**********************************************
1448 informations of the central region
1449 **********************************************/
1450
1451 WhatSpecies = (int*)malloc(sizeof(int)*(atomnum+1));
1452 fread(WhatSpecies, sizeof(int), atomnum+1, fp);
1453
1454 if (PrintLevel){
1455 for (i=0; i<=atomnum; i++) {
1456 printf("i=%2d WhatSpecies=%2d\n",i,WhatSpecies[i]);
1457 }
1458 }
1459
1460 Spe_Total_CNO = (int*)malloc(sizeof(int)*(SpeciesNum));
1461 fread(Spe_Total_CNO, sizeof(int), SpeciesNum, fp);
1462
1463 FNAN = (int*)malloc(sizeof(int)*(atomnum+1));
1464 fread(FNAN,sizeof(int), atomnum+1,fp);
1465
1466 if (PrintLevel){
1467 for (i=0; i<=atomnum; i++) {
1468 printf("i=%2d FNAN=%2d\n",i,FNAN[i]);
1469 }
1470 }
1471
1472 fread(&Max_FSNAN,sizeof(int),1,fp);
1473 fread(&ScaleSize,sizeof(double),1,fp);
1474
1475 if (PrintLevel){
1476 printf("Max_FSNAN=%2d ScaleSize=%10.5f\n",Max_FSNAN,ScaleSize);
1477 }
1478
1479 size1=(int)(Max_FSNAN*ScaleSize) + 1;
1480
1481 natn = (int**)malloc(sizeof(int*)*(atomnum+1));
1482 for (i=0; i<=(atomnum); i++) {
1483 natn[i] = (int*)malloc(sizeof(int)*size1);
1484 }
1485 for (i=0; i<=(atomnum); i++) {
1486 fread(natn[i],sizeof(int),size1,fp);
1487 }
1488
1489 if (PrintLevel){
1490 for (i=0; i<=(atomnum); i++) {
1491 for (j=0; j<size1; j++) {
1492 printf("i=%3d j=%3d natn=%2d\n",i,j,natn[i][j]);
1493 }
1494 }
1495 }
1496
1497 ncn = (int**)malloc(sizeof(int*)*(atomnum+1));
1498 for (i=0; i<=(atomnum); i++) {
1499 ncn[i] = (int*)malloc(sizeof(int)*size1);
1500 }
1501 for (i=0; i<=(atomnum); i++) {
1502 fread(ncn[i],sizeof(int),size1,fp);
1503 }
1504
1505 if (PrintLevel){
1506 for (i=0; i<=(atomnum); i++) {
1507 for (j=0; j<size1; j++) {
1508 printf("i=%3d j=%3d ncn=%2d\n",i,j,natn[i][j]);
1509 }
1510 }
1511 }
1512
1513 size1=(TCpyCell+1)*4;
1514 ia_vec=(int*)malloc(sizeof(int)*size1);
1515 fread(ia_vec,sizeof(int),size1,fp);
1516
1517 atv_ijk = (int**)malloc(sizeof(int*)*(TCpyCell+1));
1518 for (i=0; i<(TCpyCell+1); i++) {
1519 atv_ijk[i] = (int*)malloc(sizeof(int)*4);
1520 }
1521
1522 id=0;
1523 for (i=0; i<(TCpyCell+1); i++) {
1524 for (j=0; j<=3; j++) {
1525 atv_ijk[i][j] = ia_vec[id++];
1526 }
1527
1528 if (PrintLevel){
1529 printf("atv_ijk %3d %2d %2d %2d\n",i,atv_ijk[i][1],atv_ijk[i][2],atv_ijk[i][3]);
1530 }
1531
1532 }
1533 free(ia_vec);
1534
1535 /* OLP */
1536
1537 OLP = (double*****)malloc(sizeof(double****)*4);
1538 for (k=0; k<4; k++) {
1539
1540 OLP[k] = (double****)malloc(sizeof(double***)*(atomnum+1));
1541
1542 FNAN[0] = 0;
1543 for (Gc_AN=0; Gc_AN<=(atomnum); Gc_AN++){
1544
1545 if (Gc_AN==0){
1546 tno0 = 1;
1547 }
1548 else{
1549 Cwan = WhatSpecies[Gc_AN];
1550 tno0 = Spe_Total_CNO[Cwan];
1551 }
1552
1553 OLP[k][Gc_AN] = (double***)malloc(sizeof(double**)*(FNAN[Gc_AN]+1));
1554
1555 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1556
1557 if (Gc_AN==0){
1558 tno1 = 1;
1559 }
1560 else {
1561 Gh_AN = natn[Gc_AN][h_AN];
1562 Hwan = WhatSpecies[Gh_AN];
1563 tno1 = Spe_Total_CNO[Hwan];
1564 }
1565
1566 OLP[k][Gc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
1567
1568 for (i=0; i<tno0; i++){
1569 OLP[k][Gc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
1570
1571 if (Gc_AN!=0){
1572 fread(OLP[k][Gc_AN][h_AN][i],sizeof(double), tno1, fp);
1573
1574 /*
1575 */
1576
1577 }
1578 }
1579 }
1580 }
1581 }
1582
1583 /* H */
1584
1585 H = (double*****)malloc(sizeof(double****)*(SpinP_switch+1));
1586 for (k=0; k<(SpinP_switch+1); k++) {
1587
1588 H[k] = (double****)malloc(sizeof(double***)*(atomnum+1));
1589
1590 FNAN[0] = 0;
1591 for (Gc_AN=0; Gc_AN<=atomnum; Gc_AN++){
1592
1593 if (Gc_AN==0){
1594 tno0 = 1;
1595 }
1596 else{
1597 Cwan = WhatSpecies[Gc_AN];
1598 tno0 = Spe_Total_CNO[Cwan];
1599 }
1600
1601 H[k][Gc_AN] = (double***)malloc(sizeof(double**)*(FNAN[Gc_AN]+1));
1602
1603 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1604
1605 if (Gc_AN==0){
1606 tno1 = 1;
1607 }
1608 else {
1609 Gh_AN = natn[Gc_AN][h_AN];
1610 Hwan = WhatSpecies[Gh_AN];
1611 tno1 = Spe_Total_CNO[Hwan];
1612 }
1613
1614 H[k][Gc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
1615
1616 for (i=0; i<tno0; i++){
1617 H[k][Gc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
1618 if (Gc_AN!=0){
1619 fread(H[k][Gc_AN][h_AN][i],sizeof(double),tno1,fp);
1620 }
1621 }
1622 }
1623 }
1624 }
1625
1626 /* revised by Y. Xiao for Noncollinear NEGF calculations */
1627 if(SpinP_switch == 3) {
1628 iHNL = (double*****)malloc(sizeof(double****)*(2+1));
1629 for (k=0; k<(2+1); k++) {
1630
1631 iHNL[k] = (double****)malloc(sizeof(double***)*(atomnum+1));
1632
1633 FNAN[0] = 0;
1634 for (Gc_AN=0; Gc_AN<=atomnum; Gc_AN++){
1635
1636 if (Gc_AN==0){
1637 tno0 = 1;
1638 }
1639 else{
1640 Cwan = WhatSpecies[Gc_AN];
1641 tno0 = Spe_Total_CNO[Cwan];
1642 }
1643
1644 iHNL[k][Gc_AN] = (double***)malloc(sizeof(double**)*(FNAN[Gc_AN]+1));
1645
1646 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1647
1648 if (Gc_AN==0){
1649 tno1 = 1;
1650 }
1651 else {
1652 Gh_AN = natn[Gc_AN][h_AN];
1653 Hwan = WhatSpecies[Gh_AN];
1654 tno1 = Spe_Total_CNO[Hwan];
1655 }
1656
1657 iHNL[k][Gc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
1658
1659 for (i=0; i<tno0; i++){
1660 iHNL[k][Gc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
1661 if (Gc_AN!=0){
1662 fread(iHNL[k][Gc_AN][h_AN][i],sizeof(double),tno1,fp);
1663 }
1664 }
1665 }
1666 }
1667 }
1668 } /* if (SpinP_switch == 3) */
1669
1670 /* H0 (matrix elements for kinetic energy) */
1671
1672 H0 = (double*****)malloc(sizeof(double****)*4);
1673 for (k=0; k<4; k++) {
1674
1675 H0[k] = (double****)malloc(sizeof(double***)*(atomnum+1));
1676
1677 FNAN[0] = 0;
1678 for (Gc_AN=0; Gc_AN<=(atomnum); Gc_AN++){
1679
1680 if (Gc_AN==0){
1681 tno0 = 1;
1682 }
1683 else{
1684 Cwan = WhatSpecies[Gc_AN];
1685 tno0 = Spe_Total_CNO[Cwan];
1686 }
1687
1688 H0[k][Gc_AN] = (double***)malloc(sizeof(double**)*(FNAN[Gc_AN]+1));
1689
1690 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
1691
1692 if (Gc_AN==0){
1693 tno1 = 1;
1694 }
1695 else {
1696 Gh_AN = natn[Gc_AN][h_AN];
1697 Hwan = WhatSpecies[Gh_AN];
1698 tno1 = Spe_Total_CNO[Hwan];
1699 }
1700
1701 H0[k][Gc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
1702
1703 for (i=0; i<tno0; i++){
1704 H0[k][Gc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
1705
1706 if (Gc_AN!=0){
1707 fread(H0[k][Gc_AN][h_AN][i],sizeof(double), tno1, fp);
1708
1709 /*
1710 for (j=0; j<tno1; j++){
1711 printf("k=%2d Gc_AN=%2d h_AN=%2d i=%2d j=%2d H0=%15.12f\n",
1712 k,Gc_AN,h_AN,i,j,H0[k][Gc_AN][h_AN][i][j]);
1713 }
1714 */
1715
1716 }
1717 }
1718 }
1719 }
1720 }
1721
1722 /**********************************************
1723 informations of leads
1724 **********************************************/
1725
1726 for (iside=0; iside<=1; iside++) {
1727
1728 WhatSpecies_e[iside] = (int*)malloc(sizeof(int)*(atomnum_e[iside]+1));
1729 fread(WhatSpecies_e[iside], sizeof(int), atomnum_e[iside]+1, fp);
1730
1731 if (PrintLevel){
1732 for (i=0; i<=atomnum_e[iside]; i++) {
1733 printf("iside=%2d i=%2d WhatSpecies_e=%2d\n",iside,i,WhatSpecies_e[iside][i]);
1734 }
1735 }
1736
1737 Spe_Total_CNO_e[iside] = (int*)malloc(sizeof(int)*SpeciesNum_e[iside]);
1738 fread(Spe_Total_CNO_e[iside], sizeof(int), SpeciesNum_e[iside], fp);
1739
1740 FNAN_e[iside] = (int*)malloc(sizeof(int)*(atomnum_e[iside]+1));
1741 fread(FNAN_e[iside], sizeof(int), atomnum_e[iside]+1, fp);
1742
1743 fread(&Max_FSNAN_e[iside], sizeof(int), 1, fp);
1744 fread(&ScaleSize_e[iside], sizeof(double),1, fp);
1745
1746 size1=(int)(Max_FSNAN_e[iside]*ScaleSize_e[iside]) + 1;
1747
1748 natn_e[iside] = (int**)malloc(sizeof(int*)*(atomnum_e[iside]+1));
1749 for (i=0; i<=atomnum_e[iside]; i++) {
1750 natn_e[iside][i] = (int*)malloc(sizeof(int)*size1);
1751 }
1752 for (i=0; i<=atomnum_e[iside]; i++) {
1753 fread(natn_e[iside][i],sizeof(int),size1,fp);
1754 }
1755
1756 ncn_e[iside] = (int**)malloc(sizeof(int*)*(atomnum_e[iside]+1));
1757 for (i=0; i<=atomnum_e[iside]; i++) {
1758 ncn_e[iside][i] = (int*)malloc(sizeof(int)*size1);
1759 }
1760 for (i=0; i<=atomnum_e[iside]; i++) {
1761 fread(ncn_e[iside][i],sizeof(int),size1,fp);
1762 }
1763
1764 size1=(TCpyCell_e[iside]+1)*4;
1765 ia_vec=(int*)malloc(sizeof(int)*size1);
1766 fread(ia_vec,sizeof(int),size1,fp);
1767
1768 atv_ijk_e[iside] = (int**)malloc(sizeof(int*)*(TCpyCell_e[iside]+1));
1769 for (i=0; i<(TCpyCell_e[iside]+1); i++) {
1770 atv_ijk_e[iside][i] = (int*)malloc(sizeof(int)*4);
1771 }
1772
1773 id=0;
1774 for (i=0; i<(TCpyCell_e[iside]+1); i++) {
1775 for (j=0; j<=3; j++) {
1776 atv_ijk_e[iside][i][j] = ia_vec[id++];
1777 }
1778 }
1779 free(ia_vec);
1780
1781 /* overlap matrix */
1782
1783 OLP_e[iside] = (double*****)malloc(sizeof(double****)*4);
1784 for (k=0; k<4; k++) {
1785
1786 OLP_e[iside][k] = (double****)malloc(sizeof(double***)*(atomnum_e[iside]+1));
1787
1788 FNAN_e[iside][0] = 0;
1789 for (Gc_AN=0; Gc_AN<=atomnum_e[iside]; Gc_AN++){
1790
1791 if (Gc_AN==0){
1792 tno0 = 1;
1793 }
1794 else{
1795 Cwan = WhatSpecies_e[iside][Gc_AN];
1796 Cwan = WhatSpecies_e[iside][Gc_AN];
1797 tno0 = Spe_Total_CNO_e[iside][Cwan];
1798 }
1799
1800 OLP_e[iside][k][Gc_AN] = (double***)malloc(sizeof(double**)*(FNAN_e[iside][Gc_AN]+1));
1801
1802 for (h_AN=0; h_AN<=FNAN_e[iside][Gc_AN]; h_AN++){
1803
1804 if (Gc_AN==0){
1805 tno1 = 1;
1806 }
1807 else{
1808 Gh_AN = natn_e[iside][Gc_AN][h_AN];
1809 Hwan = WhatSpecies_e[iside][Gh_AN];
1810 tno1 = Spe_Total_CNO_e[iside][Hwan];
1811 }
1812
1813 OLP_e[iside][k][Gc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
1814
1815 for (i=0; i<tno0; i++){
1816 OLP_e[iside][k][Gc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
1817 if (Gc_AN!=0){
1818 fread(OLP_e[iside][k][Gc_AN][h_AN][i],sizeof(double),tno1,fp);
1819 }
1820 }
1821 }
1822 }
1823 }
1824
1825 /* Hamiltonian matrix */
1826
1827 H_e[iside] = (double*****)malloc(sizeof(double****)*(SpinP_switch+1));
1828 for (k=0; k<(SpinP_switch+1); k++) {
1829
1830 H_e[iside][k] = (double****)malloc(sizeof(double***)*(atomnum_e[iside]+1));
1831
1832 FNAN_e[iside][0] = 0;
1833 for (Gc_AN=0; Gc_AN<=atomnum_e[iside]; Gc_AN++){
1834
1835 if (Gc_AN==0){
1836 tno0 = 1;
1837 }
1838 else{
1839 Cwan = WhatSpecies_e[iside][Gc_AN];
1840 tno0 = Spe_Total_CNO_e[iside][Cwan];
1841 }
1842
1843 H_e[iside][k][Gc_AN] = (double***)malloc(sizeof(double**)*(FNAN_e[iside][Gc_AN]+1));
1844
1845 for (h_AN=0; h_AN<=FNAN_e[iside][Gc_AN]; h_AN++){
1846
1847 if (Gc_AN==0){
1848 tno1 = 1;
1849 }
1850 else{
1851 Gh_AN = natn_e[iside][Gc_AN][h_AN];
1852 Hwan = WhatSpecies_e[iside][Gh_AN];
1853 tno1 = Spe_Total_CNO_e[iside][Hwan];
1854 }
1855
1856 H_e[iside][k][Gc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
1857
1858 for (i=0; i<tno0; i++){
1859 H_e[iside][k][Gc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
1860 if (Gc_AN!=0){
1861 fread(H_e[iside][k][Gc_AN][h_AN][i],sizeof(double),tno1,fp);
1862 }
1863 }
1864 }
1865 }
1866 }
1867
1868 /* revised by Y. Xiao for Noncollinear NEGF calculations */
1869 if (SpinP_switch == 3) {
1870 iHNL_e[iside] = (double*****)malloc(sizeof(double****)*(2+1));
1871 for (k=0; k<(2+1); k++) {
1872
1873 iHNL_e[iside][k] = (double****)malloc(sizeof(double***)*(atomnum_e[iside]+1));
1874
1875 FNAN_e[iside][0] = 0;
1876 for (Gc_AN=0; Gc_AN<=atomnum_e[iside]; Gc_AN++){
1877
1878 if (Gc_AN==0){
1879 tno0 = 1;
1880 }
1881 else{
1882 Cwan = WhatSpecies_e[iside][Gc_AN];
1883 tno0 = Spe_Total_CNO_e[iside][Cwan];
1884 }
1885
1886 iHNL_e[iside][k][Gc_AN] = (double***)malloc(sizeof(double**)*(FNAN_e[iside][Gc_AN]+1));
1887
1888 for (h_AN=0; h_AN<=FNAN_e[iside][Gc_AN]; h_AN++){
1889
1890 if (Gc_AN==0){
1891 tno1 = 1;
1892 }
1893 else{
1894 Gh_AN = natn_e[iside][Gc_AN][h_AN];
1895 Hwan = WhatSpecies_e[iside][Gh_AN];
1896 tno1 = Spe_Total_CNO_e[iside][Hwan];
1897 }
1898
1899 iHNL_e[iside][k][Gc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
1900
1901 for (i=0; i<tno0; i++){
1902 iHNL_e[iside][k][Gc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
1903 if (Gc_AN!=0){
1904 fread(iHNL_e[iside][k][Gc_AN][h_AN][i],sizeof(double),tno1,fp);
1905 }
1906 }
1907 }
1908 }
1909 }
1910 }
1911
1912 }
1913
1914 /**********************************************
1915 close the file pointer
1916 **********************************************/
1917
1918 fclose(fp);
1919
1920 /************************************************************
1921 interpolation of the bias voltage effect
1922 ************************************************************/
1923
1924 if (tran_interpolate==1){
1925
1926 inconsistency = 0;
1927
1928 sprintf(fname,"%s%s",filepath,interpolate_filename2);
1929
1930 if ( (fp=fopen(fname,"r"))==NULL ) {
1931 printf("cannot open %s\n",fname);
1932 printf("in MTRAN_Read_Tran_HS\n");
1933 exit(0);
1934 }
1935
1936 /* SpinP_switch, NUM_c, and NUM_e */
1937
1938 i=0;
1939 fread(iv, sizeof(int),4,fp);
1940
1941 if (iv[i++]!=SpinP_switch) inconsistency++;
1942 if (iv[i++]!=NUM_c) inconsistency++;
1943 if (iv[i++]!=NUM_e[0]) inconsistency++;
1944 if (iv[i++]!=NUM_e[1]) inconsistency++;
1945
1946 /* chemical potential */
1947
1948 i=0;
1949 fread(v,sizeof(double),3,fp);
1950
1951 if (myid0==0){
1952 printf("The interpolation method is used to take account of the bias voltage.\n");
1953 printf(" Coefficients for the interpolation %10.6f %10.6f\n\n",interpolate_c1,interpolate_c2); fflush(stdout);
1954
1955 printf(" First SCF calc. Second SCF calc. Interporated\n");fflush(stdout);
1956 printf(" Left lead: %15.12f %15.12f %15.12f (eV)\n",
1957 ChemP_e[0]*eV2Hartree,
1958 v[1]*eV2Hartree,
1959 (interpolate_c1*ChemP_e[0] + interpolate_c2*v[1])*eV2Hartree);fflush(stdout);
1960
1961 printf(" Right lead: %15.12f %15.12f %15.12f (eV)\n",
1962 ChemP_e[1]*eV2Hartree,
1963 v[2]*eV2Hartree,
1964 (interpolate_c1*ChemP_e[1] + interpolate_c2*v[2])*eV2Hartree);fflush(stdout);
1965
1966 }
1967
1968 ChemP = interpolate_c1*ChemP + interpolate_c2*v[i++];
1969 ChemP_e[0] = interpolate_c1*ChemP_e[0] + interpolate_c2*v[i++];
1970 ChemP_e[1] = interpolate_c1*ChemP_e[1] + interpolate_c2*v[i++];
1971
1972 /* SCF_tran_bias_apply */
1973
1974 fread(iv, sizeof(int),1,fp);
1975 SCF_tran_bias_apply = iv[0];
1976
1977 /* the number of atoms */
1978
1979 i=0;
1980 fread(iv, sizeof(int),3,fp);
1981
1982 if (iv[i++]!=atomnum) inconsistency++;
1983 if (iv[i++]!=atomnum_e[0]) inconsistency++;
1984 if (iv[i++]!=atomnum_e[1]) inconsistency++;
1985
1986 /* the number of species */
1987
1988 i=0;
1989 fread(iv, sizeof(int),3,fp);
1990
1991 if (iv[i++]!=SpeciesNum) inconsistency++;
1992 if (iv[i++]!=SpeciesNum_e[0]) inconsistency++;
1993 if (iv[i++]!=SpeciesNum_e[1]) inconsistency++;
1994
1995 /* TCpyCell */
1996
1997 i=0;
1998 fread(iv, sizeof(int),3,fp);
1999
2000 if (iv[i++]!=TCpyCell) inconsistency++;
2001 if (iv[i++]!=TCpyCell_e[0]) inconsistency++;
2002 if (iv[i++]!=TCpyCell_e[1]) inconsistency++;
2003
2004 /* TRAN_region */
2005
2006 TRAN_region = (int*)malloc(sizeof(int)*(atomnum+1));
2007 fread(TRAN_region,sizeof(int),atomnum+1,fp);
2008
2009 /* TRAN_Original_Id */
2010
2011 TRAN_Original_Id = (int*)malloc(sizeof(int)*(atomnum+1));
2012 fread(TRAN_Original_Id,sizeof(int),atomnum+1,fp);
2013
2014 /**********************************************
2015 informations of the central region
2016 **********************************************/
2017
2018 WhatSpecies = (int*)malloc(sizeof(int)*(atomnum+1));
2019 fread(WhatSpecies, sizeof(int), atomnum+1, fp);
2020
2021 Spe_Total_CNO = (int*)malloc(sizeof(int)*(SpeciesNum));
2022 fread(Spe_Total_CNO, sizeof(int), SpeciesNum, fp);
2023
2024 FNAN = (int*)malloc(sizeof(int)*(atomnum+1));
2025 fread(FNAN,sizeof(int), atomnum+1,fp);
2026
2027 fread(&Max_FSNAN,sizeof(int),1,fp);
2028 fread(&ScaleSize,sizeof(double),1,fp);
2029
2030 size1=(int)(Max_FSNAN*ScaleSize) + 1;
2031
2032 natn = (int**)malloc(sizeof(int*)*(atomnum+1));
2033 for (i=0; i<=(atomnum); i++) {
2034 natn[i] = (int*)malloc(sizeof(int)*size1);
2035 }
2036 for (i=0; i<=(atomnum); i++) {
2037 fread(natn[i],sizeof(int),size1,fp);
2038 }
2039
2040 ncn = (int**)malloc(sizeof(int*)*(atomnum+1));
2041 for (i=0; i<=(atomnum); i++) {
2042 ncn[i] = (int*)malloc(sizeof(int)*size1);
2043 }
2044 for (i=0; i<=(atomnum); i++) {
2045 fread(ncn[i],sizeof(int),size1,fp);
2046 }
2047
2048 size1=(TCpyCell+1)*4;
2049 ia_vec=(int*)malloc(sizeof(int)*size1);
2050 fread(ia_vec,sizeof(int),size1,fp);
2051
2052 atv_ijk = (int**)malloc(sizeof(int*)*(TCpyCell+1));
2053 for (i=0; i<(TCpyCell+1); i++) {
2054 atv_ijk[i] = (int*)malloc(sizeof(int)*4);
2055 }
2056
2057 id=0;
2058 for (i=0; i<(TCpyCell+1); i++) {
2059 for (j=0; j<=3; j++) {
2060 atv_ijk[i][j] = ia_vec[id++];
2061 }
2062 }
2063 free(ia_vec);
2064
2065 /* OLP */
2066
2067 for (k=0; k<4; k++) {
2068
2069 for (Gc_AN=1; Gc_AN<=(atomnum); Gc_AN++){
2070
2071 Cwan = WhatSpecies[Gc_AN];
2072 tno0 = Spe_Total_CNO[Cwan];
2073
2074 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2075
2076 Gh_AN = natn[Gc_AN][h_AN];
2077 Hwan = WhatSpecies[Gh_AN];
2078 tno1 = Spe_Total_CNO[Hwan];
2079
2080 for (i=0; i<tno0; i++){
2081 fread(OLP[k][Gc_AN][h_AN][i],sizeof(double), tno1, fp);
2082 }
2083 }
2084 }
2085 }
2086
2087 /* H */
2088
2089 for (k=0; k<(SpinP_switch+1); k++) {
2090
2091 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
2092
2093 Cwan = WhatSpecies[Gc_AN];
2094 tno0 = Spe_Total_CNO[Cwan];
2095
2096 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2097
2098 Gh_AN = natn[Gc_AN][h_AN];
2099 Hwan = WhatSpecies[Gh_AN];
2100 tno1 = Spe_Total_CNO[Hwan];
2101
2102 for (i=0; i<tno0; i++){
2103
2104 fread(&vtmp[0],sizeof(double),tno1,fp);
2105
2106 for (j=0; j<tno1; j++){
2107 H[k][Gc_AN][h_AN][i][j] = interpolate_c1*H[k][Gc_AN][h_AN][i][j] + interpolate_c2*vtmp[j];
2108 }
2109 }
2110 }
2111 }
2112 }
2113
2114 /* iHNL */
2115 /*revised by Y. Xiao for Noncollinear NEGF calculations */
2116 if( SpinP_switch == 3) {
2117 for (k=0; k<(2+1); k++) {
2118
2119 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
2120
2121 Cwan = WhatSpecies[Gc_AN];
2122 tno0 = Spe_Total_CNO[Cwan];
2123
2124 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2125
2126 Gh_AN = natn[Gc_AN][h_AN];
2127 Hwan = WhatSpecies[Gh_AN];
2128 tno1 = Spe_Total_CNO[Hwan];
2129
2130 for (i=0; i<tno0; i++){
2131
2132 fread(&vtmp[0],sizeof(double),tno1,fp);
2133
2134 for (j=0; j<tno1; j++){
2135 iHNL[k][Gc_AN][h_AN][i][j] = interpolate_c1*iHNL[k][Gc_AN][h_AN][i][j] + interpolate_c2*vtmp[j];
2136 }
2137 }
2138 }
2139 }
2140 }
2141 } /* if( SpinP_switch == 3) */
2142 /* until here by Y. Xiao for Noncollinear NEGF calculations */
2143
2144 /**********************************************
2145 informations of leads
2146 **********************************************/
2147
2148 for (iside=0; iside<=1; iside++) {
2149
2150 WhatSpecies_e[iside] = (int*)malloc(sizeof(int)*(atomnum_e[iside]+1));
2151 fread(WhatSpecies_e[iside], sizeof(int), atomnum_e[iside]+1, fp);
2152
2153 Spe_Total_CNO_e[iside] = (int*)malloc(sizeof(int)*SpeciesNum_e[iside]);
2154 fread(Spe_Total_CNO_e[iside], sizeof(int), SpeciesNum_e[iside], fp);
2155
2156 FNAN_e[iside] = (int*)malloc(sizeof(int)*(atomnum_e[iside]+1));
2157 fread(FNAN_e[iside], sizeof(int), atomnum_e[iside]+1, fp);
2158
2159 fread(&Max_FSNAN_e[iside], sizeof(int), 1, fp);
2160 fread(&ScaleSize_e[iside], sizeof(double),1, fp);
2161
2162 size1=(int)(Max_FSNAN_e[iside]*ScaleSize_e[iside]) + 1;
2163
2164 natn_e[iside] = (int**)malloc(sizeof(int*)*(atomnum_e[iside]+1));
2165 for (i=0; i<=atomnum_e[iside]; i++) {
2166 natn_e[iside][i] = (int*)malloc(sizeof(int)*size1);
2167 }
2168 for (i=0; i<=atomnum_e[iside]; i++) {
2169 fread(natn_e[iside][i],sizeof(int),size1,fp);
2170 }
2171
2172 ncn_e[iside] = (int**)malloc(sizeof(int*)*(atomnum_e[iside]+1));
2173 for (i=0; i<=atomnum_e[iside]; i++) {
2174 ncn_e[iside][i] = (int*)malloc(sizeof(int)*size1);
2175 }
2176 for (i=0; i<=atomnum_e[iside]; i++) {
2177 fread(ncn_e[iside][i],sizeof(int),size1,fp);
2178 }
2179
2180 size1=(TCpyCell_e[iside]+1)*4;
2181 ia_vec=(int*)malloc(sizeof(int)*size1);
2182 fread(ia_vec,sizeof(int),size1,fp);
2183
2184 atv_ijk_e[iside] = (int**)malloc(sizeof(int*)*(TCpyCell_e[iside]+1));
2185 for (i=0; i<(TCpyCell_e[iside]+1); i++) {
2186 atv_ijk_e[iside][i] = (int*)malloc(sizeof(int)*4);
2187 }
2188
2189 id=0;
2190 for (i=0; i<(TCpyCell_e[iside]+1); i++) {
2191 for (j=0; j<=3; j++) {
2192 atv_ijk_e[iside][i][j] = ia_vec[id++];
2193 }
2194 }
2195 free(ia_vec);
2196
2197 /* overlap matrix */
2198
2199 for (k=0; k<4; k++) {
2200 for (Gc_AN=1; Gc_AN<=atomnum_e[iside]; Gc_AN++){
2201
2202 Cwan = WhatSpecies_e[iside][Gc_AN];
2203 Cwan = WhatSpecies_e[iside][Gc_AN];
2204 tno0 = Spe_Total_CNO_e[iside][Cwan];
2205
2206 for (h_AN=0; h_AN<=FNAN_e[iside][Gc_AN]; h_AN++){
2207
2208 Gh_AN = natn_e[iside][Gc_AN][h_AN];
2209 Hwan = WhatSpecies_e[iside][Gh_AN];
2210 tno1 = Spe_Total_CNO_e[iside][Hwan];
2211
2212 for (i=0; i<tno0; i++){
2213 fread(OLP_e[iside][k][Gc_AN][h_AN][i],sizeof(double),tno1,fp);
2214 }
2215 }
2216 }
2217 }
2218
2219 /* Hamiltonian matrix */
2220
2221 for (k=0; k<(SpinP_switch+1); k++) {
2222
2223 for (Gc_AN=1; Gc_AN<=atomnum_e[iside]; Gc_AN++){
2224
2225 Cwan = WhatSpecies_e[iside][Gc_AN];
2226 tno0 = Spe_Total_CNO_e[iside][Cwan];
2227
2228 for (h_AN=0; h_AN<=FNAN_e[iside][Gc_AN]; h_AN++){
2229
2230 Gh_AN = natn_e[iside][Gc_AN][h_AN];
2231 Hwan = WhatSpecies_e[iside][Gh_AN];
2232 tno1 = Spe_Total_CNO_e[iside][Hwan];
2233
2234 for (i=0; i<tno0; i++){
2235
2236 fread(&vtmp[0],sizeof(double),tno1,fp);
2237
2238 for (j=0; j<tno1; j++){
2239 H_e[iside][k][Gc_AN][h_AN][i][j] = interpolate_c1*H_e[iside][k][Gc_AN][h_AN][i][j] + interpolate_c2*vtmp[j];
2240 }
2241 }
2242 }
2243 }
2244 }
2245
2246 /* iHNL */
2247 /* revised by Y. Xiao for Noncollinear NEGF calculations */
2248 if ( SpinP_switch == 3 ) {
2249 for (k=0; k<(2+1); k++) {
2250
2251 for (Gc_AN=1; Gc_AN<=atomnum_e[iside]; Gc_AN++){
2252
2253 Cwan = WhatSpecies_e[iside][Gc_AN];
2254 tno0 = Spe_Total_CNO_e[iside][Cwan];
2255
2256 for (h_AN=0; h_AN<=FNAN_e[iside][Gc_AN]; h_AN++){
2257
2258 Gh_AN = natn_e[iside][Gc_AN][h_AN];
2259 Hwan = WhatSpecies_e[iside][Gh_AN];
2260 tno1 = Spe_Total_CNO_e[iside][Hwan];
2261
2262 for (i=0; i<tno0; i++){
2263
2264 fread(&vtmp[0],sizeof(double),tno1,fp);
2265
2266 for (j=0; j<tno1; j++){
2267 iHNL_e[iside][k][Gc_AN][h_AN][i][j] = interpolate_c1*iHNL_e[iside][k][Gc_AN][h_AN][i][j] + interpolate_c2*vtmp[j];
2268 }
2269 }
2270 }
2271 }
2272 }
2273 } /* if ( SpinP_switch == 3 ) */
2274 /* until here by Y. Xiao for Noncollinear NEGF calculations */
2275
2276 } /* iside */
2277
2278
2279 if (inconsistency!=0){
2280 printf("found some inconsistency in two files %s %s\n",interpolate_filename1,interpolate_filename2);
2281 MPI_Finalize();
2282 exit(0);
2283 }
2284
2285 } /* if (tran_interpolate==1) */
2286
2287
2288 }
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
MTRAN_Transmission(MPI_Comm comm1,int numprocs,int myid,int SpinP_switch,double ChemP_e[2],int NUM_c,int NUM_e[2],dcomplex * H00_e[2],dcomplex * S00_e[2],dcomplex * H01_e[2],dcomplex * S01_e[2],dcomplex * HCC,dcomplex * HCL,dcomplex * HCR,dcomplex * SCC,dcomplex * SCL,dcomplex * SCR,double tran_surfgreen_iteration_max,double tran_surfgreen_eps,double tran_transmission_energyrange[3],int tran_transmission_energydiv,dcomplex ** tran_transmission)2303 void MTRAN_Transmission(MPI_Comm comm1,
2304 int numprocs,
2305 int myid,
2306 int SpinP_switch,
2307 double ChemP_e[2],
2308 int NUM_c,
2309 int NUM_e[2],
2310 dcomplex *H00_e[2],
2311 dcomplex *S00_e[2],
2312 dcomplex *H01_e[2],
2313 dcomplex *S01_e[2],
2314 dcomplex *HCC,
2315 dcomplex *HCL,
2316 dcomplex *HCR,
2317 dcomplex *SCC,
2318 dcomplex *SCL,
2319 dcomplex *SCR,
2320 double tran_surfgreen_iteration_max,
2321 double tran_surfgreen_eps,
2322 double tran_transmission_energyrange[3],
2323 int tran_transmission_energydiv,
2324 dcomplex **tran_transmission)
2325 {
2326 dcomplex w;
2327 dcomplex *GRL,*GRR;
2328 dcomplex *GC_R,*GC_A;
2329 dcomplex *SigmaL_R,*SigmaL_A;
2330 dcomplex *SigmaR_R,*SigmaR_A;
2331 dcomplex *v1,*v2;
2332 dcomplex value;
2333
2334 int iw,k,iside;
2335 int ID;
2336 int tag=99;
2337
2338 int **iwIdx;
2339 int Miw,Miwmax ;
2340 int i,j;
2341 MPI_Status status;
2342
2343 v1 = (dcomplex*) malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2344 v2 = (dcomplex*) malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2345
2346 /* allocate */
2347 GRL = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[0]* NUM_e[0]);
2348 GRR = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[1]* NUM_e[1]);
2349
2350 GC_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
2351 GC_A = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
2352 SigmaL_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
2353 SigmaL_A = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
2354 SigmaR_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
2355 SigmaR_A = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
2356
2357 /* initialize */
2358
2359 for (k=0; k<=1; k++) {
2360 for (iw=0; iw<tran_transmission_energydiv ; iw++) {
2361 tran_transmission[k][iw].r = 0.0;
2362 tran_transmission[k][iw].i = 0.0;
2363 }
2364 }
2365
2366 /*parallel setup*/
2367
2368 iwIdx=(int**)malloc(sizeof(int*)*numprocs);
2369 Miwmax = (tran_transmission_energydiv)/numprocs+1;
2370 for (i=0;i<numprocs;i++) {
2371 iwIdx[i]=(int*)malloc(sizeof(int)*Miwmax);
2372 }
2373 TRAN_Distribute_Node_Idx(0, tran_transmission_energydiv-1, numprocs, Miwmax,
2374 iwIdx); /* output */
2375
2376 /* parallel global iw 0:tran_transmission_energydiv-1 */
2377 /* parallel local Miw 0:Miwmax-1 */
2378 /* parallel variable iw=iwIdx[myid][Miw] */
2379
2380 for (Miw=0; Miw<Miwmax ; Miw++) {
2381
2382 iw = iwIdx[myid][Miw];
2383
2384 if ( iw>=0 ) {
2385
2386 w.r = tran_transmission_energyrange[0] + ChemP_e[0]
2387 +
2388 (tran_transmission_energyrange[1]-tran_transmission_energyrange[0])*
2389 (double)iw/(tran_transmission_energydiv-1);
2390
2391 w.i = tran_transmission_energyrange[2];
2392
2393 /*
2394 printf("iw=%d of %d w= % 9.6e % 9.6e \n" ,iw, tran_transmission_energydiv, w.r,w.i);
2395 */
2396
2397 /* for (k=0; k<=SpinP_switch; k++) { */
2398 for (k=0; k<=0; k++) {
2399
2400 /*****************************************************************
2401 Note that retarded and advanced Green functions and self energies
2402 are not conjugate comlex in case of the k-dependent case.
2403 **************************************************************/
2404
2405 /* in case of retarded ones */
2406
2407 iside=0;
2408 TRAN_Calc_SurfGreen_direct(w,NUM_e[iside], H00_e[iside],H01_e[iside],
2409 S00_e[iside], S01_e[iside],
2410 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRL);
2411
2412 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRL, NUM_c, HCL, SCL, SigmaL_R);
2413
2414 iside=1;
2415 TRAN_Calc_SurfGreen_direct(w,NUM_e[iside], H00_e[iside],H01_e[iside],
2416 S00_e[iside], S01_e[iside],
2417 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRR);
2418
2419 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRR, NUM_c, HCR, SCR, SigmaR_R);
2420
2421 TRAN_Calc_CentGreen(w, NUM_c, SigmaL_R, SigmaR_R, HCC, SCC, GC_R);
2422
2423 /* in case of advanced ones */
2424
2425 w.i = -w.i;
2426
2427 iside=0;
2428 TRAN_Calc_SurfGreen_direct(w,NUM_e[iside], H00_e[iside],H01_e[iside],
2429 S00_e[iside], S01_e[iside],
2430 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRL);
2431
2432 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRL, NUM_c, HCL, SCL, SigmaL_A);
2433
2434 iside=1;
2435 TRAN_Calc_SurfGreen_direct(w,NUM_e[iside], H00_e[iside],H01_e[iside],
2436 S00_e[iside], S01_e[iside],
2437 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRR);
2438
2439 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRR, NUM_c, HCR, SCR, SigmaR_A);
2440
2441 TRAN_Calc_CentGreen(w, NUM_c, SigmaL_A, SigmaR_A, HCC, SCC, GC_A);
2442
2443 w.i = -w.i;
2444
2445 /* calculation of transmission */
2446
2447 TRAN_Calc_OneTransmission(NUM_c, SigmaL_R, SigmaL_A, SigmaR_R, SigmaR_A, GC_R, GC_A, v1, v2 ,&value);
2448
2449 tran_transmission[k][iw].r = value.r;
2450 tran_transmission[k][iw].i = value.i;
2451
2452 if (PrintLevel){
2453 printf("k=%2d w.r=%6.3f w.i=%6.3f value.r=%15.12f value.i=%15.12f\n",k,w.r,w.i,value.r,value.i);
2454 }
2455
2456 if (SpinP_switch==3){
2457 tran_transmission[1][iw].r = value.r;
2458 tran_transmission[1][iw].i = value.i;
2459 }
2460
2461 } /* for k */
2462 } /* if ( iw>=0 ) */
2463 } /* iw */
2464
2465 /* MPI communication */
2466
2467 for (k=0; k<=1; k++) {
2468
2469 for (ID=0; ID<numprocs; ID++) {
2470 for (Miw=0; Miw<Miwmax ; Miw++) {
2471
2472 double v[2];
2473
2474 iw = iwIdx[ID][Miw];
2475
2476 if (2<=numprocs) MPI_Barrier(comm1);
2477
2478 v[0] = tran_transmission[k][iw].r;
2479 v[1] = tran_transmission[k][iw].i;
2480
2481 if (iw>0 && 2<=numprocs) {
2482 MPI_Bcast(v, 2, MPI_DOUBLE, ID, comm1);
2483 }
2484
2485 tran_transmission[k][iw].r = v[0];
2486 tran_transmission[k][iw].i = v[1];
2487
2488 }
2489 }
2490 }
2491
2492 /* freeing of arrays */
2493
2494 free(GC_R);
2495 free(GC_A);
2496 free(SigmaL_R);
2497 free(SigmaL_A);
2498 free(SigmaR_R);
2499 free(SigmaR_A);
2500
2501 free(GRR);
2502 free(GRL);
2503 free(v2);
2504 free(v1);
2505
2506 for (i=0;i<numprocs;i++) {
2507 free(iwIdx[i]);
2508 }
2509 free(iwIdx);
2510 }
2511
2512
2513
2514
2515
2516
2517
2518
2519
MTRAN_Current(MPI_Comm comm1,int numprocs,int myid,int SpinP_switch,double ChemP_e[2],double E_Temp,int NUM_c,int NUM_e[2],dcomplex * H00_e[2],dcomplex * S00_e[2],dcomplex * H01_e[2],dcomplex * S01_e[2],dcomplex * HCC,dcomplex * HCL,dcomplex * HCR,dcomplex * SCC,dcomplex * SCL,dcomplex * SCR,double tran_surfgreen_iteration_max,double tran_surfgreen_eps,double * current,double k2,double k3)2520 void MTRAN_Current(MPI_Comm comm1,
2521 int numprocs,
2522 int myid,
2523 int SpinP_switch,
2524 double ChemP_e[2],
2525 double E_Temp,
2526 int NUM_c,
2527 int NUM_e[2],
2528 dcomplex *H00_e[2],
2529 dcomplex *S00_e[2],
2530 dcomplex *H01_e[2],
2531 dcomplex *S01_e[2],
2532 dcomplex *HCC,
2533 dcomplex *HCL,
2534 dcomplex *HCR,
2535 dcomplex *SCC,
2536 dcomplex *SCL,
2537 dcomplex *SCR,
2538 double tran_surfgreen_iteration_max,
2539 double tran_surfgreen_eps,
2540 double *current,
2541 double k2, double k3)
2542 {
2543 dcomplex w;
2544 dcomplex *GRL,*GRR;
2545 dcomplex *GC_R,*GC_A;
2546 dcomplex *SigmaL_R,*SigmaL_A;
2547 dcomplex *SigmaR_R,*SigmaR_A;
2548 /*S MitsuakiKAWAMURA2*/
2549 dcomplex *v1,*v2, *Sinv;
2550 double kvec[2];
2551 /*E MitsuakiKAWAMURA2*/
2552 dcomplex value;
2553
2554 double Beta,xL,xR,fL,fR;
2555 double my_current[4];
2556 int iw,k,iside;
2557 int ID;
2558 int tag=99;
2559
2560 int **iwIdx;
2561 int Miw,Miwmax ;
2562 int i,j;
2563 MPI_Status status;
2564
2565 Beta = 1.0/kB/E_Temp;
2566
2567 /* allocation of arrays */
2568
2569 v1 = (dcomplex*) malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2570 v2 = (dcomplex*) malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2571
2572 GRL = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[0]*NUM_e[0]);
2573 GRR = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[1]*NUM_e[1]);
2574
2575 GC_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2576 GC_A = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2577 SigmaL_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2578 SigmaL_A = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2579 SigmaR_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2580 SigmaR_A = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2581
2582 /*S MitsuakiKAWAMURA2*/
2583 Sinv = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c);
2584 TRAN_Calc_Sinv(NUM_c, SCC, Sinv);
2585 kvec[0] = k2;
2586 kvec[1] = k3;
2587 /*E MitsuakiKAWAMURA2*/
2588
2589 /* parallel setup */
2590
2591 iwIdx=(int**)malloc(sizeof(int*)*numprocs);
2592 Miwmax = (Tran_current_num_step)/numprocs+1;
2593
2594 for (i=0; i<numprocs; i++) {
2595 iwIdx[i]=(int*)malloc(sizeof(int)*Miwmax);
2596 }
2597
2598 TRAN_Distribute_Node_Idx(0, Tran_current_num_step-1, numprocs, Miwmax,
2599 iwIdx); /* output */
2600
2601 /* parallel global iw 0:tran_transmission_energydiv-1 */
2602 /* parallel local Miw 0:Miwmax-1 */
2603 /* parallel variable iw=iwIdx[myid][Miw] */
2604
2605 for (k=0; k<=1; k++) my_current[k] = 0.0;
2606
2607 for (Miw=0; Miw<Miwmax ; Miw++) {
2608
2609 iw = iwIdx[myid][Miw];
2610
2611 if ( iw>=0 ) {
2612
2613 w.r = Tran_current_lower_bound + (double)iw*Tran_current_energy_step;
2614 w.i = Tran_current_im_energy;
2615
2616 /*****************************************************************
2617 Note that retarded and advanced Green functions and self energies
2618 are not conjugate comlex in case of the k-dependent case.
2619 **************************************************************/
2620
2621 /* in case of retarded ones */
2622
2623 iside=0;
2624 TRAN_Calc_SurfGreen_direct(w,NUM_e[iside], H00_e[iside],H01_e[iside],
2625 S00_e[iside], S01_e[iside],
2626 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRL);
2627
2628 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRL, NUM_c, HCL, SCL, SigmaL_R);
2629
2630 iside=1;
2631 TRAN_Calc_SurfGreen_direct(w,NUM_e[iside], H00_e[iside],H01_e[iside],
2632 S00_e[iside], S01_e[iside],
2633 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRR);
2634
2635 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRR, NUM_c, HCR, SCR, SigmaR_R);
2636
2637 TRAN_Calc_CentGreen(w, NUM_c, SigmaL_R, SigmaR_R, HCC, SCC, GC_R);
2638
2639 /* in case of advanced ones */
2640
2641 w.i = -w.i;
2642
2643 iside=0;
2644 TRAN_Calc_SurfGreen_direct(w,NUM_e[iside], H00_e[iside],H01_e[iside],
2645 S00_e[iside], S01_e[iside],
2646 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRL);
2647
2648 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRL, NUM_c, HCL, SCL, SigmaL_A);
2649
2650 iside=1;
2651 TRAN_Calc_SurfGreen_direct(w,NUM_e[iside], H00_e[iside],H01_e[iside],
2652 S00_e[iside], S01_e[iside],
2653 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRR);
2654
2655 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRR, NUM_c, HCR, SCR, SigmaR_A);
2656
2657 TRAN_Calc_CentGreen(w, NUM_c, SigmaL_A, SigmaR_A, HCC, SCC, GC_A);
2658
2659 w.i = -w.i;
2660
2661 /* calculation of transmission */
2662
2663 TRAN_Calc_OneTransmission(NUM_c, SigmaL_R, SigmaL_A, SigmaR_R, SigmaR_A, GC_R, GC_A, v1, v2 ,&value);
2664
2665 /* add the contribution */
2666
2667 xL = (w.r - ChemP_e[0])*Beta;
2668 fL = 1.0/(1.0 + exp(xL));
2669
2670 xR = (w.r - ChemP_e[1])*Beta;
2671 fR = 1.0/(1.0 + exp(xR));
2672
2673 my_current[0] -= (fL-fR)*value.r*Tran_current_energy_step/(2.0*PI); /* in atomic unit */
2674 my_current[1] -= (fL-fR)*value.r*Tran_current_energy_step/(2.0*PI);
2675
2676 /*S MitsuakiKAWAMURA2*/
2677 if (TRAN_CurrentDensity == 1)
2678 TRAN_Calc_CurrentDensity_NC(NUM_c, GC_R, SigmaL_R, SigmaR_R, VCC, Sinv, kvec, fL, fR,
2679 Tran_current_energy_step, JLocSym, JLocASym, RhoNL, Jmat);
2680 /*E MitsuakiKAWAMURA2*/
2681 } /* if ( iw>=0 ) */
2682 } /* iw */
2683
2684
2685 /* parallel communication */
2686
2687 for (k=0; k<=1; k++) {
2688
2689 if (2<=numprocs){
2690 MPI_Allreduce(&my_current[k], ¤t[k], 1, MPI_DOUBLE, MPI_SUM, comm1);
2691 }
2692 else{
2693 current[k] = my_current[k];
2694 }
2695
2696 /**********************************************************
2697
2698 Current:
2699
2700 convert the unit from a.u to ampere
2701
2702 The unit of current is given by eEh/bar{h} in
2703 atomic unit, where e is the elementary charge,
2704 Eh Hartree, and bar{h} h/{2pi} given by
2705
2706 e = 1.60217653 * 10^{-19} C
2707 Eh = 4.35974417 * 10^{-18} J
2708 bar{h} = 1.05457168 * 10^{-34} J s
2709
2710 Therefore,
2711
2712 1 a.u.
2713 = 1.60217653 * 10^{-19} C * 4.35974417 * 10^{-18} J
2714 / (1.05457168 * 10^{-34} J s )
2715 = 6.6236178 * 10^{-3} [Cs^{-1}=ampere]
2716
2717 Electric potential:
2718
2719 convert the unit from a.u to volt
2720
2721 The unit of electric potential is given by Eh/e in
2722 atomic unit.
2723
2724 Therefore,
2725
2726 1 a.u.
2727 = 4.35974417 * 10^{-18} J/ (1.60217653 * 10^{-19} C)
2728 = 27.21138 [JC^{-1}=volt]
2729
2730 This allows us to consider that the difference of 1 eV in
2731 the chemical potential of leads can be regarded as 1 volt.
2732 *********************************************************/
2733
2734 current[k] *= 0.0066236178;
2735
2736 }
2737
2738
2739 /* freeing of arrays */
2740
2741 free(GC_R);
2742 free(GC_A);
2743 free(SigmaL_R);
2744 free(SigmaL_A);
2745 free(SigmaR_R);
2746 free(SigmaR_A);
2747
2748 free(GRR);
2749 free(GRL);
2750 free(v2);
2751 free(v1);
2752
2753 for (i=0;i<numprocs;i++) {
2754 free(iwIdx[i]);
2755 }
2756 free(iwIdx);
2757 }
2758
2759
2760
2761
2762
2763
2764
2765
2766
2767
2768
2769
2770
2771
2772
2773
2774
2775
2776
2777
2778
MTRAN_Output_Transmission(MPI_Comm comm1,int ID,char * fname,double k2,double k3,int SpinP_switch,int tran_transmission_energydiv,double tran_transmission_energyrange[3],dcomplex ** tran_transmission)2779 void MTRAN_Output_Transmission(
2780 MPI_Comm comm1,
2781 int ID,
2782 char *fname,
2783 double k2,
2784 double k3,
2785 int SpinP_switch,
2786 int tran_transmission_energydiv,
2787 double tran_transmission_energyrange[3],
2788 dcomplex **tran_transmission
2789 )
2790 {
2791 int iw,k;
2792 dcomplex w;
2793 int myid;
2794 FILE *fp;
2795
2796 MPI_Comm_rank(comm1,&myid);
2797 if (myid!=ID) { return; }
2798
2799 printf(" %s\n",fname);fflush(stdout);
2800
2801 if ( ( fp =fopen(fname,"w") )== NULL ) {
2802 printf("\ncannot open file to write transmission\n");fflush(stdout);
2803 exit(0);
2804 }
2805
2806 fprintf(fp,"#/************************************************************/\n");
2807 fprintf(fp,"#\n");
2808 fprintf(fp,"# Current:\n");
2809 fprintf(fp,"#\n");
2810 fprintf(fp,"# convert the unit from a.u to ampere\n");
2811 fprintf(fp,"#\n");
2812 fprintf(fp,"# The unit of current is given by eEh/bar{h} in\n");
2813 fprintf(fp,"# atomic unit, where e is the elementary charge,\n");
2814 fprintf(fp,"# Eh Hartree, and bar{h} h/{2pi} given by\n");
2815 fprintf(fp,"#\n");
2816 fprintf(fp,"# e = 1.60217653 * 10^{-19} C\n");
2817 fprintf(fp,"# Eh = 4.35974417 * 10^{-18} J\n");
2818 fprintf(fp,"# bar{h} = 1.05457168 * 10^{-34} J s\n");
2819 fprintf(fp,"#\n");
2820 fprintf(fp,"# Therefore,\n");
2821 fprintf(fp,"#\n");
2822 fprintf(fp,"# 1 a.u.\n");
2823 fprintf(fp,"# = 1.60217653 * 10^{-19} C * 4.35974417 * 10^{-18} J\n");
2824 fprintf(fp,"# / (1.05457168 * 10^{-34} J s )\n");
2825 fprintf(fp,"# = 6.6236178 * 10^{-3} [Cs^{-1}=ampere]\n");
2826 fprintf(fp,"#\n");
2827 fprintf(fp,"# Electric potential:\n");
2828 fprintf(fp,"#\n");
2829 fprintf(fp,"# convert the unit from a.u to volt\n");
2830 fprintf(fp,"#\n");
2831 fprintf(fp,"# The unit of electric potential is given by Eh/e in\n");
2832 fprintf(fp,"# atomic unit.\n");
2833 fprintf(fp,"#\n");
2834 fprintf(fp,"# Therefore,\n");
2835 fprintf(fp,"#\n");
2836 fprintf(fp,"# 1 a.u.\n");
2837 fprintf(fp,"# = 4.35974417 * 10^{-18} J/ (1.60217653 * 10^{-19} C)\n");
2838 fprintf(fp,"# = 27.21138 [JC^{-1}=volt]\n");
2839 fprintf(fp,"#\n");
2840 fprintf(fp,"# This allows us to consider that the difference of 1 eV in\n");
2841 fprintf(fp,"# the chemical potential of leads can be regarded as 1 volt.\n");
2842 fprintf(fp,"#\n");
2843 fprintf(fp,"#***********************************************************/\n");
2844 fprintf(fp,"#\n");
2845 fprintf(fp,"# k2= %8.4f k3= %8.4f\n",k2,k3);
2846 fprintf(fp,"# SpinP_switch= %d\n",SpinP_switch);
2847 fprintf(fp,"# Chemical potential (eV) Left, Right Leads=% 9.6e % 9.6e\n",
2848 ChemP_e[0]*eV2Hartree,ChemP_e[1]*eV2Hartree);
2849 fprintf(fp,"# Chemical potential (Hartree) Left, Right Leads=% 9.6e % 9.6e\n",
2850 ChemP_e[0],ChemP_e[1]);
2851 fprintf(fp,"# diff Chemical potential (eV) =% 9.6e\n",
2852 (ChemP_e[0]-ChemP_e[1])*eV2Hartree);
2853 fprintf(fp,"# diff Chemical potential (Hartree)=% 9.6e\n",
2854 ChemP_e[0]-ChemP_e[1]);
2855 fprintf(fp,"# tran_transmission_energydiv= %d\n", tran_transmission_energydiv);
2856 fprintf(fp,"#\n");
2857 fprintf(fp,"# iw w.real(au) w.imag(au) w.real(eV) w.imag(eV) trans.real trans.imag\n");
2858 fprintf(fp,"#\n");
2859
2860 for (iw=0; iw<tran_transmission_energydiv; iw++){
2861
2862 /* The zero energy is the chemical potential of the left lead. */
2863
2864 w.r = tran_transmission_energyrange[0] + ChemP_e[0] - ChemP_e[0]
2865 +
2866 (tran_transmission_energyrange[1]-tran_transmission_energyrange[0])*
2867 (double)iw/(tran_transmission_energydiv-1);
2868
2869 w.i = tran_transmission_energyrange[2];
2870
2871 fprintf(fp,"%3d % 9.6e % 9.6e % 9.6e % 9.6e % 9.6e % 9.6e\n",
2872 iw,
2873 w.r,
2874 w.i,
2875 w.r*eV2Hartree, w.i*eV2Hartree,
2876 tran_transmission[0][iw].r,
2877 tran_transmission[0][iw].i);
2878
2879 } /* iw */
2880
2881 fclose(fp);
2882 }
2883
2884
MTRAN_Output_Current(MPI_Comm comm1,char * fname,int TRAN_TKspace_grid2,int TRAN_TKspace_grid3,int SpinP_switch,double *** current)2885 void MTRAN_Output_Current(
2886 MPI_Comm comm1,
2887 char *fname,
2888 int TRAN_TKspace_grid2,
2889 int TRAN_TKspace_grid3,
2890 int SpinP_switch,
2891 double ***current
2892 )
2893 {
2894 int iw,i2,i3;
2895 double k2,k3;
2896 double crt0,crt1;
2897 dcomplex w;
2898 int myid;
2899 FILE *fp;
2900
2901 MPI_Comm_rank(comm1,&myid);
2902 if (myid!=Host_ID) { return; }
2903
2904 printf(" %s\n",fname);
2905
2906 if ( ( fp =fopen(fname,"w") )== NULL ) {
2907 printf("\ncannot open file to write current\n");
2908 printf("write current to stdout\n");
2909 exit(0);
2910 }
2911
2912 /* total current */
2913
2914 crt0 = 0.0;
2915 crt1 = 0.0;
2916 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
2917 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
2918 crt0 += current[i2][i3][0];
2919 crt1 += current[i2][i3][1];
2920 }
2921 }
2922
2923 crt0 /= (double)(TRAN_TKspace_grid2*TRAN_TKspace_grid3);
2924 crt1 /= (double)(TRAN_TKspace_grid2*TRAN_TKspace_grid3);
2925
2926 fprintf(fp,"#/************************************************************/\n");
2927 fprintf(fp,"#\n");
2928 fprintf(fp,"# Current:\n");
2929 fprintf(fp,"#\n");
2930 fprintf(fp,"# convert the unit from a.u to ampere\n");
2931 fprintf(fp,"#\n");
2932 fprintf(fp,"# The unit of current is given by eEh/bar{h} in\n");
2933 fprintf(fp,"# atomic unit, where e is the elementary charge,\n");
2934 fprintf(fp,"# Eh Hartree, and bar{h} h/{2pi} given by\n");
2935 fprintf(fp,"#\n");
2936 fprintf(fp,"# e = 1.60217653 * 10^{-19} C\n");
2937 fprintf(fp,"# Eh = 4.35974417 * 10^{-18} J\n");
2938 fprintf(fp,"# bar{h} = 1.05457168 * 10^{-34} J s\n");
2939 fprintf(fp,"#\n");
2940 fprintf(fp,"# Therefore,\n");
2941 fprintf(fp,"#\n");
2942 fprintf(fp,"# 1 a.u.\n");
2943 fprintf(fp,"# = 1.60217653 * 10^{-19} C * 4.35974417 * 10^{-18} J\n");
2944 fprintf(fp,"# / (1.05457168 * 10^{-34} J s )\n");
2945 fprintf(fp,"# = 6.6236178 * 10^{-3} [Cs^{-1}=ampere]\n");
2946 fprintf(fp,"#\n");
2947 fprintf(fp,"# Electric potential:\n");
2948 fprintf(fp,"#\n");
2949 fprintf(fp,"# convert the unit from a.u to volt\n");
2950 fprintf(fp,"#\n");
2951 fprintf(fp,"# The unit of electric potential is given by Eh/e in\n");
2952 fprintf(fp,"# atomic unit.\n");
2953 fprintf(fp,"#\n");
2954 fprintf(fp,"# Therefore,\n");
2955 fprintf(fp,"#\n");
2956 fprintf(fp,"# 1 a.u.\n");
2957 fprintf(fp,"# = 4.35974417 * 10^{-18} J/ (1.60217653 * 10^{-19} C)\n");
2958 fprintf(fp,"# = 27.21138 [JC^{-1}=volt]\n");
2959 fprintf(fp,"#\n");
2960 fprintf(fp,"# This allows us to consider that the difference of 1 eV in\n");
2961 fprintf(fp,"# the chemical potential of leads can be regarded as 1 volt.\n");
2962 fprintf(fp,"#\n");
2963 fprintf(fp,"#***********************************************************/\n");
2964 fprintf(fp,"#\n");
2965 fprintf(fp,"# SpinP_switch= %d\n",SpinP_switch);
2966 fprintf(fp,"# Chemical potential (eV) Left, Right Leads=% 9.6e % 9.6e\n",
2967 ChemP_e[0]*eV2Hartree,ChemP_e[1]*eV2Hartree);
2968 fprintf(fp,"# Chemical potential (Hartree) Left, Right Leads=% 9.6e % 9.6e\n",
2969 ChemP_e[0],ChemP_e[1]);
2970 fprintf(fp,"# diff Chemical potential (eV) =% 9.6e\n",
2971 (ChemP_e[0]-ChemP_e[1])*eV2Hartree);
2972 fprintf(fp,"# diff Chemical potential (Hartree)=% 9.6e\n",
2973 ChemP_e[0]-ChemP_e[1]);
2974 fprintf(fp,"# Corresponding bias voltage (Volt)=% 9.6e\n",
2975 (ChemP_e[0]-ChemP_e[1])*27.21138 );
2976 fprintf(fp,"# average current (ampere)=% 9.6e \n",crt0);
2977 fprintf(fp,"#\n");
2978 fprintf(fp,"# i2 i3 k2 k3 current (ampere) \n");
2979 fprintf(fp,"#\n");
2980
2981 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
2982 k2 = -0.5 + (2.0*(double)i2+1.0)/(2.0*(double)TRAN_TKspace_grid2);
2983 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
2984 k3 = -0.5 + (2.0*(double)i3+1.0)/(2.0*(double)TRAN_TKspace_grid3);
2985
2986 fprintf(fp,"%3d %3d %8.4f %8.4f % 9.6e \n",
2987 i2,i3,k2,k3,current[i2][i3][0]);
2988 }
2989 }
2990
2991 fprintf(fp,"\n\n");
2992
2993 fclose(fp);
2994 }
2995
2996
2997
MTRAN_Output_Conductance(MPI_Comm comm1,char * fname,int * T_k_ID,int T_knum,int TRAN_TKspace_grid2,int TRAN_TKspace_grid3,int * T_IGrids2,int * T_IGrids3,double * T_KGrids2,double * T_KGrids3,int SpinP_switch,int tran_transmission_energydiv,double tran_transmission_energyrange[3],dcomplex **** tran_transmission)2998 void MTRAN_Output_Conductance(
2999 MPI_Comm comm1,
3000 char *fname,
3001 int *T_k_ID,
3002 int T_knum,
3003 int TRAN_TKspace_grid2,
3004 int TRAN_TKspace_grid3,
3005 int *T_IGrids2,
3006 int *T_IGrids3,
3007 double *T_KGrids2,
3008 double *T_KGrids3,
3009 int SpinP_switch,
3010 int tran_transmission_energydiv,
3011 double tran_transmission_energyrange[3],
3012 dcomplex ****tran_transmission
3013 )
3014 {
3015 int spin,k,i2,i3,iw,iw0,iw1;
3016 int po,ii2,ii3,ID;
3017 double k2,k3,Av_ChemP;
3018 double ***conductance;
3019 double w0,w1,e0,e1;
3020 double wm,wn,cond[2];
3021 int myid;
3022 FILE *fp;
3023
3024 MPI_Comm_rank(comm1,&myid);
3025 if (myid==Host_ID) printf(" %s\n",fname);
3026
3027 /***********************************************
3028 allocation of arrays
3029 ***********************************************/
3030
3031 conductance = (double***)malloc(sizeof(double**)*TRAN_TKspace_grid2);
3032 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
3033 conductance[i2] = (double**)malloc(sizeof(double*)*TRAN_TKspace_grid3);
3034 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
3035 conductance[i2][i3] = (double*)malloc(sizeof(double)*2);
3036 for (spin=0; spin<2; spin++){
3037 conductance[i2][i3][spin] = 0.0;
3038 }
3039 }
3040 }
3041
3042 /***********************************************
3043 find the conductance at each spin and k-point
3044 ***********************************************/
3045
3046 Av_ChemP = 0.50*(ChemP_e[0] + ChemP_e[1]);
3047
3048 po = 0;
3049
3050 for (iw=0; iw<tran_transmission_energydiv; iw++){
3051
3052 /* The zero energy is the chemical potential of the left lead. */
3053
3054 w0 = w1;
3055 w1 = tran_transmission_energyrange[0] + ChemP_e[0] - ChemP_e[0]
3056 +
3057 (tran_transmission_energyrange[1]-tran_transmission_energyrange[0])*
3058 (double)iw/(tran_transmission_energydiv-1);
3059
3060 /* find iw0 which is the index near the average chemical potential */
3061
3062 if (Av_ChemP<(w1+ChemP_e[0]) && po==0){
3063
3064 po = 1;
3065 iw0 = iw - 1;
3066 iw1 = iw;
3067 e0 = w0 + ChemP_e[0];
3068 e1 = w1 + ChemP_e[0];
3069 }
3070 }
3071
3072 wm = (Av_ChemP - e0)/(e1-e0);
3073 wn = (e1 - Av_ChemP)/(e1-e0);
3074
3075 /* calculate the conductance by a linear interpolation */
3076
3077 for (spin=0; spin<=1; spin++){
3078 for (k=0; k<T_knum; k++){
3079
3080 ID = T_k_ID[k];
3081
3082 k2 = T_KGrids2[k];
3083 k3 = T_KGrids3[k];
3084 i2 = T_IGrids2[k];
3085 i3 = T_IGrids3[k];
3086
3087 if (myid==ID){
3088
3089 conductance[i2][i3][spin] = tran_transmission[i2][i3][spin][iw0].r*wn
3090 + tran_transmission[i2][i3][spin][iw1].r*wm;
3091
3092 }
3093
3094 MPI_Bcast(conductance[i2][i3], 2, MPI_DOUBLE, ID, comm1);
3095
3096 }
3097 }
3098
3099 /***********************************************
3100 calculate average conductance
3101 ***********************************************/
3102
3103 for (spin=0; spin<=1; spin++){
3104
3105 cond[spin] = 0.0;
3106 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
3107 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
3108 cond[spin] += conductance[i2][i3][spin];
3109 }
3110 }
3111
3112 cond[spin] = cond[spin]/(double)(TRAN_TKspace_grid2*TRAN_TKspace_grid3);
3113 }
3114
3115 /***********************************************
3116 save conductance
3117 ***********************************************/
3118
3119 if (myid==Host_ID){
3120
3121 if ( ( fp =fopen(fname,"w") )== NULL ) {
3122 printf("\ncannot open file to write conductance\n");fflush(stdout);
3123 exit(0);
3124 }
3125
3126 fprintf(fp,"#/************************************************************/\n");
3127 fprintf(fp,"#\n");
3128 fprintf(fp,"# Conductance:\n");
3129 fprintf(fp,"#\n");
3130 fprintf(fp,"# G0 = e^2/h\n");
3131 fprintf(fp,"# = (1.60217653 * 10^{-19})^2/(6.626076 * 10^{-34})\n");
3132 fprintf(fp,"# = 3.87404194169 * 10^{-5} [C^2 J^{-1} s^{-1}]\n");
3133 fprintf(fp,"# note that\n");
3134 fprintf(fp,"# e = 1.60217653 * 10^{-19} C\n");
3135 fprintf(fp,"# h = 6.626076 * 10^{-34} J s\n");
3136 fprintf(fp,"#***********************************************************/\n");
3137 fprintf(fp,"#\n");
3138
3139 fprintf(fp,"# SpinP_switch= %d\n",SpinP_switch);
3140 fprintf(fp,"# Chemical potential (eV) Left, Right Leads=% 9.6e % 9.6e\n",
3141 ChemP_e[0]*eV2Hartree,ChemP_e[1]*eV2Hartree);
3142 fprintf(fp,"# Chemical potential (Hartree) Left, Right Leads=% 9.6e % 9.6e\n",
3143 ChemP_e[0],ChemP_e[1]);
3144 fprintf(fp,"# diff Chemical potential (eV) =% 9.6e\n",
3145 (ChemP_e[0]-ChemP_e[1])*eV2Hartree);
3146 fprintf(fp,"# diff Chemical potential (Hartree)=% 9.6e\n",
3147 ChemP_e[0]-ChemP_e[1]);
3148 fprintf(fp,"# Corresponding bias voltage (Volt)=% 9.6e\n",
3149 (ChemP_e[0]-ChemP_e[1])*27.21138 );
3150 fprintf(fp,"# average conductance (G0)=% 9.6e \n",cond[0]);
3151 fprintf(fp,"#\n");
3152 fprintf(fp,"# i2 i3 k2 k3 conductance (G0)\n");
3153 fprintf(fp,"#\n");
3154
3155 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
3156 k2 = -0.5 + (2.0*(double)i2+1.0)/(2.0*(double)TRAN_TKspace_grid2);
3157 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
3158 k3 = -0.5 + (2.0*(double)i3+1.0)/(2.0*(double)TRAN_TKspace_grid3);
3159
3160 fprintf(fp,"%3d %3d %8.4f %8.4f % 9.6e \n",
3161 i2,i3,k2,k3,conductance[i2][i3][0]);
3162 }
3163 }
3164
3165 fclose(fp);
3166 }
3167
3168 /***********************************************
3169 freeing of arrays
3170 ***********************************************/
3171
3172 for (i2=0; i2<TRAN_TKspace_grid2; i2++){
3173 for (i3=0; i3<TRAN_TKspace_grid3; i3++){
3174 free(conductance[i2][i3]);
3175 }
3176 free(conductance[i2]);
3177 }
3178 free(conductance);
3179 }
3180
3181
3182
3183
3184
3185
3186
3187
3188
3189
3190
MTRAN_Input(MPI_Comm comm1,int argc,char * fname,double ChemP_e[2],int * TRAN_TKspace_grid2,int * TRAN_TKspace_grid3,int * SpinP_switch,double * E_Temp,int * tran_surfgreen_iteration_max,double * tran_surfgreen_eps,int * tran_transmission_on,double tran_transmission_energyrange[3],int * tran_transmission_energydiv,dcomplex **** (* tran_transmission))3191 void MTRAN_Input(MPI_Comm comm1,
3192 int argc,
3193 char *fname,
3194 double ChemP_e[2],
3195 int *TRAN_TKspace_grid2,
3196 int *TRAN_TKspace_grid3,
3197 int *SpinP_switch,
3198 double *E_Temp,
3199 int *tran_surfgreen_iteration_max,
3200 double *tran_surfgreen_eps,
3201 int *tran_transmission_on,
3202 double tran_transmission_energyrange[3],
3203 int *tran_transmission_energydiv,
3204 dcomplex ****(*tran_transmission)
3205 )
3206 {
3207 int i,i2,i3,k,po;
3208 int side0,side1;
3209 double r_vec[20];
3210 int i_vec[20];
3211 int i_vec2[20];
3212 char *s_vec[20];
3213 double Beta,f0,f1;
3214 double x,x0,x1,tmpx0,tmpx1;
3215 double Av_ChemP;
3216 int myid;
3217 /* S MitsuakiKAWAMURA */
3218 char buf[BUFSIZE];
3219 FILE *fp;
3220 /* E MitsuakiKAWAMURA*/
3221
3222 MPI_Comm_rank(comm1,&myid);
3223
3224 /* open the input file */
3225
3226 if (input_open(fname)==0){
3227 MPI_Finalize();
3228 exit(0);
3229 }
3230
3231 s_vec[0]="Off"; s_vec[1]="On"; s_vec[2]="NC";
3232 i_vec[0]=0 ; i_vec[1]=1 ; i_vec[2]=3;
3233 input_string2int("scf.SpinPolarization", SpinP_switch, 3, s_vec,i_vec);
3234 input_double("scf.ElectronicTemperature", E_Temp, 300);
3235 /* chage the unit from K to a.u. */
3236 *E_Temp = *E_Temp/eV2Hartree;
3237
3238 Beta = 1.0/kB/(*E_Temp);
3239
3240 input_int( "NEGF.Surfgreen.iterationmax", tran_surfgreen_iteration_max, 600);
3241 input_double("NEGF.Surfgreen.convergeeps", tran_surfgreen_eps, 1.0e-12);
3242
3243 /**** k-points parallel to the layer, which are used for the transmission calc. ****/
3244
3245 i_vec2[0]=1;
3246 i_vec2[1]=1;
3247 input_intv("NEGF.tran.Kgrid",2,i_vec,i_vec2);
3248 *TRAN_TKspace_grid2 = i_vec[0];
3249 *TRAN_TKspace_grid3 = i_vec[1];
3250
3251 if (*TRAN_TKspace_grid2<=0){
3252 printf("NEGF.tran.Kgrid should be over 1\n");
3253 MPI_Finalize();
3254 exit(1);
3255 }
3256
3257 if (*TRAN_TKspace_grid3<=0){
3258 printf("NEGF.tran.Kgrid should be over 1\n");
3259 MPI_Finalize();
3260 exit(1);
3261 }
3262
3263 input_logical("NEGF.tran.on",tran_transmission_on,1); /* default=on */
3264
3265 if (tran_transmission_on) {
3266 i=0;
3267 r_vec[i++] = -10.0;
3268 r_vec[i++] = 10.0;
3269 r_vec[i++] = 1.0e-3;
3270 input_doublev("NEGF.tran.energyrange",i, tran_transmission_energyrange, r_vec); /* in eV */
3271 /* change the unit from eV to Hartree */
3272 tran_transmission_energyrange[0] /= eV2Hartree;
3273 tran_transmission_energyrange[1] /= eV2Hartree;
3274 tran_transmission_energyrange[2] /= eV2Hartree;
3275
3276 input_int("NEGF.tran.energydiv",tran_transmission_energydiv,200);
3277
3278 *tran_transmission = (dcomplex****)malloc(sizeof(dcomplex***)*(*TRAN_TKspace_grid2));
3279 for (i2=0; i2<(*TRAN_TKspace_grid2); i2++){
3280 (*tran_transmission)[i2] = (dcomplex***)malloc(sizeof(dcomplex**)*(*TRAN_TKspace_grid3));
3281 for (i3=0; i3<(*TRAN_TKspace_grid3); i3++){
3282 (*tran_transmission)[i2][i3] = (dcomplex**)malloc(sizeof(dcomplex*)*3);
3283 for (i=0; i<3; i++) {
3284 (*tran_transmission)[i2][i3][i] = (dcomplex*)malloc(sizeof(dcomplex)*(*tran_transmission_energydiv));
3285 }
3286 }
3287 }
3288 }
3289 else {
3290 tran_transmission_energydiv=0;
3291 *tran_transmission=NULL;
3292 }
3293
3294 /* set Order_Lead_Side */
3295
3296 if (ChemP_e[0]<ChemP_e[1]){
3297 Order_Lead_Side[0] = 0;
3298 Order_Lead_Side[1] = 1;
3299 }
3300 else {
3301 Order_Lead_Side[0] = 1;
3302 Order_Lead_Side[1] = 0;
3303 }
3304
3305 input_double("NEGF.current.energy.step", &Tran_current_energy_step, 0.01); /* in eV */
3306 if (Tran_current_energy_step<0.0) {
3307 printf("NEGF.bias.neq.energy.step should be positive.\n");
3308 MPI_Finalize();
3309 exit(1);
3310 }
3311
3312 /* change the unit from eV to Hartree */
3313 Tran_current_energy_step /= eV2Hartree;
3314
3315 input_double("NEGF.current.im.energy", &Tran_current_im_energy, tran_transmission_energyrange[2]*eV2Hartree); /* in eV */
3316
3317 /* change the unit from eV to Hartree */
3318 Tran_current_im_energy /= eV2Hartree;
3319
3320 input_double("NEGF.current.cutoff", &Tran_current_cutoff, 1.0e-8); /* dimensionless */
3321
3322 side0 = Order_Lead_Side[0];
3323 side1 = Order_Lead_Side[1];
3324
3325 Av_ChemP = 0.5*(ChemP_e[side0] + ChemP_e[side1]);
3326 x = Av_ChemP;
3327
3328 po = 0;
3329 do {
3330
3331 f0 = 1.0/(1.0+exp((x-ChemP_e[side0])*Beta));
3332 f1 = 1.0/(1.0+exp((x-ChemP_e[side1])*Beta));
3333
3334 if ( fabs(f1-f0)<Tran_current_cutoff ){
3335
3336 po = 1;
3337 x1 = x;
3338 };
3339
3340 x += Tran_current_energy_step*0.1;
3341
3342 } while (po==0);
3343
3344 x0 = Av_ChemP - (x1 - Av_ChemP);
3345 Tran_current_lower_bound = x0;
3346 Tran_current_num_step = (int)((x1-x0)/(double)Tran_current_energy_step);
3347
3348 if (Tran_current_num_step<10 && myid==0){
3349 printf("NEGF.current.energy.step %8.4e seems to be large for the calculation of current in the bias voltage %8.4e\n",
3350 Tran_current_energy_step*eV2Hartree,(ChemP_e[side1]-ChemP_e[side0])*eV2Hartree);
3351 printf("The recommended Tran.current.energy.step is %8.4e (eV).\n",
3352 (x1-x0)/50.0*eV2Hartree);
3353 }
3354
3355 /* S MitsuakiKAWAMURA */
3356
3357 /********************************************************
3358 Parameters for the Transmission EigenChannel
3359 *********************************************************/
3360
3361 /* Check Whether eigenchannel is computed or not */
3362
3363 input_logical("NEGF.tran.channel", &TRAN_Channel, 1);
3364
3365 if (TRAN_Channel == 1){
3366
3367 /* The number of k points where eigchannels are computed */
3368
3369 input_int("NEGF.Channel.Nkpoint", &TRAN_Channel_Nkpoint, 1);
3370 if (TRAN_Channel_Nkpoint < 1){
3371 printf("NEGF.Channel.Nkpoint must be >= 1 \n");
3372 MPI_Finalize();
3373 exit(1);
3374 }
3375
3376 /* Allocate and read k points where eigenchannel are computed */
3377
3378 TRAN_Channel_kpoint = (double**)malloc(sizeof(double*)*TRAN_Channel_Nkpoint);
3379 for (i = 0; i < TRAN_Channel_Nkpoint; i++){
3380 TRAN_Channel_kpoint[i] = (double*)malloc(sizeof(double) * 2);
3381 } /* for (i = 0; i<(TRAN_Channel_Nkpoint + 1); i++) */
3382
3383 if (fp = input_find("<NEGF.Channel.kpoint")) {
3384 for (i = 0; i < TRAN_Channel_Nkpoint - 1; i++){
3385 fgets(buf, BUFSIZE, fp);
3386 sscanf(buf, "%lf %lf",
3387 &TRAN_Channel_kpoint[i][0], &TRAN_Channel_kpoint[i][1]);
3388 } /* for (i = 0; i<TRAN_Channel_Nkpoint; i++) */
3389 fscanf(fp, "%lf %lf",
3390 &TRAN_Channel_kpoint[i][0], &TRAN_Channel_kpoint[i][1]);
3391
3392 if (!input_last("NEGF.Channel.kpoint>")) {
3393 /* format error */
3394 printf("Format error for NEGF.Channel.kpoint\n");
3395 MPI_Finalize();
3396 exit(1);
3397 } /* if (!input_last("NEGE.Channel.kpoint>")) */
3398
3399 /* k points are read as the fractional coordinate */
3400 } /* if (fp = input_find("<NEGF.Channel.kpoint")) */
3401 else if (TRAN_Channel_Nkpoint == 1){
3402 /* In default, the first k point is 0.0 0.0 */
3403 TRAN_Channel_kpoint[0][0] = 0.0;
3404 TRAN_Channel_kpoint[0][1] = 0.0;
3405 } /* (TRAN_Channel_Nkpoint == 1) */
3406 else if (TRAN_Channel_Nkpoint > 1){
3407 printf("NEGF.Channel.Nkpoint > 1 and NEGF.Channel.kpoint is NOT found \n");
3408 MPI_Finalize();
3409 exit(1);
3410 } /* else if (TRAN_Channel_Nkpoint > 1) */
3411
3412 if (myid == Host_ID) for (i = 0; i < TRAN_Channel_Nkpoint; i++)
3413 printf(" TRAN_Channel_kpoint %2d %10.6f %10.6f\n",
3414 i, TRAN_Channel_kpoint[i][0], TRAN_Channel_kpoint[i][1]);
3415
3416 /* The number of energy where eigenchannels are computed */
3417
3418 input_int("NEGF.Channel.Nenergy", &TRAN_Channel_Nenergy, 1);
3419 if (TRAN_Channel_Nenergy < 1){
3420 printf("NEGF.Channel.Nenergy must be >= 1 \n");
3421 MPI_Finalize();
3422 exit(1);
3423 }
3424
3425 /* Allocate and read energies where eigenchannels are computed */
3426
3427 TRAN_Channel_energy = (double*)malloc(sizeof(double)*(TRAN_Channel_Nenergy));
3428
3429 if (fp = input_find("<NEGF.Channel.Energy")){
3430 for (i = 0; i<TRAN_Channel_Nenergy - 1; i++){
3431 fgets(buf, BUFSIZE, fp);
3432 sscanf(buf, "%lf", &TRAN_Channel_energy[i]);
3433 } /* for (i = 0; i<TRAN_Channel_Nenergy; i++) */
3434 fscanf(fp, "%lf", &TRAN_Channel_energy[i]);
3435
3436 if (!input_last("NEGF.Channel.Energy>")){
3437 /* format error */
3438 printf("Format error for NEGF.Channel.energy\n");
3439 MPI_Finalize();
3440 exit(1);
3441 } /* if (!input_last("NEGE.Channel.energy>")) */
3442
3443 } /* if (fp = input_find("<NEGF.Channel.Energy")) */
3444 else if (TRAN_Channel_Nenergy == 1){
3445 /* In default, the first energy is 0.0 (E_F of the left lead) */
3446 TRAN_Channel_energy[0] = 0.0;
3447 } /* else if (TRAN_Channel_Nenergy == 1) */
3448 else if (TRAN_Channel_Nenergy > 1){
3449 printf("NEGF.Channel.Nenergy > 1 and NEGF.Channel.energy is NOT found \n");
3450 MPI_Finalize();
3451 exit(1);
3452 } /* else if (TRAN_Channel_Nenergy > 1) */
3453
3454 for (i = 0; i < TRAN_Channel_Nenergy; i++){
3455 if (myid == Host_ID)
3456 printf(" TRAN_Channel_energy %2d %10.6f eV\n", i, TRAN_Channel_energy[i]);
3457 TRAN_Channel_energy[i] /= eV2Hartree;
3458 }
3459
3460 /* The number of eigenchannel-cube file */
3461
3462 input_int("NEGF.Channel.Num", &TRAN_Channel_Num, 10);
3463 if (myid == Host_ID)
3464 printf(" TRAN_Channel_Num %d \n", TRAN_Channel_Num);
3465 if (TRAN_Channel_Num < 0){
3466 printf("NEGF.Channel.Num must be >= 0 \n");
3467 MPI_Finalize();
3468 exit(1);
3469 }
3470 }
3471 /***********************************************************************
3472 Current density calculation flag
3473 ***********************************************************************/
3474
3475 input_logical("NEGF.tran.CurrentDensity", &TRAN_CurrentDensity, 1);
3476 input_logical("NEGF.OffDiagonalCurrent", &TRAN_OffDiagonalCurrent, 0);
3477 if (TRAN_OffDiagonalCurrent == 1 && TRAN_CurrentDensity != 1) {
3478 TRAN_CurrentDensity = 1;
3479 if (myid == Host_ID)
3480 printf("Since NEGF.OffDiagonalCurrent is on, NEGF.tran.CurrentDensity is automatically set to on.\n");
3481 }
3482
3483 /* E MitsuakiKAWAMURA */
3484 /* print information */
3485
3486 if (myid==0){
3487 printf("\n");
3488 printf("Parameters for the calculation of the current\n");
3489 printf(" lower bound: %15.12f (eV)\n",x0*eV2Hartree);
3490 printf(" upper bound: %15.12f (eV)\n",x1*eV2Hartree);
3491 printf(" energy step: %15.12f (eV)\n",Tran_current_energy_step*eV2Hartree);
3492 printf(" imginary energy %15.12f (eV)\n",Tran_current_im_energy*eV2Hartree);
3493 printf(" number of steps: %3d \n",Tran_current_num_step);
3494 }
3495
3496 input_close();
3497
3498 }
3499
3500
MTRAN_Input_Sys(int argc,char * file,char * filepath,char * filename,int * tran_interpolate,char * interpolate_filename1,char * interpolate_filename2)3501 void MTRAN_Input_Sys(int argc,char *file, char *filepath, char *filename,
3502 int *tran_interpolate,
3503 char *interpolate_filename1,
3504 char *interpolate_filename2)
3505 {
3506 double r_vec[20];
3507 double r_vec1[20];
3508
3509 if (input_open(file)==0){
3510 MPI_Finalize();
3511 exit(0);
3512 }
3513
3514 input_string("System.CurrrentDirectory",filepath,"./");
3515 input_string("System.Name",filename,"default");
3516
3517 /* check whether the interpolation is made or not */
3518 input_logical("NEGF.tran.interpolate",tran_interpolate,0); /* default=off */
3519 input_string("NEGF.tran.interpolate.file1",interpolate_filename1,"file1.tranb1");
3520 input_string("NEGF.tran.interpolate.file2",interpolate_filename2,"file2.tranb2");
3521
3522 r_vec[0] = 1.0; r_vec[1] = 0.0;
3523 input_doublev("NEGF.tran.interpolate.coes",2, r_vec1, r_vec);
3524 interpolate_c1 = r_vec1[0];
3525 interpolate_c2 = r_vec1[1];
3526
3527 if (1.0e-9<fabs(r_vec1[0]+r_vec1[1]-1.0)){
3528 printf("The sum of coefficients for the interpolation should be unity.\n");
3529 MPI_Finalize();
3530 exit(0);
3531 }
3532
3533 input_close();
3534 }
3535
3536
MTRAN_Set_MP(int job,int anum,int * WhatSpecies,int * Spe_Total_CNO,int * NUM,int * MP)3537 void MTRAN_Set_MP(
3538 int job,
3539 int anum, int *WhatSpecies, int *Spe_Total_CNO,
3540 int *NUM, /* output */
3541 int *MP /* output */
3542 )
3543 {
3544 int Anum, i, wanA, tnoA;
3545
3546 /* setup MP */
3547 Anum = 1;
3548 for (i=1; i<=anum; i++){
3549 if (job) MP[i]=Anum;
3550 wanA = WhatSpecies[i];
3551 tnoA = Spe_Total_CNO[wanA];
3552 Anum += tnoA;
3553 }
3554
3555 *NUM=Anum-1;
3556 }
3557
3558
3559
3560
3561
3562
MTRAN_Set_SurfOverlap(char * position,double k2,double k3,int SpinP_switch,int atomnum_e[2],double ***** OLP_e[2],double ***** H_e[2],double ***** iHNL_e[2],int SpeciesNum_e[2],int * WhatSpecies_e[2],int * Spe_Total_CNO_e[2],int * FNAN_e[2],int ** natn_e[2],int ** ncn_e[2],int ** atv_ijk_e[2],dcomplex * S00_e[2],dcomplex * S01_e[2],dcomplex * H00_e[2],dcomplex * H01_e[2])3563 void MTRAN_Set_SurfOverlap(char *position,
3564 double k2,
3565 double k3,
3566 int SpinP_switch,
3567 int atomnum_e[2],
3568 double *****OLP_e[2],
3569 double *****H_e[2],
3570 double *****iHNL_e[2],
3571 int SpeciesNum_e[2],
3572 int *WhatSpecies_e[2],
3573 int *Spe_Total_CNO_e[2],
3574 int *FNAN_e[2],
3575 int **natn_e[2],
3576 int **ncn_e[2],
3577 int **atv_ijk_e[2],
3578 dcomplex *S00_e[2],
3579 dcomplex *S01_e[2],
3580 dcomplex *H00_e[2],
3581 dcomplex *H01_e[2]
3582 )
3583 /* #define S00_ref(i,j) ( ((j)-1)*NUM+(i)-1 ) */
3584 #define S00_ref(i,j) ( ((j)-1)*NUM+(i)-1 )
3585 #define S01_ref(i,j) ( ((j+NUM0)-1)*NUM+(i)-1 )
3586 #define S10_ref(i,j) ( ((j)-1)*NUM+(i+NUM0)-1 )
3587 #define S11_ref(i,j) ( ((j+NUM0)-1)*NUM+(i+NUM0)-1 )
3588 {
3589 int NUM,n2;
3590 int Anum, Bnum, wanA, tnoA;
3591 int i,j,k;
3592 int GA_AN;
3593 int GB_AN, LB_AN,wanB, tnoB,Rn;
3594 int l1,l2,l3;
3595 int spin,MN;
3596 int SpinP_switch_e[2];
3597 int direction,iside;
3598 double si,co,kRn;
3599 double s,h[10];
3600 static double epscutoff=1.0e-12;
3601 double epscutoff2;
3602 int *MP;
3603 double a00,a11,a01r,a01i,b00,b11,b01;
3604 int NUM0;
3605 /*debug */
3606 char msg[100];
3607 /*end debug */
3608
3609 SpinP_switch_e[0] = SpinP_switch;
3610 SpinP_switch_e[1] = SpinP_switch;
3611
3612 /* position -> direction */
3613
3614 if ( strcasecmp(position,"left")==0) {
3615 direction=-1;
3616 iside=0;
3617 }
3618 else if ( strcasecmp(position,"right")==0) {
3619 direction= 1;
3620 iside=1;
3621 }
3622
3623 /* set MP */
3624 /* MTRAN_Set_MP(0, atomnum_e[iside], WhatSpecies_e[iside], Spe_Total_CNO_e[iside], &NUM, &i);*/
3625 MP = (int*)malloc(sizeof(int)*(atomnum_e[iside]+1));
3626 MTRAN_Set_MP(1, atomnum_e[iside], WhatSpecies_e[iside], Spe_Total_CNO_e[iside], &NUM, MP);
3627 NUM0=NUM;
3628 NUM=2*NUM;
3629
3630 n2 = NUM ;
3631
3632 for (i=0; i<n2*n2; i++){
3633 S00_e[iside][i].r = 0.0;
3634 S00_e[iside][i].i = 0.0;
3635 S01_e[iside][i].r = 0.0;
3636 S01_e[iside][i].i = 0.0;
3637 }
3638
3639 for (i=0; i<n2*n2; i++){
3640 H00_e[iside][i].r = 0.0;
3641 H00_e[iside][i].i = 0.0;
3642 H01_e[iside][i].r = 0.0;
3643 H01_e[iside][i].i = 0.0;
3644 }
3645
3646 for (GA_AN=1; GA_AN<=atomnum_e[iside]; GA_AN++){
3647
3648 wanA = WhatSpecies_e[iside][GA_AN];
3649 tnoA = Spe_Total_CNO_e[iside][wanA];
3650 Anum = MP[GA_AN];
3651
3652 for (LB_AN=0; LB_AN<=FNAN_e[iside][GA_AN]; LB_AN++){
3653
3654 GB_AN = natn_e[iside][GA_AN][LB_AN];
3655 Rn = ncn_e[iside][GA_AN][LB_AN];
3656 wanB = WhatSpecies_e[iside][GB_AN];
3657 tnoB = Spe_Total_CNO_e[iside][wanB];
3658 Bnum = MP[GB_AN];
3659
3660 l1 = atv_ijk_e[iside][Rn][1];
3661 l2 = atv_ijk_e[iside][Rn][2];
3662 l3 = atv_ijk_e[iside][Rn][3];
3663
3664 kRn = k2*(double)l2 + k3*(double)l3;
3665 si = sin(2.0*PI*kRn);
3666 co = cos(2.0*PI*kRn);
3667
3668 if (l1==0) {
3669
3670 for (i=0; i<tnoA; i++){
3671 for (j=0; j<tnoB; j++){
3672
3673 /* S_alpha-alpha */
3674 S00_e[iside][S00_ref(Anum+i,Bnum+j)].r += co*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3675 S00_e[iside][S00_ref(Anum+i,Bnum+j)].i += si*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3676 /* S_beta-beta */
3677 S00_e[iside][S11_ref(Anum+i,Bnum+j)].r += co*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3678 S00_e[iside][S11_ref(Anum+i,Bnum+j)].i += si*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3679
3680 a00=H_e[iside][0][GA_AN][LB_AN][i][j];
3681 a11=H_e[iside][1][GA_AN][LB_AN][i][j];
3682 a01r=H_e[iside][2][GA_AN][LB_AN][i][j];
3683 a01i=H_e[iside][3][GA_AN][LB_AN][i][j];
3684
3685 b00=iHNL_e[iside][0][GA_AN][LB_AN][i][j];
3686 b11=iHNL_e[iside][1][GA_AN][LB_AN][i][j];
3687 b01=iHNL_e[iside][2][GA_AN][LB_AN][i][j];
3688 /* H_alpha-alpha */
3689 H00_e[iside][S00_ref(Anum+i,Bnum+j)].r += co*a00 - si*b00;
3690 H00_e[iside][S00_ref(Anum+i,Bnum+j)].i += co*b00 + si*a00;
3691 /* H_alpha-beta */
3692 H00_e[iside][S01_ref(Anum+i,Bnum+j)].r += co*a01r - si*(a01i+b01);
3693 H00_e[iside][S01_ref(Anum+i,Bnum+j)].i += si*a01r + co*(a01i+b01);
3694 /* H_beta-beta */
3695 H00_e[iside][S11_ref(Anum+i,Bnum+j)].r += co*a11 - si*b11;
3696 H00_e[iside][S11_ref(Anum+i,Bnum+j)].i += co*b11 + si*a11;
3697
3698 }
3699 }
3700 }
3701
3702 if (l1==direction) {
3703
3704 for (i=0; i<tnoA; i++){
3705 for (j=0; j<tnoB; j++){
3706 /* S_alpha-alpha */
3707 S01_e[iside][S00_ref(Anum+i,Bnum+j)].r += co*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3708 S01_e[iside][S00_ref(Anum+i,Bnum+j)].i += si*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3709 /* S_beta-beta */
3710 S01_e[iside][S11_ref(Anum+i,Bnum+j)].r += co*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3711 S01_e[iside][S11_ref(Anum+i,Bnum+j)].i += si*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3712
3713 a00=H_e[iside][0][GA_AN][LB_AN][i][j];
3714 a11=H_e[iside][1][GA_AN][LB_AN][i][j];
3715 a01r=H_e[iside][2][GA_AN][LB_AN][i][j];
3716 a01i=H_e[iside][3][GA_AN][LB_AN][i][j];
3717
3718 b00=iHNL_e[iside][0][GA_AN][LB_AN][i][j];
3719 b11=iHNL_e[iside][1][GA_AN][LB_AN][i][j];
3720 b01=iHNL_e[iside][2][GA_AN][LB_AN][i][j];
3721 /* H_alpha-alpha */
3722 H01_e[iside][S00_ref(Anum+i,Bnum+j)].r += co*a00 - si*b00;
3723 H01_e[iside][S00_ref(Anum+i,Bnum+j)].i += co*b00 + si*a00;
3724 /* H_alpha-beta */
3725 H01_e[iside][S01_ref(Anum+i,Bnum+j)].r += co*a01r - si*(a01i+b01);
3726 H01_e[iside][S01_ref(Anum+i,Bnum+j)].i += si*a01r + co*(a01i+b01);
3727 /* H_beta-beta */
3728 H01_e[iside][S11_ref(Anum+i,Bnum+j)].r += co*a11 - si*b11;
3729 H01_e[iside][S11_ref(Anum+i,Bnum+j)].i += co*b11 + si*a11;
3730 /*
3731 S01_e[iside][S00_ref(Anum+i,Bnum+j)].r += co*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3732 S01_e[iside][S00_ref(Anum+i,Bnum+j)].i += si*OLP_e[iside][0][GA_AN][LB_AN][i][j];
3733
3734 for (k=0; k<=SpinP_switch_e[iside]; k++ ){
3735 H01_e[iside][k][S00_ref(Anum+i,Bnum+j)].r += co*H_e[iside][k][GA_AN][LB_AN][i][j];
3736 H01_e[iside][k][S00_ref(Anum+i,Bnum+j)].i += si*H_e[iside][k][GA_AN][LB_AN][i][j];
3737 }
3738 */
3739 }
3740 }
3741 }
3742 }
3743 } /* GA_AN */
3744
3745 for (i=1; i<=NUM0; i++) {
3746 for (j=1; j<=NUM0; j++) {
3747
3748 H00_e[iside][(j-1)*NUM+i+NUM0-1].r = H00_e[iside][(i+NUM0-1)*NUM+j-1].r;
3749 H00_e[iside][(j-1)*NUM+i+NUM0-1].i = -H00_e[iside][(i+NUM0-1)*NUM+j-1].i;
3750
3751 H01_e[iside][(j-1)*NUM+i+NUM0-1].r = H01_e[iside][(i+NUM0-1)*NUM+j-1].r;
3752 H01_e[iside][(j-1)*NUM+i+NUM0-1].i = -H01_e[iside][(i+NUM0-1)*NUM+j-1].i;
3753
3754 }
3755 }
3756
3757
3758 /* to set some very small value to zero to avoid numerical instability*/
3759 /*
3760 for (i=0; i<n2*n2; i++){
3761
3762 if ( (fabs(S00_e[iside][i].r)+fabs(S00_e[iside][i].r)) < epscutoff ) {
3763 S00_e[iside][i].r = 0.0;
3764 S00_e[iside][i].i = 0.0;
3765 }
3766
3767 if ( (fabs(S01_e[iside][i].r)+fabs(S01_e[iside][i].r)) < epscutoff ) {
3768 S01_e[iside][i].r = 0.0;
3769 S01_e[iside][i].i = 0.0;
3770 }
3771
3772 if ( (fabs(H00_e[iside][i].r)+fabs(H00_e[iside][i].r)) < epscutoff ) {
3773 H00_e[iside][i].r = 0.0;
3774 H00_e[iside][i].i = 0.0;
3775 }
3776
3777 if ( (fabs(H01_e[iside][i].r)+fabs(H01_e[iside][i].r)) < epscutoff ) {
3778 H01_e[iside][i].r = 0.0;
3779 H01_e[iside][i].i = 0.0;
3780 }
3781
3782 }
3783 */
3784
3785 /* free arrays */
3786
3787 free(MP);
3788 }
3789
3790
3791
3792
MTRAN_Set_CentOverlap(int job,int SpinP_switch,double k2,double k3,int NUM_c,int NUM_e[2],double ***** H,double ***** iHNL,double ***** OLP,int atomnum,int atomnum_e[2],int * WhatSpecies,int * WhatSpecies_e[2],int * Spe_Total_CNO,int * Spe_Total_CNO_e[2],int * FNAN,int ** natn,int ** ncn,int ** atv_ijk,int * TRAN_region,int * TRAN_Original_Id)3793 void MTRAN_Set_CentOverlap(
3794 int job,
3795 int SpinP_switch,
3796 double k2,
3797 double k3,
3798 int NUM_c,
3799 int NUM_e[2],
3800 double *****H,
3801 double *****iHNL,
3802 double *****OLP,
3803 int atomnum,
3804 int atomnum_e[2],
3805 int *WhatSpecies,
3806 int *WhatSpecies_e[2],
3807 int *Spe_Total_CNO,
3808 int *Spe_Total_CNO_e[2],
3809 int *FNAN,
3810 int **natn,
3811 int **ncn,
3812 int **atv_ijk,
3813 int *TRAN_region,
3814 int *TRAN_Original_Id
3815 )
3816
3817 #define SCC00_ref(i,j) ( ((j)-1)*nc+(i)-1 )
3818 #define SCC01_ref(i,j) ( ((j+nc0)-1)*nc+(i)-1 )
3819 #define SCC10_ref(i,j) ( ((j)-1)*nc+(i+nc0)-1 )
3820 #define SCC11_ref(i,j) ( ((j+nc0)-1)*nc+(i+nc0)-1 )
3821
3822 #define SCL00_ref(i,j) ( ((j)-1)*nc+(i)-1 )
3823 #define SCL01_ref(i,j) ( ((j+na0)-1)*nc+(i)-1 )
3824 #define SCL10_ref(i,j) ( ((j)-1)*nc+(i+nc0)-1 )
3825 #define SCL11_ref(i,j) ( ((j+na0)-1)*nc+(i+nc0)-1 )
3826
3827 #define SCR00_ref(i,j) ( ((j)-1)*nc+(i)-1 )
3828 #define SCR01_ref(i,j) ( ((j+nb0)-1)*nc+(i)-1 )
3829 #define SCR10_ref(i,j) ( ((j)-1)*nc+(i+nc0)-1 )
3830 #define SCR11_ref(i,j) ( ((j+nb0)-1)*nc+(i+nc0)-1 )
3831
3832 #define S00l00_ref(i,j) ( ((j)-1)*na+(i)-1 )
3833 #define S00l01_ref(i,j) ( ((j+na0)-1)*na+(i)-1 )
3834 #define S00l10_ref(i,j) ( ((j)-1)*na+(i+na0)-1 )
3835 #define S00l11_ref(i,j) ( ((j+na0)-1)*na+(i+na0)-1 )
3836
3837 #define S00r00_ref(i,j) ( ((j)-1)*nb+(i)-1 )
3838 #define S00r01_ref(i,j) ( ((j+nb0)-1)*nb+(i)-1 )
3839 #define S00r10_ref(i,j) ( ((j)-1)*nb+(i+nb0)-1 )
3840 #define S00r11_ref(i,j) ( ((j+nb0)-1)*nb+(i+nb0)-1 )
3841
3842 {
3843 int *MP, *MP_e[2];
3844 int i;
3845 int nc0,nc,na0,na,nb0,nb;
3846
3847 /* setup MP */
3848 MP = (int*)malloc(sizeof(int)*(atomnum+1));
3849 MTRAN_Set_MP( 1, atomnum, WhatSpecies, Spe_Total_CNO, &i, MP);
3850 nc0=i;
3851 nc=2*i;
3852
3853 MP_e[0] = (int*)malloc(sizeof(int)*(atomnum_e[0]+1));
3854 MTRAN_Set_MP( 1, atomnum_e[0], WhatSpecies_e[0], Spe_Total_CNO_e[0], &i, MP_e[0]);
3855 na0=i;
3856 na=2*i;
3857
3858 MP_e[1] = (int*)malloc(sizeof(int)*(atomnum_e[1]+1));
3859 MTRAN_Set_MP( 1, atomnum_e[1], WhatSpecies_e[1], Spe_Total_CNO_e[1], &i, MP_e[1]);
3860 nb0=i;
3861 nb=2*i;
3862
3863
3864 for (i=0; i<NUM_c*NUM_c; i++) {
3865 SCC[i].r = 0.0;
3866 SCC[i].i = 0.0;
3867 }
3868 for (i=0; i<NUM_c*NUM_e[0]; i++) {
3869 SCL[i].r = 0.0;
3870 SCL[i].i = 0.0;
3871 }
3872 for (i=0; i<NUM_c*NUM_e[1]; i++) {
3873 SCR[i].r = 0.0;
3874 SCR[i].i = 0.0;
3875 }
3876
3877 for (i=0; i<NUM_c*NUM_c; i++) {
3878 HCC[i].r = 0.0;
3879 HCC[i].i = 0.0;
3880 /* S MitsuakiKAWAMURA */
3881 VCC[i].r = 0.0;
3882 VCC[i].i = 0.0;
3883 /* E MitsuakiKAWAMURA */
3884 }
3885
3886 for (i=0; i<NUM_c*NUM_e[0]; i++) {
3887 HCL[i].r = 0.0;
3888 HCL[i].i = 0.0;
3889 }
3890
3891 for (i=0; i<NUM_c*NUM_e[1]; i++) {
3892 HCR[i].r = 0.0;
3893 HCR[i].i = 0.0;
3894 }
3895
3896
3897 if ((job&1)==1) {
3898
3899 int MA_AN,GA_AN, wanA, tnoA, Anum;
3900 int LB_AN, GB_AN, wanB, tnoB, l1,l2,l3, Bnum;
3901 int i,j,k,q,AN;
3902 int Rn;
3903 double kRn,si,co;
3904 double tmp,s0,h0;
3905 double a00,a11,a01r,a01i,b00,b11,b01;
3906
3907 /* make Overlap , HCC, SCC */
3908 /* parallel global GA_AN 1:atomnum */
3909
3910 for (GA_AN=1; GA_AN<=atomnum; GA_AN++){
3911
3912 wanA = WhatSpecies[GA_AN];
3913 tnoA = Spe_Total_CNO[wanA];
3914 Anum = MP[GA_AN];
3915
3916 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
3917
3918 GB_AN = natn[GA_AN][LB_AN];
3919 Rn = ncn[GA_AN][LB_AN];
3920 wanB = WhatSpecies[GB_AN];
3921 tnoB = Spe_Total_CNO[wanB];
3922 Bnum = MP[GB_AN];
3923
3924 l1 = atv_ijk[Rn][1];
3925 l2 = atv_ijk[Rn][2];
3926 l3 = atv_ijk[Rn][3];
3927
3928 kRn = k2*(double)l2 + k3*(double)l3;
3929 si = sin(2.0*PI*kRn);
3930 co = cos(2.0*PI*kRn);
3931
3932 for (i=0; i<tnoA; i++){
3933 for (j=0; j<tnoB; j++){
3934
3935 /* l1 is the direction to the electrode */
3936
3937 if (l1==0){
3938
3939 /* S_alpha-alpha */
3940 SCC[SCC00_ref(Anum+i,Bnum+j)].r += co*OLP[0][GA_AN][LB_AN][i][j];
3941 SCC[SCC00_ref(Anum+i,Bnum+j)].i += si*OLP[0][GA_AN][LB_AN][i][j];
3942 /* S_beta-beta */
3943 SCC[SCC11_ref(Anum+i,Bnum+j)].r += co*OLP[0][GA_AN][LB_AN][i][j];
3944 SCC[SCC11_ref(Anum+i,Bnum+j)].i += si*OLP[0][GA_AN][LB_AN][i][j];
3945
3946 a00=H[0][GA_AN][LB_AN][i][j];
3947 a11=H[1][GA_AN][LB_AN][i][j];
3948 a01r=H[2][GA_AN][LB_AN][i][j];
3949 a01i=H[3][GA_AN][LB_AN][i][j];
3950
3951 b00=iHNL[0][GA_AN][LB_AN][i][j];
3952 b11=iHNL[1][GA_AN][LB_AN][i][j];
3953 b01=iHNL[2][GA_AN][LB_AN][i][j];
3954
3955 /* H_alpha-alpha */
3956 HCC[SCC00_ref(Anum+i,Bnum+j)].r += co*a00 - si*b00;
3957 HCC[SCC00_ref(Anum+i,Bnum+j)].i += co*b00 + si*a00;
3958 /* H_alpha-beta */
3959 HCC[SCC01_ref(Anum+i,Bnum+j)].r += co*a01r - si*(a01i+b01);
3960 HCC[SCC01_ref(Anum+i,Bnum+j)].i += si*a01r + co*(a01i+b01);
3961 /* H_beta-beta */
3962 HCC[SCC11_ref(Anum+i,Bnum+j)].r += co*a11 - si*b11;
3963 HCC[SCC11_ref(Anum+i,Bnum+j)].i += co*b11 + si*a11;
3964
3965 a00 = (H[0][GA_AN][LB_AN][i][j] - H0[0][GA_AN][LB_AN][i][j]);
3966 a11 = (H[1][GA_AN][LB_AN][i][j] - H0[0][GA_AN][LB_AN][i][j]);
3967
3968 /* H_alpha-alpha */
3969 VCC[SCC00_ref(Anum + i, Bnum + j)].r += co*a00 - si*b00;
3970 VCC[SCC00_ref(Anum + i, Bnum + j)].i += co*b00 + si*a00;
3971 /* H_alpha-beta */
3972 VCC[SCC01_ref(Anum + i, Bnum + j)].r += co*a01r - si*(a01i + b01);
3973 VCC[SCC01_ref(Anum + i, Bnum + j)].i += si*a01r + co*(a01i + b01);
3974 /* H_beta-beta */
3975 VCC[SCC11_ref(Anum + i, Bnum + j)].r += co*a11 - si*b11;
3976 VCC[SCC11_ref(Anum + i, Bnum + j)].i += co*b11 + si*a11;
3977
3978 }
3979
3980 }
3981 }
3982
3983 } /* LB_AN */
3984 } /* MA_AN */
3985
3986 for (i=1; i<=nc0; i++) {
3987 for (j=1; j<=nc0; j++) {
3988
3989 HCC[(j - 1)*nc + i + nc0 - 1].r = HCC[(i + nc0 - 1)*nc + j - 1].r;
3990 HCC[(j - 1)*nc + i + nc0 - 1].i = -HCC[(i + nc0 - 1)*nc + j - 1].i;
3991
3992 VCC[(j - 1)*nc + i + nc0 - 1].r = VCC[(i + nc0 - 1)*nc + j - 1].r;
3993 VCC[(j - 1)*nc + i + nc0 - 1].i = -VCC[(i + nc0 - 1)*nc + j - 1].i;
3994
3995 }
3996 }
3997
3998
3999 } /* job&1 */
4000
4001
4002 if ( (job&2) == 2 ) {
4003
4004 {
4005 int MA_AN, GA_AN, wanA, tnoA, Anum;
4006 int GA_AN_e, Anum_e;
4007 int GB_AN, wanB, tnoB, Bnum;
4008 int GB_AN_e, Bnum_e;
4009 int i,j,k;
4010 int iside;
4011
4012 /* overwrite CL1 region */
4013
4014 iside = 0;
4015
4016 /* parallel global GA_AN 1:atomnum */
4017
4018 for (GA_AN=1; GA_AN<=atomnum; GA_AN++){
4019
4020 wanA = WhatSpecies[GA_AN];
4021 tnoA = Spe_Total_CNO[wanA];
4022 Anum = MP[GA_AN];
4023
4024 for (GB_AN=1; GB_AN<=atomnum; GB_AN++){
4025
4026 if ( TRAN_region[GA_AN]==12 && TRAN_region[GB_AN]==12 ) {
4027
4028 GA_AN_e = TRAN_Original_Id[GA_AN];
4029 Anum_e = MP_e[iside][GA_AN_e];
4030
4031 wanB = WhatSpecies[GB_AN];
4032 tnoB = Spe_Total_CNO[wanB];
4033 Bnum = MP[GB_AN];
4034
4035 GB_AN_e = TRAN_Original_Id[GB_AN];
4036 Bnum_e = MP_e[iside][GB_AN_e];
4037
4038 for (i=0; i<tnoA; i++){
4039 for (j=0; j<tnoB; j++){
4040
4041 /* S_alpha-alpha */
4042 SCC[SCC00_ref(Anum + i, Bnum + j)].r = S00_e[iside][S00l00_ref(Anum_e + i, Bnum_e + j)].r;
4043 /* S_beta-beta */
4044 SCC[SCC11_ref(Anum + i, Bnum + j)].r = S00_e[iside][S00l11_ref(Anum_e + i, Bnum_e + j)].r;
4045
4046 /* V_alpha-alpha */
4047 VCC[SCC00_ref(Anum + i, Bnum + j)].r = VCC[SCC00_ref(Anum + i, Bnum + j)].r
4048 + H00_e[iside][S00l00_ref(Anum_e + i, Bnum_e + j)].r - HCC[SCC00_ref(Anum + i, Bnum + j)].r;
4049 /* V_alpha-beta */
4050 VCC[SCC01_ref(Anum + i, Bnum + j)].r = VCC[SCC01_ref(Anum + i, Bnum + j)].r
4051 + H00_e[iside][S00l01_ref(Anum_e + i, Bnum_e + j)].r - HCC[SCC01_ref(Anum + i, Bnum + j)].r;
4052 /* V_beta-alpha */
4053 VCC[SCC10_ref(Anum + i, Bnum + j)].r = VCC[SCC10_ref(Anum + i, Bnum + j)].r
4054 + H00_e[iside][S00l10_ref(Anum_e + i, Bnum_e + j)].r - HCC[SCC10_ref(Anum + i, Bnum + j)].r;
4055 /* V_beta-beta */
4056 VCC[SCC11_ref(Anum + i, Bnum + j)].r = VCC[SCC11_ref(Anum + i, Bnum + j)].r
4057 + H00_e[iside][S00l11_ref(Anum_e + i, Bnum_e + j)].r - HCC[SCC11_ref(Anum + i, Bnum + j)].r;
4058
4059 /* H_alpha-alpha */
4060 HCC[SCC00_ref(Anum + i, Bnum + j)].r = H00_e[iside][S00l00_ref(Anum_e + i, Bnum_e + j)].r;
4061 /* H_alpha-beta */
4062 HCC[SCC01_ref(Anum + i, Bnum + j)].r = H00_e[iside][S00l01_ref(Anum_e + i, Bnum_e + j)].r;
4063 /* H_beta-alpha */
4064 HCC[SCC10_ref(Anum + i, Bnum + j)].r = H00_e[iside][S00l10_ref(Anum_e + i, Bnum_e + j)].r;
4065 /* H_beta-beta */
4066 HCC[SCC11_ref(Anum + i, Bnum + j)].r = H00_e[iside][S00l11_ref(Anum_e + i, Bnum_e + j)].r;
4067
4068 /* S_alpha-alpha */
4069 SCC[SCC00_ref(Anum + i, Bnum + j)].i = S00_e[iside][S00l00_ref(Anum_e + i, Bnum_e + j)].i;
4070 /* S_beta-beta */
4071 SCC[SCC11_ref(Anum + i, Bnum + j)].i = S00_e[iside][S00l11_ref(Anum_e + i, Bnum_e + j)].i;
4072
4073 /* V_alpha-alpha */
4074 VCC[SCC00_ref(Anum + i, Bnum + j)].i = VCC[SCC00_ref(Anum + i, Bnum + j)].i
4075 + H00_e[iside][S00l00_ref(Anum_e + i, Bnum_e + j)].i - HCC[SCC00_ref(Anum + i, Bnum + j)].i;
4076 /* V_alpha-beta */
4077 VCC[SCC01_ref(Anum + i, Bnum + j)].i = VCC[SCC01_ref(Anum + i, Bnum + j)].i
4078 + H00_e[iside][S00l01_ref(Anum_e + i, Bnum_e + j)].i - HCC[SCC01_ref(Anum + i, Bnum + j)].i;
4079 /* V_beta-alpha */
4080 VCC[SCC10_ref(Anum + i, Bnum + j)].i = VCC[SCC10_ref(Anum + i, Bnum + j)].i
4081 + H00_e[iside][S00l10_ref(Anum_e + i, Bnum_e + j)].i - HCC[SCC10_ref(Anum + i, Bnum + j)].i;
4082 /* V_beta-beta */
4083 VCC[SCC11_ref(Anum + i, Bnum + j)].i = VCC[SCC11_ref(Anum + i, Bnum + j)].i
4084 + H00_e[iside][S00l11_ref(Anum_e + i, Bnum_e + j)].i - HCC[SCC11_ref(Anum + i, Bnum + j)].i;
4085
4086 /* H_alpha-alpha */
4087 HCC[SCC00_ref(Anum + i, Bnum + j)].i = H00_e[iside][S00l00_ref(Anum_e + i, Bnum_e + j)].i;
4088 /* H_alpha-beta */
4089 HCC[SCC01_ref(Anum + i, Bnum + j)].i = H00_e[iside][S00l01_ref(Anum_e + i, Bnum_e + j)].i;
4090 /* H_beta-alpha */
4091 HCC[SCC10_ref(Anum + i, Bnum + j)].i = H00_e[iside][S00l10_ref(Anum_e + i, Bnum_e + j)].i;
4092 /* H_beta-beta */
4093 HCC[SCC11_ref(Anum + i, Bnum + j)].i = H00_e[iside][S00l11_ref(Anum_e + i, Bnum_e + j)].i;
4094 /*
4095 SCC[SCC_ref(Anum+i,Bnum+j)] = S00_e[iside][S00l_ref(Anum_e+i,Bnum_e+j)];
4096
4097 for (k=0; k<=SpinP_switch; k++) {
4098 HCC[k][SCC_ref(Anum+i,Bnum+j)] = H00_e[iside][k][S00l_ref(Anum_e+i,Bnum_e+j)];
4099 } */
4100 }
4101 }
4102
4103 }
4104 }
4105 }
4106 }
4107
4108 {
4109 int MA_AN, GA_AN, wanA, tnoA, Anum;
4110 int GA_AN_e, Anum_e;
4111 int GB_AN, wanB, tnoB, Bnum;
4112 int GB_AN_e, Bnum_e;
4113 int i,j,k;
4114 int iside;
4115
4116 /* overwrite CR1 region */
4117
4118 iside = 1;
4119
4120 /*parallel global GA_AN 1:atomnum */
4121 /*parallel local MA_AN 1:Matomnum */
4122 /*parallel variable GA_AN = M2G[MA_AN] */
4123
4124 for (GA_AN=1; GA_AN<=atomnum; GA_AN++){
4125
4126 wanA = WhatSpecies[GA_AN];
4127 tnoA = Spe_Total_CNO[wanA];
4128 Anum = MP[GA_AN];
4129
4130 for (GB_AN=1; GB_AN<=atomnum; GB_AN++){
4131
4132 if ( TRAN_region[GA_AN]==13 && TRAN_region[GB_AN]==13 ) {
4133
4134 GA_AN_e = TRAN_Original_Id[GA_AN];
4135 Anum_e = MP_e[iside][GA_AN_e]; /* = Anum */
4136 wanB = WhatSpecies[GB_AN];
4137 tnoB = Spe_Total_CNO[wanB];
4138 Bnum = MP[GB_AN];
4139
4140 GB_AN_e = TRAN_Original_Id[GB_AN];
4141 Bnum_e = MP_e[iside][GB_AN_e];
4142
4143 for (i=0; i<tnoA; i++){
4144 for (j=0; j<tnoB; j++){
4145
4146 /* S_alpha-alpha */
4147 SCC[SCC00_ref(Anum + i, Bnum + j)].r = S00_e[iside][S00r00_ref(Anum_e + i, Bnum_e + j)].r;
4148 /* S_beta-beta */
4149 SCC[SCC11_ref(Anum + i, Bnum + j)].r = S00_e[iside][S00r11_ref(Anum_e + i, Bnum_e + j)].r;
4150
4151 /* H_alpha-alpha */
4152 VCC[SCC00_ref(Anum + i, Bnum + j)].r = VCC[SCC00_ref(Anum + i, Bnum + j)].r
4153 +H00_e[iside][S00r00_ref(Anum_e + i, Bnum_e + j)].r - HCC[SCC00_ref(Anum + i, Bnum + j)].r;
4154 /* H_alpha-beta */
4155 VCC[SCC01_ref(Anum + i, Bnum + j)].r = VCC[SCC01_ref(Anum + i, Bnum + j)].r
4156 +H00_e[iside][S00r01_ref(Anum_e + i, Bnum_e + j)].r - HCC[SCC01_ref(Anum + i, Bnum + j)].r;
4157 /* H_beta-alpha */
4158 VCC[SCC10_ref(Anum + i, Bnum + j)].r = VCC[SCC10_ref(Anum + i, Bnum + j)].r
4159 +H00_e[iside][S00r10_ref(Anum_e + i, Bnum_e + j)].r - HCC[SCC10_ref(Anum + i, Bnum + j)].r;
4160 /* H_beta-beta */
4161 VCC[SCC11_ref(Anum + i, Bnum + j)].r = VCC[SCC11_ref(Anum + i, Bnum + j)].r
4162 +H00_e[iside][S00r11_ref(Anum_e + i, Bnum_e + j)].r - HCC[SCC11_ref(Anum + i, Bnum + j)].r;
4163
4164 /* H_alpha-alpha */
4165 HCC[SCC00_ref(Anum + i, Bnum + j)].r = H00_e[iside][S00r00_ref(Anum_e + i, Bnum_e + j)].r;
4166 /* H_alpha-beta */
4167 HCC[SCC01_ref(Anum + i, Bnum + j)].r = H00_e[iside][S00r01_ref(Anum_e + i, Bnum_e + j)].r;
4168 /* H_beta-alpha */
4169 HCC[SCC10_ref(Anum + i, Bnum + j)].r = H00_e[iside][S00r10_ref(Anum_e + i, Bnum_e + j)].r;
4170 /* H_beta-beta */
4171 HCC[SCC11_ref(Anum + i, Bnum + j)].r = H00_e[iside][S00r11_ref(Anum_e + i, Bnum_e + j)].r;
4172
4173 /* S_alpha-alpha */
4174 SCC[SCC00_ref(Anum + i, Bnum + j)].i = S00_e[iside][S00r00_ref(Anum_e + i, Bnum_e + j)].i;
4175 /* S_beta-beta */
4176 SCC[SCC11_ref(Anum + i, Bnum + j)].i = S00_e[iside][S00r11_ref(Anum_e + i, Bnum_e + j)].i;
4177
4178 /* H_alpha-alpha */
4179 VCC[SCC00_ref(Anum + i, Bnum + j)].i = VCC[SCC00_ref(Anum + i, Bnum + j)].i
4180 + H00_e[iside][S00r00_ref(Anum_e + i, Bnum_e + j)].i - HCC[SCC00_ref(Anum + i, Bnum + j)].i;
4181 /* H_alpha-beta */
4182 VCC[SCC01_ref(Anum + i, Bnum + j)].i = VCC[SCC01_ref(Anum + i, Bnum + j)].i
4183 + H00_e[iside][S00r01_ref(Anum_e + i, Bnum_e + j)].i - HCC[SCC01_ref(Anum + i, Bnum + j)].i;
4184 /* H_beta-alpha */
4185 VCC[SCC10_ref(Anum + i, Bnum + j)].i = VCC[SCC10_ref(Anum + i, Bnum + j)].i
4186 + H00_e[iside][S00r10_ref(Anum_e + i, Bnum_e + j)].i - HCC[SCC10_ref(Anum + i, Bnum + j)].i;
4187 /* H_beta-beta */
4188 VCC[SCC11_ref(Anum + i, Bnum + j)].i = VCC[SCC11_ref(Anum + i, Bnum + j)].i
4189 + H00_e[iside][S00r11_ref(Anum_e + i, Bnum_e + j)].i - HCC[SCC11_ref(Anum + i, Bnum + j)].i;
4190
4191 /* H_alpha-alpha */
4192 HCC[SCC00_ref(Anum + i, Bnum + j)].i = H00_e[iside][S00r00_ref(Anum_e + i, Bnum_e + j)].i;
4193 /* H_alpha-beta */
4194 HCC[SCC01_ref(Anum + i, Bnum + j)].i = H00_e[iside][S00r01_ref(Anum_e + i, Bnum_e + j)].i;
4195 /* H_beta-alpha */
4196 HCC[SCC10_ref(Anum + i, Bnum + j)].i = H00_e[iside][S00r10_ref(Anum_e + i, Bnum_e + j)].i;
4197 /* H_beta-beta */
4198 HCC[SCC11_ref(Anum + i, Bnum + j)].i = H00_e[iside][S00r11_ref(Anum_e + i, Bnum_e + j)].i;
4199 /*
4200 SCC[SCC_ref(Anum+i,Bnum+j)] = S00_e[iside][S00r_ref(Anum_e+i,Bnum_e+j)];
4201
4202 for (k=0; k<=SpinP_switch; k++) {
4203 HCC[k][SCC_ref(Anum+i,Bnum+j)] = H00_e[iside][k][S00r_ref(Anum_e+i,Bnum_e+j)];
4204 } */
4205 }
4206 }
4207
4208 }
4209
4210 }
4211 }
4212 }
4213
4214 {
4215 int iside;
4216 int MA_AN, GA_AN, wanA, tnoA, Anum, GA_AN_e, Anum_e;
4217 int GB_AN_e, wanB_e, tnoB_e, Bnum_e;
4218 int i,j,k;
4219
4220 /* make Overlap, HCL, SCL from OLP_e, and H_e*/
4221
4222 iside = 0;
4223
4224 /*parallel global GA_AN 1:atomnum */
4225 /*parallel local MA_AN 1:Matomnum */
4226 /*parallel variable GA_AN = M2G[MA_AN] */
4227
4228 for (GA_AN=1; GA_AN<=atomnum; GA_AN++){
4229
4230 if (TRAN_region[GA_AN]%10==2 || TRAN_region[GA_AN]==6) {
4231
4232 wanA = WhatSpecies[GA_AN];
4233 tnoA = Spe_Total_CNO[wanA];
4234 Anum = MP[GA_AN]; /* GA_AN is in C */
4235
4236 GA_AN_e = TRAN_Original_Id[GA_AN];
4237 Anum_e = MP_e[iside][GA_AN_e];
4238
4239 for (GB_AN_e=1; GB_AN_e<=atomnum_e[iside]; GB_AN_e++) {
4240
4241 wanB_e = WhatSpecies_e[iside][GB_AN_e];
4242 tnoB_e = Spe_Total_CNO_e[iside][wanB_e];
4243 Bnum_e = MP_e[iside][GB_AN_e];
4244
4245 for (i=0; i<tnoA; i++){
4246 for (j=0; j<tnoB_e; j++){
4247
4248 /* S_alpha-alpha */
4249 SCL[SCL00_ref(Anum+i,Bnum_e+j)].r = S01_e[iside][S00l00_ref(Anum_e+i,Bnum_e+j)].r;
4250 /* S_beta-beta */
4251 SCL[SCL11_ref(Anum+i,Bnum_e+j)].r = S01_e[iside][S00l11_ref(Anum_e+i,Bnum_e+j)].r;
4252
4253 /* H_alpha-alpha */
4254 HCL[SCL00_ref(Anum+i,Bnum_e+j)].r = H01_e[iside][S00l00_ref(Anum_e+i,Bnum_e+j)].r;
4255 /* H_alpha-beta */
4256 HCL[SCL01_ref(Anum+i,Bnum_e+j)].r = H01_e[iside][S00l01_ref(Anum_e+i,Bnum_e+j)].r;
4257 /* H_beta-alpha */
4258 HCL[SCL10_ref(Anum+i,Bnum_e+j)].r = H01_e[iside][S00l10_ref(Anum_e+i,Bnum_e+j)].r;
4259 /* H_beta-beta */
4260 HCL[SCL11_ref(Anum+i,Bnum_e+j)].r = H01_e[iside][S00l11_ref(Anum_e+i,Bnum_e+j)].r;
4261 /* S_alpha-alpha */
4262 SCL[SCL00_ref(Anum+i,Bnum_e+j)].i = S01_e[iside][S00l00_ref(Anum_e+i,Bnum_e+j)].i;
4263 /* S_beta-beta */
4264 SCL[SCL11_ref(Anum+i,Bnum_e+j)].i = S01_e[iside][S00l11_ref(Anum_e+i,Bnum_e+j)].i;
4265
4266 /* H_alpha-alpha */
4267 HCL[SCL00_ref(Anum+i,Bnum_e+j)].i = H01_e[iside][S00l00_ref(Anum_e+i,Bnum_e+j)].i;
4268 /* H_alpha-beta */
4269 HCL[SCL01_ref(Anum+i,Bnum_e+j)].i = H01_e[iside][S00l01_ref(Anum_e+i,Bnum_e+j)].i;
4270 /* H_beta-alpha */
4271 HCL[SCL10_ref(Anum+i,Bnum_e+j)].i = H01_e[iside][S00l10_ref(Anum_e+i,Bnum_e+j)].i;
4272 /* H_beta-beta */
4273 HCL[SCL11_ref(Anum+i,Bnum_e+j)].i = H01_e[iside][S00l11_ref(Anum_e+i,Bnum_e+j)].i;
4274 /*
4275 SCL[SCL_ref(Anum+i,Bnum_e+j)] = S01_e[iside][ S00l_ref(Anum_e+i, Bnum_e+j)];
4276
4277 for (k=0; k<=SpinP_switch; k++) {
4278 HCL[k][SCL_ref(Anum+i,Bnum_e+j)] = H01_e[iside][k][ S00l_ref(Anum_e+i, Bnum_e+j)];
4279 } */
4280 }
4281 }
4282 }
4283 }
4284 }
4285 }
4286
4287 {
4288 int iside;
4289 int MA_AN, GA_AN, wanA, tnoA, Anum, GA_AN_e, Anum_e;
4290 int GB_AN_e, wanB_e, tnoB_e, Bnum_e;
4291 int i,j,k;
4292
4293 /* make Overlap , HCR, SCR from OLP_e, and H_e*/
4294
4295 iside = 1;
4296
4297 /* for (i=0; i<NUM_c*NUM_e[iside]; i++) {
4298 SCR_nc[i].r = 0.0;
4299 SCR_nc[i].i = 0.0;
4300 HCR_nc[i].r = 0.0;
4301 HCR_nc[i].i = 0.0;
4302 } */
4303
4304 for (GA_AN=1; GA_AN<=atomnum; GA_AN++){
4305
4306 if (TRAN_region[GA_AN]%10==3 || TRAN_region[GA_AN]==6 ) {
4307
4308 wanA = WhatSpecies[GA_AN];
4309 tnoA = Spe_Total_CNO[wanA];
4310 Anum = MP[GA_AN]; /* GA_AN is in C */
4311 GA_AN_e = TRAN_Original_Id[GA_AN];
4312 Anum_e = MP_e[iside][GA_AN_e];
4313
4314 for (GB_AN_e=1; GB_AN_e<=atomnum_e[iside];GB_AN_e++) {
4315 wanB_e = WhatSpecies_e[iside][GB_AN_e];
4316 tnoB_e = Spe_Total_CNO_e[iside][wanB_e];
4317 Bnum_e = MP_e[iside][GB_AN_e];
4318 for (i=0; i<tnoA; i++){
4319 for (j=0; j<tnoB_e; j++){
4320
4321 /* S_alpha-alpha */
4322 SCR[SCR00_ref(Anum+i,Bnum_e+j)].r = S01_e[iside][S00r00_ref(Anum_e+i,Bnum_e+j)].r;
4323 /* S_beta-beta */
4324 SCR[SCR11_ref(Anum+i,Bnum_e+j)].r = S01_e[iside][S00r11_ref(Anum_e+i,Bnum_e+j)].r;
4325
4326 /* H_alpha-alpha */
4327 HCR[SCR00_ref(Anum+i,Bnum_e+j)].r = H01_e[iside][S00r00_ref(Anum_e+i,Bnum_e+j)].r;
4328 /* H_alpha-beta */
4329 HCR[SCR01_ref(Anum+i,Bnum_e+j)].r = H01_e[iside][S00r01_ref(Anum_e+i,Bnum_e+j)].r;
4330 /* H_beta-alpha */
4331 HCR[SCR10_ref(Anum+i,Bnum_e+j)].r = H01_e[iside][S00r10_ref(Anum_e+i,Bnum_e+j)].r;
4332 /* H_beta-beta */
4333 HCR[SCR11_ref(Anum+i,Bnum_e+j)].r = H01_e[iside][S00r11_ref(Anum_e+i,Bnum_e+j)].r;
4334
4335 /* S_alpha-alpha */
4336 SCR[SCR00_ref(Anum+i,Bnum_e+j)].i = S01_e[iside][S00r00_ref(Anum_e+i,Bnum_e+j)].i;
4337 /* S_beta-beta */
4338 SCR[SCR11_ref(Anum+i,Bnum_e+j)].i = S01_e[iside][S00r11_ref(Anum_e+i,Bnum_e+j)].i;
4339
4340 /* H_alpha-alpha */
4341 HCR[SCR00_ref(Anum+i,Bnum_e+j)].i = H01_e[iside][S00r00_ref(Anum_e+i,Bnum_e+j)].i;
4342 /* H_alpha-beta */
4343 HCR[SCR01_ref(Anum+i,Bnum_e+j)].i = H01_e[iside][S00r01_ref(Anum_e+i,Bnum_e+j)].i;
4344 /* H_beta-alpha */
4345 HCR[SCR10_ref(Anum+i,Bnum_e+j)].i = H01_e[iside][S00r10_ref(Anum_e+i,Bnum_e+j)].i;
4346 /* H_beta-beta */
4347 HCR[SCR11_ref(Anum+i,Bnum_e+j)].i = H01_e[iside][S00r11_ref(Anum_e+i,Bnum_e+j)].i;
4348 /*
4349 SCR[SCR_ref(Anum+i,Bnum_e+j)] = S01_e[iside][ S00r_ref(Anum_e+i, Bnum_e+j)];
4350
4351 for (k=0; k<=SpinP_switch; k++) {
4352 HCR[k][SCR_ref(Anum+i,Bnum_e+j)] = H01_e[iside][k][ S00r_ref(Anum_e+i, Bnum_e+j)];
4353 } */
4354 }
4355 }
4356 }
4357 }
4358 }
4359 }
4360
4361 } /* job&2 */
4362
4363
4364 free(MP);
4365 free(MP_e[1]);
4366 free(MP_e[0]);
4367 }
4368
4369
4370
MTRAN_Allocate_HS(int NUM_c,int NUM_e[2],int SpinP_switch)4371 void MTRAN_Allocate_HS(int NUM_c,
4372 int NUM_e[2],
4373 int SpinP_switch)
4374 {
4375 int i,side;
4376
4377 for (side=0; side<=1; side++) {
4378
4379 S00_e[side] = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[side]*NUM_e[side] );
4380 S01_e[side] = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[side]*NUM_e[side] );
4381 for (i=0; i<(NUM_e[side]*NUM_e[side]); i++) {
4382 S00_e[side][i].r = 0.0;
4383 S00_e[side][i].i = 0.0;
4384 S01_e[side][i].r = 0.0;
4385 S01_e[side][i].i = 0.0;
4386 }
4387
4388 H00_e[side] = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[side]*NUM_e[side] );
4389 H01_e[side] = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[side]*NUM_e[side] );
4390 for (i=0; i<(NUM_e[side]*NUM_e[side]); i++) {
4391 H00_e[side][i].r = 0.0;
4392 H00_e[side][i].i = 0.0;
4393 H01_e[side][i].r = 0.0;
4394 H01_e[side][i].i = 0.0;
4395 }
4396 }
4397
4398 SCC = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c );
4399 SCL = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_e[0] );
4400 SCR = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_e[1] );
4401
4402 for (i=0; i<NUM_c*NUM_c; i++) {
4403 SCC[i].r = 0.0;
4404 SCC[i].i = 0.0;
4405 }
4406 for (i=0; i<NUM_c*NUM_e[0]; i++) {
4407 SCL[i].r = 0.0;
4408 SCL[i].i = 0.0;
4409 }
4410 for (i=0; i<NUM_c*NUM_e[1]; i++) {
4411 SCR[i].r = 0.0;
4412 SCR[i].i = 0.0;
4413 }
4414
4415 /* S MitsuakiKAWAMURA*/
4416 HCC = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c);
4417 VCC = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_c);
4418 for (i = 0; i<NUM_c*NUM_c; i++) {
4419 HCC[i].r = 0.0;
4420 HCC[i].i = 0.0;
4421 VCC[i].r = 0.0;
4422 VCC[i].i = 0.0;
4423 }
4424 /* S MitsuakiKAWAMURA*/
4425
4426 HCL = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_e[0]);
4427 for (i=0; i<NUM_c*NUM_e[0]; i++) {
4428 HCL[i].r = 0.0;
4429 HCL[i].i = 0.0;
4430 }
4431
4432 HCR = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c*NUM_e[1]);
4433 for (i=0; i<NUM_c*NUM_e[1]; i++) {
4434 HCR[i].r = 0.0;
4435 HCR[i].i = 0.0;
4436 }
4437
4438 }
4439
4440 /*S MitsuakiKawaMura*/
MTRAN_Free_All_NC()4441 static void MTRAN_Free_All_NC()
4442 {
4443 int i, i2, i3, k, iside;
4444 int Gc_AN, Cwan, tno0;
4445 int h_AN;
4446
4447 for (iside = 0; iside <= 1; iside++) {
4448 free(S00_e[iside]);
4449 free(S01_e[iside]);
4450 free(H00_e[iside]);
4451 free(H01_e[iside]);
4452 }
4453
4454 free(SCC);
4455 free(HCC);
4456 free(VCC);
4457 free(SCL);
4458 free(SCR);
4459 free(HCL);
4460 free(HCR);
4461 /*
4462 Malloc in MTRAN_Input
4463 */
4464 if (tran_transmission_on) {
4465
4466 for (i2 = 0; i2 < TRAN_TKspace_grid2; i2++) {
4467 for (i3 = 0; i3 < TRAN_TKspace_grid3; i3++) {
4468 for (i = 0; i<3; i++) {
4469 free(tran_transmission[i2][i3][i]);
4470 }
4471 free(tran_transmission[i2][i3]);
4472 }
4473 free(tran_transmission[i2]);
4474 }
4475 free(tran_transmission);
4476 }
4477
4478 if (TRAN_Channel == 1) {
4479
4480 for (i = 0; i < TRAN_Channel_Nkpoint; i++) {
4481 free(TRAN_Channel_kpoint[i]);
4482 } /* for (i = 0; i<(TRAN_Channel_Nkpoint + 1); i++) */
4483 free(TRAN_Channel_kpoint);
4484
4485 free(TRAN_Channel_energy);
4486 }
4487 /*
4488 Malloc in MTRAN_Read_Tran_HS
4489 */
4490 free(TRAN_region);
4491 free(TRAN_Original_Id);
4492 free(WhatSpecies);
4493 free(Spe_Total_CNO);
4494 free(FNAN);
4495
4496 for (i = 0; i <= (atomnum); i++) {
4497 free(natn[i]);
4498 free(ncn[i]);
4499 free(atv_ijk[i]);
4500 }
4501 free(natn);
4502 free(ncn);
4503 free(atv_ijk);
4504
4505 for (k = 0; k<4; k++) {
4506
4507 FNAN[0] = 0;
4508 for (Gc_AN = 0; Gc_AN <= (atomnum); Gc_AN++) {
4509
4510 if (Gc_AN == 0) {
4511 tno0 = 1;
4512 }
4513 else {
4514 Cwan = WhatSpecies[Gc_AN];
4515 tno0 = Spe_Total_CNO[Cwan];
4516 }
4517
4518 for (h_AN = 0; h_AN <= FNAN[Gc_AN]; h_AN++) {
4519
4520 for (i = 0; i < tno0; i++) {
4521 free(OLP[k][Gc_AN][h_AN][i]);
4522 free(H[k][Gc_AN][h_AN][i]);
4523 free(H0[k][Gc_AN][h_AN][i]);
4524 free(iHNL[k][Gc_AN][h_AN][i]);
4525 }
4526 free(OLP[k][Gc_AN][h_AN]);
4527 free(H[k][Gc_AN][h_AN]);
4528 free(H0[k][Gc_AN][h_AN]);
4529 free(iHNL[k][Gc_AN][h_AN]);
4530 }
4531 free(OLP[k][Gc_AN]);
4532 free(H[k][Gc_AN]);
4533 free(H0[k][Gc_AN]);
4534 free(iHNL[k][Gc_AN]);
4535 }
4536 free(OLP[k]);
4537 free(H[k]);
4538 free(H0[k]);
4539 free(iHNL[k]);
4540 }
4541 free(OLP);
4542 free(H);
4543 free(H0);
4544 free(iHNL);
4545
4546 /**********************************************
4547 informations of leads
4548 **********************************************/
4549
4550 for (iside = 0; iside <= 1; iside++) {
4551
4552 for (k = 0; k<4; k++) {
4553
4554 FNAN_e[iside][0] = 0;
4555 for (Gc_AN = 0; Gc_AN <= atomnum_e[iside]; Gc_AN++) {
4556
4557 if (Gc_AN == 0) {
4558 tno0 = 1;
4559 }
4560 else {
4561 Cwan = WhatSpecies_e[iside][Gc_AN];
4562 Cwan = WhatSpecies_e[iside][Gc_AN];
4563 tno0 = Spe_Total_CNO_e[iside][Cwan];
4564 }
4565
4566 for (h_AN = 0; h_AN <= FNAN_e[iside][Gc_AN]; h_AN++) {
4567
4568 for (i = 0; i < tno0; i++) {
4569 free(OLP_e[iside][k][Gc_AN][h_AN][i]);
4570 free(H_e[iside][k][Gc_AN][h_AN][i]);
4571 free(iHNL_e[iside][k][Gc_AN][h_AN][i]);
4572 }
4573 free(OLP_e[iside][k][Gc_AN][h_AN]);
4574 free(H_e[iside][k][Gc_AN][h_AN]);
4575 free(iHNL_e[iside][k][Gc_AN][h_AN]);
4576 }
4577 free(OLP_e[iside][k][Gc_AN]);
4578 free(H_e[iside][k][Gc_AN]);
4579 free(iHNL_e[iside][k][Gc_AN]);
4580 }
4581 free(OLP_e[iside][k]);
4582 free(H_e[iside][k]);
4583 free(iHNL_e[iside][k]);
4584 }
4585 free(OLP_e[iside]);
4586 free(H_e[iside]);
4587 free(iHNL_e[iside]);
4588 }
4589
4590 for (iside = 0; iside <= 1; iside++) {
4591
4592 free(WhatSpecies_e[iside]);
4593 free(Spe_Total_CNO_e[iside]);
4594 free(FNAN_e[iside]);
4595 for (i = 0; i <= atomnum_e[iside]; i++) {
4596 free(natn_e[iside][i]);
4597 free(ncn_e[iside][i]);
4598 }
4599 free(natn_e[iside]);
4600 free(ncn_e[iside]);
4601
4602 for (i=0; i<(TCpyCell+1); i++) {
4603 free(atv_ijk_e[iside][i]);
4604 }
4605 free(atv_ijk_e[iside]);
4606 }
4607
4608 }
4609 /*E MitsuakiKawamura*/
4610
4611
4612
4613
4614
4615
4616
4617