1 /**********************************************************************
2 TRAN_Input_std_Atoms.c:
3
4 TRAN_Input_std_Atoms.c is a subroutine to read the input data.
5
6 Log of TRAN_Input_std_Atoms.c:
7
8 24/July/2008 Released by H.Kino and T.Ozaki
9
10 ***********************************************************************/
11 /* revised by Y. Xiao for Noncollinear NEGF calculations */
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <string.h>
16 #include <math.h>
17
18 #include "Inputtools.h"
19 #include "openmx_common.h"
20 #include <mpi.h>
21 #include "tran_prototypes.h"
22 #include "tran_variables.h"
23
24 #define MAXBUF 1024
25
26 int OrbPol2int(char OrbPol[YOUSO10]);
27
28
TRAN_Input_std_Atoms(MPI_Comm comm1,int Solver)29 void TRAN_Input_std_Atoms( MPI_Comm comm1, int Solver )
30 {
31 int po=0;
32 int idim=1;
33 int myid;
34
35 FILE *fp;
36 char *s_vec[20];
37 int i_vec[20];
38 double r_vec[20];
39
40 int i,j,spe,spe_e;
41 char Species[YOUSO10];
42 double GO_XR,GO_YR,GO_ZR;
43 double GO_XL,GO_YL,GO_ZL;
44 char OrbPol[YOUSO10];
45 char buf[MAXBUF];
46 /* revised by Y. Xiao for Noncollinear NEGF calculations */
47 double mx,my,mz,tmp,S_coordinate[3];
48 /*until here by Y. Xiao for Noncollinear NEGF calculations */
49
50 if (Solver!=4) return;
51
52 MPI_Comm_rank(comm1,&myid);
53
54 /* center */
55 input_int("Atoms.Number",&Catomnum,0);
56 if (Catomnum<=0){
57
58 if (myid==Host_ID){
59 printf("Atoms.Number may be wrong.\n");
60 }
61 MPI_Finalize();
62 exit(1);
63 }
64
65 /* left */
66 input_int("LeftLeadAtoms.Number",&Latomnum,0);
67 if (Latomnum<=0){
68
69 if (myid==Host_ID){
70 printf("LeftLeadAtoms.Number may be wrong.\n");
71 }
72 MPI_Finalize();
73 exit(1);
74 }
75
76 /* right */
77 input_int("RightLeadAtoms.Number",&Ratomnum,0);
78 if (Ratomnum<=0){
79
80 if (myid==Host_ID){
81 printf("RightLeadAtoms.Number may be wrong.\n");
82 }
83 MPI_Finalize();
84 exit(1);
85 }
86
87 atomnum = Catomnum + Latomnum + Ratomnum;
88 List_YOUSO[1] = atomnum + 1;
89
90 /* memory allocation */
91
92 Allocate_Arrays(1);
93
94 /* memory allocation for TRAN_* */
95 TRAN_Allocate_Atoms(atomnum);
96
97 s_vec[0]="Ang"; s_vec[1]="AU";
98 i_vec[0]= 0; i_vec[1]=1;
99 input_string2int("Atoms.SpeciesAndCoordinates.Unit",
100 &coordinates_unit,2,s_vec,i_vec);
101
102 /* left lead */
103
104 if (fp=input_find("<LeftLeadAtoms.SpeciesAndCoordinates") ) {
105
106 for (i=1; i<=Latomnum; i++){
107
108 fgets(buf,MAXBUF,fp);
109 /* revised by Y. Xiao for Noncollinear NEGF calculations */
110
111 /* spin non-collinear */
112 if (SpinP_switch==3){
113
114 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %s",
115 &j, Species,
116 &Gxyz[i][1],&Gxyz[i][2],&Gxyz[i][3],
117 &InitN_USpin[i],&InitN_DSpin[i],
118 &Angle0_Spin[i], &Angle1_Spin[i],
119 &Angle0_Orbital[i], &Angle1_Orbital[i],
120 &Constraint_SpinAngle[i],
121 OrbPol );
122
123 if (fabs(Angle0_Spin[i])<1.0e-10){
124 Angle0_Spin[i] = Angle0_Spin[i] + rnd(1.0e-5);
125 }
126
127 Angle0_Spin[i] = PI*Angle0_Spin[i]/180.0;
128 Angle1_Spin[i] = PI*Angle1_Spin[i]/180.0;
129 InitAngle0_Spin[i] = Angle0_Spin[i];
130 InitAngle1_Spin[i] = Angle1_Spin[i];
131
132 if (fabs(Angle0_Orbital[i])<1.0e-10){
133 Angle0_Orbital[i] = Angle0_Orbital[i] + rnd(1.0e-5);
134 }
135
136 Angle0_Orbital[i] = PI*Angle0_Orbital[i]/180.0;
137 Angle1_Orbital[i] = PI*Angle1_Orbital[i]/180.0;
138 InitAngle0_Orbital[i] = Angle0_Orbital[i];
139 InitAngle1_Orbital[i] = Angle1_Orbital[i];
140
141 /*check whether the Euler angle measured from the direction (1,0) is used*/
142 if ( (InitN_USpin[i]-InitN_DSpin[i]) < 0.0){
143
144 mx = -sin(InitAngle0_Spin[i])*cos(InitAngle1_Spin[i]);
145 my = -sin(InitAngle0_Spin[i])*sin(InitAngle1_Spin[i]);
146 mz = -cos(InitAngle0_Spin[i]);
147
148 xyz2spherical(mx,my,mz,0.0,0.0,0.0,S_coordinate);
149
150 Angle0_Spin[i] = S_coordinate[1];
151 Angle1_Spin[i] = S_coordinate[2];
152 InitAngle0_Spin[i] = Angle0_Spin[i];
153 InitAngle1_Spin[i] = Angle1_Spin[i];
154
155 tmp = InitN_USpin[i];
156 InitN_USpin[i] = InitN_DSpin[i];
157 InitN_DSpin[i] = tmp;
158
159 mx = -sin(InitAngle0_Orbital[i])*cos(InitAngle1_Orbital[i]);
160 my = -sin(InitAngle0_Orbital[i])*sin(InitAngle1_Orbital[i]);
161 mz = -cos(InitAngle0_Orbital[i]);
162
163 xyz2spherical(mx,my,mz,0.0,0.0,0.0,S_coordinate);
164
165 Angle0_Orbital[i] = S_coordinate[1];
166 Angle1_Orbital[i] = S_coordinate[2];
167 InitAngle0_Orbital[i] = Angle0_Orbital[i];
168 InitAngle1_Orbital[i] = Angle1_Orbital[i];
169
170 } /* if ( (InitN_USpin[i]-InitN_DSpin[i]) < 0.0) */
171 } /* if (SpinP_switch == 3) */
172
173 /* spin collinear */
174 else {
175
176 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %s",&j,Species,
177 &Gxyz[i][1],&Gxyz[i][2],&Gxyz[i][3],
178 &InitN_USpin[i],&InitN_DSpin[i], OrbPol);
179
180 }
181
182 WhatSpecies[i] = Species2int(Species);
183 TRAN_region[i]=2;
184 TRAN_Original_Id[i]=j;
185
186 if (Hub_U_switch==1) OrbPol_flag[i] = OrbPol2int(OrbPol);
187
188 /* check the consistency of basis set */
189
190 spe = WhatSpecies[i];
191 spe_e = WhatSpecies_e[0][i];
192
193 if (i!=j){
194
195 if (myid==Host_ID){
196 printf("Error of sequential number %i in <LeftLeadAtoms.SpeciesAndCoordinates\n",j);
197 }
198 MPI_Finalize();
199 exit(1);
200 }
201
202 if (2<=level_stdout && myid==Host_ID){
203
204 printf("<Input_std> L_AN=%2d T_AN=%2d WhatSpecies=%2d USpin=%7.4f DSpin=%7.4f\n",
205 i,i,
206 WhatSpecies[i],
207 InitN_USpin[i],
208 InitN_DSpin[i]);
209 }
210 }
211
212 ungetc('\n',fp);
213
214 if (!input_last("LeftLeadAtoms.SpeciesAndCoordinates>")) {
215 /* format error */
216
217 if (myid==Host_ID){
218 printf("Format error for LeftLeadAtoms.SpeciesAndCoordinates\n");
219 }
220 MPI_Finalize();
221 exit(1);
222 }
223 }
224
225 /* center */
226 if (fp=input_find("<Atoms.SpeciesAndCoordinates") ) {
227 for (i=1; i<=Catomnum; i++){
228
229 fgets(buf,MAXBUF,fp);
230 /* spin non-collinear */
231 if (SpinP_switch==3){
232
233 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %s",
234 &j, Species,
235 &Gxyz[Latomnum+i][1],&Gxyz[Latomnum+i][2],&Gxyz[Latomnum+i][3],
236 &InitN_USpin[Latomnum+i],&InitN_DSpin[Latomnum+i],
237 &Angle0_Spin[Latomnum+i], &Angle1_Spin[Latomnum+i],
238 &Angle0_Orbital[Latomnum+i], &Angle1_Orbital[Latomnum+i],
239 &Constraint_SpinAngle[Latomnum+i],
240 OrbPol );
241
242 if (fabs(Angle0_Spin[Latomnum+i])<1.0e-10){
243 Angle0_Spin[Latomnum+i] = Angle0_Spin[Latomnum+i] + rnd(1.0e-5);
244 }
245
246 Angle0_Spin[Latomnum+i] = PI*Angle0_Spin[Latomnum+i]/180.0;
247 Angle1_Spin[Latomnum+i] = PI*Angle1_Spin[Latomnum+i]/180.0;
248 InitAngle0_Spin[Latomnum+i] = Angle0_Spin[Latomnum+i];
249 InitAngle1_Spin[Latomnum+i] = Angle1_Spin[Latomnum+i];
250
251 if (fabs(Angle0_Orbital[Latomnum+i])<1.0e-10){
252 Angle0_Orbital[Latomnum+i] = Angle0_Orbital[Latomnum+i] + rnd(1.0e-5);
253 }
254
255 Angle0_Orbital[Latomnum+i] = PI*Angle0_Orbital[Latomnum+i]/180.0;
256 Angle1_Orbital[Latomnum+i] = PI*Angle1_Orbital[Latomnum+i]/180.0;
257 InitAngle0_Orbital[Latomnum+i] = Angle0_Orbital[Latomnum+i];
258 InitAngle1_Orbital[Latomnum+i] = Angle1_Orbital[Latomnum+i];
259
260 /*check whether the Euler angle measured from the direction (1,0) is used*/
261 if ( (InitN_USpin[Latomnum+i]-InitN_DSpin[Latomnum+i]) < 0.0){
262
263 mx = -sin(InitAngle0_Spin[Latomnum+i])*cos(InitAngle1_Spin[Latomnum+i]);
264 my = -sin(InitAngle0_Spin[Latomnum+i])*sin(InitAngle1_Spin[Latomnum+i]);
265 mz = -cos(InitAngle0_Spin[Latomnum+i]);
266
267 xyz2spherical(mx,my,mz,0.0,0.0,0.0,S_coordinate);
268
269 Angle0_Spin[Latomnum+i] = S_coordinate[1];
270 Angle1_Spin[Latomnum+i] = S_coordinate[2];
271 InitAngle0_Spin[Latomnum+i] = Angle0_Spin[Latomnum+i];
272 InitAngle1_Spin[Latomnum+i] = Angle1_Spin[Latomnum+i];
273
274 tmp = InitN_USpin[Latomnum+i];
275 InitN_USpin[Latomnum+i] = InitN_DSpin[Latomnum+i];
276 InitN_DSpin[Latomnum+i] = tmp;
277
278 mx = -sin(InitAngle0_Orbital[Latomnum+i])*cos(InitAngle1_Orbital[Latomnum+i]);
279 my = -sin(InitAngle0_Orbital[Latomnum+i])*sin(InitAngle1_Orbital[Latomnum+i]);
280 mz = -cos(InitAngle0_Orbital[Latomnum+i]);
281
282 xyz2spherical(mx,my,mz,0.0,0.0,0.0,S_coordinate);
283
284 Angle0_Orbital[Latomnum+i] = S_coordinate[1];
285 Angle1_Orbital[Latomnum+i] = S_coordinate[2];
286 InitAngle0_Orbital[Latomnum+i] = Angle0_Orbital[Latomnum+i];
287 InitAngle1_Orbital[Latomnum+i] = Angle1_Orbital[Latomnum+i];
288
289 } /* if ( (InitN_USpin[i]-InitN_DSpin[i]) < 0.0) */
290 } /* if (SpinP_switch == 3) */
291
292 /* spin collinear */
293 else {
294
295 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %s",
296 &j,Species,
297 &Gxyz[Latomnum+i][1],
298 &Gxyz[Latomnum+i][2],
299 &Gxyz[Latomnum+i][3],
300 &InitN_USpin[Latomnum+i],
301 &InitN_DSpin[Latomnum+i],
302 OrbPol);
303
304 }
305
306 WhatSpecies[Latomnum+i] = Species2int(Species);
307 TRAN_region[Latomnum+i]= 1;
308 TRAN_Original_Id[Latomnum+i]= j;
309
310 if (Hub_U_switch==1) OrbPol_flag[Latomnum+i] = OrbPol2int(OrbPol);
311
312 if (i!=j){
313
314 if (myid==Host_ID){
315 printf("Error of sequential number %i in <Atoms.SpeciesAndCoordinates\n",j);
316 }
317 MPI_Finalize();
318 exit(1);
319 }
320
321 if (2<=level_stdout && myid==Host_ID){
322 printf("<Input_std> ct_AN=%2d WhatSpecies=%2d USpin=%7.4f DSpin=%7.4f\n",
323 Latomnum+i,WhatSpecies[Latomnum+i],InitN_USpin[Latomnum+i],InitN_DSpin[Latomnum+i]);
324 }
325 }
326
327 ungetc('\n',fp);
328
329 if (!input_last("Atoms.SpeciesAndCoordinates>")) {
330 /* format error */
331 if (myid==Host_ID){
332 printf("Format error for Atoms.SpeciesAndCoordinates\n");
333 }
334 MPI_Finalize();
335 exit(1);
336
337 }
338 }
339
340 /* right */
341
342 if (fp=input_find("<RightLeadAtoms.SpeciesAndCoordinates") ) {
343
344 for (i=1; i<=Ratomnum; i++){
345
346 fgets(buf,MAXBUF,fp);
347 /* spin non-collinear */
348 if (SpinP_switch==3){
349
350 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %s",
351 &j, Species,
352 &Gxyz[Catomnum+Latomnum+i][1],&Gxyz[Catomnum+Latomnum+i][2],&Gxyz[Catomnum+Latomnum+i][3],
353 &InitN_USpin[Catomnum+Latomnum+i],&InitN_DSpin[Catomnum+Latomnum+i],
354 &Angle0_Spin[Catomnum+Latomnum+i], &Angle1_Spin[Catomnum+Latomnum+i],
355 &Angle0_Orbital[Catomnum+Latomnum+i], &Angle1_Orbital[Catomnum+Latomnum+i],
356 &Constraint_SpinAngle[Catomnum+Latomnum+i],
357 OrbPol );
358
359 if (fabs(Angle0_Spin[Catomnum+Latomnum+i])<1.0e-10){
360 Angle0_Spin[Catomnum+Latomnum+i] = Angle0_Spin[Catomnum+Latomnum+i] + rnd(1.0e-5);
361 }
362
363 Angle0_Spin[Catomnum+Latomnum+i] = PI*Angle0_Spin[Catomnum+Latomnum+i]/180.0;
364 Angle1_Spin[Catomnum+Latomnum+i] = PI*Angle1_Spin[Catomnum+Latomnum+i]/180.0;
365 InitAngle0_Spin[Catomnum+Latomnum+i] = Angle0_Spin[Catomnum+Latomnum+i];
366 InitAngle1_Spin[Catomnum+Latomnum+i] = Angle1_Spin[Catomnum+Latomnum+i];
367
368 if (fabs(Angle0_Orbital[Catomnum+Latomnum+i])<1.0e-10){
369 Angle0_Orbital[Catomnum+Latomnum+i] = Angle0_Orbital[Catomnum+Latomnum+i] + rnd(1.0e-5);
370 }
371
372 Angle0_Orbital[Catomnum+Latomnum+i] = PI*Angle0_Orbital[Catomnum+Latomnum+i]/180.0;
373 Angle1_Orbital[Catomnum+Latomnum+i] = PI*Angle1_Orbital[Catomnum+Latomnum+i]/180.0;
374 InitAngle0_Orbital[Catomnum+Latomnum+i] = Angle0_Orbital[Catomnum+Latomnum+i];
375 InitAngle1_Orbital[Catomnum+Latomnum+i] = Angle1_Orbital[Catomnum+Latomnum+i];
376
377 /*check whether the Euler angle measured from the direction (1,0) is used*/
378 if ( (InitN_USpin[Catomnum+Latomnum+i]-InitN_DSpin[Catomnum+Latomnum+i]) < 0.0){
379
380 mx = -sin(InitAngle0_Spin[Catomnum+Latomnum+i])*cos(InitAngle1_Spin[Catomnum+Latomnum+i]);
381 my = -sin(InitAngle0_Spin[Catomnum+Latomnum+i])*sin(InitAngle1_Spin[Catomnum+Latomnum+i]);
382 mz = -cos(InitAngle0_Spin[Catomnum+Latomnum+i]);
383
384 xyz2spherical(mx,my,mz,0.0,0.0,0.0,S_coordinate);
385
386 Angle0_Spin[Catomnum+Latomnum+i] = S_coordinate[1];
387 Angle1_Spin[Catomnum+Latomnum+i] = S_coordinate[2];
388 InitAngle0_Spin[Catomnum+Latomnum+i] = Angle0_Spin[Catomnum+Latomnum+i];
389 InitAngle1_Spin[Catomnum+Latomnum+i] = Angle1_Spin[Catomnum+Latomnum+i];
390
391 tmp = InitN_USpin[Catomnum+Latomnum+i];
392 InitN_USpin[Catomnum+Latomnum+i] = InitN_DSpin[Catomnum+Latomnum+i];
393 InitN_DSpin[Catomnum+Latomnum+i] = tmp;
394
395 mx = -sin(InitAngle0_Orbital[Catomnum+Latomnum+i])*cos(InitAngle1_Orbital[Catomnum+Latomnum+i]);
396 my = -sin(InitAngle0_Orbital[Catomnum+Latomnum+i])*sin(InitAngle1_Orbital[Catomnum+Latomnum+i]);
397 mz = -cos(InitAngle0_Orbital[Catomnum+Latomnum+i]);
398
399 xyz2spherical(mx,my,mz,0.0,0.0,0.0,S_coordinate);
400
401 Angle0_Orbital[Catomnum+Latomnum+i] = S_coordinate[1];
402 Angle1_Orbital[Catomnum+Latomnum+i] = S_coordinate[2];
403 InitAngle0_Orbital[Catomnum+Latomnum+i] = Angle0_Orbital[Catomnum+Latomnum+i];
404 InitAngle1_Orbital[Catomnum+Latomnum+i] = Angle1_Orbital[Catomnum+Latomnum+i];
405
406 } /* if ( (InitN_USpin[i]-InitN_DSpin[i]) < 0.0) */
407 } /* if (SpinP_switch == 3) */
408
409 /* spin collinear */
410 else {
411 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %s",
412 &j,Species,
413 &Gxyz[Catomnum+Latomnum+i][1],
414 &Gxyz[Catomnum+Latomnum+i][2],
415 &Gxyz[Catomnum+Latomnum+i][3],
416 &InitN_USpin[Catomnum+Latomnum+i],
417 &InitN_DSpin[Catomnum+Latomnum+i], OrbPol);
418
419 }
420
421 WhatSpecies[Catomnum+Latomnum+i] = Species2int(Species);
422
423 TRAN_region[Catomnum+Latomnum+i]= 3;
424 TRAN_Original_Id[Catomnum+Latomnum+i]= j;
425
426 if (Hub_U_switch==1) OrbPol_flag[Catomnum+Latomnum+i] = OrbPol2int(OrbPol);
427
428 if (i!=j){
429 if (myid==Host_ID){
430 printf("Error of sequential number %i in <RightLeadAtoms.SpeciesAndCoordinates\n",j);
431 }
432 MPI_Finalize();
433 exit(1);
434 }
435
436 if (2<=level_stdout && myid==Host_ID){
437 printf("<Input_std> R_AN=%2d T_AN=%2d WhatSpecies=%2d USpin=%7.4f DSpin=%7.4f\n",
438 i,Catomnum+Latomnum+i,
439 WhatSpecies[Catomnum+Latomnum+i],
440 InitN_USpin[Catomnum+Latomnum+i],
441 InitN_DSpin[Catomnum+Latomnum+i]);
442 }
443 }
444
445 ungetc('\n',fp);
446
447 if (!input_last("RightLeadAtoms.SpeciesAndCoordinates>")) {
448 /* format error */
449
450 if (myid==Host_ID){
451 printf("Format error for RightLeadAtoms.SpeciesAndCoordinates\n");
452 }
453 MPI_Finalize();
454 exit(1);
455
456 }
457 }
458
459 if (coordinates_unit==0){
460 for (i=1; i<=atomnum; i++){
461 Gxyz[i][1] = Gxyz[i][1]/BohrR;
462 Gxyz[i][2] = Gxyz[i][2]/BohrR;
463 Gxyz[i][3] = Gxyz[i][3]/BohrR;
464 }
465 }
466
467 /* compare the coordinates with those used for the band calculation of the left lead */
468
469 for (i=1; i<=Latomnum; i++){
470 printf("ABC1 X i=%2d %15.12f %15.12f %15.12f %15.12f\n",i,Gxyz_e[0][i][1],Gxyz[i][1],Gxyz[1][1],Gxyz_e[0][1][1]);
471 printf("ABC1 Y i=%2d %15.12f %15.12f %15.12f %15.12f\n",i,Gxyz_e[0][i][2],Gxyz[i][2],Gxyz[1][2],Gxyz_e[0][1][2]);
472 printf("ABC1 Z i=%2d %15.12f %15.12f %15.12f %15.12f\n",i,Gxyz_e[0][i][3],Gxyz[i][3],Gxyz[1][3],Gxyz_e[0][1][3]);
473 }
474
475 MPI_Finalize();
476 exit(1);
477
478 for (i=1; i<=Latomnum; i++){
479
480 if ( 1.0e-12<fabs(Gxyz_e[0][i][1]-Gxyz[i][1]+Gxyz[1][1]-Gxyz_e[0][1][1])
481 || 1.0e-12<fabs(Gxyz_e[0][i][2]-Gxyz[i][2]+Gxyz[1][2]-Gxyz_e[0][1][2])
482 || 1.0e-12<fabs(Gxyz_e[0][i][3]-Gxyz[i][3]+Gxyz[1][3]-Gxyz_e[0][1][3]) ){
483
484 if (myid==Host_ID){
485 printf("\n\n");
486 printf("The LEFT lead cannot be superposed on the original cell even after the translation.\n");
487 printf("Check your atomic coordinates of the LEFT lead.\n\n");
488 }
489
490 MPI_Finalize();
491 exit(1);
492 }
493 }
494
495 /* compare the coordinates with those used for the band calculation of the right lead */
496
497 for (i=1; i<=Ratomnum; i++){
498
499 if ( 1.0e-12<fabs((Gxyz_e[1][i][1]+(Gxyz[Catomnum+Latomnum+1][1]-Gxyz_e[1][1][1]))-Gxyz[Catomnum+Latomnum+i][1])
500 || 1.0e-12<fabs((Gxyz_e[1][i][2]+(Gxyz[Catomnum+Latomnum+1][2]-Gxyz_e[1][1][2]))-Gxyz[Catomnum+Latomnum+i][2])
501 || 1.0e-12<fabs((Gxyz_e[1][i][3]+(Gxyz[Catomnum+Latomnum+1][3]-Gxyz_e[1][1][3]))-Gxyz[Catomnum+Latomnum+i][3]) ){
502
503 if (myid==Host_ID){
504 printf("\n\n");
505 printf("The RIGHT lead cannot be superposed on the original cell even after the translation.\n");
506 printf("Check your atomic coordinates of the RIGHT lead.\n\n");
507 }
508
509 MPI_Finalize();
510 exit(1);
511 }
512
513 }
514
515 /****************************************************
516 Unit cell
517 ****************************************************/
518
519 for (i=1; i<=3; i++){
520
521 Left_tv[i][1] = tv_e[0][i][1];
522 Left_tv[i][2] = tv_e[0][i][2];
523 Left_tv[i][3] = tv_e[0][i][3];
524
525 Right_tv[i][1] = tv_e[1][i][1];
526 Right_tv[i][2] = tv_e[1][i][2];
527 Right_tv[i][3] = tv_e[1][i][3];
528 }
529
530 for (i=2; i<=3; i++){
531 tv[i][1] = tv_e[0][i][1];
532 tv[i][2] = tv_e[0][i][2];
533 tv[i][3] = tv_e[0][i][3];
534 }
535
536 /*****************************************************
537 set the grid origin and the unit vectors for the NEGF
538 calculations so that the boundaries between leads and
539 the extended central region match with those in the
540 band calculations for the leads.
541 *****************************************************/
542
543 GO_XL = Gxyz[1][1] - Gxyz_e[0][1][1] + Grid_Origin_e[0][1];
544 GO_YL = Gxyz[1][2] - Gxyz_e[0][1][2] + Grid_Origin_e[0][2];
545 GO_ZL = Gxyz[1][3] - Gxyz_e[0][1][3] + Grid_Origin_e[0][3];
546
547 GO_XR = Gxyz[Catomnum+Latomnum+1][1] - Gxyz_e[1][1][1] + Grid_Origin_e[1][1];
548 GO_YR = Gxyz[Catomnum+Latomnum+1][2] - Gxyz_e[1][1][2] + Grid_Origin_e[1][2];
549 GO_ZR = Gxyz[Catomnum+Latomnum+1][3] - Gxyz_e[1][1][3] + Grid_Origin_e[1][3];
550
551 /* large cell = left lead + central region + right lead */
552
553 tv[idim][1] = GO_XR + Right_tv[idim][1] - GO_XL;
554 tv[idim][2] = GO_YR + Right_tv[idim][2] - GO_YL;
555 tv[idim][3] = GO_ZR + Right_tv[idim][3] - GO_ZL;
556
557 scf_fixed_origin[0] = GO_XL;
558 scf_fixed_origin[1] = GO_YL;
559 scf_fixed_origin[2] = GO_ZL;
560
561 /*
562 if (myid==0){
563
564 printf("Right_tv\n");
565 printf("%15.12f %15.12f %15.12f\n",Right_tv[1][1],Right_tv[1][2],Right_tv[1][3]);
566
567 printf("GO_XR=%15.12f GO_XL=%15.12f\n",GO_XR,GO_XL);
568 printf("GO_YR=%15.12f GO_YL=%15.12f\n",GO_YR,GO_YL);
569 printf("GO_ZR=%15.12f GO_ZL=%15.12f\n",GO_ZR,GO_ZL);
570
571 printf("tv\n");
572 printf("%15.12f %15.12f %15.12f\n",tv[1][1],tv[1][2],tv[1][3]);
573 printf("%15.12f %15.12f %15.12f\n",tv[2][1],tv[2][2],tv[2][3]);
574 printf("%15.12f %15.12f %15.12f\n",tv[3][1],tv[3][2],tv[3][3]);
575
576 }
577
578 MPI_Finalize();
579 exit(0);
580 */
581
582
583 }
584