1 /**********************************************************************
2 MTRAN_EigenChannel.c:
3
4 TRAN_Calc_EChannel.c is a subroutine to calculate the Eigen Channel
5 See PRB 76, 115117 (2007).
6 It is made by modifing TRAN_Calc_OneTransmission.c.
7
8 input: SigmaL, SigmaR, G_CC^R(w)
9 assuming z= w+i delta
10
11 Gamma(z) = i (Sigma_R(z) - Sigma_R^+(z))
12
13 A(z) = G_CC_R(z) Gamma(z) G_CC_R^+(z)
14
15 Log of MTRAN_EigenChannel.c:
16
17 xx/Xxx/2015 Released by M. Kawamura
18
19 ***********************************************************************/
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <math.h>
24 #include <mpi.h>
25 #include "tran_prototypes.h"
26
27 int TRAN_Calc_OrtSpace(
28 int NUM_c,
29 dcomplex *SCC,
30 dcomplex *rtS,
31 dcomplex *rtSinv
32 ); /* int TRAN_Calc_OrtSpace */
33
34 void TRAN_Calc_Linewidth(
35 int NUM_c,
36 dcomplex *SigmaL_R,
37 dcomplex *SigmaR_R
38 ); /* void TRAN_Calc_Linewidth */
39
40 /* G Gamma_L G^+ -> Sigma_L */
41
42 void TRAN_Calc_MatTrans(
43 int NUM_c,
44 dcomplex *SigmaL_R,
45 dcomplex *GC_R,
46 char *trans1,
47 char *trans2
48 ); /* void TRAN_Calc_MatTrans */
49
50 /* L\"owdin orthogonalization */
51
52 void TRAN_Calc_LowdinOrt(
53 int NUM_c,
54 dcomplex *SigmaL_R,
55 dcomplex *SigmaR_R,
56 int NUM_cs,
57 dcomplex *rtS,
58 dcomplex *rtSinv,
59 dcomplex *ALbar,
60 dcomplex *GamRbar
61 ); /* void TRAN_Calc_LowdinOrt */
62
63 /* Diagonalize ALbar -> its eigenvector */
64 /* Then scale it : Albar_{ml} -> ALbar * sqrt(eig_l) */
65
66 void TRAN_Calc_Diagonalize(
67 int NUM_cs,
68 dcomplex *ALbar,
69 double *eval,
70 int lscale
71 ); /* void TRAN_Calc_Diagonalize */
72
73 /* Transform eigenchannel into non-orthogonal basis space */
74
75 void TRAN_Calc_ChannelLCAO(
76 int NUM_c,
77 int NUM_cs,
78 dcomplex *rtSinv,
79 dcomplex *ALbar,
80 dcomplex *GamRbar,
81 double *eval,
82 dcomplex *GC_R,
83 int TRAN_Channel_Num,
84 dcomplex **EChannel,
85 double *eigentrans
86 ); /* void TRAN_Calc_ChannelLCAO */
87
88 /* Output EigenChannel in the basis space */
89
90 void TRAN_Output_ChannelLCAO(
91 int myid0,
92 int kloop,
93 int iw,
94 int ispin,
95 int NUM_c,
96 double *TRAN_Channel_kpoint,
97 double TRAN_Channel_energy,
98 double *eval,
99 dcomplex *GC_R,
100 double **eigentrans_sum
101 ); /* void TRAN_Output_ChannelLCAO */
102
MTRAN_EigenChannel(MPI_Comm comm1,int numprocs,int myid,int myid0,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_Channel_Nenergy,double * TRAN_Channel_energy,int TRAN_Channel_Num,int kloop,double * TRAN_Channel_kpoint,dcomplex **** EChannel,double *** eigentrans,double ** eigentrans_sum)103 void MTRAN_EigenChannel(
104 MPI_Comm comm1,
105 int numprocs,
106 int myid,
107 int myid0,
108 int SpinP_switch,
109 double ChemP_e[2],
110 int NUM_c,
111 int NUM_e[2],
112 dcomplex **H00_e[2],
113 dcomplex *S00_e[2],
114 dcomplex **H01_e[2],
115 dcomplex *S01_e[2],
116 dcomplex **HCC,
117 dcomplex **HCL,
118 dcomplex **HCR,
119 dcomplex *SCC,
120 dcomplex *SCL,
121 dcomplex *SCR,
122 double tran_surfgreen_iteration_max,
123 double tran_surfgreen_eps,
124 double tran_transmission_energyrange[3],
125 int TRAN_Channel_Nenergy,
126 double *TRAN_Channel_energy,
127 int TRAN_Channel_Num,
128 int kloop,
129 double *TRAN_Channel_kpoint,
130 dcomplex ****EChannel,
131 double ***eigentrans,
132 double **eigentrans_sum)
133 {
134 dcomplex w;
135 dcomplex *GRL, *GRR, *GC_R;
136 dcomplex *SigmaL_R, *SigmaR_R;
137 dcomplex *rtS, *rtSinv;
138 double *eval;
139 dcomplex *ALbar, *GamRbar;
140
141 int iw, ispin, iside;
142 int **iwIdx;
143 int Miw, Miwmax;
144 int NUM_cs;
145
146 /* allocate */
147 GRL = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[0] * NUM_e[0]);
148 GRR = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[1] * NUM_e[1]);
149
150 GC_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
151 SigmaL_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
152 SigmaR_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
153 rtS = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
154 rtSinv = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
155
156 /*parallel setup*/
157
158 iwIdx = (int**)malloc(sizeof(int*)*numprocs);
159 Miwmax = (TRAN_Channel_Nenergy) / numprocs + 1;
160 for (iw = 0; iw<numprocs; iw++) {
161 iwIdx[iw] = (int*)malloc(sizeof(int)*Miwmax);
162 }
163 TRAN_Distribute_Node_Idx(0, TRAN_Channel_Nenergy - 1, numprocs, Miwmax,
164 iwIdx); /* output */
165
166 /* L\owdin orthogonalization : obtain S^{1/2} & S^{-1/2} */
167
168 NUM_cs = TRAN_Calc_OrtSpace(NUM_c, SCC, rtS, rtSinv);
169 printf(" myid0 = %3d, #k : %4d, N_{ort} / N_{nonort} : %d / %d \n",
170 myid0, kloop, NUM_cs, NUM_c);
171
172 ALbar = (dcomplex*)malloc(sizeof(dcomplex)*NUM_cs* NUM_cs);
173 GamRbar = (dcomplex*)malloc(sizeof(dcomplex)*NUM_cs* NUM_cs);
174 eval = (double*)malloc(sizeof(double)*NUM_c);
175
176 /* parallel global iw 0:tran_transmission_energydiv-1 */
177 /* parallel local Miw 0:Miwmax-1 */
178 /* parallel variable iw=iwIdx[myid][Miw] */
179
180 for (Miw = 0; Miw<Miwmax; Miw++) {
181
182 iw = iwIdx[myid][Miw];
183
184 if (iw >= 0) {
185
186 w.r = TRAN_Channel_energy[iw] + ChemP_e[0];
187 w.i = tran_transmission_energyrange[2];
188
189 /*
190 printf("iw=%d of %d w= % 9.6e % 9.6e \n" ,iw, tran_transmission_energydiv, w.r,w.i);
191 */
192
193 for (ispin = 0; ispin <= SpinP_switch; ispin++) {
194
195 /*****************************************************************
196 Note that retarded and advanced Green functions and self energies
197 are not conjugate comlex in case of the k-dependent case.
198 **************************************************************/
199
200 /* in case of retarded ones */
201
202 iside = 0;
203 TRAN_Calc_SurfGreen_direct(w, NUM_e[iside], H00_e[iside][ispin], H01_e[iside][ispin],
204 S00_e[iside], S01_e[iside],
205 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRL);
206
207 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRL, NUM_c, HCL[ispin], SCL, SigmaL_R);
208
209 iside = 1;
210 TRAN_Calc_SurfGreen_direct(w, NUM_e[iside], H00_e[iside][ispin], H01_e[iside][ispin],
211 S00_e[iside], S01_e[iside],
212 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRR);
213
214 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRR, NUM_c, HCR[ispin], SCR, SigmaR_R);
215
216 TRAN_Calc_CentGreen(w, NUM_c, SigmaL_R, SigmaR_R, HCC[ispin], SCC, GC_R);
217
218 /* in case of advanced ones */
219 /* Advanced Green's function is Hermit conjugate of the retare one */
220
221 /* Gamma = i(Sigma - Sigma^+) -> Sigma */
222
223 TRAN_Calc_Linewidth(NUM_c, SigmaL_R, SigmaR_R);
224
225 /* G Gamma_L G^+ -> Sigma_L */
226
227 TRAN_Calc_MatTrans(NUM_c, SigmaL_R, GC_R, "N", "C");
228
229 /* L\"owdin orthogonalization */
230
231 TRAN_Calc_LowdinOrt(NUM_c, SigmaL_R, SigmaR_R, NUM_cs, rtS, rtSinv, ALbar, GamRbar);
232
233 /* Diagonalize ALbar -> its eigenvector */
234 /* Then scale it : Albar_{ml} -> ALbar * sqrt(eig_l) */
235
236 TRAN_Calc_Diagonalize(NUM_cs, ALbar, eval, 1);
237
238 /* Albar^+ Gambar Albar -> Gambar */
239
240 TRAN_Calc_MatTrans(NUM_cs, GamRbar, ALbar, "C", "N");
241
242 /* Diagonalize transformed GamRbar */
243
244 TRAN_Calc_Diagonalize(NUM_cs, GamRbar, eval, 0);
245
246 /* Transform eigenchannel into non-orthogonal basis space */
247
248 TRAN_Calc_ChannelLCAO(NUM_c, NUM_cs, rtSinv, ALbar, GamRbar, eval, GC_R, TRAN_Channel_Num,
249 EChannel[iw][ispin],eigentrans[iw][ispin]);
250
251 /* Output EigenChannel in the basis space */
252
253 TRAN_Output_ChannelLCAO(myid0,kloop, iw, ispin, NUM_c,
254 TRAN_Channel_kpoint, TRAN_Channel_energy[iw], eval, GC_R, eigentrans_sum);
255
256 } /* for ispin */
257 } /* if ( iw>=0 ) */
258 } /* iw */
259
260 /* freeing of arrays */
261 free(GC_R);
262 free(SigmaR_R);
263 free(SigmaL_R);
264
265 free(GRR);
266 free(GRL);
267
268 free(rtS);
269 free(rtSinv);
270 free(eval);
271
272 free(ALbar);
273 free(GamRbar);
274
275 for (iw = 0; iw<numprocs; iw++) {
276 free(iwIdx[iw]);
277 }
278 free(iwIdx);
279 } /* void MTRAN_EigenChannel */
280
MTRAN_EigenChannel_NC(MPI_Comm comm1,int numprocs,int myid,int myid0,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_Channel_Nenergy,double * TRAN_Channel_energy,int TRAN_Channel_Num,int kloop,double * TRAN_Channel_kpoint,dcomplex *** EChannel,double ** eigentrans,double ** eigentrans_sum)281 void MTRAN_EigenChannel_NC(
282 MPI_Comm comm1,
283 int numprocs,
284 int myid,
285 int myid0,
286 int SpinP_switch,
287 double ChemP_e[2],
288 int NUM_c,
289 int NUM_e[2],
290 dcomplex *H00_e[2],
291 dcomplex *S00_e[2],
292 dcomplex *H01_e[2],
293 dcomplex *S01_e[2],
294 dcomplex *HCC,
295 dcomplex *HCL,
296 dcomplex *HCR,
297 dcomplex *SCC,
298 dcomplex *SCL,
299 dcomplex *SCR,
300 double tran_surfgreen_iteration_max,
301 double tran_surfgreen_eps,
302 double tran_transmission_energyrange[3],
303 int TRAN_Channel_Nenergy,
304 double *TRAN_Channel_energy,
305 int TRAN_Channel_Num,
306 int kloop,
307 double *TRAN_Channel_kpoint,
308 dcomplex ***EChannel,
309 double **eigentrans,
310 double **eigentrans_sum)
311 {
312 dcomplex w;
313 dcomplex *GRL, *GRR, *GC_R;
314 dcomplex *SigmaL_R, *SigmaR_R;
315 dcomplex *rtS, *rtSinv;
316 double *eval;
317 dcomplex *ALbar, *GamRbar;
318
319 int iw, iside;
320 int **iwIdx;
321 int Miw, Miwmax;
322 int NUM_cs;
323
324 /* allocate */
325 GRL = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[0] * NUM_e[0]);
326 GRR = (dcomplex*)malloc(sizeof(dcomplex)*NUM_e[1] * NUM_e[1]);
327
328 GC_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
329 SigmaL_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
330 SigmaR_R = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
331 rtS = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
332 rtSinv = (dcomplex*)malloc(sizeof(dcomplex)*NUM_c* NUM_c);
333
334 /*parallel setup*/
335
336 iwIdx = (int**)malloc(sizeof(int*)*numprocs);
337 Miwmax = (TRAN_Channel_Nenergy) / numprocs + 1;
338 for (iw = 0; iw<numprocs; iw++) {
339 iwIdx[iw] = (int*)malloc(sizeof(int)*Miwmax);
340 }
341 TRAN_Distribute_Node_Idx(0, TRAN_Channel_Nenergy - 1, numprocs, Miwmax,
342 iwIdx); /* output */
343
344 /* L\owdin orthogonalization : obtain S^{1/2} & S^{-1/2} */
345
346 NUM_cs = TRAN_Calc_OrtSpace(NUM_c, SCC, rtS, rtSinv);
347 printf(" myid0 = %3d, #k : %4d, N_{ort} / N_{nonort} : %d / %d \n",
348 myid0, kloop, NUM_cs, NUM_c);
349
350 ALbar = (dcomplex*)malloc(sizeof(dcomplex)*NUM_cs* NUM_cs);
351 GamRbar = (dcomplex*)malloc(sizeof(dcomplex)*NUM_cs* NUM_cs);
352 eval = (double*)malloc(sizeof(double)*NUM_c);
353
354 /* parallel global iw 0:tran_transmission_energydiv-1 */
355 /* parallel local Miw 0:Miwmax-1 */
356 /* parallel variable iw=iwIdx[myid][Miw] */
357
358 for (Miw = 0; Miw<Miwmax; Miw++) {
359
360 iw = iwIdx[myid][Miw];
361
362 if (iw >= 0) {
363
364 w.r = TRAN_Channel_energy[iw] + ChemP_e[0];
365 w.i = tran_transmission_energyrange[2];
366
367 /*
368 printf("iw=%d of %d w= % 9.6e % 9.6e \n" ,iw, tran_transmission_energydiv, w.r,w.i);
369 */
370
371 /*****************************************************************
372 Note that retarded and advanced Green functions and self energies
373 are not conjugate comlex in case of the k-dependent case.
374 **************************************************************/
375
376 /* in case of retarded ones */
377
378 iside = 0;
379 TRAN_Calc_SurfGreen_direct(w, NUM_e[iside], H00_e[iside], H01_e[iside],
380 S00_e[iside], S01_e[iside],
381 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRL);
382
383 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRL, NUM_c, HCL, SCL, SigmaL_R);
384
385 iside = 1;
386 TRAN_Calc_SurfGreen_direct(w, NUM_e[iside], H00_e[iside], H01_e[iside],
387 S00_e[iside], S01_e[iside],
388 tran_surfgreen_iteration_max, tran_surfgreen_eps, GRR);
389
390 TRAN_Calc_SelfEnergy(w, NUM_e[iside], GRR, NUM_c, HCR, SCR, SigmaR_R);
391
392 TRAN_Calc_CentGreen(w, NUM_c, SigmaL_R, SigmaR_R, HCC, SCC, GC_R);
393
394 /* in case of advanced ones */
395 /* Advanced Green's function is Hermit conjugate of the retare one */
396
397 /* Gamma = i(Sigma - Sigma^+) -> Sigma */
398
399 TRAN_Calc_Linewidth(NUM_c, SigmaL_R, SigmaR_R);
400
401 /* G Gamma_L G^+ -> Sigma_L */
402
403 TRAN_Calc_MatTrans(NUM_c, SigmaL_R, GC_R, "N", "C");
404
405 /* L\"owdin orthogonalization */
406
407 TRAN_Calc_LowdinOrt(NUM_c, SigmaL_R, SigmaR_R, NUM_cs, rtS, rtSinv, ALbar, GamRbar);
408
409 /* Diagonalize ALbar -> its eigenvector */
410 /* Then scale it : Albar_{ml} -> ALbar * sqrt(eig_l) */
411
412 TRAN_Calc_Diagonalize(NUM_cs, ALbar, eval, 1);
413
414 /* Albar^+ Gambar Albar -> Gambar */
415
416 TRAN_Calc_MatTrans(NUM_cs, GamRbar, ALbar, "C", "N");
417
418 /* Diagonalize transformed GamRbar */
419
420 TRAN_Calc_Diagonalize(NUM_cs, GamRbar, eval, 0);
421
422 /* Transform eigenchannel into non-orthogonal basis space */
423
424 TRAN_Calc_ChannelLCAO(NUM_c, NUM_cs, rtSinv, ALbar, GamRbar, eval, GC_R, TRAN_Channel_Num,
425 EChannel[iw],eigentrans[iw]);
426
427 /* Output EigenChannel in the basis space */
428
429 TRAN_Output_ChannelLCAO(myid0, kloop, iw, 0, NUM_c,
430 TRAN_Channel_kpoint, TRAN_Channel_energy[iw], eval, GC_R, eigentrans_sum);
431
432 } /* if ( iw>=0 ) */
433 } /* iw */
434
435 /* freeing of arrays */
436 free(GC_R);
437 free(SigmaR_R);
438 free(SigmaL_R);
439
440 free(GRR);
441 free(GRL);
442
443 free(rtS);
444 free(rtSinv);
445 free(eval);
446
447 free(ALbar);
448 free(GamRbar);
449
450 for (iw = 0; iw<numprocs; iw++) {
451 free(iwIdx[iw]);
452 }
453 free(iwIdx);
454 } /* void MTRAN_EigenChannel_NC */
455
456