1 #include <math.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <ctype.h>
6 #include "openmx_common.h"
7 #include "Inputtools.h"
8 #include "mpi.h"
9 #include <omp.h>
10
11
12
13 #include "tran_prototypes.h"
14
15 #define MAXBUF 1024
16 #define Max_Num_WF_Projs 15
17
18
19 void SpeciesString2int(int p);
20 void kpath_changeunit(double tv[4][4],double tv0[4][4],int Band_Nkpath,
21 double ***Band_kpath);
22 void kpoint_changeunit(double tv[4][4],double tv0[4][4],int MO_Nkpoint,
23 double **MO_kpoint);
24 void Set_Cluster_UnitCell(double tv[4][4],int unitflag);
25
26 int OrbPol2int(char OrbPol[YOUSO10]);
27 char *ToCapital(char *s);
28 int divisible_cheker(int N);
29 static void Set_In_First_Cell();
30 static void Remake_RestartFile(int numprocs_new, int numprocs_old, int N1, int N2, int N3);
31
32 /* hmweng */
33 int Analyze_Wannier_Projectors(int p, char ctmp[YOUSO10],
34 int **tmp_Wannier_Pro_SelMat,
35 double ***tmp_Wannier_Projector_Hybridize_Matrix);
36 void Get_Rotational_Matrix(double alpha, double beta, double gamma, int L, double tmpRotMat[7][7]);
37 int Calc_Factorial(int arg);
38 void Get_Euler_Rotation_Angle(
39 double zx, double zy, double zz,
40 double xx, double xy, double xz,
41 double *alpha_r, double *beta_r, double *gamma_r);
42
43
44
Input_std(char * file)45 void Input_std(char *file)
46 {
47 FILE *fp,*fp_check;
48 int i,j,k,itmp;
49 int num_wannier_total_projectors;
50 int l,mul; /* added by MJ */
51 int po=0; /* error count */
52 double r_vec[40];
53 int i_vec[40],i_vec2[40];
54 char *s_vec[40],Species[YOUSO10];
55 char OrbPol[YOUSO10];
56 double ecutoff1dfft;
57 double mx,my,mz,tmp;
58 double tmpx,tmpy,tmpz;
59 double S_coordinate[3];
60 double length,x,y,z;
61 int orbitalopt;
62 char buf[MAXBUF];
63 char file_check[YOUSO10];
64 int numprocs,myid;
65 int output_hks;
66 int numprocs1;
67
68 MPI_Comm_size(MPI_COMM_WORLD1,&numprocs);
69 MPI_Comm_rank(MPI_COMM_WORLD1,&myid);
70
71 if (myid==Host_ID && 0<level_stdout){
72 printf("*******************************************************\n");
73 printf(" read the input file and initializing \n");
74 printf("*******************************************************\n\n");
75 }
76
77 /****************************************************
78 open a file
79 ****************************************************/
80
81 if (input_open(file)==0){
82 MPI_Finalize();
83 exit(0);
84 }
85
86 input_string("System.CurrrentDirectory",filepath,"./");
87 input_string("System.Name",filename,"default");
88 //input_string("DATA.PATH","/usr/local/share/openmx/DFT_DATA13","../DFT_DATA13");
89 input_int("level.of.stdout", &level_stdout,1);
90 input_int("level.of.fileout",&level_fileout,1);
91 input_logical("memory.usage.fileout",&memoryusage_fileout,0); /* default=off */
92
93 if (level_stdout<0 || 3<level_stdout){
94 printf("Invalid value of level.of.stdout\n");
95 po++;
96 }
97
98 if (level_fileout<0 || 3<level_fileout){
99 printf("Invalid value of level.of.fileout\n");
100 po++;
101 }
102
103 /****************************************************
104 projector expansion of VNA
105 ****************************************************/
106
107 input_logical("scf.ProExpn.VNA",&ProExpn_VNA,1); /* default=on */
108 input_int("scf.BufferL.VNA", &BufferL_ProVNA,6);
109 input_int("scf.RadialF.VNA", &List_YOUSO[34],12);
110
111 /****************************************************
112 cutoff energy
113 ****************************************************/
114
115 /* for cutoff energy */
116
117 input_double("scf.energycutoff",&Grid_Ecut,(double)150.0);
118 input_logical("scf.MPI.tuned.grids",&MPI_tunedgrid_flag,0);
119
120 /* for fixed Ngrids */
121
122 i_vec2[0]=0;
123 i_vec2[1]=0;
124 i_vec2[2]=0;
125 input_intv("scf.Ngrid",3,i_vec,i_vec2);
126 Ngrid1 = i_vec[0];
127 Ngrid2 = i_vec[1];
128 Ngrid3 = i_vec[2];
129
130 if (Ngrid1==0 && Ngrid2==0 && Ngrid3==0)
131 Ngrid_fixed_flag = 0;
132 else
133 Ngrid_fixed_flag = 1;
134
135 if (Ngrid_fixed_flag==1){
136 i = divisible_cheker(Ngrid1);
137 j = divisible_cheker(Ngrid2);
138 k = divisible_cheker(Ngrid3);
139
140 if ( (i*j*k)==0 ) {
141
142 printf("scf.Ngrid must be divisible by \n");
143
144 printf(" ");
145 for (i=0; i<NfundamentalNum; i++){
146 printf("%3d ",fundamentalNum[i]);
147 }
148 printf("\n");
149
150 MPI_Finalize();
151 exit(0);
152 }
153 }
154
155 /****************************************************
156 definition of atomic species
157 ****************************************************/
158
159 input_int("Species.Number",&SpeciesNum,0);
160 real_SpeciesNum = SpeciesNum;
161
162 if (SpeciesNum<=0){
163 printf("Species.Number may be wrong.\n");
164 po++;
165 }
166 List_YOUSO[18] = SpeciesNum;
167
168 /* memory allocation */
169 Allocate_Arrays(0);
170
171 /*************************************************************
172 for LDA+U
173 Hub_U_switch should be called before Allocate_Arrays(1);
174 *************************************************************/
175
176 input_logical("scf.Hubbard.U",&Hub_U_switch, 0); /* --- MJ */
177
178 /* default Hub_U_occupation = 2; */
179
180 s_vec[0]="DUAL"; i_vec[0]=2;
181 s_vec[1]="ONSITE"; i_vec[1]=0;
182 s_vec[2]="FULL" ; i_vec[2]=1;
183
184 input_string2int("scf.Hubbard.Occupation",&Hub_U_occupation, 3, s_vec,i_vec);
185
186 /****************************************************
187 Orbital optimization
188 ****************************************************/
189
190 s_vec[0]="OFF"; s_vec[1]="Atoms"; s_vec[2]="Species"; s_vec[3]="Atoms2"; s_vec[4]="Species2";
191 i_vec[0]=0; i_vec[1]=1; i_vec[2]=2; i_vec[3]=3; i_vec[4]=4;
192 input_string2int("orbitalOpt.Method",&orbitalopt,5,s_vec,i_vec);
193
194 switch (orbitalopt) {
195 case 0: { Cnt_switch=0; } break;
196 case 1: { Cnt_switch=1; RCnt_switch=1; SCnt_switch=0; } break;
197 case 2: { Cnt_switch=1; RCnt_switch=1; ACnt_switch=1; SCnt_switch=0; } break;
198 case 3: { Cnt_switch=1; RCnt_switch=1; SCnt_switch=1; } break;
199 case 4: { Cnt_switch=1; RCnt_switch=1; ACnt_switch=1; SCnt_switch=1; } break;
200 }
201
202 /*************************************************************
203 read species
204 *************************************************************/
205
206 if (fp=input_find("<Definition.of.Atomic.Species")) {
207
208 for (i=0; i<SpeciesNum; i++){
209 fgets(buf,MAXBUF,fp);
210 sscanf(buf,"%s %s %s %lf",SpeName[i],SpeBasis[i],SpeVPS[i],&Spe_AtomicMass[i]);
211 SpeciesString2int(i);
212 }
213
214 ungetc('\n',fp);
215
216 if (! input_last("Definition.of.Atomic.Species>")) {
217 /* format error */
218
219 po++;
220
221 if (myid==Host_ID){
222 printf("Format error for Definition.of.Atomic.Species\n");
223 }
224 MPI_Finalize();
225 exit(0);
226 }
227 }
228
229 if (2<=level_stdout){
230 for (i=0; i<SpeciesNum; i++){
231 printf("<Input_std> %i Name %s\n",i,SpeName[i]);
232 printf("<Input_std> %i Basis %s\n",i,SpeBasis[i]);
233 printf("<Input_std> %i VPS %s\n",i,SpeVPS[i]);
234 }
235 }
236
237 List_YOUSO[35] = 0;
238 for (i=0; i<SpeciesNum; i++){
239 if (List_YOUSO[35]<Spe_MaxL_Basis[i]) List_YOUSO[35] = Spe_MaxL_Basis[i];
240 }
241 List_YOUSO[35] = List_YOUSO[35] + BufferL_ProVNA;
242
243 /****************************************************
244 change of atomic mass
245 ****************************************************/
246
247 /*
248 if (fp=input_find("<Mass.of.Atomic.Species")) {
249
250 for (i=0; i<SpeciesNum; i++){
251 fscanf(fp,"%s %s %s",SpeName[i],SpeBasis[i],SpeVPS[i]);
252
253 }
254
255 if (! input_last("Mass.of.Atomic.Species>")) {
256
257 po++;
258 if (myid==Host_ID){
259 printf("Format error for Mass.of.Atomic.Species\n");
260 }
261 MPI_Finalize();
262 exit(0);
263 }
264 }
265 */
266
267 /******************************************************
268 optimization for cubic or tetragonal cell
269 by a fitting procedure
270 ******************************************************/
271
272 i=0;
273 s_vec[i]="off"; i_vec[i]=0; i++; /* off */
274 s_vec[i]="CellCub"; i_vec[i]=1; i++; /* opt. of cubic cell */
275 s_vec[i]="CellTet"; i_vec[i]=2; i++; /* opt. of tetragonal cell */
276
277 j = input_string2int("Cell.Opt",&CellOpt_switch, i, s_vec,i_vec);
278 if (j==-1){
279 MPI_Finalize();
280 exit(0);
281 }
282
283 /****************************************************
284 Molecular dynamics or geometry optimization
285 ****************************************************/
286
287 i=0;
288 s_vec[i]="NOMD"; i_vec[i]=0; i++;
289 s_vec[i]="NVE" ; i_vec[i]=1; i++;
290 s_vec[i]="NVT_VS"; i_vec[i]=2; i++; /* modified by mari */
291 s_vec[i]="Opt"; i_vec[i]=3; i++;
292 s_vec[i]="EF"; i_vec[i]=4; i++;
293 s_vec[i]="BFGS"; i_vec[i]=5; i++;
294 s_vec[i]="RF"; i_vec[i]=6; i++; /* RF method by hmweng */
295 s_vec[i]="DIIS"; i_vec[i]=7; i++;
296 s_vec[i]="Constraint_DIIS"; i_vec[i]=8; i++; /* not used */
297 s_vec[i]="NVT_NH"; i_vec[i]=9; i++;
298 s_vec[i]="Opt_LBFGS"; i_vec[i]=10; i++;
299 s_vec[i]="NVT_VS2"; i_vec[i]=11; i++; /* modified by Ohwaki */
300 s_vec[i]="EvsLC"; i_vec[i]=12; i++;
301 s_vec[i]="NEB"; i_vec[i]=13; i++;
302 s_vec[i]="NVT_VS4"; i_vec[i]=14; i++; /* modified by Ohwaki */
303 s_vec[i]="NVT_Langevin"; i_vec[i]=15; i++; /* modified by Ohwaki */
304 s_vec[i]="DF"; i_vec[i]=16; i++; /* delta-factor */
305 s_vec[i]="OptC1"; i_vec[i]=17; i++; /* cell opt with fixed fractional coordinates by SD */
306 s_vec[i]="OptC2"; i_vec[i]=18; i++; /* cell opt with fixed fractional coordinates and angles fixed by SD */
307 s_vec[i]="OptC3"; i_vec[i]=19; i++; /* cell opt with fixed fractional coordinates, angles fixed and |a1|=|a2|=|a3| by SD */
308 s_vec[i]="OptC4"; i_vec[i]=20; i++; /* cell opt with fixed fractional coordinates, angles fixed and |a1|=|a2|!=|a3| by SD */
309 s_vec[i]="OptC5"; i_vec[i]=21; i++; /* cell opt with no constraint for cell and coordinates by SD */
310 s_vec[i]="RFC1"; i_vec[i]=22; i++; /* cell opt with fixed fractional coordinates by RF */
311 s_vec[i]="RFC2"; i_vec[i]=23; i++; /* cell opt with fixed fractional coordinates and angles fixed by RF */
312 s_vec[i]="RFC3"; i_vec[i]=24; i++; /* cell opt with fixed fractional coordinates, angles fixed and |a1|=|a2|=|a3| by RF */
313 s_vec[i]="RFC4"; i_vec[i]=25; i++; /* cell opt with fixed fractional coordinates, angles fixed and |a1|=|a2|!=|a3| by RF */
314 s_vec[i]="RFC5"; i_vec[i]=26; i++; /* cell opt with no constraint for cell and coordinates by RF */
315
316 j = input_string2int("MD.Type",&MD_switch, i, s_vec,i_vec);
317
318 if (j==-1){
319 MPI_Finalize();
320 exit(0);
321 }
322
323 /* cell optimization by SD */
324
325 if (MD_switch==17){
326 MD_switch = 17;
327 cellopt_swtich = 1; /* OptC1: no constraint for cell vectors */
328 }
329
330 else if (MD_switch==18){
331 MD_switch = 17;
332 cellopt_swtich = 2; /* OptC2: angles fixed for cell vectors */
333 }
334
335 else if (MD_switch==19){
336
337 MD_switch = 17;
338 cellopt_swtich = 3; /* OptC3: angles fixed and |a1|=|a2|=|a3| for cell vectors */
339
340 double a1,a2,a3;
341
342 a1 = tv[1][1]*tv[1][1]+tv[1][2]*tv[1][2]+tv[1][3]*tv[1][3];
343 a2 = tv[2][1]*tv[2][1]+tv[2][2]*tv[2][2]+tv[2][3]*tv[2][3];
344 a3 = tv[3][1]*tv[3][1]+tv[3][2]*tv[3][2]+tv[3][3]*tv[3][3];
345
346 if ( 0.000000001<fabs(a1-a2) || 0.000000001<fabs(a1-a3) ){
347
348 if (myid==Host_ID){
349 printf("\nThe condition |a1|=|a2|=|a3| should be satisfied for the use of OptC3\n\n");
350 }
351
352 MPI_Finalize();
353 exit(0);
354 }
355 }
356
357 else if (MD_switch==20){
358
359 MD_switch = 17;
360 cellopt_swtich = 4; /* OptC4: angles fixed and |a1|=|a2|!=|a3| for cell vectors */
361
362 double a1,a2,a3;
363
364 a1 = tv[1][1]*tv[1][1]+tv[1][2]*tv[1][2]+tv[1][3]*tv[1][3];
365 a2 = tv[2][1]*tv[2][1]+tv[2][2]*tv[2][2]+tv[2][3]*tv[2][3];
366 a3 = tv[3][1]*tv[3][1]+tv[3][2]*tv[3][2]+tv[3][3]*tv[3][3];
367
368 if ( 0.000000001<fabs(a1-a2) ){
369
370 if (myid==Host_ID){
371 printf("\nThe condition |a1|=|a2| should be satisfied for the use of OptC4\n\n");
372 }
373
374 MPI_Finalize();
375 exit(0);
376 }
377 }
378
379 else if (MD_switch==21){
380 MD_switch = 17;
381 cellopt_swtich = 5; /* OptC5: no constraint for cell and coordinates */
382 }
383
384 /* cell optimization by RF */
385
386 if (MD_switch==22){
387 MD_switch = 18;
388 cellopt_swtich = 1; /* RFC1: no constraint for cell vectors */
389
390 if (myid==Host_ID){
391 printf("The optimizer is not supported.\n");
392 }
393
394 MPI_Finalize();
395 exit(1);
396 }
397
398 else if (MD_switch==23){
399 MD_switch = 18;
400 cellopt_swtich = 2; /* RFC2: angles fixed for cell vectors */
401
402 if (myid==Host_ID){
403 printf("The optimizer is not supported.\n");
404 }
405
406 MPI_Finalize();
407 exit(1);
408 }
409
410 else if (MD_switch==24){
411
412 MD_switch = 18;
413 cellopt_swtich = 3; /* RFC3: angles fixed and |a1|=|a2|=|a3| for cell vectors */
414
415 double a1,a2,a3;
416
417 a1 = tv[1][1]*tv[1][1]+tv[1][2]*tv[1][2]+tv[1][3]*tv[1][3];
418 a2 = tv[2][1]*tv[2][1]+tv[2][2]*tv[2][2]+tv[2][3]*tv[2][3];
419 a3 = tv[3][1]*tv[3][1]+tv[3][2]*tv[3][2]+tv[3][3]*tv[3][3];
420
421 if ( 0.000000001<fabs(a1-a2) || 0.000000001<fabs(a1-a3) ){
422
423 if (myid==Host_ID){
424 printf("\nThe condition |a1|=|a2|=|a3| should be satisfied for the use of OptC3\n\n");
425 }
426
427 MPI_Finalize();
428 exit(0);
429 }
430
431 if (myid==Host_ID){
432 printf("The optimizer is not supported.\n");
433 }
434
435 MPI_Finalize();
436 exit(1);
437 }
438
439 else if (MD_switch==25){
440
441 MD_switch = 18;
442 cellopt_swtich = 4; /* RFC4: angles fixed and |a1|=|a2|!=|a3| for cell vectors */
443
444 double a1,a2,a3;
445
446 a1 = tv[1][1]*tv[1][1]+tv[1][2]*tv[1][2]+tv[1][3]*tv[1][3];
447 a2 = tv[2][1]*tv[2][1]+tv[2][2]*tv[2][2]+tv[2][3]*tv[2][3];
448 a3 = tv[3][1]*tv[3][1]+tv[3][2]*tv[3][2]+tv[3][3]*tv[3][3];
449
450 if ( 0.000000001<fabs(a1-a2) ){
451
452 if (myid==Host_ID){
453 printf("\nThe condition |a1|=|a2| should be satisfied for the use of OptC4\n\n");
454 }
455
456 MPI_Finalize();
457 exit(0);
458 }
459
460 if (myid==Host_ID){
461 printf("The optimizer is not supported.\n");
462 }
463
464 MPI_Finalize();
465 exit(1);
466 }
467
468 else if (MD_switch==26){
469 MD_switch = 18;
470 cellopt_swtich = 5; /* RFC5: no constraint for cell and coordinates */
471 }
472
473 /* MD.maxIter */
474
475 input_int("MD.maxIter",&MD_IterNumber,1);
476 if (MD_IterNumber<1){
477 printf("MD_IterNumber=%i should be over 0.\n",MD_IterNumber);
478 po++;
479 }
480
481 if (MD_switch==16 && MD_IterNumber!=7){ /* delta-factor */
482 printf("MD_IterNumber=%i should be 7 for the delta-factor calculation.\n",MD_IterNumber);
483 po++;
484 }
485
486 input_int("MD.Current.Iter",&MD_Current_Iter,0);
487
488 input_double("MD.TimeStep",&MD_TimeStep,(double)0.5);
489 if (MD_TimeStep<0.0){
490 printf("MD.TimeStep=%lf should be over 0.\n",MD_TimeStep);
491 po++;
492 }
493
494 input_double("MD.Opt.criterion",&MD_Opt_criterion,(double)0.0003);
495 input_int("MD.Opt.DIIS.History",&M_GDIIS_HISTORY,3);
496 input_int("MD.Opt.StartDIIS",&OptStartDIIS,5);
497 input_int("MD.Opt.EveryDIIS",&OptEveryDIIS,200);
498
499 input_double("MD.EvsLC.Step",&MD_EvsLattice_Step,(double)0.4);
500
501 i_vec[0]=1; i_vec[1]=1, i_vec[2]=1;
502 input_intv("MD.EvsLC.flag",3,MD_EvsLattice_flag,i_vec);
503
504 input_logical("MD.Out.ABC",&MD_OutABC,0); /* default=off */
505
506 /*
507 input_double("MD.Initial.MaxStep",&SD_scaling_user,(double)0.02);
508 */
509
510 i=0;
511 s_vec[i]="Schlegel" ; i_vec[i]=1; i++;
512 s_vec[i]="iden"; i_vec[i]=0; i++;
513 s_vec[i]="FF" ; i_vec[i]=2; i++;
514 j = input_string2int("MD.Opt.Init.Hessian",&Initial_Hessian_flag, i, s_vec,i_vec);
515
516 /* Ang -> a.u. */
517 /*
518 SD_scaling_user /= BohrR;
519 */
520
521 /*
522 input_double("MD.Opt.DIIS.Mixing",&Gdiis_Mixing,(double)0.1);
523 */
524
525 if (19<M_GDIIS_HISTORY){
526 printf("MD.Opt.DIIS.History should be lower than 19.\n");
527 MPI_Finalize();
528 exit(0);
529 }
530
531 if (MD_switch==2 || MD_switch==9 || MD_switch==11 || MD_switch==14 || MD_switch==15){
532
533 if (fp=input_find("<MD.TempControl")) {
534
535 fscanf(fp,"%i",&TempNum);
536
537 /* added by mari */
538 if (MD_switch==2 || MD_switch==11 || MD_switch==14) { /* NVT_VS or NVT_VS2 or NVT_VS4 */
539 NumScale[0] = 0;
540
541 for (i=1; i<=TempNum; i++){
542 fscanf(fp,"%d %d %lf %lf",&NumScale[i],&IntScale[i],
543 &TempScale[i],&RatScale[i]);
544 TempPara[i][1] = NumScale[i];
545 TempPara[i][2] = TempScale[i];
546 }
547
548 TempPara[0][1] = 0;
549 TempPara[0][2] = TempPara[1][2];
550 }
551
552 /* added by mari */
553 else if (MD_switch==9 || MD_switch==15) { /* NVT_NH or NVT_Langevin*/
554 for (i=1; i<=TempNum; i++){
555 fscanf(fp,"%lf %lf",&TempPara[i][1],&TempPara[i][2]);
556 }
557
558 TempPara[0][1] = 0;
559 TempPara[0][2] = TempPara[1][2];
560 }
561
562 if ( ! input_last("MD.TempControl>") ) {
563 /* format error */
564 printf("Format error for MD.TempControl\n");
565 po++;
566 }
567
568 }
569 }
570
571 if (fp=input_find("<MD.CellPressureControl")) {
572 fscanf(fp,"%i",&PreNum);
573 for (i=1; i<=TempNum; i++){
574 fscanf(fp,"%lf %lf",&PrePara[i][1],&PrePara[i][2]);
575 }
576 if ( ! input_last("MD.CellPressureControl>") ) {
577 /* format error */
578 printf("Format error for MD.CellPressureControl\n");
579 po++;
580 }
581 }
582
583 input_double("NH.Mass.HeatBath",&TempQ,(double)20.0);
584
585 input_double("Langevin.Friction.Factor",&FricFac,(double)0.001);
586
587 /****************************************************
588 solver of the eigenvalue problem
589 ****************************************************/
590
591 s_vec[0]="Recursion"; s_vec[1]="Cluster"; s_vec[2]="Band";
592 s_vec[3]="NEGF"; s_vec[4]="DC"; s_vec[5]="GDC";
593 s_vec[6]="Cluster-DIIS"; s_vec[7]="Krylov"; s_vec[8]="Cluster2";
594 s_vec[9]="EC";
595
596 i_vec[0]=1; i_vec[1]=2; i_vec[2]=3;
597 i_vec[3]=4; i_vec[4]=5; i_vec[5]=6;
598 i_vec[6]=7; i_vec[7]=8; i_vec[8]=9;
599 i_vec[9]=10;
600
601 input_string2int("scf.EigenvalueSolver", &Solver, 10, s_vec,i_vec);
602
603 if (Solver==1){
604 if (myid==Host_ID){
605 printf("The recursion method is not supported anymore.\n");
606 }
607 MPI_Finalize();
608 exit(0);
609 }
610
611 if (Solver==6){
612 if (myid==Host_ID){
613 printf("The GDC method is not supported anymore.\n");
614 }
615 MPI_Finalize();
616 exit(0);
617 }
618
619 /* default=dstevx */
620 s_vec[0]="dstevx"; s_vec[1]="dstegr"; s_vec[2]="dstedc"; s_vec[3]="dsteqr";
621 i_vec[0]=2; i_vec[1]=0; i_vec[2]=1; i_vec[3]=3;
622 input_string2int("scf.lapack.dste", &dste_flag, 4, s_vec,i_vec);
623
624 s_vec[0]="elpa1"; s_vec[1]="lapack";
625 i_vec[0]=1; i_vec[1]=0;
626 input_string2int("scf.eigen.lib", &scf_eigen_lib_flag, 3, s_vec,i_vec);
627
628 if (Solver==1){
629 if (myid==Host_ID){
630 printf("Recursion method is not supported in this version.\n");
631 printf("Check your input file.\n\n");
632 }
633 MPI_Finalize();
634 exit(1);
635 }
636
637 /****************************************************
638 for generation of Monkhorst-Pack k-points
639 ****************************************************/
640
641 input_double("scf.MP.criterion",&Criterion_MP_Special_Kpt,(double)1.0e-5);
642
643 s_vec[0]="REGULAR"; s_vec[1]="MP";
644 i_vec[0]=1 ; i_vec[1]=2 ;
645 input_string2int("scf.Generation.Kpoint", &way_of_kpoint, 2, s_vec,i_vec);
646
647 if (Solver==4 && way_of_kpoint==2){
648 if (myid==Host_ID){
649 printf("The Monkhorst-Pack is not supported for NEGF.\n");
650 printf("Check your input file.\n\n");
651 }
652 MPI_Finalize();
653 exit(1);
654 }
655
656 /****************************************************
657 flags for saving and reading Fourier transformed
658 quantities generated in FT_*.c
659 ****************************************************/
660
661 input_logical("FT.files.save",&FT_files_save,0);
662 input_logical("FT.files.read",&FT_files_read,0);
663
664 /****************************************************
665 SCF or electronic system
666 ****************************************************/
667
668 s_vec[0]="LDA"; s_vec[1]="LSDA-CA"; s_vec[2]="LSDA-PW"; s_vec[3]="GGA-PBE";s_vec[4]="EXX-TEST";
669 i_vec[0]=1; i_vec[1]=2; i_vec[2]= 3; i_vec[3]= 4; i_vec[4]=5;
670 input_string2int("scf.XcType", &XC_switch, 5, s_vec,i_vec);
671
672 s_vec[0]="Off"; s_vec[1]="On"; s_vec[2]="NC";
673 i_vec[0]=0 ; i_vec[1]=1 ; i_vec[2]=3;
674 input_string2int("scf.SpinPolarization", &SpinP_switch, 3, s_vec,i_vec);
675 if (SpinP_switch==0) List_YOUSO[23] = 1;
676 else if (SpinP_switch==1) List_YOUSO[23] = 2;
677 else if (SpinP_switch==3) List_YOUSO[23] = 4;
678
679 if (XC_switch==3 && Solver==8){
680 if (myid==Host_ID){
681 printf("Krylov subspace method is not supported for non-collinear calculations.\n");
682 }
683 MPI_Finalize();
684 exit(1);
685 }
686
687 if (XC_switch==1 && 1<=SpinP_switch){
688 if (myid==Host_ID){
689 printf("SpinP_switch should be OFF for this exchange functional.\n");
690 }
691 MPI_Finalize();
692 exit(1);
693 }
694
695 /* scf.Constraint.NC.Spin */
696
697 s_vec[0]="off"; s_vec[1]="on"; s_vec[2]="on2";
698 i_vec[0]=0 ; i_vec[1]=1 ; i_vec[2]=2;
699 input_string2int("scf.Constraint.NC.Spin", &Constraint_NCS_switch, 3, s_vec,i_vec);
700
701 if (SpinP_switch!=3 && 1<=Constraint_NCS_switch){
702 if (myid==Host_ID){
703 printf("The constraint scheme is not supported for a collinear DFT calculation.\n");
704 printf("Check your input file.\n\n");
705 }
706 MPI_Finalize();
707 exit(1);
708 }
709
710 if (1<=Constraint_NCS_switch && Hub_U_occupation!=2){
711 if (myid==Host_ID){
712 printf("The constraint scheme is supported in case of scf.Hubbard.Occupation=dual.\n");
713 printf("Check your input file.\n\n");
714 }
715 MPI_Finalize();
716 exit(1);
717 }
718
719 input_double("scf.Constraint.NC.Spin.V",&Constraint_NCS_V,(double)0.0); /* in eV */
720 /* eV to Hartree */
721 Constraint_NCS_V = Constraint_NCS_V/eV2Hartree;
722
723 /* scf.SpinOrbit.Coupling */
724
725 input_logical("scf.SpinOrbit.Coupling",&SO_switch,0);
726
727 if (SpinP_switch!=3 && SO_switch==1){
728 if (myid==Host_ID){
729 printf("Spin-orbit coupling is not supported for collinear DFT calculations.\n");
730 printf("Check your input file.\n\n");
731 }
732 MPI_Finalize();
733 exit(1);
734 }
735
736 if (SpinP_switch==0 && SO_switch==1){
737 if (myid==Host_ID){
738 printf("scf.SpinOrbit.Coupling must be OFF when scf.SpinPolarization=OFF\n");
739 }
740 MPI_Finalize();
741 exit(1);
742 }
743
744 /* scf.SO.factor */
745
746 SO_factor_flag = 0;
747
748 if (SpinP_switch==3 && SO_switch==1){
749
750 if (fp=input_find("<scf.SO.factor")) {
751
752 /* switch on SO_factor_flag */
753
754 SO_factor_flag = 1;
755
756 /* initialize the SO_factor */
757 for (i=0; i<SpeciesNum; i++){
758 for (l=0; l<=3; l++){
759 SO_factor[i][l] = 1.0 ;
760 }
761 }
762
763 /* read SO_factor from the '.dat' file */
764 for (i=0; i<SpeciesNum; i++){
765 fscanf(fp,"%s",Species);
766 j = Species2int(Species);
767
768 for (l=0; l<=3; l++){
769 fscanf(fp,"%s %lf", buf, &SO_factor[j][l]) ;
770 }
771 }
772
773 if (! input_last("scf.SO.factor>") ) {
774 /* format error */
775 printf("Format error for scf.SO.factor\n");
776 po++;
777 }
778
779 } /* if (fp=input_find("<scf.SO.factor")) */
780
781 } /* if (SpinP_switch==3 && SO_switch==1) */
782
783 /* scf.partialCoreCorrection */
784
785 input_logical("scf.partialCoreCorrection",&PCC_switch,1);
786
787 if (PCC_switch==0){
788 if (myid==Host_ID){
789 printf("scf.partialCoreCorrection should be always switched on.\n");
790 printf("Check your input file.\n\n");
791 }
792 MPI_Finalize();
793 exit(1);
794 }
795
796 /* scf.pcc.opencore */
797
798 if (fp=input_find("<scf.pcc.opencore")) {
799
800 for (i=0; i<SpeciesNum; i++){
801 fscanf(fp,"%s",Species);
802 j = Species2int(Species);
803 fscanf(fp,"%d", &Spe_OpenCore_flag[j]) ;
804
805 if (Spe_OpenCore_flag[j]==0 || Spe_OpenCore_flag[j]==1 || Spe_OpenCore_flag[j]==-1){
806 }
807 else{
808 /* format error */
809 printf("Valid input for scf.pcc.opencore is 0, -1, or 1.\n");
810 po++;
811 }
812 }
813
814 if (! input_last("scf.pcc.opencore>") ) {
815 /* format error */
816 printf("Format error for scf.pcc.opencore\n");
817 po++;
818 }
819
820 } /* if (fp=input_find("<scf.pcc.opencore")) */
821
822 /* scf.NC.Zeeman.Spin */
823
824 input_logical("scf.NC.Zeeman.Spin",&Zeeman_NCS_switch,0);
825
826 if (SpinP_switch!=3 && Zeeman_NCS_switch==1){
827 if (myid==Host_ID){
828 printf("The Zeeman term is not supported for a collinear DFT calculation.\n");
829 printf("Check your input file.\n\n");
830 }
831 MPI_Finalize();
832 exit(1);
833 }
834
835 if (Zeeman_NCS_switch==1 && Hub_U_occupation!=2){
836 if (myid==Host_ID){
837 printf("The Zeeman term for spin is supported in case of scf.Hubbard.Occupation=dual.\n");
838 printf("Check your input file.\n\n");
839 }
840 MPI_Finalize();
841 exit(1);
842 }
843
844 if (1<=Constraint_NCS_switch && Zeeman_NCS_switch==1){
845 if (myid==Host_ID){
846 printf("For spin magnetic moment, the constraint scheme and the Zeeman term\n");
847 printf("are mutually exclusive. Check your input file.\n\n");
848 }
849 MPI_Finalize();
850 exit(1);
851 }
852
853 input_double("scf.NC.Mag.Field.Spin",&Mag_Field_Spin,(double)0.0);
854
855 /**************************************
856 Tesla to a.u.
857 1 Tesla = 1/(2.35051742*10^5) a.u.
858 ***************************************/
859
860 Mag_Field_Spin = Mag_Field_Spin/(2.35051742*100000.0);
861
862 /* scf.NC.Zeeman.Orbital */
863
864 input_logical("scf.NC.Zeeman.Orbital",&Zeeman_NCO_switch,0);
865
866 if (SpinP_switch!=3 && Zeeman_NCO_switch==1){
867 if (myid==Host_ID){
868 printf("The Zeeman term is not supported for a collinear DFT calculation.\n");
869 printf("Check your input file.\n\n");
870 }
871 MPI_Finalize();
872 exit(1);
873 }
874
875 if (Zeeman_NCO_switch==1 && Hub_U_occupation!=2){
876 if (myid==Host_ID){
877 printf("The Zeeman term for orbital is supported in case of scf.Hubbard.Occupation=dual.\n");
878 printf("Check your input file.\n\n");
879 }
880 MPI_Finalize();
881 exit(1);
882 }
883
884 if (Zeeman_NCO_switch==1 && SO_switch==0){
885 if (myid==Host_ID){
886 printf("The Zeeman term for orbital moment is not supported without the SO term.\n");
887 printf("Check your input file.\n\n");
888 }
889 MPI_Finalize();
890 exit(1);
891 }
892
893 input_double("scf.NC.Mag.Field.Orbital",&Mag_Field_Orbital,(double)0.0);
894
895 /**************************************
896 Tesla to a.u.
897 1 Tesla = 1/(2.35051742*10^5) a.u.
898 ***************************************/
899
900 Mag_Field_Orbital = Mag_Field_Orbital/(2.35051742*100000.0);
901
902 if (SpinP_switch==0) List_YOUSO[5] = 1;
903 else if (SpinP_switch==1) List_YOUSO[5] = 2;
904 else if (SpinP_switch==3 && SO_switch==0) List_YOUSO[5] = 3;
905 else if (SpinP_switch==3 && SO_switch==1) List_YOUSO[5] = 3;
906
907 i_vec2[0]=4;
908 i_vec2[1]=4;
909 i_vec2[2]=4;
910 input_intv("scf.Kgrid",3,i_vec,i_vec2);
911 Kspace_grid1 = i_vec[0];
912 Kspace_grid2 = i_vec[1];
913 Kspace_grid3 = i_vec[2];
914
915 if (Kspace_grid1<=0){
916 printf("Kspace_grid1 should be over 1\n");
917 MPI_Finalize();
918 exit(1);
919 }
920 if (Kspace_grid2<=0){
921 printf("Kspace_grid2 should be over 1\n");
922 MPI_Finalize();
923 exit(1);
924 }
925 if (Kspace_grid3<=0){
926 printf("Kspace_grid3 should be over 1\n");
927 MPI_Finalize();
928 exit(1);
929 }
930
931 if (Solver!=3 && Solver!=4){
932 List_YOUSO[27] = 1;
933 List_YOUSO[28] = 1;
934 List_YOUSO[29] = 1;
935 }
936 else{
937 List_YOUSO[27] = Kspace_grid1;
938 List_YOUSO[28] = Kspace_grid2;
939 List_YOUSO[29] = Kspace_grid3;
940 }
941
942 /* set PeriodicGamma_flag in 1 in the band calc. with only the gamma point */
943 PeriodicGamma_flag = 0;
944 if (Solver==3 && Kspace_grid1==1 && Kspace_grid2==1 && Kspace_grid3==1){
945
946 if (myid==Host_ID){
947 printf("When only the gamma point is considered, the eigenvalue solver is changed to 'Cluster' with the periodic boundary condition.\n");fflush(stdout);
948 }
949 PeriodicGamma_flag = 1;
950 Solver = 2;
951 }
952
953 input_double("scf.ElectronicTemperature",&E_Temp,(double)300.0);
954 E_Temp = E_Temp/eV2Hartree;
955 Original_E_Temp = E_Temp;
956
957 s_vec[0]="Simple"; i_vec[0]=0;
958 s_vec[1]="RMM-DIIS"; i_vec[1]=1;
959 s_vec[2]="GR-Pulay"; i_vec[2]=2;
960 s_vec[3]="Kerker"; i_vec[3]=3;
961 s_vec[4]="RMM-DIISK"; i_vec[4]=4;
962 s_vec[5]="RMM-DIISH"; i_vec[5]=5;
963 s_vec[6]="RMM-ADIIS"; i_vec[6]=6;
964
965 input_string2int("scf.Mixing.Type",&Mixing_switch,7,s_vec,i_vec);
966 if (Mixing_switch==3 || Mixing_switch==4) Kmixing_flag = 1;
967 else Kmixing_flag = 0;
968
969 input_double("scf.Init.Mixing.Weight",&Mixing_weight,(double)0.3);
970 input_double("scf.Min.Mixing.Weight",&Min_Mixing_weight,(double)0.001);
971 input_double("scf.Max.Mixing.Weight",&Max_Mixing_weight,(double)0.4);
972 /* if Kerker_factor is not set here, later Kerker factor is automatically determined. */
973 input_double("scf.Kerker.factor",&Kerker_factor,(double)-1.0);
974 input_int("scf.Mixing.History",&Num_Mixing_pDM,5);
975 input_int("scf.Mixing.StartPulay",&Pulay_SCF,6);
976 Pulay_SCF_original = Pulay_SCF;
977 input_int("scf.Mixing.EveryPulay",&EveryPulay_SCF,1);
978 input_int("scf.ExtCharge.History",&Extrapolated_Charge_History,3);
979
980 if (Pulay_SCF<4 && Mixing_switch==1){
981 if (myid==Host_ID){
982 printf("For RMM-DIIS, scf.Mixing.StartPulay should be set to larger than 3.\n");
983 }
984 MPI_Finalize();
985 exit(1);
986 }
987
988 /* increase electric temperature in case of SCF oscillation, default=off */
989 input_logical("scf.Mixing.Control.Temp", &SCF_Control_Temp, 0);
990
991 if (Mixing_switch==0){
992 List_YOUSO[16] = 3;
993 List_YOUSO[38] = 1;
994 }
995 else if (Mixing_switch==1 || Mixing_switch==2 || Mixing_switch==6) {
996 List_YOUSO[16] = Num_Mixing_pDM + 2;
997 List_YOUSO[38] = 1;
998 List_YOUSO[39] = Num_Mixing_pDM + 3;
999 }
1000 else if (Mixing_switch==3) {
1001 List_YOUSO[16] = 3;
1002 List_YOUSO[38] = 3;
1003 }
1004 else if (Mixing_switch==4) {
1005 List_YOUSO[16] = 3;
1006 List_YOUSO[38] = Num_Mixing_pDM + 2;
1007 }
1008 else if (Mixing_switch==5){
1009 List_YOUSO[16] = 3;
1010 List_YOUSO[38] = 1;
1011 List_YOUSO[39] = Num_Mixing_pDM + 3;
1012 }
1013
1014 input_double("scf.criterion",&SCF_Criterion,(double)1.0e-6);
1015 if (SCF_Criterion<0.0){
1016 printf("SCF_Criterion=%10.9f should be larger than 0.\n",SCF_Criterion);
1017 po++;
1018 }
1019
1020 input_double("scf.system.charge",&system_charge,(double)0.0);
1021
1022 /* scf.fixed.grid */
1023
1024 r_vec[0]=1.0e+9; r_vec[1]=1.0e+9; r_vec[2]=1.0e+9;
1025 input_doublev("scf.fixed.grid",3,scf_fixed_origin,r_vec);
1026
1027 TRAN_Input_std(MPI_COMM_WORLD1, Solver, SpinP_switch, filepath, kB,
1028 eV2Hartree, E_Temp, &output_hks);
1029
1030 /********************************************************
1031 Effective Screening Medium (ESM) method Calculation
1032 added by T.Ohwaki
1033 *********************************************************/
1034
1035 s_vec[0]="off"; s_vec[1]="on1"; s_vec[2]="on2"; s_vec[3]="on3"; s_vec[4]="on4";
1036 i_vec[0]=0 ; i_vec[1]=1 ; i_vec[2]=2 ; i_vec[3]=3 ; i_vec[4]=4 ;
1037 input_string2int("ESM.switch", &ESM_switch, 5, s_vec,i_vec);
1038
1039 s_vec[0]="off"; s_vec[1]="on";
1040 i_vec[0]=0 ; i_vec[1]=1 ;
1041 input_string2int("ESM.wall.switch", &ESM_wall_switch, 2, s_vec,i_vec);
1042
1043 if (myid==Host_ID && 0<level_stdout){
1044
1045 if (ESM_switch==1){
1046 printf("\n");
1047 printf("********************************************************** \n");
1048 printf(" Effective Screening Medium (ESM) method calculation \n");
1049 printf(" \n");
1050 printf(" The following calc. is implemented with ESM method. \n");
1051 printf(" Boundary condition = Vacuum|Vacuum|Vacuum \n");
1052 printf(" \n");
1053 printf(" Copyright (C), 2011-2013, T.Ohwaki and M.Otani \n");
1054 printf("********************************************************** \n");
1055 printf("\n");
1056 }
1057
1058 else if (ESM_switch==2){
1059 printf("\n");
1060 printf("********************************************************** \n");
1061 printf(" Effective Screening Medium (ESM) method calculation \n");
1062 printf(" \n");
1063 printf(" The following calc. is implemented with ESM method. \n");
1064 printf(" Boundary condition = Metal|Vacuum|Metal \n");
1065 printf(" \n");
1066 printf(" Copyright (C), 2011-2013, T.Ohwaki and M.Otani \n");
1067 printf("********************************************************** \n");
1068 printf("\n");
1069 }
1070
1071 else if (ESM_switch==3){
1072 printf("\n");
1073 printf("********************************************************** \n");
1074 printf(" Effective Screening Medium (ESM) method calculation \n");
1075 printf(" \n");
1076 printf(" The following calc. is implemented with ESM method. \n");
1077 printf(" Boundary condition = Vacuum|Vacuum|Metal \n");
1078 printf(" \n");
1079 printf(" Copyright (C), 2011-2013, T.Ohwaki and M.Otani \n");
1080 printf("********************************************************** \n");
1081 printf("\n");
1082 }
1083
1084 else if (ESM_switch==4){
1085 printf("\n");
1086 printf("********************************************************** \n");
1087 printf(" Effective Screening Medium (ESM) method calculation \n");
1088 printf(" \n");
1089 printf(" The following calc. is implemented with ESM method. \n");
1090 printf(" Boundary condition = Metal|Vacuum|Metal \n");
1091 printf(" plus Uniform electric field \n");
1092 printf(" \n");
1093 printf(" Copyright (C), 2011-2013, T.Ohwaki and M.Otani \n");
1094 printf("********************************************************** \n");
1095 printf("\n");
1096 }
1097
1098 }
1099
1100 input_double("ESM.potential.diff",&V_ESM,(double)0.0); /* eV */
1101 /* change the unit from eV to Hartree */
1102 V_ESM = V_ESM/eV2Hartree;
1103
1104 if (ESM_switch!=0 && ESM_switch!=4 && 1.0e-13<fabs(V_ESM)){
1105 if (myid==Host_ID){
1106 printf("<ESM:Warning> Non-zero ESM.ElectricField is not valid except for on4.\n\n");
1107 }
1108
1109 MPI_Finalize();
1110 exit(1);
1111 }
1112
1113 /* This means the distance from the upper edge along the a-axis (x-coordinate). */
1114 input_double("ESM.wall.position",&ESM_wall_position,(double)1.0); /* Angstrom */
1115 /* change the unit from ang. to a.u. */
1116 ESM_wall_position /= BohrR;
1117
1118 if (ESM_wall_switch==1){
1119
1120 if (ESM_wall_position<0.0){
1121 if (myid==Host_ID){
1122 printf("<ESM:Warning> ESM.wall.position must be positive.\n\n");
1123 }
1124
1125 MPI_Finalize();
1126 exit(1);
1127 }
1128
1129 input_double("ESM.wall.height",&ESM_wall_height,(double)100.0); /* eV */
1130 /* change the unit from eV to Hartree */
1131 ESM_wall_height /= eV2Hartree;
1132
1133 input_double("ESM.buffer.range",&ESM_buffer_range,(double)10.0); /* Angstrom */
1134 /* change the unit from ang. to a.u. */
1135 ESM_buffer_range /= BohrR;
1136
1137 } /* ESM_wall_switch */
1138
1139 /* Artificially imposed force */
1140 /* added by T.Ohwaki */
1141
1142 s_vec[0]="off"; s_vec[1]="on";
1143 i_vec[0]=0 ; i_vec[1]=1 ;
1144 Arti_Force = 0; /* default = off */
1145 input_string2int("MD.Artificial_Force", &Arti_Force, 2, s_vec,i_vec);
1146 input_double("MD.Artificial_Grad",&Arti_Grad,(double)0.0); /* Hartree/Bohr */
1147
1148 if (myid==Host_ID && Arti_Force==1 && 0<level_stdout){
1149 printf("\n");
1150 printf("##################################################### \n");
1151 printf("## ## \n");
1152 printf("## An artificial force is imposed on the 1st atom. ## \n");
1153 printf("## * Gradient = %12.9f (Hartree/Bohr) \n",Arti_Grad);
1154 printf("## ## \n");
1155 printf("##################################################### \n");
1156 printf("\n");
1157 }
1158
1159 /*****************************************************
1160 if restart files for geometry optimization exist,
1161 read data from them.
1162 default = off
1163 *****************************************************/
1164
1165 input_logical("geoopt.restart",&GeoOpt_RestartFromFile, 0);
1166
1167 /****************************************************
1168 atoms
1169 ****************************************************/
1170
1171 /* except for NEGF */
1172
1173 if (Solver!=4){
1174
1175 /* atom */
1176
1177 input_int("Atoms.Number",&atomnum,0);
1178
1179 if (atomnum<=0){
1180 printf("Atoms.Number may be wrong.\n");
1181 po++;
1182 }
1183 List_YOUSO[1] = atomnum + 1;
1184
1185 /* memory allocation */
1186 Allocate_Arrays(1);
1187
1188 /* initialize */
1189
1190 s_vec[0]="Ang"; s_vec[1]="AU"; s_vec[2]="FRAC";
1191 i_vec[0]= 0; i_vec[1]= 1; i_vec[2]= 2;
1192 input_string2int("Atoms.SpeciesAndCoordinates.Unit",
1193 &coordinates_unit,3,s_vec,i_vec);
1194
1195 if (fp=input_find("<Atoms.SpeciesAndCoordinates") ) {
1196
1197 for (i=1; i<=atomnum; i++){
1198 fgets(buf,MAXBUF,fp);
1199
1200 /* spin non-collinear */
1201 if (SpinP_switch==3){
1202
1203 /*******************************************************
1204 (1) spin non-collinear
1205 *******************************************************/
1206
1207 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %s",
1208 &j, Species,
1209 &Gxyz[i][1],&Gxyz[i][2],&Gxyz[i][3],
1210 &InitN_USpin[i],&InitN_DSpin[i],
1211 &Angle0_Spin[i], &Angle1_Spin[i],
1212 &Angle0_Orbital[i], &Angle1_Orbital[i],
1213 &Constraint_SpinAngle[i],
1214 OrbPol );
1215
1216 /* atomMove is initialzed as 1 */
1217
1218 if (fabs(Angle0_Spin[i])<1.0e-10){
1219 Angle0_Spin[i] = Angle0_Spin[i] + rnd(1.0e-5);
1220 }
1221
1222 Angle0_Spin[i] = PI*Angle0_Spin[i]/180.0;
1223 Angle1_Spin[i] = PI*Angle1_Spin[i]/180.0;
1224 InitAngle0_Spin[i] = Angle0_Spin[i];
1225 InitAngle1_Spin[i] = Angle1_Spin[i];
1226
1227 if (fabs(Angle0_Orbital[i])<1.0e-10){
1228 Angle0_Orbital[i] = Angle0_Orbital[i] + rnd(1.0e-5);
1229 }
1230
1231 Angle0_Orbital[i] = PI*Angle0_Orbital[i]/180.0;
1232 Angle1_Orbital[i] = PI*Angle1_Orbital[i]/180.0;
1233 InitAngle0_Orbital[i] = Angle0_Orbital[i];
1234 InitAngle1_Orbital[i] = Angle1_Orbital[i];
1235
1236
1237 /*************************************************************************
1238 check whether the Euler angle measured from the direction (1,0) is used,
1239 if not, change the Euler angle
1240 *************************************************************************/
1241
1242 if ( (InitN_USpin[i]-InitN_DSpin[i])<0.0 ){
1243
1244 mx = -sin(InitAngle0_Spin[i])*cos(InitAngle1_Spin[i]);
1245 my = -sin(InitAngle0_Spin[i])*sin(InitAngle1_Spin[i]);
1246 mz = -cos(InitAngle0_Spin[i]);
1247
1248 xyz2spherical(mx,my,mz,0.0,0.0,0.0,S_coordinate);
1249
1250 Angle0_Spin[i] = S_coordinate[1];
1251 Angle1_Spin[i] = S_coordinate[2];
1252 InitAngle0_Spin[i] = Angle0_Spin[i];
1253 InitAngle1_Spin[i] = Angle1_Spin[i];
1254
1255 tmp = InitN_USpin[i];
1256 InitN_USpin[i] = InitN_DSpin[i];
1257 InitN_DSpin[i] = tmp;
1258
1259 mx = -sin(InitAngle0_Orbital[i])*cos(InitAngle1_Orbital[i]);
1260 my = -sin(InitAngle0_Orbital[i])*sin(InitAngle1_Orbital[i]);
1261 mz = -cos(InitAngle0_Orbital[i]);
1262
1263 xyz2spherical(mx,my,mz,0.0,0.0,0.0,S_coordinate);
1264
1265 Angle0_Orbital[i] = S_coordinate[1];
1266 Angle1_Orbital[i] = S_coordinate[2];
1267 InitAngle0_Orbital[i] = Angle0_Orbital[i];
1268 InitAngle1_Orbital[i] = Angle1_Orbital[i];
1269 }
1270 }
1271
1272 /**************************************************
1273 (2) spin collinear
1274 **************************************************/
1275
1276 else{
1277
1278 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %s",
1279 &j, Species,
1280 &Gxyz[i][1],&Gxyz[i][2],&Gxyz[i][3],
1281 &InitN_USpin[i],&InitN_DSpin[i], OrbPol );
1282 }
1283
1284 WhatSpecies[i] = Species2int(Species);
1285
1286 if (Hub_U_switch==1) OrbPol_flag[i] = OrbPol2int(OrbPol);
1287
1288 if (i!=j){
1289 printf("Format error of the sequential number %i in <Atoms.SpeciesAndCoordinates\n",j);
1290 po++;
1291 }
1292
1293 if (2<=level_stdout){
1294 printf("<Input_std> ct_AN=%2d WhatSpecies=%2d USpin=%7.4f DSpin=%7.4f\n",
1295 i,WhatSpecies[i],InitN_USpin[i],InitN_DSpin[i]);
1296 }
1297 }
1298
1299 ungetc('\n',fp);
1300 if (!input_last("Atoms.SpeciesAndCoordinates>")) {
1301 /* format error */
1302 printf("Format error for Atoms.SpeciesAndCoordinates\n");
1303 po++;
1304 }
1305
1306 }
1307
1308 /****************************************************
1309 the unit of atomic coordinates is transformed
1310 ****************************************************/
1311
1312 /* Ang to AU */
1313 if (coordinates_unit==0){
1314 for (i=1; i<=atomnum; i++){
1315 Gxyz[i][1] = Gxyz[i][1]/BohrR;
1316 Gxyz[i][2] = Gxyz[i][2]/BohrR;
1317 Gxyz[i][3] = Gxyz[i][3]/BohrR;
1318 }
1319 }
1320
1321 /****************************************************
1322 unit cell
1323 ****************************************************/
1324
1325 s_vec[0]="Ang"; s_vec[1]="AU";
1326 i_vec[0]=0; i_vec[1]=1;
1327 input_string2int("Atoms.UnitVectors.Unit",&unitvector_unit,2,s_vec,i_vec);
1328
1329 if (fp=input_find("<Atoms.Unitvectors")) {
1330
1331 for (i=1; i<=3; i++){
1332 fscanf(fp,"%lf %lf %lf",&tv[i][1],&tv[i][2],&tv[i][3]);
1333 }
1334 if ( ! input_last("Atoms.Unitvectors>") ) {
1335 /* format error */
1336 printf("Format error for Atoms.Unitvectors\n");
1337 po++;
1338 }
1339
1340 /* Ang to AU */
1341 if (unitvector_unit==0){
1342 for (i=1; i<=3; i++){
1343 tv[i][1] = tv[i][1]/BohrR;
1344 tv[i][2] = tv[i][2]/BohrR;
1345 tv[i][3] = tv[i][3]/BohrR;
1346 }
1347 }
1348
1349 /* check the cell vectors for the band calculation of the lead */
1350
1351 if (output_hks==1){
1352
1353 if ( tv[1][1]<0.0 || tv[1][2]<0.0 || tv[1][3]<0.0 ||
1354 1.0e-12<fabs(tv[1][1]*tv[1][2]) ||
1355 1.0e-12<fabs(tv[1][2]*tv[1][3]) ||
1356 1.0e-12<fabs(tv[1][3]*tv[1][1])
1357 ){
1358
1359 if (myid==Host_ID){
1360 printf("The a-axis of the unit cell of leads used for the NEGF calculation\n");
1361 printf("must be x, y, or, z-axis in the cartesian coordinate.\n");
1362 printf("The following a-axis will be accepted:\n\n");
1363
1364 printf(" ax 0.0 0.0\n");
1365 printf(" or\n");
1366 printf(" 0.0 ay 0.0\n");
1367 printf(" or\n");
1368 printf(" 0.0 0.0 az\n\n");
1369 printf("where ax, ay, and az must be POSITIVE.\n\n");
1370 }
1371
1372 MPI_Finalize();
1373 exit(1);
1374 }
1375 }
1376
1377 /* Effective Screening Medium (ESM) method Calculation */
1378 /* added by T.Ohwaki */
1379
1380 if (ESM_switch!=0){
1381
1382 if ( tv[1][1]<0.0 || 1.0e-13<fabs(tv[1][2]) || 1.0e-13<fabs(tv[1][3]) ||
1383 1.0e-13<fabs(tv[1][1]*tv[1][2]) ||
1384 1.0e-13<fabs(tv[1][2]*tv[1][3]) ||
1385 1.0e-13<fabs(tv[1][3]*tv[1][1])
1386 ){
1387
1388 if (myid==Host_ID){
1389 printf("<ESM:Warning>\n");
1390 printf("The a-axis of the unit cell used for the ESM calculation\n");
1391 printf("must be (x, 0, 0) in the cartesian coordinate,\n");
1392 printf("where x must be POSITIVE.\n\n");
1393 }
1394
1395 MPI_Finalize();
1396 exit(1);
1397 }
1398
1399 if ( 1.0e-13<fabs(tv[2][1]) || 1.0e-13<fabs(tv[3][1]) ){
1400
1401 if (myid==Host_ID){
1402 printf("<ESM:Warning>\n");
1403 printf("The b- and c-axes of the unit cell used for the ESM calculation\n");
1404 printf("must be orthogonal to the a-axis.\n\n");
1405 }
1406
1407 MPI_Finalize();
1408 exit(1);
1409 }
1410 }
1411
1412 }
1413
1414 else {
1415 /* automatically set up the unit cell */
1416
1417 /* Effective Screening Medium (ESM) method Calculation */
1418 /* added by T.Ohwaki */
1419
1420 if (ESM_switch!=0){
1421 if (myid==Host_ID){
1422 printf("<ESM:Warning> A unit cell must be specified in the ESM calculation.\n\n");
1423 }
1424
1425 MPI_Finalize();
1426 exit(1);
1427 }
1428
1429 /* In case of Solver==3 */
1430
1431 if (Solver==3){
1432 if (myid==Host_ID){
1433 printf("You have to give a unit cell for the band calc.\n");
1434 }
1435 MPI_Finalize();
1436 exit(1);
1437 }
1438
1439 /* The other case */
1440
1441 if (coordinates_unit==2){
1442 if (myid==Host_ID){
1443 printf("You have to give a unit cell in case of use of fractional coordinates.\n");
1444 }
1445 MPI_Finalize();
1446 exit(1);
1447 }
1448
1449 /* Set a unit cell */
1450
1451 Set_Cluster_UnitCell(tv,unitvector_unit);
1452 Determine_Cell_from_ECutoff(tv,Grid_Ecut+0.001);
1453 }
1454
1455 /* FRAC to AU */
1456 if (coordinates_unit==2){
1457
1458 /* The fractional coordinates should be kept within 0 to 1. */
1459
1460 for (i=1; i<=atomnum; i++){
1461 for (k=1; k<=3; k++){
1462
1463 itmp = (int)Gxyz[i][k];
1464
1465 if (1.0<Gxyz[i][k]){
1466
1467 /* Effective Screening Medium (ESM) method Calculation */
1468 /* added by T.Ohwaki */
1469
1470 if (ESM_switch!=0 && k==1){
1471 if (myid==Host_ID){
1472 printf("<ESM:Warning>\n");
1473 printf("The fractional coordinate of a-axis for atom %d = %16.9f \n",i,Gxyz[i][k]);
1474 printf("The fractional coordinate of a-axis should be kept within 0 to 1 in the ESM calculation.\n\n");
1475 }
1476
1477 MPI_Finalize();
1478 exit(1);
1479 }
1480
1481 /* ended by T.Ohwaki */
1482
1483 if (GeoOpt_RestartFromFile==0){
1484
1485 Gxyz[i][k] = Gxyz[i][k] - (double)itmp;
1486
1487 if (myid==Host_ID){
1488 if (k==1) printf("The fractional coordinate of a-axis for atom %d was translated within the range (0 to 1).\n",i);
1489 if (k==2) printf("The fractional coordinate of b-axis for atom %d was translated within the range (0 to 1).\n",i);
1490 if (k==3) printf("The fractional coordinate of c-axis for atom %d was translated within the range (0 to 1).\n",i);
1491 }
1492 }
1493 }
1494 else if (Gxyz[i][k]<0.0){
1495
1496 /* Effective Screening Medium (ESM) method Calculation */
1497 /* added by T.Ohwaki */
1498
1499 if (ESM_switch!=0 && k==1){
1500 if (myid==Host_ID){
1501 printf("<ESM:Warning>\n");
1502 printf("The fractional coordinate of a-axis for atom %d = %16.9f \n",i,Gxyz[i][k]);
1503 printf("The fractional coordinate of a-axis should be kept within 0 to 1 in the ESM calculation.\n\n");
1504 }
1505
1506 MPI_Finalize();
1507 exit(1);
1508 }
1509
1510 /* ended by T.Ohwaki */
1511
1512 if (GeoOpt_RestartFromFile==0){
1513
1514 Gxyz[i][k] = Gxyz[i][k] + (double)(abs(itmp)+1);
1515
1516 if (myid==Host_ID){
1517 if (k==1) printf("The fractional coordinate of a-axis for atom %d was translated within the range (0 to 1).\n",i);
1518 if (k==2) printf("The fractional coordinate of b-axis for atom %d was translated within the range (0 to 1).\n",i);
1519 if (k==3) printf("The fractional coordinate of c-axis for atom %d was translated within the range (0 to 1).\n",i);
1520 }
1521 }
1522 }
1523
1524 }
1525 }
1526
1527 /* calculation of xyz-coordinate in A.U. The grid origin is Grid_Origin or zero. */
1528
1529 if ( 1.0e+8<scf_fixed_origin[0] &&
1530 1.0e+8<scf_fixed_origin[1] &&
1531 1.0e+8<scf_fixed_origin[2] ){
1532
1533 scf_fixed_origin[0] = 0.0;
1534 scf_fixed_origin[1] = 0.0;
1535 scf_fixed_origin[2] = 0.0;
1536
1537 tmpx = 0.0;
1538 tmpy = 0.0;
1539 tmpz = 0.0;
1540 }
1541 else {
1542 tmpx = scf_fixed_origin[0];
1543 tmpy = scf_fixed_origin[1];
1544 tmpz = scf_fixed_origin[2];
1545 }
1546
1547 for (i=1; i<=atomnum; i++){
1548 x = Gxyz[i][1]*tv[1][1] + Gxyz[i][2]*tv[2][1] + Gxyz[i][3]*tv[3][1] + tmpx;
1549 y = Gxyz[i][1]*tv[1][2] + Gxyz[i][2]*tv[2][2] + Gxyz[i][3]*tv[3][2] + tmpy;
1550 z = Gxyz[i][1]*tv[1][3] + Gxyz[i][2]*tv[2][3] + Gxyz[i][3]*tv[3][3] + tmpz;
1551 Gxyz[i][1] = x;
1552 Gxyz[i][2] = y;
1553 Gxyz[i][3] = z;
1554 }
1555 }
1556
1557 } /* if (Solver!=4){ */
1558
1559 /*******************************
1560 NEGF
1561 *******************************/
1562
1563 else{
1564
1565 TRAN_Input_std_Atoms(MPI_COMM_WORLD1, Solver);
1566 }
1567
1568 /**************************************************
1569 store the initial magnetic moment
1570 **************************************************/
1571
1572 if (Constraint_NCS_switch==2){
1573 for (i=1; i<=atomnum; i++){
1574 InitMagneticMoment[i] = fabs(InitN_USpin[i] - InitN_DSpin[i]);
1575 }
1576 }
1577
1578 /********************************************************
1579 Automatic determination of Kerker factor
1580
1581 Kerker_factor = 1/G0 inplies
1582 Kerker_factor * G0 -> constant
1583
1584 Furthermore, to take account of anisotropy of system,
1585 an enhancement factor is introduced by
1586
1587 y = a*(DG/AG) + 1
1588
1589 where,
1590 AG = (G12+G22+G32)/3;
1591 DG = (fabs(G22-G12)+fabs(G32-G22)+fabs(G12-G32))/3;
1592
1593 In an extreme case that G22<<G12 and G32<<G12,
1594 AG approaches to G12/3, and DG approaches to 2/3 G12.
1595 So, DG/AG becomes 2 in the case.
1596
1597 On the other hand, if G12=G22=G32, y becomes 1.
1598
1599 Taking appropriate prefactors, we define as
1600 Kerker_factor = 0.5/G0*y = 0.5/G0*(4.0*DG/AG+1.0);
1601 *******************************************************/
1602
1603 if ( Kerker_factor<0.0 && (Mixing_switch==3 || Mixing_switch==4) ){
1604
1605 double tmp[4];
1606 double CellV,G0,G12,G22,G32,DG,AG;
1607
1608 Cross_Product(tv[2],tv[3],tmp);
1609 CellV = Dot_Product(tv[1],tmp);
1610 Cell_Volume = fabs(CellV);
1611
1612 Cross_Product(tv[2],tv[3],tmp);
1613 rtv[1][1] = 2.0*PI*tmp[1]/CellV;
1614 rtv[1][2] = 2.0*PI*tmp[2]/CellV;
1615 rtv[1][3] = 2.0*PI*tmp[3]/CellV;
1616
1617 Cross_Product(tv[3],tv[1],tmp);
1618 rtv[2][1] = 2.0*PI*tmp[1]/CellV;
1619 rtv[2][2] = 2.0*PI*tmp[2]/CellV;
1620 rtv[2][3] = 2.0*PI*tmp[3]/CellV;
1621
1622 Cross_Product(tv[1],tv[2],tmp);
1623 rtv[3][1] = 2.0*PI*tmp[1]/CellV;
1624 rtv[3][2] = 2.0*PI*tmp[2]/CellV;
1625 rtv[3][3] = 2.0*PI*tmp[3]/CellV;
1626
1627 G12 = rtv[1][1]*rtv[1][1] + rtv[1][2]*rtv[1][2] + rtv[1][3]*rtv[1][3];
1628 G22 = rtv[2][1]*rtv[2][1] + rtv[2][2]*rtv[2][2] + rtv[2][3]*rtv[2][3];
1629 G32 = rtv[3][1]*rtv[3][1] + rtv[3][2]*rtv[3][2] + rtv[3][3]*rtv[3][3];
1630
1631 G12 = sqrt(G12);
1632 G22 = sqrt(G22);
1633 G32 = sqrt(G32);
1634
1635 if (G12<G22) G0 = G12;
1636 else G0 = G22;
1637 if (G32<G0) G0 = G32;
1638
1639 AG = (G12+G22+G32)/3.0;
1640 DG = (fabs(G22-G12)+fabs(G32-G22)+fabs(G12-G32))/3.0;
1641
1642 Kerker_factor = 0.5/G0*(2.0*DG/AG+1.0);
1643
1644 if (myid==Host_ID && 0<level_stdout){
1645 printf("Automatic determination of Kerker_factor: %15.12f\n",Kerker_factor);
1646 }
1647 }
1648
1649 /****************************************************
1650 DFT-D:
1651 Added by Okuno in March 2011
1652 ****************************************************/
1653
1654 /* for vdW okuno */
1655 input_logical("scf.dftD",&dftD_switch,0);
1656 input_int("version.dftD",&version_dftD,2); /* Ellner */
1657
1658 if (version_dftD!=2 && version_dftD!=3){
1659
1660 if (myid==Host_ID){
1661 printf("version.dftD should be 1 or 2.\n");
1662 }
1663
1664 MPI_Finalize();
1665 exit(1);
1666 }
1667
1668
1669 /* DFT D okuno */
1670 if(dftD_switch){
1671 s_vec[1]="Ang"; s_vec[0]="AU"; /* Ellner CHANGED ORDER */
1672 i_vec[1]=0; i_vec[0]=1;
1673 input_string2int("DFTD.Unit",&unit_dftD,2,s_vec,i_vec);
1674 input_double("DFTD.rcut_dftD",&rcut_dftD,100.0);
1675
1676 input_double("DFTD.d",&beta_dftD,20.0);
1677 if(unit_dftD==0){ /* change Ang to AU */
1678 rcut_dftD = rcut_dftD/BohrR;
1679 }
1680
1681 input_double("DFTD.scale6",&scal6_dftD,0.75);
1682
1683 i_vec2[0]=1;
1684 i_vec2[1]=1;
1685 i_vec2[2]=1;
1686 input_intv("DFTD.IntDirection",3,i_vec,i_vec2);
1687 DFTD_IntDir1 = i_vec[0];
1688 DFTD_IntDir2 = i_vec[1];
1689 DFTD_IntDir3 = i_vec[2];
1690
1691 /*
1692 if(myid==Host_ID){
1693 printf("DFTD.Unit = %i \n",unit_dftD);
1694 printf("DFTD.rcut_dftD =%f\n",rcut_dftD);
1695 printf("DFTD.beta =%f\n",beta_dftD);
1696 printf("DFTD.scale6 =%f\n",scal6_dftD);
1697 }
1698 */
1699
1700 for (i=1; i<=atomnum; i++){
1701 Gxyz[i][60] = 1.0;
1702 }
1703
1704 if (fp=input_find("<DFTD.periodicity")) {
1705
1706 for (i=1; i<=atomnum; i++){
1707 fscanf(fp,"%d %lf",&j,&Gxyz[i][60]);
1708 }
1709
1710 if ( ! input_last("DFTD.periodicity>") ) {
1711 /* format error */
1712 printf("Format error for DFTD.periodicity\n");
1713 po++;
1714 }
1715 }
1716
1717 /* Ellner DFT-D3 */
1718 s_vec[0]="ZERO"; s_vec[1]="BJ";
1719 i_vec[0]=1; i_vec[1]=2;
1720 input_string2int("DFTD3.damp",&DFTD3_damp_dftD,2,s_vec,i_vec);
1721 if(version_dftD==3) input_double("DFTD.scale6",&s6_dftD,1.0); /* FIXED */
1722
1723 if ( XC_switch==4 ){
1724 input_double("DFTD.sr6",&sr6_dftD,1.217);
1725 input_double("DFTD.a1",&a1_dftD,0.4289);
1726 input_double("DFTD.a2",&a2_dftD,4.4407);
1727 if (DFTD3_damp_dftD==1) input_double("DFTD.scale8",&s8_dftD,0.722); /* GGA-PBE ZERO DAMP */
1728 if (DFTD3_damp_dftD==2) input_double("DFTD.scale8",&s8_dftD,0.7875);/* GGA-PBE BJ DAMP */
1729 }
1730 else{
1731 input_double("DFTD.scale8",&s8_dftD,0.5);
1732 input_double("DFTD.sr6",&sr6_dftD,1.0);
1733 input_double("DFTD.a1",&a1_dftD,0.5);
1734 input_double("DFTD.a2",&a2_dftD,5.0);
1735 }
1736
1737 input_double("DFTD.cncut_dftD",&cncut_dftD,40.0);
1738 if(unit_dftD==0){ /* change Ang to AU */
1739 cncut_dftD = cncut_dftD/BohrR;
1740 }
1741
1742 } /* if (dftD_switch) */
1743
1744 /************************************************************
1745 set fixed components of cell vectors in cell optimization
1746
1747 1: fixed
1748 0: relaxed
1749 ************************************************************/
1750
1751 for (i=1; i<=3; i++){
1752 for (j=1; j<=3; j++){
1753 Cell_Fixed_XYZ[i][j] = 0;
1754 }
1755 }
1756
1757 if (fp=input_find("<MD.Fixed.Cell.Vectors")) {
1758
1759 for (i=1; i<=3; i++){
1760 fscanf(fp,"%d %d %d",&Cell_Fixed_XYZ[i][1],&Cell_Fixed_XYZ[i][2],&Cell_Fixed_XYZ[i][3]);
1761 }
1762
1763 if ( ! input_last("MD.Fixed.Cell.Vectors>") ) {
1764 /* format error */
1765 printf("Format error for MD.Fixed.Cell.Vectors\n");
1766 po++;
1767 }
1768 }
1769
1770 /* RFC5 */
1771 if (MD_switch==18 && cellopt_swtich==5){
1772 for (i=1; i<=3; i++){
1773 for (j=1; j<=3; j++){
1774 Cell_Fixed_XYZ[i][j] = 0;
1775 }
1776 }
1777 }
1778
1779 /****************************************************
1780 set fixed atomic position in geometry optimization
1781 and MD:
1782
1783 1: fixed
1784 0: relaxed
1785 ****************************************************/
1786
1787 if (fp=input_find("<MD.Fixed.XYZ")) {
1788
1789 for (i=1; i<=atomnum; i++){
1790 fscanf(fp,"%d %d %d %d",
1791 &j,&atom_Fixed_XYZ[i][1],&atom_Fixed_XYZ[i][2],&atom_Fixed_XYZ[i][3]);
1792 }
1793
1794 if ( ! input_last("MD.Fixed.XYZ>") ) {
1795 /* format error */
1796 printf("Format error for MD.Fixed.XYZ\n");
1797 po++;
1798 }
1799 }
1800
1801 /****************************************************
1802 set initial velocities for MD
1803 ****************************************************/
1804
1805 MD_Init_Velocity = 0;
1806
1807 if (fp=input_find("<MD.Init.Velocity")) {
1808
1809 MD_Init_Velocity = 1;
1810
1811 for (i=1; i<=atomnum; i++){
1812
1813 fscanf(fp,"%d %lf %lf %lf",&j,&Gxyz[i][24],&Gxyz[i][25],&Gxyz[i][26]);
1814
1815 /***********************************************
1816 5.291772083*10^{-11} m / 2.418884*10^{-17} s
1817 = 2.1876917*10^6 m/s
1818 = 1 a.u. for velocity
1819
1820 1 m/s = 0.4571028 * 10^{-6} a.u.
1821 ***********************************************/
1822
1823 for (j=1; j<=3; j++){
1824 if (atom_Fixed_XYZ[i][j]==0){
1825 Gxyz[i][23+j] = Gxyz[i][23+j]*0.4571028*0.000001;
1826 Gxyz[i][26+j] = Gxyz[i][23+j];
1827 }
1828
1829 else{
1830 Gxyz[i][23+j] = 0.0;
1831 Gxyz[i][26+j] = 0.0;
1832 }
1833 }
1834 }
1835
1836 if ( ! input_last("MD.Init.Velocity>") ) {
1837 /* format error */
1838 printf("Format error for MD.Init.Velocity\n");
1839 po++;
1840 }
1841 }
1842
1843 /*************************************************************
1844 set atom groups for multi heat-bath thermostat for MD (VS4):
1845 added by T. Ohwaki
1846 *************************************************************/
1847
1848 input_int("MD.num.AtomGroup",&num_AtGr,1);
1849
1850 Allocate_Arrays(10);
1851 int chk_vs4;
1852
1853 for (k=1; k<=num_AtGr; k++){
1854 atnum_AtGr[k] = 0;
1855 }
1856
1857 if (fp=input_find("<MD.AtomGroup")){
1858
1859 for (i=1; i<=atomnum; i++){
1860 fscanf(fp,"%d %d", &j, &AtomGr[i]);
1861 chk_vs4 = 0;
1862
1863 for (k=1; k<=num_AtGr; k++){
1864 if (AtomGr[i] == k){
1865 atnum_AtGr[k]+=1;
1866 chk_vs4 = 1;
1867 }
1868 }
1869
1870 if (chk_vs4 ==0){
1871 printf("Please check your input for atom group!!\n");
1872 po++;
1873 }
1874 }
1875
1876 if ( ! input_last("MD.AtomGroup>") ) {
1877 /* format error */
1878 printf("Format error for MD.AtomGroup\n");
1879 po++;
1880 }
1881
1882 if (myid==Host_ID && 0<level_stdout){
1883 printf("\n");
1884 printf("*************************************************************** \n");
1885 printf(" Multi heat-bath MD calculation with velocity scaling method \n");
1886 printf(" \n");
1887 printf(" Number of atom groups = %d \n", num_AtGr);
1888 for (k=1; k<=num_AtGr; k++){
1889 printf(" Number of atoms in group %d = %d \n", k, atnum_AtGr[k]);
1890 }
1891 for (i=1; i<=atomnum; i++){
1892 printf(" Atom %d : group %d \n", i, AtomGr[i]);
1893 }
1894 printf(" \n");
1895 printf("*************************************************************** \n");
1896 printf("\n");
1897 }
1898 }
1899
1900 /****************************************************
1901 Starting point of LDA+U --- by MJ
1902 ****************************************************/
1903
1904 /* for LDA+U */
1905
1906 if (Hub_U_switch == 1){ /* --- MJ */
1907
1908 if (fp=input_find("<Hubbard.U.values")) {
1909
1910 /* initialize the U-values */
1911 for (i=0; i<SpeciesNum; i++){
1912 for (l=0; l<=Spe_MaxL_Basis[i]; l++){
1913 for (mul=0; mul<Spe_Num_Basis[i][l]; mul++){
1914 Hub_U_Basis[i][l][mul]=0.0 ;
1915 }
1916 }
1917 }
1918
1919 /* read the U-values from the '.dat' file */ /* --- MJ */
1920 for (i=0; i<SpeciesNum; i++){
1921 fscanf(fp,"%s",Species);
1922 j = Species2int(Species);
1923
1924 for (l=0; l<=Spe_MaxL_Basis[j]; l++){
1925 for (mul=0; mul<Spe_Num_Basis[j][l]; mul++){
1926 fscanf(fp,"%s %lf", buf, &Hub_U_Basis[j][l][mul]) ;
1927 }
1928 }
1929 }
1930
1931 if (! input_last("Hubbard.U.values>") ) {
1932 /* format error */
1933 printf("Format error for Hubbard.U.values\n");
1934 po++;
1935 }
1936
1937 } /* if (fp=input_find("<Hubbard.U.values")) */
1938
1939 for (i=0; i<SpeciesNum; i++){
1940 for (l=0; l<=Spe_MaxL_Basis[i]; l++){
1941 for (mul=0; mul<Spe_Num_Basis[i][l]; mul++){
1942 Hub_U_Basis[i][l][mul] = Hub_U_Basis[i][l][mul]/eV2Hartree;
1943 }
1944 }
1945 }
1946
1947 } /* if (Hub_U_switch == 1) */
1948
1949 /****************************************************
1950 Ending point of LDA+U --- by MJ
1951 ****************************************************/
1952
1953 /****************************************************
1954 the maximum number of SCF iterations
1955 ****************************************************/
1956
1957 input_int("scf.maxIter",&DFTSCF_loop,40);
1958 if (DFTSCF_loop<0) {
1959 printf("scf.maxIter should be over 0.\n");
1960 po++;
1961 }
1962
1963 /****************************************************
1964 parameters for O(N^2) method
1965 ****************************************************/
1966
1967 input_int("scf.Npoles.ON2",&ON2_Npoles,100);
1968 if (ON2_Npoles<0) {
1969 printf("scf.Npoles.ON2 should be over 0.\n");
1970 po++;
1971 }
1972
1973 /****************************************************
1974 Net charge of the system
1975 ****************************************************/
1976
1977 input_double("Atoms.NetCharge",&Given_Total_Charge,(double)0.0);
1978
1979 /*****************************************************
1980 if restart files exist, read data from them.
1981 default = off
1982 *****************************************************/
1983
1984 s_vec[0]="off"; i_vec[0] = 0; /* default, off for first MD, later on */
1985 s_vec[1]="alloff"; i_vec[1] = -1; /* switch off explicitly */
1986 s_vec[2]="on"; i_vec[2] = 1; /* switch on explicitly */
1987
1988 input_string2int("scf.restart", &Scf_RestartFromFile, 3, s_vec,i_vec);
1989
1990 /* check the number of processors */
1991
1992 if (Scf_RestartFromFile==1){
1993
1994 MPI_Comm_size(mpi_comm_level1,&numprocs1);
1995 sprintf(file_check,"%s%s_rst/%s.crst_check",filepath,filename,filename);
1996
1997 if ((fp = fopen(file_check,"r")) != NULL){
1998 fscanf(fp,"%d %d %d %d",&i_vec[0],&i_vec[1],&i_vec[2],&i_vec[3]);
1999 fclose(fp);
2000
2001 if (i_vec[0]!=numprocs1){
2002 Remake_RestartFile(numprocs1,i_vec[0],i_vec[1],i_vec[2],i_vec[3]);
2003 }
2004 }
2005 else{
2006 printf("Failure of saving %s\n",file_check);
2007 }
2008 }
2009
2010 /****************************************************
2011 Band dispersion
2012 ****************************************************/
2013
2014 input_logical("Band.dispersion",&Band_disp_switch, 0);
2015
2016 Band_kpath=NULL;
2017 Band_kname=NULL;
2018 Band_N_perpath=NULL;
2019 input_int("Band.Nkpath",&Band_Nkpath,0);
2020 if (2<=level_stdout) printf("Band.Nkpath=%d\n",Band_Nkpath);
2021
2022 if (Band_Nkpath>0) {
2023
2024 Band_kPathUnit=0;
2025 if (fp=input_find("<Band.kpath.UnitCell") ) {
2026 Band_kPathUnit=1;
2027 for (i=1;i<=3;i++) {
2028 fscanf(fp,"%lf %lf %lf",&Band_UnitCell[i][1],
2029 &Band_UnitCell[i][2],&Band_UnitCell[i][3]);
2030 }
2031 if ( ! input_last("Band.kpath.UnitCell>") ) {
2032 /* format error */
2033 printf("Format error near Band.kpath.UnitCell>\n");
2034 po++;
2035 }
2036
2037 if (myid==Host_ID){
2038 for (i=1;i<=3;i++) {
2039 printf("%lf %lf %lf\n",Band_UnitCell[i][1],
2040 Band_UnitCell[i][2],Band_UnitCell[i][3]);
2041 }
2042 }
2043
2044 if (unitvector_unit==0) {
2045 for (i=1;i<=3;i++)
2046 for (j=1;j<=3;j++)
2047 Band_UnitCell[i][j]=Band_UnitCell[i][j]/BohrR;
2048 }
2049 }
2050 else {
2051 for (i=1;i<=3;i++)
2052 for (j=1;j<=3;j++)
2053 Band_UnitCell[i][j]=tv[i][j];
2054 }
2055
2056 /* allocate */
2057
2058 Band_N_perpath=(int*)malloc(sizeof(int)*(Band_Nkpath+1));
2059 for (i=0; i<(Band_Nkpath+1); i++) Band_N_perpath[i] = 0;
2060
2061 Band_kpath = (double***)malloc(sizeof(double**)*(Band_Nkpath+1));
2062 for (i=0; i<(Band_Nkpath+1); i++){
2063 Band_kpath[i] = (double**)malloc(sizeof(double*)*3);
2064 for (j=0; j<3; j++){
2065 Band_kpath[i][j] = (double*)malloc(sizeof(double)*4);
2066 for (k=0; k<4; k++) Band_kpath[i][j][k] = 0.0;
2067 }
2068 }
2069
2070 Band_kname = (char***)malloc(sizeof(char**)*(Band_Nkpath+1));
2071 for (i=0; i<(Band_Nkpath+1); i++){
2072 Band_kname[i] = (char**)malloc(sizeof(char*)*3);
2073 for (j=0; j<3; j++){
2074 Band_kname[i][j] = (char*)malloc(sizeof(char)*YOUSO10);
2075 }
2076 }
2077
2078 /* end of allocation */
2079
2080 if (myid==Host_ID && 2<=level_stdout) printf("kpath\n");
2081
2082 if (fp=input_find("<Band.kpath") ) {
2083 for (i=1; i<=Band_Nkpath; i++) {
2084 fscanf(fp,"%d %lf %lf %lf %lf %lf %lf %s %s",
2085 &Band_N_perpath[i] ,
2086 &Band_kpath[i][1][1], &Band_kpath[i][1][2],&Band_kpath[i][1][3],
2087 &Band_kpath[i][2][1], &Band_kpath[i][2][2],&Band_kpath[i][2][3],
2088 Band_kname[i][1],Band_kname[i][2]);
2089
2090 if (myid==Host_ID && 2<=level_stdout){
2091 printf("%d (%lf %lf %lf) (%lf %lf %lf) %s %s\n",
2092 Band_N_perpath[i] ,
2093 Band_kpath[i][1][1], Band_kpath[i][1][2],Band_kpath[i][1][3],
2094 Band_kpath[i][2][1], Band_kpath[i][2][2],Band_kpath[i][2][3],
2095 Band_kname[i][1],Band_kname[i][2]);
2096 }
2097
2098 }
2099 if ( ! input_last("Band.kpath>") ) {
2100 /* format error */
2101 printf("Format error near Band.kpath>\n");
2102 po++;
2103 }
2104 }
2105 else {
2106 /* format error */
2107 printf("<Band.kpath is necessary.\n");
2108 po++;
2109 }
2110 if (Band_kPathUnit){
2111 kpath_changeunit( tv, Band_UnitCell, Band_Nkpath, Band_kpath );
2112 }
2113
2114 }
2115
2116
2117
2118
2119 /****************************************************
2120 One dimentional FFT
2121 ****************************************************/
2122
2123 input_int("1DFFT.NumGridK",&Ngrid_NormK,900);
2124
2125 if (Ngrid_NormK<GL_Mesh)
2126 List_YOUSO[15] = GL_Mesh + 10;
2127 else
2128 List_YOUSO[15] = Ngrid_NormK + 2;
2129
2130 input_double("1DFFT.EnergyCutoff",&ecutoff1dfft,(double)3600.0);
2131 PAO_Nkmax=sqrt(ecutoff1dfft);
2132
2133 input_int("1DFFT.NumGridR",&OneD_Grid,900);
2134
2135 /************************************************************
2136 it it not allowed to make the simultaneuos specification
2137 of spin non-collinear and orbital optimization
2138 ************************************************************/
2139
2140 if (SpinP_switch==3 && Cnt_switch==1){
2141 if (myid==Host_ID){
2142 printf("unsupported:\n");
2143 printf("simultaneuos specification of spin non-collinear\n");
2144 printf("and orbital optimization\n");
2145 }
2146
2147 MPI_Finalize();
2148 exit(1);
2149 }
2150
2151 /************************************************************
2152 it it not allowed to make the simultaneuos specification
2153 of NEGF and orbital optimization
2154 ************************************************************/
2155
2156 if (Solver==4 && Cnt_switch==1){
2157 if (myid==Host_ID){
2158 printf("unsupported:\n");
2159 printf("simultaneuos specification of NEGF and orbital optimization\n");
2160 }
2161
2162 MPI_Finalize();
2163 exit(1);
2164 }
2165
2166 s_vec[0]="Symmetrical"; s_vec[1]="Free"; s_vec[2]="Simple";
2167 i_vec[0]=1; i_vec[1]=0; i_vec[2]=2;
2168 input_string2int("orbitalOpt.InitCoes",&SICnt_switch,3,s_vec,i_vec);
2169 input_int("orbitalOpt.scf.maxIter",&orbitalOpt_SCF,40);
2170 input_int("orbitalOpt.Opt.maxIter",&orbitalOpt_MD,100);
2171 input_int("orbitalOpt.per.MDIter",&orbitalOpt_per_MDIter,1000000);
2172 input_double("orbitalOpt.criterion",&orbitalOpt_criterion,(double)1.0e-4);
2173 input_double("orbitalOpt.SD.step",&orbitalOpt_SD_step,(double)0.001);
2174 input_int("orbitalOpt.HistoryPulay",&orbitalOpt_History,20);
2175 input_int("orbitalOpt.StartPulay",&orbitalOpt_StartPulay,10);
2176 input_logical("orbitalOpt.Force.Skip" ,&orbitalOpt_Force_Skip,0);
2177
2178 s_vec[0]="DIIS"; s_vec[1]="EF";
2179 i_vec[0]=1; i_vec[1]=2;
2180 input_string2int("orbitalOpt.Opt.Method",&OrbOpt_OptMethod,2,s_vec,i_vec);
2181
2182 /****************************************************
2183 order-N method for SCF
2184 ****************************************************/
2185
2186 input_double("orderN.HoppingRanges",&BCR,(double)7.0);
2187 BCR=BCR/BohrR;
2188
2189 input_int("orderN.NumHoppings",&NOHS_L,2);
2190 if (Solver==2 || Solver==3 || Solver==4 || Solver==7 || Solver==9){
2191 NOHS_L = 1;
2192 BCR = 1.0;
2193 }
2194
2195 NOHS_C = NOHS_L;
2196
2197 /* start EC */
2198
2199 input_int("orderN.EC.Sub.Dim",&EC_Sub_Dim,400);
2200
2201 /* end EC */
2202
2203 /* start Krylov */
2204
2205 /* input_int("orderN.Kgrid",&orderN_Kgrid,5); */
2206
2207 input_int("orderN.KrylovH.order",&KrylovH_order,400);
2208 input_int("orderN.KrylovS.order",&KrylovS_order,4*KrylovH_order);
2209 input_logical("orderN.Recalc.Buffer",&recalc_EM,1);
2210 input_logical("orderN.Exact.Inverse.S",&EKC_Exact_invS_flag,1);
2211 input_logical("orderN.Expand.Core",&EKC_expand_core_flag,1);
2212 input_logical("orderN.Inverse.S",&EKC_invS_flag,0); /* ?? */
2213
2214 if (EKC_Exact_invS_flag==0){
2215 EKC_invS_flag = 1;
2216 }
2217 else {
2218 EKC_invS_flag = 0;
2219 }
2220
2221 orderN_FNAN_SNAN_flag = 0;
2222
2223 if (fp=input_find("<orderN.FNAN+SNAN")) {
2224
2225 orderN_FNAN_SNAN_flag = 1;
2226
2227 for (i=1; i<=atomnum; i++){
2228 fscanf(fp,"%i %i",&j,&orderN_FNAN_SNAN[i]);
2229
2230 if (i!=j) {
2231 if (myid==Host_ID){
2232 printf("Format error for orderN.FNAN+SNAN\n");
2233 }
2234
2235 MPI_Finalize();
2236 exit(0);
2237 }
2238 }
2239
2240 if (! input_last("orderN.FNAN+SNAN>")) {
2241 /* format error */
2242
2243 po++;
2244
2245 if (myid==Host_ID){
2246 printf("Format error for orderN.FNAN+SNAN\n");
2247 }
2248 MPI_Finalize();
2249 exit(0);
2250 }
2251 }
2252
2253 /* end Krylov */
2254
2255 input_int("orderN.RecursiveLevels",&rlmax,10);
2256 List_YOUSO[3] = rlmax + 3;
2257
2258 s_vec[0]="NO"; s_vec[1]="SQRT";
2259 i_vec[0]=1 ; i_vec[1]=2;
2260 input_string2int("orderN.TerminatorType",&T_switch,2, s_vec,i_vec);
2261
2262 s_vec[0]="Recursion"; s_vec[1]="Divide";
2263 i_vec[0]=1; i_vec[1]=2;
2264 input_string2int("orderN.InverseSolver",&IS_switch,2,s_vec,i_vec);
2265
2266 input_int("orderN.InvRecursiveLevels",&rlmax_IS,30);
2267 List_YOUSO[9] = rlmax_IS + 3;
2268
2269 input_double("orderN.ChargeDeviation",&CN_Error,(double)1.0e-7);
2270 input_double("orderN.InitChemPot",&ChemP,(double)-3.0); /* in eV */
2271 ChemP = ChemP/eV2Hartree;
2272
2273 input_int("orderN.AvNumTerminater",&Av_num,1);
2274 input_int("orderN.NumPoles",&POLES,10);
2275
2276 /****************************************************
2277 control patameters for outputting wave functions
2278 ****************************************************/
2279
2280 input_logical("MO.fileout",&MO_fileout,0);
2281 input_int("num.HOMOs",&num_HOMOs,1);
2282 input_int("num.LUMOs",&num_LUMOs,1);
2283
2284 if (MO_fileout==0){
2285 num_HOMOs = 1;
2286 num_LUMOs = 1;
2287 }
2288
2289 if ((Solver!=2 && Solver!=3 && Solver!=7) && MO_fileout==1){
2290
2291 s_vec[0]="Recursion"; s_vec[1]="Cluster"; s_vec[2]="Band";
2292 s_vec[3]="NEGF"; s_vec[4]="DC"; s_vec[5]="GDC";
2293 s_vec[6]="Cluster-DIIS"; s_vec[7]="Krylov"; s_vec[8]="Cluster2";
2294
2295 printf("MO.fileout=ON is not supported in case of scf.EigenvalueSolver=%s\n",
2296 s_vec[Solver-1]);
2297 printf("MO.fileout is changed to OFF\n");
2298 MO_fileout = 0;
2299 }
2300
2301 if ( (Solver==2 || Solver==3 || Solver==7) && MO_fileout==1 ){
2302 List_YOUSO[31] = num_HOMOs;
2303 List_YOUSO[32] = num_LUMOs;
2304 }
2305 else{
2306 List_YOUSO[31] = 1;
2307 List_YOUSO[32] = 1;
2308 }
2309
2310 /* for bulk */
2311
2312 input_int("MO.Nkpoint",&MO_Nkpoint,1);
2313 if (MO_Nkpoint<0){
2314 printf("MO_Nkpoint must be positive.\n");
2315 po++;
2316 }
2317
2318 if (MO_fileout && Solver==2 && MO_Nkpoint!=1){
2319 if (myid==Host_ID){
2320 printf("MO.Nkpoint should be 1 for the cluster calculation.\n");
2321 }
2322 MPI_Finalize();
2323 exit(0);
2324 }
2325
2326 if ( (Solver==2 || Solver==3 || Solver==7) && MO_fileout==1 ){
2327 List_YOUSO[33] = MO_Nkpoint;
2328 }
2329 else{
2330 List_YOUSO[33] = 1;
2331 }
2332
2333 /* memory allocation */
2334 Allocate_Arrays(5);
2335
2336 if (fp=input_find("<MO.kpoint")) {
2337 for (i=0; i<MO_Nkpoint; i++){
2338 fscanf(fp,"%lf %lf %lf",&MO_kpoint[i][1],&MO_kpoint[i][2],&MO_kpoint[i][3]);
2339
2340 if (2<=level_stdout){
2341 printf("<Input_std> MO_kpoint %2d %10.6f %10.6f %10.6f\n",
2342 i,MO_kpoint[i][1],MO_kpoint[i][2],MO_kpoint[i][3]);
2343 }
2344 }
2345 if (!input_last("MO.kpoint>")) {
2346 /* format error */
2347 printf("Format error for MO.kpoint\n");
2348 po++;
2349 }
2350
2351 if (Band_kPathUnit){
2352 kpoint_changeunit(tv,Band_UnitCell,MO_Nkpoint,MO_kpoint);
2353 }
2354 }
2355
2356 /****************************************************
2357 Natural Bond Orbital Analysis
2358 ****************************************************/
2359
2360 s_vec[0]="off"; s_vec[1]="on1"; s_vec[2]="on2"; s_vec[3]="on3"; s_vec[4]="on4";
2361 i_vec[0]=0 ; i_vec[1]=1 ; i_vec[2]=2 ; i_vec[3]=3 ; i_vec[4]=4 ;
2362 input_string2int("NBO.switch", &NBO_switch, 5, s_vec,i_vec);
2363
2364 if (NBO_switch!=0){
2365
2366 input_logical("NAO.only",&NAO_only,1);
2367
2368 input_double("NAO.threshold",&NAO_Occ_or_Ryd,0.85);
2369
2370 input_int("NBO.Num.CenterAtoms",&Num_NBO_FCenter,1);
2371 if (Num_NBO_FCenter<0){
2372 printf("NBO.Num.CenterAtoms must be positive.\n");
2373 po++;
2374 }
2375
2376 if (Num_NBO_FCenter>atomnum){
2377 printf("NBO ERROR: NBO.Num.CenterAtoms should be less than the total number of atoms.");
2378 po++;
2379 }
2380
2381 Allocate_Arrays(11);
2382
2383 if (fp=input_find("<NBO.CenterAtoms")) {
2384 for (i=0; i<Num_NBO_FCenter; i++){
2385 fscanf(fp,"%d",&NBO_FCenter[i]);
2386 if (NBO_FCenter[i] <= 0 || NBO_FCenter[i] > atomnum){
2387 printf ("NBO ERROR: Index of center atom is wrong.");
2388 po++;
2389 }
2390
2391 if (0<=level_stdout && myid == Host_ID){
2392 printf("<Input_std> NBO_Center %d : %d \n",i,NBO_FCenter[i]);
2393 }
2394 }
2395 }
2396
2397 input_logical("NHO.fileout",&NHO_fileout,0);
2398 input_logical("NBO.fileout",&NBO_fileout,0);
2399
2400 input_logical("NBO.SmallCell.Switch",&NBO_SmallCell_Switch,0);
2401
2402 if (fp=input_find("<NBO.SmallCell")) {
2403 for (i=1; i<=3; i++){
2404 fscanf(fp,"%lf %lf",&NBO_SmallCellFrac[i][1],&NBO_SmallCellFrac[i][2]);
2405
2406 if (1<=level_stdout && myid == Host_ID){
2407 printf("<Input_std> NBO_SmallCell %10.6f %10.6f \n",
2408 NBO_SmallCellFrac[i][1],NBO_SmallCellFrac[i][2]);
2409 }
2410 }
2411
2412 if (!input_last("NBO.SmallCell>")) {
2413 printf("Format error for NBO.SmallCell\n");
2414 po++;
2415 }
2416 }
2417
2418 /*
2419 input_int("NAO.Nkpoint",&NAO_Nkpoint,1);
2420 if (NAO_Nkpoint<0){
2421 printf("NAO_Nkpoint must be positive.\n");
2422 po++;
2423 }
2424
2425 if (fp=input_find("<NAO.kpoint")) {
2426 for (i=0; i<NAO_Nkpoint; i++){
2427 fscanf(fp,"%lf %lf %lf",&NAO_kpoint[i][1],&NAO_kpoint[i][2],&NAO_kpoint[i][3]);
2428
2429 if (0<=level_stdout && myid == Host_ID){
2430 printf("<Input_std> NAO_kpoint %2d %10.6f %10.6f %10.6f\n",
2431 i,NAO_kpoint[i][1],NAO_kpoint[i][2],NAO_kpoint[i][3]);
2432 }
2433 }
2434 if (!input_last("NAO.kpoint>")) {
2435 printf("Format error for MO.kpoint\n");
2436 po++;
2437 }
2438
2439 if (Band_kPathUnit){
2440 kpoint_changeunit2(tv,Band_UnitCell,NAO_Nkpoint,NAO_kpoint);
2441 }
2442
2443 }
2444 */
2445
2446 }
2447
2448 /* added by Chi-Cheng for unfolding */
2449 /****************************************************
2450 Unfolding Band Structure
2451 ****************************************************/
2452
2453 input_logical("Unfolding.Electronic.Band",&unfold_electronic_band,0);
2454
2455 if (unfold_electronic_band==1) {
2456
2457 unfold_abc=(double**)malloc(sizeof(double*)*3);
2458 for (i=0; i<3; i++) unfold_abc[i]=(double*)malloc(sizeof(double)*3);
2459 if (fp=input_find("<Unfolding.ReferenceVectors")) {
2460 for (i=0; i<3; i++){
2461 fscanf(fp,"%lf %lf %lf",&unfold_abc[i][0],&unfold_abc[i][1],&unfold_abc[i][2]);
2462 }
2463 if (!input_last("Unfolding.ReferenceVectors>")) {
2464 /* format error */
2465 printf("Format error in Unfolding.ReferenceVectors\n");
2466 po++;
2467 }
2468 /* Ang to AU */
2469 if (unitvector_unit==0) {
2470 for (i=0; i<3; i++){
2471 unfold_abc[i][0] = unfold_abc[i][0]/BohrR;
2472 unfold_abc[i][1] = unfold_abc[i][1]/BohrR;
2473 unfold_abc[i][2] = unfold_abc[i][2]/BohrR;
2474 }
2475 }
2476 }
2477 else {
2478 for (i=0; i<3; i++) {
2479 unfold_abc[i][0] = tv[i+1][1];
2480 unfold_abc[i][1] = tv[i+1][2];
2481 unfold_abc[i][2] = tv[i+1][3];
2482 }
2483 }
2484
2485 unfold_origin=(double*)malloc(sizeof(double)*3);
2486 if (fp=input_find("<Unfolding.Referenceorigin")) {
2487 fscanf(fp,"%lf %lf %lf",&unfold_origin[0],&unfold_origin[1],&unfold_origin[2]);
2488 if (!input_last("Unfolding.Referenceorigin>")) {
2489 /* format error */
2490 printf("Format error in Unfolding.Referenceorigin\n");
2491 po++;
2492 }
2493 /* Ang to AU */
2494 if (unitvector_unit==0) {
2495 for (i=0; i<3; i++) unfold_origin[i] = unfold_origin[i]/BohrR;
2496 }
2497 }
2498 else {
2499 /* default value will be calculated later on */
2500 unfold_origin[0]=-0.9999900023114;
2501 unfold_origin[1]=-0.9999900047614;
2502 unfold_origin[2]=-0.9999900058914;
2503 }
2504
2505 input_double("Unfolding.LowerBound",&unfold_lbound,-10.);
2506 input_double("Unfolding.UpperBound",&unfold_ubound,10.); /* in eV */
2507 if (unfold_lbound>=unfold_ubound){
2508 printf("Unfolding: Lower bound must be lower than the upper bound.\n");
2509 po++;
2510 }
2511 unfold_lbound=unfold_lbound/eV2Hartree;
2512 unfold_ubound=unfold_ubound/eV2Hartree;
2513
2514 unfold_mapN2n=(int*)malloc(sizeof(int)*atomnum);
2515 if (fp=input_find("<Unfolding.Map")) {
2516 for (i=0; i<atomnum; i++){
2517 fscanf(fp,"%i %i",&j,&unfold_mapN2n[i]);
2518 if ((j!=i+1)||(unfold_mapN2n[i]<0)) { printf("Format error in Unfolding.Map! (Values cannot be negative.)\n"); po++;}
2519 }
2520 if (!input_last("Unfolding.Map>")) {
2521 /* format error */
2522 printf("Format error in Unfolding.Map\n");
2523 po++;
2524 }
2525 }
2526 else {
2527 for (i=0; i<atomnum; i++) unfold_mapN2n[i]=i;
2528 }
2529
2530 if (Solver!=2 && Solver!=3){
2531
2532 s_vec[0]="Recursion"; s_vec[1]="Cluster"; s_vec[2]="Band";
2533 s_vec[3]="NEGF"; s_vec[4]="DC"; s_vec[5]="GDC";
2534 s_vec[6]="Cluster-DIIS"; s_vec[7]="Krylov"; s_vec[8]="Cluster2";
2535
2536 printf("Unfolding.fileout=ON is not supported in case of scf.EigenvalueSolver=%s\n",
2537 s_vec[Solver-1]);
2538 printf("Unfolding.fileout is changed to OFF\n");
2539 unfold_electronic_band = 0;
2540 }
2541
2542 /* for bulk */
2543
2544 input_int("Unfolding.Nkpoint",&unfold_Nkpoint,0);
2545 if (unfold_Nkpoint<0){
2546 printf("Unfolding.Nkpoint must be positive.\n");
2547 po++;
2548 }
2549
2550 input_int("Unfolding.desired_totalnkpt",&unfold_nkpts,0);
2551 if (unfold_nkpts<0){
2552 printf("Unfolding.desired_totalnkpt must be positive.\n");
2553 po++;
2554 }
2555
2556 /* memory allocation */
2557 unfold_kpoint = (double**)malloc(sizeof(double*)*(unfold_Nkpoint+1));
2558 for (i=0; i<(unfold_Nkpoint+1); i++){
2559 unfold_kpoint[i] = (double*)malloc(sizeof(double)*4);
2560 }
2561
2562 unfold_kpoint_name = (char**)malloc(sizeof(char*)*(unfold_Nkpoint+1));
2563 for (i=0; i<(unfold_Nkpoint+1); i++){
2564 unfold_kpoint_name[i] = (char*)malloc(sizeof(char)*YOUSO10);
2565 }
2566
2567 if (fp=input_find("<Unfolding.kpoint")) {
2568 for (i=0; i<unfold_Nkpoint; i++){
2569 fscanf(fp,"%s %lf %lf %lf",
2570 unfold_kpoint_name[i],
2571 &unfold_kpoint[i][1],
2572 &unfold_kpoint[i][2],
2573 &unfold_kpoint[i][3]);
2574
2575 if (2<=level_stdout){
2576 printf("<Input_std> Unfolding.kpoint %2d %10.6f %10.6f %10.6f\n",
2577 i,unfold_kpoint[i][1],unfold_kpoint[i][2],unfold_kpoint[i][3]);
2578 }
2579 }
2580 if (!input_last("Unfolding.kpoint>")) {
2581 /* format error */
2582 printf("Format error for Unfolding.kpoint\n");
2583 po++;
2584 }
2585 }
2586 }
2587 /* end unfolding */
2588
2589 /****************************************************
2590 OutData_bin_flag
2591 ****************************************************/
2592
2593 input_logical("OutData.bin.flag",&OutData_bin_flag,0); /* default=off */
2594
2595 /****************************************************
2596 for output of contracted orbitals
2597 ****************************************************/
2598
2599 input_logical("CntOrb.fileout",&CntOrb_fileout,0);
2600 if ((!(Cnt_switch==1 && RCnt_switch==1)) && CntOrb_fileout==1){
2601 printf("CntOrb.fileout=on is valid in case of orbitalOpt.Method=Restricted or Speciese\n");
2602 po++;
2603 }
2604
2605 if (CntOrb_fileout==1){
2606 input_int("Num.CntOrb.Atoms",&Num_CntOrb_Atoms,1);
2607 CntOrb_Atoms = (int*)malloc(sizeof(int)*Num_CntOrb_Atoms);
2608
2609 if (fp=input_find("<Atoms.Cont.Orbitals")) {
2610 for (i=0; i<Num_CntOrb_Atoms; i++){
2611 fscanf(fp,"%i",&CntOrb_Atoms[i]);
2612 if (CntOrb_Atoms[i]<1 || atomnum<CntOrb_Atoms[i]){
2613 printf("Invalid atom in <Atoms.Cont.Orbitals\n" );
2614 po++;
2615 }
2616 }
2617
2618 if (!input_last("Atoms.Cont.Orbitals>")) {
2619 /* format error */
2620 printf("Format error for Atoms.Cont.Orbitals\n");
2621 po++;
2622 }
2623 }
2624 }
2625
2626 /****************************************************
2627 external electric field
2628 ****************************************************/
2629
2630 r_vec[0]=0.0; r_vec[1]=0.0; r_vec[2]=0.0;
2631 input_doublev("scf.Electric.Field",3,E_Field,r_vec);
2632
2633 if ( 1.0e-50<fabs(E_Field[0]) || 1.0e-50<fabs(E_Field[1]) || 1.0e-50<fabs(E_Field[2]) ){
2634
2635 E_Field_switch = 1;
2636
2637 /*******************************************
2638 unit transformation
2639
2640 V/m = J/C/m
2641 = 1/(4.35975*10^{-18}) Hatree
2642 *(1/(1.602177*10^{-19}) e )^{-1}
2643 *(1/(0.5291772*10^{-10}) a0 )^{-1}
2644 = 0.1944688 * 10^{-11} Hartree/e/a0
2645
2646 input unit: GV/m = 10^9 V/m
2647 used unit: Hartree/e/a0
2648
2649 GV/m = 0.1944688 * 10^{-2} Hartree/e/a0
2650 *******************************************/
2651
2652 length = sqrt( Dot_Product(tv[1], tv[1]) );
2653 x = E_Field[0]*tv[1][1]/length;
2654 y = E_Field[0]*tv[1][2]/length;
2655 z = E_Field[0]*tv[1][3]/length;
2656
2657 length = sqrt( Dot_Product(tv[2], tv[2]) );
2658 x += E_Field[1]*tv[2][1]/length;
2659 y += E_Field[1]*tv[2][2]/length;
2660 z += E_Field[1]*tv[2][3]/length;
2661
2662 length = sqrt( Dot_Product(tv[3], tv[3]) );
2663 x += E_Field[2]*tv[3][1]/length;
2664 y += E_Field[2]*tv[3][2]/length;
2665 z += E_Field[2]*tv[3][3]/length;
2666
2667 length = sqrt( x*x + y*y + z*z );
2668 x = x/length;
2669 y = y/length;
2670 z = z/length;
2671
2672 if (myid==Host_ID && 0<level_stdout){
2673 printf("\n");
2674 printf("<Applied External Electric Field>\n");
2675 printf(" direction (x, y, z) = %10.5f %10.5f %10.5f\n",x,y,z);
2676 printf(" magnitude (10^9 V/m) = %10.5f\n\n",length);
2677 }
2678
2679 E_Field[0] = 0.1944688*0.01*E_Field[0];
2680 E_Field[1] = 0.1944688*0.01*E_Field[1];
2681 E_Field[2] = 0.1944688*0.01*E_Field[2];
2682 }
2683 else {
2684 E_Field_switch = 0;
2685 }
2686
2687 /****************************************************
2688 DOS, PDOS
2689 ****************************************************/
2690
2691 input_logical("OpticalConductivity.fileout",&Opticalconductivity_fileout,0);
2692
2693 input_logical("Dos.fileout",&Dos_fileout,0);
2694 input_logical("DosGauss.fileout",&DosGauss_fileout,0);
2695 input_int("DosGauss.Num.Mesh",&DosGauss_Num_Mesh,200);
2696 input_double("DosGauss.Width",&DosGauss_Width,0.2); /* in eV */
2697
2698 /* change the unit from eV to Hartree */
2699 DosGauss_Width = DosGauss_Width/eV2Hartree;
2700
2701 if (Dos_fileout && DosGauss_fileout){
2702
2703 if (myid==Host_ID){
2704 printf("Dos.fileout and DosGauss.fileout are mutually exclusive.\n");
2705 }
2706 MPI_Finalize();
2707 exit(1);
2708 }
2709
2710 if ( DosGauss_fileout && (Solver!=3 && PeriodicGamma_flag!=1) ){ /* band */
2711
2712 if (myid==Host_ID){
2713 printf("DosGauss.fileout is supported for only band calculation.\n");
2714 }
2715 MPI_Finalize();
2716 exit(1);
2717 }
2718
2719 r_vec[0]=-20.0; r_vec[1]=20.0;
2720 input_doublev("Dos.Erange",2,Dos_Erange,r_vec);
2721 /* change the unit from eV to Hartree */
2722 Dos_Erange[0]= Dos_Erange[0]/eV2Hartree;
2723 Dos_Erange[1]= Dos_Erange[1]/eV2Hartree;
2724
2725 i_vec[0]=Kspace_grid1; i_vec[1]=Kspace_grid2, i_vec[2]=Kspace_grid3;
2726 input_intv("Dos.Kgrid",3,Dos_Kgrid,i_vec);
2727
2728 /**********************************************************
2729 calculation of partial charge for scanning tunneling
2730 microscopy (STM) simulations by the Tersoff-Hamann scheme
2731 ***********************************************************/
2732
2733 input_logical("partial.charge",&cal_partial_charge,0); /* default=off */
2734
2735 if (cal_partial_charge && (Dos_fileout==0 && DosGauss_fileout==0)){
2736 if (myid==Host_ID){
2737 printf("partial.charge can be switched on in case that\n");
2738 printf("either Dos.fileout or DosGauss.fileout is switched on.\n");
2739
2740 }
2741 MPI_Finalize();
2742 exit(1);
2743 }
2744
2745 if ( cal_partial_charge && (Solver!=2 && Solver!=3) ){
2746 if (myid==Host_ID){
2747 printf("The calculation of partial charge is supported only for the Cluster and Band methods.\n");
2748
2749 }
2750 MPI_Finalize();
2751 exit(1);
2752 }
2753
2754 input_double("partial.charge.energy.window",&ene_win_partial_charge,0.0); /* in eV */
2755 ene_win_partial_charge /= eV2Hartree; /* eV to Hartee */
2756
2757 /****************************************************
2758 write a binary file, filename.scfout, which includes
2759 Hamiltonian, overlap, and density matrices.
2760 ****************************************************/
2761
2762 input_logical("HS.fileout",&HS_fileout,0);
2763
2764 /****************************************************
2765 Energy decomposition
2766 ****************************************************/
2767
2768 input_logical("Energy.Decomposition",&Energy_Decomposition_flag,0);
2769
2770 if (Energy_Decomposition_flag==1 && Cnt_switch==1){
2771 if (myid==Host_ID){
2772 printf("Energy decomposition is not supported for orbital optimization.\n");
2773 }
2774 MPI_Finalize();
2775 exit(0);
2776 }
2777
2778 /****************************************************
2779 Voronoi charge
2780 ****************************************************/
2781
2782 input_logical("Voronoi.charge",&Voronoi_Charge_flag,0);
2783
2784 /****************************************************
2785 Voronoi orbital moment
2786 ****************************************************/
2787
2788 input_logical("Voronoi.orbital.moment",&Voronoi_OrbM_flag,0);
2789
2790 /****************************************************
2791 parameters on Wannier funtions by hmweng
2792 ****************************************************/
2793
2794 input_logical("Wannier.Func.Calc",&Wannier_Func_Calc,0); /* default=off */
2795 input_logical("Wannier90.fileout",&Wannier90_fileout,0); /* default=off */
2796
2797 if (Wannier_Func_Calc){
2798
2799 int **tmp_Wannier_Pro_SelMat;
2800 double ***tmp_Wannier_Projector_Hybridize_Matrix;
2801
2802 input_int("Wannier.Func.Num", &Wannier_Func_Num,0);
2803 input_double("Wannier.Outer.Window.Bottom",&Wannier_Outer_Window_Bottom,(double)-10.0); /* eV */
2804 input_double("Wannier.Outer.Window.Top", &Wannier_Outer_Window_Top,(double)10.0); /* eV */
2805
2806 input_double("Wannier.Inner.Window.Bottom",&Wannier_Inner_Window_Bottom,(double)0.0); /* eV */
2807 input_double("Wannier.Inner.Window.Top", &Wannier_Inner_Window_Top,(double)0.0); /* eV */
2808 input_logical("Wannier.Initial.Guess", &Wannier_Initial_Guess,1); /* on|off */
2809
2810 if( (Wannier_Outer_Window_Top-Wannier_Outer_Window_Bottom)<=0.0 ){
2811
2812 if (myid==Host_ID){
2813 printf("Error:WF The top of the OUTER window should be higher than the bottom.\n");
2814 }
2815
2816 MPI_Finalize();
2817 exit(0);
2818 }
2819
2820 if( (Wannier_Inner_Window_Top-Wannier_Inner_Window_Bottom)<0.0){
2821
2822 if (myid==Host_ID){
2823 printf("Error:WF The top of the INNER window should be higher than the bottom.\n");
2824 }
2825
2826 MPI_Finalize();
2827 exit(0);
2828 }
2829
2830 if ( 0.0<(Wannier_Inner_Window_Top-Wannier_Inner_Window_Bottom) ){
2831
2832 if( Wannier_Inner_Window_Bottom<Wannier_Outer_Window_Bottom
2833 || Wannier_Inner_Window_Top>Wannier_Outer_Window_Top ){
2834
2835 if (myid==Host_ID){
2836 printf("Error:WF INNER window (%10.5f,%10.5f) is not inside of OUTER window (%10.5f,%10.5f).\n",
2837 Wannier_Inner_Window_Bottom,Wannier_Inner_Window_Top,
2838 Wannier_Outer_Window_Bottom,Wannier_Outer_Window_Top);
2839 }
2840
2841 MPI_Finalize();
2842 exit(0);
2843 }
2844 }
2845
2846
2847 if(fabs(Wannier_Inner_Window_Bottom-Wannier_Inner_Window_Top)<1e-6){
2848
2849 if (myid==Host_ID){
2850 printf("Message:WF The top and bottom of INNER window is the same.\n");
2851 printf("Message:WF We assume that you DO NOT want to use an inner window and no states will be frozen.\n");
2852 }
2853
2854 Wannier_Inner_Window_Top=Wannier_Inner_Window_Bottom;
2855 }
2856
2857 s_vec[0]="Ang"; s_vec[1]="AU"; s_vec[2]="FRAC";
2858 i_vec[0]= 0; i_vec[1]= 1; i_vec[2]= 2;
2859 input_string2int("Wannier.Initial.Projectors.Unit",
2860 &Wannier_unit,3,s_vec,i_vec);
2861
2862
2863 if (fp=input_find("<Wannier.Initial.Projectors")) {
2864
2865 {
2866 int po,num_lines;
2867 char ctmp[200];
2868 double tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9;
2869 po = 0;
2870 num_lines = 0;
2871 do{
2872 fscanf(fp,"%s %lf %lf %lf %lf %lf %lf %lf %lf %lf",
2873 ctmp,&tmp1,&tmp2,&tmp3,&tmp4,&tmp5,&tmp6,&tmp7,&tmp8,&tmp9);
2874
2875
2876 if (strcmp(ctmp,"Wannier.Initial.Projectors>")==0){
2877 po = 1;
2878 }
2879 else{
2880 num_lines++;
2881 }
2882
2883 } while (po==0 && num_lines<10000);
2884
2885 if (strcmp(ctmp,"Wannier.Initial.Projectors>")!=0){
2886 if (myid==Host_ID){
2887 printf("Format error for Wannier.Initial.Projectors\n");
2888 }
2889 MPI_Finalize();
2890 exit(0);
2891 }
2892
2893 Wannier_Num_Kinds_Projectors = num_lines;
2894 }
2895 }
2896
2897 MPI_Barrier(mpi_comm_level1);
2898
2899 Allocate_Arrays(8);
2900
2901 /* allocation of temporal arrays */
2902 tmp_Wannier_Pro_SelMat=(int**)malloc(sizeof(int*)*Wannier_Num_Kinds_Projectors);
2903 for(i=0; i<Wannier_Num_Kinds_Projectors; i++){
2904 tmp_Wannier_Pro_SelMat[i]=(int*)malloc(sizeof(int)*Max_Num_WF_Projs);
2905 for(j=0;j<Max_Num_WF_Projs;j++){
2906 tmp_Wannier_Pro_SelMat[i][j]=-1;
2907 }
2908 }
2909 tmp_Wannier_Projector_Hybridize_Matrix=(double***)malloc(sizeof(double**)*Wannier_Num_Kinds_Projectors);
2910 for(i=0; i<Wannier_Num_Kinds_Projectors; i++){
2911 tmp_Wannier_Projector_Hybridize_Matrix[i]=(double**)malloc(sizeof(double*)*Max_Num_WF_Projs);
2912 for(j=0;j<Max_Num_WF_Projs;j++){
2913 tmp_Wannier_Projector_Hybridize_Matrix[i][j]=(double*)malloc(sizeof(double)*Max_Num_WF_Projs);
2914 }
2915 }
2916
2917 num_wannier_total_projectors = 0;
2918
2919 if (fp=input_find("<Wannier.Initial.Projectors")) {
2920
2921 {
2922 char ctmp[YOUSO10];
2923
2924 for (i=0; i<Wannier_Num_Kinds_Projectors; i++){
2925
2926 fscanf(fp,"%s %lf %lf %lf %lf %lf %lf %lf %lf %lf",
2927 ctmp,
2928 &Wannier_Pos[i][1],
2929 &Wannier_Pos[i][2],
2930 &Wannier_Pos[i][3],
2931
2932 &Wannier_Z_Direction[i][1],
2933 &Wannier_Z_Direction[i][2],
2934 &Wannier_Z_Direction[i][3],
2935 &Wannier_X_Direction[i][1],
2936 &Wannier_X_Direction[i][2],
2937 &Wannier_X_Direction[i][3]);
2938
2939 num_wannier_total_projectors += Analyze_Wannier_Projectors(i,ctmp,tmp_Wannier_Pro_SelMat,tmp_Wannier_Projector_Hybridize_Matrix);
2940
2941 if (Wannier_unit==2){
2942
2943 /* adjust position of Wannier funtion within 0 to 1 */
2944
2945 for (k=1; k<=3; k++){
2946
2947 /*adjust it to -0.5 to 0.5 */
2948
2949 if (2<=level_stdout) printf(" %10.5f -> ",Wannier_Pos[i][k]);
2950
2951 itmp = (int)Wannier_Pos[i][k];
2952 if (1.0<Wannier_Pos[i][k]){
2953 Wannier_Pos[i][k] = Wannier_Pos[i][k] - (double)itmp;
2954 }
2955 else if (Wannier_Pos[i][k]<0.0){
2956 Wannier_Pos[i][k] = Wannier_Pos[i][k] + (double)(abs(itmp)+1);
2957 }
2958
2959 if (2<=level_stdout) printf(" %10.5f ",Wannier_Pos[i][k]);
2960 }
2961
2962 /* calculate Cartesian coordiantes */
2963
2964 x = Wannier_Pos[i][1]*tv[1][1] + Wannier_Pos[i][2]*tv[2][1] + Wannier_Pos[i][3]*tv[3][1];
2965 y = Wannier_Pos[i][1]*tv[1][2] + Wannier_Pos[i][2]*tv[2][2] + Wannier_Pos[i][3]*tv[3][2];
2966 z = Wannier_Pos[i][1]*tv[1][3] + Wannier_Pos[i][2]*tv[2][3] + Wannier_Pos[i][3]*tv[3][3];
2967
2968 Wannier_Pos[i][1] = x;
2969 Wannier_Pos[i][2] = y;
2970 Wannier_Pos[i][3] = z;
2971
2972 }
2973
2974 if(Wannier_unit==0){
2975
2976 /* transfer the unit to AU since all the atomic site are experessed in AU */
2977
2978 Wannier_Pos[i][1]=Wannier_Pos[i][1]/BohrR;
2979 Wannier_Pos[i][2]=Wannier_Pos[i][2]/BohrR;
2980 Wannier_Pos[i][3]=Wannier_Pos[i][3]/BohrR;
2981 }
2982
2983 if (myid==Host_ID){
2984 printf("Projector Center(Ang.): (%10.5f, %10.5f, %10.5f)\n",
2985 Wannier_Pos[i][1]*BohrR,
2986 Wannier_Pos[i][2]*BohrR,
2987 Wannier_Pos[i][3]*BohrR);
2988 }
2989
2990 /* hmweng */
2991
2992 Get_Euler_Rotation_Angle(Wannier_Z_Direction[i][1],
2993 Wannier_Z_Direction[i][2],
2994 Wannier_Z_Direction[i][3],
2995 Wannier_X_Direction[i][1],
2996 Wannier_X_Direction[i][2],
2997 Wannier_X_Direction[i][3],
2998 &Wannier_Euler_Rotation_Angle[i][0],
2999 &Wannier_Euler_Rotation_Angle[i][1],
3000 &Wannier_Euler_Rotation_Angle[i][2]);
3001
3002 /* for each kind of Projector, find rotation matrix to the new coordinate system */
3003
3004 for(k=1; k<=3; k++){ /* L=0 is not included since s orbital is rotationally invariant */
3005
3006 double tmpRotMat[7][7];
3007 int jmp, jm;
3008
3009 if(Wannier_NumL_Pro[i][k]!=0){ /* if this l component is involved */
3010
3011 Get_Rotational_Matrix(Wannier_Euler_Rotation_Angle[i][0],
3012 Wannier_Euler_Rotation_Angle[i][1],
3013 Wannier_Euler_Rotation_Angle[i][2],
3014 k, tmpRotMat);
3015
3016 for(jmp=0;jmp<2*k+1;jmp++){
3017 for(jm=0;jm<2*k+1;jm++){
3018
3019 Wannier_RotMat_for_Real_Func[i][k][jmp][jm]=tmpRotMat[jmp][jm];
3020 }
3021 }
3022 }
3023 }
3024
3025 }
3026 }
3027 }
3028
3029 if (num_wannier_total_projectors!=Wannier_Func_Num && Wannier_Initial_Guess==1){
3030
3031 if (myid==Host_ID){
3032 printf("Error:WF The number (%d) of projectors is not same as that (%d) of Wannier functions\n",
3033 num_wannier_total_projectors,Wannier_Func_Num);
3034 }
3035
3036 MPI_Finalize();
3037 exit(0);
3038 }
3039
3040 { /* for each kind of projectors, find the selection matrix which can tell the position of
3041 projectors among the envolved basis set.
3042 For example,
3043 if projector(s) is(are) definded as dxy, which means the envolved basis set is local d
3044 orbitals, the projector is dxy which is in the third one arranged in the following order:
3045 dz2 --> 0, dx2-y2 --> 1, dxy -->2, dxz --> 3, dyz -->4.
3046 if projectors are defined as sp2, which means the envolved basis set are local s and p
3047 orbitals, the used local orbitals are s, px, py, this seletion matrix will have the values
3048 like:
3049 s --> 0, px -->1, py --> 2 (pz is in the basis set, but will not be selected in
3050 the final Amn matrix)
3051 */
3052
3053 Allocate_Arrays(9);
3054
3055 j=0;
3056
3057 for(i=0; i<Wannier_Num_Kinds_Projectors; i++){
3058
3059 for(k=0; k<Wannier_Num_Pro[i]; k++){
3060 Wannier_Select_Matrix[i][k] = tmp_Wannier_Pro_SelMat[i][k];
3061 }
3062
3063 for(k=0;k<Wannier_Num_Pro[i];k++){
3064 for(l=0;l<Wannier_Num_Pro[i];l++){
3065 Wannier_Projector_Hybridize_Matrix[i][k][l] = tmp_Wannier_Projector_Hybridize_Matrix[i][k][l];
3066 }
3067 }
3068
3069
3070 for(k=0;k<Wannier_Num_Pro[i];k++){
3071 Wannier_Guide[0][j]=Wannier_Pos[i][1];
3072 Wannier_Guide[1][j]=Wannier_Pos[i][2];
3073 Wannier_Guide[2][j]=Wannier_Pos[i][3];
3074 j++;
3075 }
3076
3077 } /*for each kind of projector*/
3078 }
3079
3080 i_vec2[0]=4;
3081 i_vec2[1]=4;
3082 i_vec2[2]=4;
3083 input_intv("Wannier.Kgrid",3,i_vec,i_vec2);
3084 Wannier_grid1 = i_vec[0];
3085 Wannier_grid2 = i_vec[1];
3086 Wannier_grid3 = i_vec[2];
3087
3088 input_int("Wannier.MaxShells", &Wannier_MaxShells,12);
3089 input_int("Wannier.Minimizing.Max.Steps", &Wannier_Minimizing_Max_Steps,200);
3090
3091 input_logical("Wannier.Interpolated.Bands", &Wannier_Draw_Int_Bands,0);
3092
3093 if (Wannier_Draw_Int_Bands){
3094 if (!Band_disp_switch){
3095
3096 if (myid==Host_ID){
3097 printf("Error:WF You are requiring for plotting the interpolated bands but without turn on key word Band.dispersion.\n");
3098 }
3099 MPI_Finalize();
3100 exit(0);
3101
3102 }
3103 }
3104
3105 input_logical("Wannier.Function.Plot", &Wannier_Draw_MLWF,0);
3106
3107 i_vec2[0]=0; i_vec2[1]=0; i_vec2[2]=0;
3108 input_intv("Wannier.Function.Plot.SuperCells",3,Wannier_Plot_SuperCells,i_vec2);
3109
3110 if ( Wannier_grid1<Wannier_Plot_SuperCells[0]
3111 || Wannier_grid2<Wannier_Plot_SuperCells[1]
3112 || Wannier_grid3<Wannier_Plot_SuperCells[2]){
3113
3114 if (myid==Host_ID){
3115 printf("Error:WF For each component, the following condition should be satisfied\n");
3116 printf(" Wannier.Function.Plot.SuperCells<=Wannier.Kgrid\n");
3117
3118 }
3119 MPI_Finalize();
3120 exit(0);
3121 }
3122
3123 input_double("Wannier.Dis.Mixing.Para",&Wannier_Dis_Mixing_Para,(double)0.5);
3124 input_double("Wannier.Dis.Conv.Criterion",&Wannier_Dis_Conv_Criterion,(double)1e-8);
3125 input_int("Wannier.Dis.SCF.Max.Steps",&Wannier_Dis_SCF_Max_Steps, 200);
3126 input_int("Wannier.Minimizing.Scheme",&Wannier_Min_Scheme, 0);
3127 input_double("Wannier.Minimizing.StepLength",&Wannier_Min_StepLength,(double)2.0);
3128 input_int("Wannier.Minimizing.Secant.Steps",&Wannier_Min_Secant_Steps, 5);
3129 input_double("Wannier.Minimizing.Secant.StepLength",&Wannier_Min_Secant_StepLength, (double)2.0);
3130 input_double("Wannier.Minimizing.Conv.Criterion",&Wannier_Min_Conv_Criterion, (double)1e-8);
3131 input_logical("Wannier.Output.kmesh",&Wannier_Output_kmesh, 0);
3132
3133 input_logical("Wannier.Readin.Overlap.Matrix",&Wannier_Readin_Overlap_Matrix,0);
3134
3135 if (Wannier_Readin_Overlap_Matrix==1){
3136 sprintf(file_check,"%s%s.mmn",filepath,filename);
3137 if ((fp_check = fopen(file_check,"r")) != NULL){
3138 fclose(fp_check);
3139 }
3140 else {
3141 if (myid==Host_ID){
3142 printf("\nError:WF not found %s\n",file_check);
3143 }
3144
3145 MPI_Finalize();
3146 exit(1);
3147 }
3148 }
3149
3150 if (Wannier_Readin_Overlap_Matrix==0)
3151 Wannier_Output_Overlap_Matrix = 1;
3152 else
3153 Wannier_Output_Overlap_Matrix = 0;
3154
3155 Wannier_Readin_Projection_Matrix = 0;
3156
3157 if (Wannier_Initial_Guess==1)
3158 Wannier_Output_Projection_Matrix = 1;
3159 else
3160 Wannier_Output_Projection_Matrix = 0;
3161
3162 if(Wannier_Readin_Overlap_Matrix==1) DFTSCF_loop=1;
3163
3164 /* freeing of temporal arrays */
3165
3166 for(i=0; i<Wannier_Num_Kinds_Projectors; i++){
3167 free(tmp_Wannier_Pro_SelMat[i]);
3168 }
3169 free(tmp_Wannier_Pro_SelMat);
3170
3171 for(i=0; i<Wannier_Num_Kinds_Projectors; i++){
3172 free(tmp_Wannier_Projector_Hybridize_Matrix[i]);
3173 }
3174 free(tmp_Wannier_Projector_Hybridize_Matrix);
3175
3176 } /* if(Wannier_Func_Calc) */
3177
3178 /****************************************************
3179 parameters for exact exchange methods by TOYODA
3180 ****************************************************/
3181
3182 /*---------- added by TOYODA 08/FEB/2010 */
3183 input_logical("exx.step1_skip", &g_exx_skip1, 0);
3184 input_logical("exx.step2_skip", &g_exx_skip2, 0);
3185 input_string("exx.cachedir", g_exx_cachedir, ".");
3186 input_int("exx.liberi.lmax", &g_exx_liberi_lmax, 20);
3187 input_int("exx.liberi.ngrid", &g_exx_liberi_ngrid, 1024);
3188 input_int("exx.liberi.ngl", &g_exx_liberi_ngl, 100);
3189 input_double("exx.rc_cut", &g_exx_rc_cut, -1.0);
3190 input_double("exx.w_scr", &g_exx_w_scr, -1.0);
3191 input_logical("exx.output_DM", &g_exx_switch_output_DM, 0);
3192 input_logical("exx.rhox", &g_exx_switch_rhox, 0);
3193 /*--------- until here */
3194
3195 /****************************************************
3196 parameters for cell optimizaiton
3197 ****************************************************/
3198
3199 /* scf.stress.tensor */
3200
3201 input_logical("scf.stress.tensor",&scf_stress_flag,0); /* default=off */
3202
3203 /****************************************************
3204 set flags for calculations of stress tensor
3205 scf_stress_flag: default=off
3206 MD_cellopt_flag: default=off
3207 ****************************************************/
3208
3209 /* In case of cell optimization, */
3210
3211 MD_cellopt_flag = 0;
3212
3213 if ( MD_switch==17
3214 || MD_switch==18
3215 )
3216 {
3217 scf_stress_flag = 1;
3218 MD_cellopt_flag = 1;
3219 }
3220
3221 /* check the compatibility with other functionlities */
3222
3223 if (scf_stress_flag==1 && Cnt_switch==1){
3224 if (myid==Host_ID){
3225 printf("Orbital optimization is not compatible with stress calculation.\n");
3226 }
3227 MPI_Finalize();
3228 exit(0);
3229 }
3230
3231 if (scf_stress_flag==1 && Hub_U_switch==1 && Hub_U_occupation==1){
3232 if (myid==Host_ID){
3233 printf("\n'full' local propjector in DFT+U is not compatible with stress calculation.\n");
3234 }
3235 MPI_Finalize();
3236 exit(0);
3237 }
3238
3239 if (scf_stress_flag==1 && SpinP_switch==3){
3240 if (myid==Host_ID){
3241 printf("\n'Non-collinear calculation is not compatible with stress calculation.\n");
3242 }
3243 MPI_Finalize();
3244 exit(0);
3245 }
3246
3247 /****************************************************
3248 input_close
3249 ****************************************************/
3250
3251 input_close();
3252
3253 if (po>0 || input_errorCount()>0) {
3254 printf("errors in the inputfile\n");
3255 MPI_Finalize();
3256 exit(1);
3257 }
3258
3259 /****************************************************
3260 adjustment of atomic position
3261 Atomic positions are not adjusted in Ver. 3.6
3262 ****************************************************/
3263
3264 /* if (Solver!=4) Set_In_First_Cell(); */
3265
3266 /****************************************************
3267 Gxyz -> His_Gxyz[0]
3268 ****************************************************/
3269
3270 k = 0;
3271 for (i=1; i<=atomnum; i++){
3272 His_Gxyz[0][k] = Gxyz[i][1]; k++;
3273 His_Gxyz[0][k] = Gxyz[i][2]; k++;
3274 His_Gxyz[0][k] = Gxyz[i][3]; k++;
3275 }
3276
3277 /****************************************************
3278 generate Monkhorst-Pack k-points
3279 ****************************************************/
3280
3281 if (Solver==3 && way_of_kpoint==2){
3282
3283 Generating_MP_Special_Kpt(/* input */
3284 atomnum, SpeciesNum,
3285 tv, Gxyz,
3286 InitN_USpin, InitN_DSpin,
3287 Criterion_MP_Special_Kpt,
3288 SpinP_switch,
3289 WhatSpecies,
3290 Kspace_grid1, Kspace_grid2, Kspace_grid3
3291 /* implicit output
3292 num_non_eq_kpt,
3293 NE_KGrids1, NE_KGrids2, NE_KGrids3,
3294 NE_T_k_op */ );
3295 }
3296
3297 /**************************************************************
3298 Effective Screening Medium (ESM) method Calculation
3299 added by T.Ohwaki
3300 ***************************************************************/
3301
3302 #if 0
3303
3304 if (ESM_switch!=0){
3305
3306 double xb0,xb1;
3307
3308 xb0 = 0.0;
3309 xb1 = tv[1][1];
3310
3311 /* check whether the wall position does not overlap with atoms. */
3312
3313 for (i=1; i<=atomnum; i++){
3314
3315 if ( ESM_wall_switch==1){
3316 if ( Gxyz[i][1]>(xb1-ESM_wall_position)){
3317
3318 if (myid==Host_ID){
3319 printf("<ESM:Warning>\n");
3320 printf("The coordinate of a-axis for atom %d = %16.9f \n",i,Gxyz[i][1]);
3321 printf("The coordinate of the wall position = %16.9f \n",xb1-ESM_wall_position);
3322 printf("The axis of the unit cell is too short.\n");
3323 printf("The wall position overlaps with atoms.\n");
3324 printf("Please increase the length of the a-axis.\n");
3325 }
3326
3327 MPI_Finalize();
3328 exit(1);
3329 }
3330 }
3331
3332 /* check the size of the upper vaccum */
3333
3334 if ( Gxyz[i][1]>(xb1-ESM_buffer_range)){
3335
3336 if (myid==Host_ID){
3337 printf("<ESM:Warning>\n");
3338 printf("The coordinate of a-axis for atom %d = %16.9f \n",i,Gxyz[i][1]);
3339 printf("The coordinate of the upper vacuum edge = %16.9f \n",xb1-ESM_buffer_range);
3340 printf("The region of the upper vaccum along the a-axis is too short.\n");
3341 printf("Please increase the length of the a-axis.\n");
3342 }
3343
3344 MPI_Finalize();
3345 exit(1);
3346 }
3347
3348 /* check the size of the lower vaccum */
3349
3350 if ( Gxyz[i][1]<(xb0+ESM_buffer_range)){
3351
3352 if (myid==Host_ID){
3353 printf("<ESM:Warning>\n");
3354 printf("The coordinate of a-axis for atom %d = %16.9f \n",i,Gxyz[i][1]);
3355 printf("The coordinate of the lower vacuum edge = %16.9f \n",xb0+ESM_buffer_range);
3356 printf("The region of the lower vaccum along the a-axis is too short.\n");
3357 printf("Please increase the length of the a-axis.\n");
3358 }
3359
3360 MPI_Finalize();
3361 exit(1);
3362 }
3363
3364 }
3365 }
3366
3367 #endif
3368
3369 /****************************************************
3370 print out to std
3371 ****************************************************/
3372
3373 if (myid==Host_ID && 0<level_stdout){
3374 printf("\n\n<Input_std> Your input file was normally read.\n");
3375 printf("<Input_std> The system includes %i species and %i atoms.\n",
3376 real_SpeciesNum,atomnum);
3377 }
3378
3379 }
3380
3381
3382
Remake_RestartFile(int numprocs_new,int numprocs_old,int N1,int N2,int N3)3383 void Remake_RestartFile(int numprocs_new, int numprocs_old, int N1, int N2, int N3)
3384 {
3385 int myid,ID,N2D,n2Ds_new,n2De_new;
3386 int n2Ds_old,n2De_old,num,spin,i,n,n3,m;
3387 int IDs_old,IDe_old,n2D;
3388 FILE *fp;
3389 char fileCD[YOUSO10];
3390 double *tmp_array0,*tmp_array1;
3391
3392 MPI_Comm_rank(mpi_comm_level1,&myid);
3393
3394 /*
3395 printf("myid=%2d numprocs_new=%2d numprocs_old=%2d\n",
3396 myid,numprocs_new,numprocs_old);
3397 */
3398
3399 N2D = N1*N2;
3400 n2Ds_new = (myid*N2D+numprocs_new-1)/numprocs_new;
3401 n2De_new = ((myid+1)*N2D+numprocs_new-1)/numprocs_new;
3402 IDs_old = n2Ds_new*numprocs_old/N2D;
3403 IDe_old = n2De_new*numprocs_old/N2D;
3404
3405 /* allocation of array */
3406 num = (n2De_new-n2Ds_new)*N3;
3407 tmp_array1 = (double*)malloc(sizeof(double)*num);
3408
3409 /*
3410 printf("n2Ds_new=%2d n2De_new=%2d\n",n2Ds_new,n2De_new);
3411 printf("IDs_old=%2d IDe_old=%2d\n",IDs_old,IDe_old);
3412 printf("num=%2d\n",num);
3413 */
3414
3415 /* restructure files */
3416
3417 for (spin=0; spin<=SpinP_switch; spin++){
3418 for (i=0; i<Extrapolated_Charge_History; i++){
3419
3420 n = 0;
3421 for (ID=IDs_old; ID<=IDe_old; ID++){
3422
3423 n2Ds_old = (ID*N2D+numprocs_old-1)/numprocs_old;
3424 n2De_old = ((ID+1)*N2D+numprocs_old-1)/numprocs_old;
3425
3426 /* allocation of array */
3427 num = (n2De_old-n2Ds_old)*N3;
3428 tmp_array0 = (double*)malloc(sizeof(double)*num);
3429
3430 /* set file name */
3431 sprintf(fileCD,"%s%s_rst/%s.crst%i_%i_%i",filepath,filename,filename,spin,ID,i);
3432
3433 /* read files */
3434 if ((fp = fopen(fileCD,"rb")) != NULL){
3435
3436 fread(tmp_array0,sizeof(double),num,fp);
3437
3438 fclose(fp);
3439
3440 for (n2D=n2Ds_old; n2D<n2De_old; n2D++){
3441 if (n2Ds_new<=n2D && n2D<n2De_new){
3442 for (n3=0; n3<N3; n3++){
3443 m = (n2D-n2Ds_old)*N3 + n3;
3444 tmp_array1[n] = tmp_array0[m];
3445 n++;
3446 }
3447 }
3448 }
3449 }
3450
3451 /* freeing of array */
3452 free(tmp_array0);
3453
3454 } /* ID */
3455
3456 /* set file name */
3457 sprintf(fileCD,"%s%s_rst/%s.crst%i_%i_%i",filepath,filename,filename,spin,myid,i);
3458
3459 if (n!=0){
3460 /* save data */
3461 if ((fp = fopen(fileCD,"wb")) != NULL){
3462 fwrite(tmp_array1,sizeof(double),n,fp);
3463 fclose(fp);
3464 }
3465 else{
3466 printf("Could not open a file %s\n",fileCD);
3467 }
3468 }
3469
3470 } /* i */
3471 } /* spin */
3472
3473 /* freeing of array */
3474 free(tmp_array1);
3475 }
3476
3477
3478
Set_In_First_Cell()3479 static void Set_In_First_Cell()
3480 {
3481 int i,Gc_AN;
3482 int itmp;
3483 double Cxyz[4];
3484 double tmp[4];
3485 double xc,yc,zc;
3486 double CellV;
3487
3488 /* calculate the reciprocal vectors */
3489
3490 Cross_Product(tv[2],tv[3],tmp);
3491 CellV = Dot_Product(tv[1],tmp);
3492 Cell_Volume = fabs(CellV);
3493
3494 Cross_Product(tv[2],tv[3],tmp);
3495 rtv[1][1] = 2.0*PI*tmp[1]/CellV;
3496 rtv[1][2] = 2.0*PI*tmp[2]/CellV;
3497 rtv[1][3] = 2.0*PI*tmp[3]/CellV;
3498
3499 Cross_Product(tv[3],tv[1],tmp);
3500 rtv[2][1] = 2.0*PI*tmp[1]/CellV;
3501 rtv[2][2] = 2.0*PI*tmp[2]/CellV;
3502 rtv[2][3] = 2.0*PI*tmp[3]/CellV;
3503
3504 Cross_Product(tv[1],tv[2],tmp);
3505 rtv[3][1] = 2.0*PI*tmp[1]/CellV;
3506 rtv[3][2] = 2.0*PI*tmp[2]/CellV;
3507 rtv[3][3] = 2.0*PI*tmp[3]/CellV;
3508
3509 /* find the center of coordinates */
3510
3511 xc = 0.0;
3512 yc = 0.0;
3513 zc = 0.0;
3514
3515 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
3516 xc += Gxyz[Gc_AN][1];
3517 yc += Gxyz[Gc_AN][2];
3518 zc += Gxyz[Gc_AN][3];
3519 }
3520
3521 xc = xc/(double)atomnum;
3522 yc = yc/(double)atomnum;
3523 zc = zc/(double)atomnum;
3524
3525 X_Center_Coordinate = xc;
3526 Y_Center_Coordinate = yc;
3527 Z_Center_Coordinate = zc;
3528
3529 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
3530
3531 Cxyz[1] = Gxyz[Gc_AN][1] - xc;
3532 Cxyz[2] = Gxyz[Gc_AN][2] - yc;
3533 Cxyz[3] = Gxyz[Gc_AN][3] - zc;
3534
3535 Cell_Gxyz[Gc_AN][1] = Dot_Product(Cxyz,rtv[1])*0.5/PI;
3536 Cell_Gxyz[Gc_AN][2] = Dot_Product(Cxyz,rtv[2])*0.5/PI;
3537 Cell_Gxyz[Gc_AN][3] = Dot_Product(Cxyz,rtv[3])*0.5/PI;
3538
3539 for (i=1; i<=3; i++){
3540 if (1.0<fabs(Cell_Gxyz[Gc_AN][i])){
3541 if (0.0<=Cell_Gxyz[Gc_AN][i]){
3542 itmp = (int)Cell_Gxyz[Gc_AN][i];
3543 Cell_Gxyz[Gc_AN][i] = Cell_Gxyz[Gc_AN][i] - (double)itmp;
3544 }
3545 else{
3546 itmp = abs((int)Cell_Gxyz[Gc_AN][i]) + 1;
3547 Cell_Gxyz[Gc_AN][i] = Cell_Gxyz[Gc_AN][i] + (double)itmp;
3548 }
3549 }
3550 }
3551
3552 Gxyz[Gc_AN][1] = Cell_Gxyz[Gc_AN][1]*tv[1][1]
3553 + Cell_Gxyz[Gc_AN][2]*tv[2][1]
3554 + Cell_Gxyz[Gc_AN][3]*tv[3][1] + xc;
3555
3556 Gxyz[Gc_AN][2] = Cell_Gxyz[Gc_AN][1]*tv[1][2]
3557 + Cell_Gxyz[Gc_AN][2]*tv[2][2]
3558 + Cell_Gxyz[Gc_AN][3]*tv[3][2] + yc;
3559
3560 Gxyz[Gc_AN][3] = Cell_Gxyz[Gc_AN][1]*tv[1][3]
3561 + Cell_Gxyz[Gc_AN][2]*tv[2][3]
3562 + Cell_Gxyz[Gc_AN][3]*tv[3][3] + zc;
3563 }
3564
3565 }
3566
3567
3568
3569
3570
3571
Analyze_Wannier_Projectors(int p,char ctmp[YOUSO10],int ** tmp_Wannier_Pro_SelMat,double *** tmp_Wannier_Projector_Hybridize_Matrix)3572 int Analyze_Wannier_Projectors(int p, char ctmp[YOUSO10],
3573 int **tmp_Wannier_Pro_SelMat,
3574 double ***tmp_Wannier_Projector_Hybridize_Matrix)
3575 {
3576 /*******************************************
3577 list of supported projectors:
3578 s(1),
3579 px(1),py(1),pz(1),
3580 dz2(1),dx2-y2(1),dxy(1),dxz(1),dyz(1),
3581 fz3(1),fxz2(1),fyz2(1),fzx2(1),fxyz(1),fx3-3xy2(1),f3yx2-y3(1),
3582 sp (s+px)(2),
3583 sp2 (s+px+py)(3),
3584 sp3 (s+px+py+pz)(4),
3585 sp3dz2(5),sp3deg(6),
3586 p (all p)(3)
3587 d (all d)(5)
3588 f (all f)(7)
3589 where the number of parenthesis is
3590 the total number of projectors.
3591 The number should be less than Max_Num_WF_Projs.
3592 *******************************************/
3593
3594 int i,po,ip,num,L,j,k,po_ok;
3595 /* Here 100 means the all the posibilitty of combination of initial
3596 guess. Total number of initial guess orbital should smaller than 100. */
3597 int Num_Projectors[100];
3598 int NumL[100][4];
3599 /* The following array has dimension 20, that means the maximume number
3600 of projectors included in one orbital combination. sp3 means 4 orbitals.
3601 sp3deg means 6 orbitals */
3602 int SelMat[100][Max_Num_WF_Projs];
3603 double HybridMat[100][Max_Num_WF_Projs][Max_Num_WF_Projs];
3604 char Name_Wannier_Projectors[100][30];
3605 char c;
3606 int myid;
3607
3608 /* get MPI ID */
3609 MPI_Comm_rank(MPI_COMM_WORLD1,&myid);
3610
3611 /* analysis of projectors used for calculation of Wannier functions */
3612
3613 i = 0;
3614 po = 0;
3615
3616 while( ((c=ctmp[i])!='\0' || po<=1) && i<YOUSO10 ){
3617
3618 if (c=='-' && po==0){
3619 po = 1;
3620 Wannier_ProSpeName[p][i] = '\0';
3621 ip = i;
3622 }
3623 else if (po==0){
3624 Wannier_ProSpeName[p][i] = ctmp[i];
3625 }
3626 else if (po==1){
3627 Wannier_ProName[p][i-ip-1] = ctmp[i];
3628 }
3629 else if (c=='\n'){
3630 po = 2;
3631 Wannier_ProName[p][i-ip-1] = '\0';
3632 }
3633
3634 i++;
3635 }
3636
3637 if (myid==Host_ID){
3638 printf("Projector name %s, orbitals %s.\n",Wannier_ProSpeName[p],Wannier_ProName[p]);
3639 }
3640
3641 /* check whether Wannier_ProSpeName is found in SpeBasisName or not */
3642
3643 po = 0;
3644 for (i=0; i<SpeciesNum; i++){
3645
3646 /*
3647 if (myid==Host_ID) printf("i=%2d %s\n",i,SpeName[i]);
3648 */
3649
3650 if (strcmp(Wannier_ProSpeName[p],SpeName[i])==0){
3651 po = 1;
3652 WannierPro2SpeciesNum[p] = i;
3653
3654 /* check the multiple basis specification of the same angular moment */
3655
3656 po_ok = 0;
3657
3658 /* s1 case */
3659 if (Spe_Num_Basis[i][0]==1 &&
3660 Spe_Num_Basis[i][1]==0 &&
3661 Spe_Num_Basis[i][2]==0 &&
3662 Spe_Num_Basis[i][3]==0 &&
3663 Spe_Num_Basis[i][4]==0) po_ok = 1;
3664
3665 /* s1p1 case */
3666 if (Spe_Num_Basis[i][0]==1 &&
3667 Spe_Num_Basis[i][1]==1 &&
3668 Spe_Num_Basis[i][2]==0 &&
3669 Spe_Num_Basis[i][3]==0 &&
3670 Spe_Num_Basis[i][4]==0) po_ok = 2;
3671
3672 /* s1p1d1 case */
3673 if (Spe_Num_Basis[i][0]==1 &&
3674 Spe_Num_Basis[i][1]==1 &&
3675 Spe_Num_Basis[i][2]==1 &&
3676 Spe_Num_Basis[i][3]==0 &&
3677 Spe_Num_Basis[i][4]==0) po_ok = 3;
3678
3679 /* s1p1d1f1 case */
3680 if (Spe_Num_Basis[i][0]==1 &&
3681 Spe_Num_Basis[i][1]==1 &&
3682 Spe_Num_Basis[i][2]==1 &&
3683 Spe_Num_Basis[i][3]==1 &&
3684 Spe_Num_Basis[i][4]==0) po_ok = 4;
3685
3686 if (po_ok==0){
3687 if (myid==Host_ID){
3688 printf("Could not assign the multiple projector for the same angular part of '%s'\n",
3689 Wannier_ProSpeName[p]);
3690 printf("The following specification can be allowed:\n");
3691 printf(" s1, s1p1, s1p1d1, or s1p1d1f1\n");
3692 }
3693
3694 MPI_Finalize();
3695 exit(0);
3696 }
3697
3698 }
3699 }
3700
3701 if (po==0){
3702 if (myid==Host_ID){
3703 printf("could not find %s among species\n",Wannier_ProSpeName[p]);
3704 }
3705 MPI_Finalize();
3706 exit(0);
3707 }
3708
3709 /* find the number of projectors */
3710
3711 Num_Wannier_Template_Projectors = 24;
3712
3713 for (i=0; i<Num_Wannier_Template_Projectors; i++){
3714 for (L=0; L<=3; L++){
3715 NumL[i][L] = 0;
3716 }
3717 }
3718
3719 /* selection matrix and Hybridize Matrix*/
3720
3721 /* s */
3722 strcpy(Name_Wannier_Projectors[ 0], "s"); Num_Projectors[0] = 1; NumL[0][0]=1;
3723 SelMat[0][0]=0;
3724 HybridMat[0][0][0]=1.0;
3725
3726 /* px */
3727 strcpy(Name_Wannier_Projectors[ 1], "px"); Num_Projectors[1] = 1; NumL[1][1]=1;
3728 SelMat[1][0]=0;
3729 HybridMat[1][0][0]=1.0;
3730
3731 /* py */
3732 strcpy(Name_Wannier_Projectors[ 2], "py"); Num_Projectors[2] = 1; NumL[2][1]=1;
3733 SelMat[2][0]=1;
3734 HybridMat[2][0][0]=1.0;
3735
3736 /* pz */
3737 strcpy(Name_Wannier_Projectors[ 3], "pz"); Num_Projectors[3] = 1; NumL[3][1]=1;
3738 SelMat[3][0]=2;
3739 HybridMat[3][0][0]=1.0;
3740
3741 /* dz2 */
3742 strcpy(Name_Wannier_Projectors[ 4], "dz2"); Num_Projectors[4] = 1; NumL[4][2]=1;
3743 SelMat[4][0]=0;
3744 HybridMat[4][0][0]=1.0;
3745
3746 /* dx2-y2 */
3747 strcpy(Name_Wannier_Projectors[ 5], "dx2-y2"); Num_Projectors[5] = 1; NumL[5][2]=1;
3748 SelMat[5][0]=1;
3749 HybridMat[5][0][0]=1.0;
3750
3751 /* dxy */
3752 strcpy(Name_Wannier_Projectors[ 6], "dxy"); Num_Projectors[6] = 1; NumL[6][2]=1;
3753 SelMat[6][0]=2;
3754 HybridMat[6][0][0]=1.0;
3755
3756 /* dxz */
3757 strcpy(Name_Wannier_Projectors[ 7], "dxz"); Num_Projectors[7] = 1; NumL[7][2]=1;
3758 SelMat[7][0]=3;
3759 HybridMat[7][0][0]=1.0;
3760
3761 /* dyz */
3762 strcpy(Name_Wannier_Projectors[ 8], "dyz"); Num_Projectors[8] = 1; NumL[8][2]=1;
3763 SelMat[8][0]=4;
3764 HybridMat[8][0][0]=1.0;
3765
3766 /* fz3 */
3767 strcpy(Name_Wannier_Projectors[ 9], "fz3"); Num_Projectors[9] = 1; NumL[9][3]=1;
3768 SelMat[9][0]=0;
3769 HybridMat[9][0][0]=1.0;
3770
3771 /* fxz2 */
3772 strcpy(Name_Wannier_Projectors[ 10], "fxz2"); Num_Projectors[10] = 1; NumL[10][3]=1;
3773 SelMat[10][0]=1;
3774 HybridMat[10][0][0]=1.0;
3775
3776 /* fyz2 */
3777 strcpy(Name_Wannier_Projectors[ 11], "fyz2"); Num_Projectors[11] = 1; NumL[11][3]=1;
3778 SelMat[11][0]=2;
3779 HybridMat[11][0][0]=1.0;
3780
3781 /* fzx2 */
3782 strcpy(Name_Wannier_Projectors[ 12], "fzx2"); Num_Projectors[12] = 1; NumL[12][3]=1;
3783 SelMat[12][0]=3;
3784 HybridMat[12][0][0]=1.0;
3785
3786 /* fxyz */
3787 strcpy(Name_Wannier_Projectors[ 13], "fxyz"); Num_Projectors[13] = 1; NumL[13][3]=1;
3788 SelMat[13][0]=4;
3789 HybridMat[13][0][0]=1.0;
3790
3791 /* fx3-3xy2 */
3792 strcpy(Name_Wannier_Projectors[ 14], "fx3-3xy2"); Num_Projectors[14] = 1; NumL[14][3]=1;
3793 SelMat[14][0]=5;
3794 HybridMat[14][0][0]=1.0;
3795
3796 /* f3yx2-y3 */
3797 strcpy(Name_Wannier_Projectors[ 15], "f3yx2-y3"); Num_Projectors[15] = 1; NumL[15][3]=1;
3798 SelMat[15][0]=6;
3799 HybridMat[15][0][0]=1.0;
3800
3801 /* sp */
3802 strcpy(Name_Wannier_Projectors[ 16], "sp"); Num_Projectors[16] = 2; NumL[16][0]=1; NumL[16][1]=1;
3803 SelMat[16][0]=0; /* this sp means the s and px hybirdization, therefore s and px should be selected. */
3804 SelMat[16][1]=1;
3805 HybridMat[16][0][0]=1.0/sqrt(2.0); HybridMat[16][0][1]=1.0/sqrt(2.0);
3806 HybridMat[16][1][0]=1.0/sqrt(2.0); HybridMat[16][1][1]=-1.0/sqrt(2.0);
3807
3808 /* sp2 */
3809 strcpy(Name_Wannier_Projectors[ 17], "sp2"); Num_Projectors[17] = 3; NumL[17][0]=1; NumL[17][1]=1;
3810 SelMat[17][0]=0; /* this sp2 means the s, px and py hybirdization, therefore s, px, py should be selected. */
3811 SelMat[17][1]=1;
3812 SelMat[17][2]=2;
3813
3814 HybridMat[17][0][0] = 1.0/sqrt(3.0);
3815 HybridMat[17][0][1] =-1.0/sqrt(6.0);
3816 HybridMat[17][0][2] = 1.0/sqrt(2.0);
3817
3818 HybridMat[17][1][0] = 1.0/sqrt(3.0);
3819 HybridMat[17][1][1] =-1.0/sqrt(6.0);
3820 HybridMat[17][1][2]= -1.0/sqrt(2.0);
3821
3822 HybridMat[17][2][0] = 1.0/sqrt(3.0);
3823 HybridMat[17][2][1] = 2.0/sqrt(6.0);
3824 HybridMat[17][2][2] = 0.0;
3825
3826 /* sp3 */
3827 strcpy(Name_Wannier_Projectors[ 18], "sp3"); Num_Projectors[18] = 4; NumL[18][0]=1; NumL[18][1]=1;
3828 SelMat[18][0]=0; /* this sp3 means the s, px, py and pz hybirdization, therefore s, px, py, pz should be selected. */
3829 SelMat[18][1]=1;
3830 SelMat[18][2]=2;
3831 SelMat[18][3]=3;
3832
3833 HybridMat[18][0][0] = 1.0/2.0;
3834 HybridMat[18][0][1] = 1.0/2.0;
3835 HybridMat[18][0][2] = 1.0/2.0;
3836 HybridMat[18][0][3] = 1.0/2.0;
3837
3838 HybridMat[18][1][0] = 1.0/2.0;
3839 HybridMat[18][1][1] = 1.0/2.0;
3840 HybridMat[18][1][2] =-1.0/2.0;
3841 HybridMat[18][1][3] =-1.0/2.0;
3842
3843 HybridMat[18][2][0] = 1.0/2.0;
3844 HybridMat[18][2][1] =-1.0/2.0;
3845 HybridMat[18][2][2] = 1.0/2.0;
3846 HybridMat[18][2][3] =-1.0/2.0;
3847
3848 HybridMat[18][3][0] = 1.0/2.0;
3849 HybridMat[18][3][1] =-1.0/2.0;
3850 HybridMat[18][3][2] =-1.0/2.0;
3851 HybridMat[18][3][3] = 1.0/2.0;
3852
3853 /* sp3dz2 */
3854 strcpy(Name_Wannier_Projectors[ 19], "sp3dz2");
3855 Num_Projectors[19] = 5; NumL[19][0]=1; NumL[19][1]=1; NumL[19][2]=1;
3856
3857 /* this sp3dz2 means the s, px, py and pz, and dz2 hybirdization,
3858 therefore s, px, py, pz and dz2 should be selected. */
3859
3860 SelMat[19][0]=0;
3861 SelMat[19][1]=1;
3862 SelMat[19][2]=2;
3863 SelMat[19][3]=3;
3864 SelMat[19][4]=4;
3865
3866 HybridMat[19][0][0] = 1.0/sqrt(3.0);
3867 HybridMat[19][0][1] =-1.0/sqrt(6.0);
3868 HybridMat[19][0][2] = 1.0/sqrt(2.0);
3869 HybridMat[19][0][3] = 0.0;
3870 HybridMat[19][0][4] = 0.0;
3871
3872 HybridMat[19][1][0] = 1.0/sqrt(3.0);
3873 HybridMat[19][1][1] =-1.0/sqrt(6.0);
3874 HybridMat[19][1][2] =-1.0/sqrt(2.0);
3875 HybridMat[19][1][3] = 0.0;
3876 HybridMat[19][1][4] = 0.0;
3877
3878 HybridMat[19][2][0] = 1.0/sqrt(3.0);
3879 HybridMat[19][2][1] = 2.0/sqrt(6.0);
3880 HybridMat[19][2][2] = 0.0;
3881 HybridMat[19][2][3] = 0.0;
3882 HybridMat[19][2][4] = 0.0;
3883
3884 HybridMat[19][3][0] = 0.0;
3885 HybridMat[19][3][1] = 0.0;
3886 HybridMat[19][3][2] = 0.0;
3887 HybridMat[19][3][3] = 1.0/sqrt(2.0);
3888 HybridMat[19][3][4] = 1.0/sqrt(2.0);
3889
3890 HybridMat[19][4][0] = 0.0;
3891 HybridMat[19][4][1] = 0.0;
3892 HybridMat[19][4][2] = 0.0;
3893 HybridMat[19][4][3] =-1.0/sqrt(2.0);
3894 HybridMat[19][4][4] = 1.0/sqrt(2.0);
3895
3896 /* sp3deg */
3897
3898 strcpy(Name_Wannier_Projectors[ 20], "sp3deg");
3899 Num_Projectors[20] = 6; NumL[20][0]=1; NumL[20][1]=1; NumL[20][2]=1;
3900
3901 /* this sp3deg means the s, px, py and pz, and dz2, dx2-y2 hybirdization,
3902 therefore these orbitals should be selected. */
3903
3904 SelMat[20][0]=0;
3905 SelMat[20][1]=1;
3906 SelMat[20][2]=2;
3907 SelMat[20][3]=3;
3908 SelMat[20][4]=4;
3909 SelMat[20][5]=5;
3910
3911 HybridMat[20][0][0] = 1.0/sqrt(6.0);
3912 HybridMat[20][0][1] =-1.0/sqrt(2.0);
3913 HybridMat[20][0][2] = 0.0;
3914 HybridMat[20][0][3] = 0.0;
3915 HybridMat[20][0][4] =-1.0/sqrt(12.0);
3916 HybridMat[20][0][5] = 0.5;
3917
3918 HybridMat[20][1][0] = 1.0/sqrt(6.0);
3919 HybridMat[20][1][1] = 1.0/sqrt(2.0);
3920 HybridMat[20][1][2] = 0.0;
3921 HybridMat[20][1][3] = 0.0;
3922 HybridMat[20][1][4] =-1.0/sqrt(12.0);
3923 HybridMat[20][1][5] = 0.5;
3924
3925 HybridMat[20][2][0] = 1.0/sqrt(6.0);
3926 HybridMat[20][2][1] = 0.0;
3927 HybridMat[20][2][2] =-1.0/sqrt(2.0);
3928 HybridMat[20][2][3] = 0.0;
3929 HybridMat[20][2][4] =-1.0/sqrt(12.0);
3930 HybridMat[20][2][5] =-0.5;
3931
3932 HybridMat[20][3][0] = 1.0/sqrt(6.0);
3933 HybridMat[20][3][1] = 0.0;
3934 HybridMat[20][3][2] = 1.0/sqrt(2.0);
3935 HybridMat[20][3][3] = 0.0;
3936 HybridMat[20][3][4] =-1.0/sqrt(12.0);
3937 HybridMat[20][3][5] =-0.5;
3938
3939 HybridMat[20][4][0] = 1.0/sqrt(6.0);
3940 HybridMat[20][4][1] = 0.0;
3941 HybridMat[20][4][2] = 0.0;
3942 HybridMat[20][4][3] =-1.0/sqrt(2.0);
3943 HybridMat[20][4][4] = 1.0/sqrt(3.0);
3944 HybridMat[20][4][5] = 0.0;
3945
3946 HybridMat[20][5][0] = 1.0/sqrt(6.0);
3947 HybridMat[20][5][1] = 0.0;
3948 HybridMat[20][5][2] = 0.0;
3949 HybridMat[20][5][3] = 1.0/sqrt(2.0);
3950 HybridMat[20][5][4] = 1.0/sqrt(3.0);
3951 HybridMat[20][5][5] = 0.0;
3952
3953 /* p, px, py and pz. no hybirdization. */
3954
3955 strcpy(Name_Wannier_Projectors[ 21], "p");
3956 Num_Projectors[21] = 3; NumL[21][0]=0; NumL[21][1]=1; NumL[21][2]=0;
3957
3958 SelMat[21][0]=0;
3959 SelMat[21][1]=1;
3960 SelMat[21][2]=2;
3961
3962 HybridMat[21][0][0]=1.0; HybridMat[21][0][1]=0.0; HybridMat[21][0][2]=0.0;
3963 HybridMat[21][1][0]=0.0; HybridMat[21][1][1]=1.0; HybridMat[21][1][2]=0.0;
3964 HybridMat[21][2][0]=0.0; HybridMat[21][2][1]=0.0; HybridMat[21][2][2]=1.0;
3965
3966 /* d, all five d orbitals.*/
3967
3968 strcpy(Name_Wannier_Projectors[ 22], "d");
3969 Num_Projectors[22] = 5; NumL[22][0]=0; NumL[22][1]=0; NumL[22][2]=1;
3970
3971 SelMat[22][0]=0;
3972 SelMat[22][1]=1;
3973 SelMat[22][2]=2;
3974 SelMat[22][3]=3;
3975 SelMat[22][4]=4;
3976
3977 for (i=0; i<=4; i++){
3978 for (j=0; j<=4; j++){
3979 HybridMat[22][i][j]=0.0;
3980 if (i==j) HybridMat[22][i][j]=1.0;
3981 }
3982 }
3983 /* f, all five f orbitals.*/
3984
3985 strcpy(Name_Wannier_Projectors[ 23], "f");
3986 Num_Projectors[23] = 7; NumL[23][0]=0; NumL[23][1]=0; NumL[23][2]=0; NumL[23][3]=1;
3987
3988 SelMat[23][0]=0;
3989 SelMat[23][1]=1;
3990 SelMat[23][2]=2;
3991 SelMat[23][3]=3;
3992 SelMat[23][4]=4;
3993 SelMat[23][5]=5;
3994 SelMat[23][6]=6;
3995
3996 for (i=0; i<=6; i++){
3997 for (j=0; j<=6; j++){
3998 HybridMat[23][i][j]=0.0;
3999 if (i==j) HybridMat[23][i][j]=1.0;
4000 }
4001 }
4002
4003 po = 0;
4004 i = 0;
4005 do {
4006 if (strcmp(Name_Wannier_Projectors[i],Wannier_ProName[p])==0){
4007 po = 1;
4008 num = Num_Projectors[i];
4009
4010 Wannier_ProName2Num[p]=i;
4011
4012 Wannier_Num_Pro[p] = Num_Projectors[i];
4013
4014 Wannier_NumL_Pro[p][0] = NumL[i][0];
4015 Wannier_NumL_Pro[p][1] = NumL[i][1];
4016 Wannier_NumL_Pro[p][2] = NumL[i][2];
4017 Wannier_NumL_Pro[p][3] = NumL[i][3];
4018
4019 /* for each kind of projectors, find the selection matrix which can tell the position of
4020 projectors among the envolved basis set.
4021 For example,
4022 if projector(s) is(are) definded as dxy, which means the envolved basis set is local d
4023 orbitals, the projector is dxy which is in the third one arranged in the following order:
4024 dz2 --> 0, dx2-y2 --> 1, dxy -->2, dxz --> 3, dyz -->4.
4025 if projectors are defined as sp2, which means the envolved basis set are local s and p
4026 orbitals, the used local orbitals are s, px, py, this seletion matrix will have the values
4027 like:
4028 s --> 0, px -->1, py --> 2 (pz is in the basis set, but will not be selected in
4029 the final Amn matrix)
4030 */
4031
4032 for(j=0;j<num;j++){
4033 tmp_Wannier_Pro_SelMat[p][j]=SelMat[i][j];
4034 for(k=0;k<num;k++){
4035 tmp_Wannier_Projector_Hybridize_Matrix[p][j][k]=HybridMat[i][j][k];
4036 }
4037 }
4038 }
4039 i++;
4040 } while(po==0 && i<Num_Wannier_Template_Projectors);
4041
4042 if (po==0){
4043
4044 if (myid==Host_ID){
4045 printf("Error:WF could not find %s among defined projectors\n",Wannier_ProName[p]);
4046 }
4047
4048 MPI_Finalize();
4049 exit(0);
4050 }
4051
4052 return num;
4053 }
4054
4055
4056
SpeciesString2int(int p)4057 void SpeciesString2int(int p)
4058 {
4059 int i,l,n,po,k,base;
4060 int tmp;
4061 int nlist[10];
4062 char c,cstr[YOUSO10*3];
4063 int numprocs,myid;
4064
4065 MPI_Comm_size(MPI_COMM_WORLD1,&numprocs);
4066 MPI_Comm_rank(MPI_COMM_WORLD1,&myid);
4067
4068 /* Get basis name */
4069
4070 sprintf(SpeBasisName[p],"");
4071
4072 i = 0;
4073 po = 0;
4074 while( ((c=SpeBasis[p][i])!='\0' || po==0) && i<YOUSO10 ){
4075 if (c=='-'){
4076 po = 1;
4077 SpeBasisName[p][i] = '\0';
4078 }
4079 if (po==0) SpeBasisName[p][i] = SpeBasis[p][i];
4080 i++;
4081 }
4082
4083 if (2<=level_stdout && myid==Host_ID){
4084 printf("<Input_std> SpeBasisName=%s\n",SpeBasisName[p]);
4085 }
4086
4087 /* Get basis type */
4088
4089 for (l=0; l<5; l++){
4090 Spe_Num_Basis[p][l] = 0;
4091 }
4092
4093 i = 0;
4094 po = 0;
4095 n = 0;
4096
4097 /* get informations of basis set */
4098
4099 while((c=SpeBasis[p][i])!='\0'){
4100
4101 if (po==1){
4102
4103 /* in case of orbital optimization, set Spe_Num_CBasis */
4104 if ( (c=='s' || c=='p' || c=='d' || c=='f' || c=='g')
4105 && n!=0 && Cnt_switch==1){
4106
4107 Spe_Num_CBasis[p][l] = 0;
4108 base = 1;
4109 for (k=(n-1); 0<=k; k--){
4110 Spe_Num_CBasis[p][l] += nlist[k]*base;
4111 base *= 10;
4112 }
4113 }
4114
4115 /* analysis of the string */
4116
4117 if (c=='s'){ l=0; n=0; }
4118 else if (c=='p'){ l=1; n=0; }
4119 else if (c=='d'){ l=2; n=0; }
4120 else if (c=='f'){ l=3; n=0; }
4121 else if (c=='g'){ l=4; n=0; }
4122 else{
4123
4124 /* no orbital optimization */
4125 if (Cnt_switch==0){
4126
4127 if (c=='>') {
4128 printf("Format error in Definition of Atomic Species\n");
4129 MPI_Finalize();
4130 exit(1);
4131 }
4132
4133 if (n==0){
4134 cstr[0] = c;
4135 cstr[1] = '\0';
4136 Spe_Num_Basis[p][l] = atoi(cstr); /* corrected by t.ohwaki */
4137 n++;
4138 }
4139 else if (n==1){
4140 cstr[0] = c;
4141 cstr[1] = '\0';
4142 tmp = Spe_Num_Basis[p][l];
4143 Spe_Num_Basis[p][l] = 10*tmp + atoi(cstr); /* corrected by t.ohwaki */
4144 n++;
4145 }
4146 else if (n==2){
4147 cstr[0] = c;
4148 cstr[1] = '\0';
4149 tmp = Spe_Num_Basis[p][l];
4150 Spe_Num_Basis[p][l] = 10*tmp + atoi(cstr); /* corrected by t.ohwaki */
4151 n++;
4152 }
4153 else {
4154 printf("Format error in Definition of Atomic Species\n");
4155 MPI_Finalize();
4156 exit(1);
4157 }
4158 }
4159
4160 /* orbital optimization */
4161 else if (Cnt_switch==1){
4162
4163 if (c=='>'){
4164
4165 /* set Spe_Num_Basis */
4166 Spe_Num_Basis[p][l] = 0;
4167 base = 1;
4168 for (k=(n-1); 0<=k; k--){
4169 Spe_Num_Basis[p][l] += nlist[k]*base;
4170 base *= 10;
4171 }
4172
4173 /* reset n */
4174 n = 0;
4175 }
4176
4177 else {
4178 cstr[0] = c;
4179 cstr[1] = '\0';
4180 nlist[n] = atoi(cstr);
4181 n++;
4182 }
4183 }
4184
4185 }
4186 }
4187
4188 if (SpeBasis[p][i]=='-') po = 1;
4189 i++;
4190 }
4191
4192 /* in case of orbital optimization, set Spe_Num_CBasis */
4193
4194 if ( n!=0 && Cnt_switch==1){
4195
4196 Spe_Num_CBasis[p][l] = 0;
4197 base = 1;
4198 for (k=(n-1); 0<=k; k--){
4199 Spe_Num_CBasis[p][l] += nlist[k]*base;
4200 base *= 10;
4201 }
4202 }
4203
4204 /* check the number of primitive and contracted basis functions */
4205
4206 for (l=0; l<5; l++){
4207 if (Spe_Num_Basis[p][l]!=0) Spe_MaxL_Basis[p] = l;
4208
4209 if (Spe_Num_Basis[p][l]<Spe_Num_CBasis[p][l]){
4210
4211 printf("# of contracted orbitals are larger than # of primitive oribitals\n");
4212 printf("Primitive=%3d Contracted=%3d\n",Spe_Num_Basis[p][l],Spe_Num_CBasis[p][l]);
4213 MPI_Finalize();
4214 exit(1);
4215 }
4216
4217 if (2<=level_stdout && Cnt_switch==0 && myid==Host_ID){
4218 printf("<Input_std> p=%2d l=%2d Primitive=%3d\n",
4219 p,l,Spe_Num_Basis[p][l]);
4220 }
4221 else if (Cnt_switch==1 && myid==Host_ID) { /* added by t.ohwaki */
4222 printf("<Input_std> p=%2d l=%2d Primitive=%3d Contracted=%3d\n",
4223 p,l,Spe_Num_Basis[p][l],Spe_Num_CBasis[p][l]);
4224 }
4225
4226 }
4227
4228 if (2<=level_stdout || Cnt_switch==1) printf("\n");
4229 }
4230
4231
4232
4233
4234
Species2int(char Species[YOUSO10])4235 int Species2int(char Species[YOUSO10])
4236 {
4237 int i,po;
4238
4239 i = 0;
4240 po = 0;
4241 while (i<SpeciesNum && po==0){
4242 if (SEQ(Species,SpeName[i])==1){
4243 po = 1;
4244 }
4245 if (po==0) i++;
4246 };
4247
4248 if (po==0){
4249 printf("Found an undefined species name %s\n",Species);
4250 printf("in Atoms.SpeciesAndCoordinates or Hubbard.U.values\n");
4251 printf("Please check your input file\n");
4252 MPI_Finalize();
4253 exit(1);
4254 }
4255
4256 return i;
4257 }
4258
4259
4260
OrbPol2int(char OrbPol[YOUSO10])4261 int OrbPol2int(char OrbPol[YOUSO10])
4262 {
4263 int i,po;
4264 char opns[3][YOUSO10]={"OFF","ON","EX"};
4265
4266 i = 0;
4267 po = 0;
4268
4269 ToCapital(OrbPol);
4270
4271 while (i<3 && po==0){
4272 if (SEQ(OrbPol,opns[i])==1){
4273 po = 1;
4274 }
4275 if (po==0) i++;
4276 };
4277
4278 if (po==0){
4279 printf("Invalid flag for LDA+U (Atoms.SpeciesAndCoordinates) %s\n",OrbPol);
4280 printf("Please check your input file\n");
4281 MPI_Finalize();
4282 exit(1);
4283 }
4284
4285 return i;
4286 }
4287
4288
4289
ToCapital(char * s)4290 char *ToCapital(char *s)
4291 {
4292 char *p;
4293 for (p=s; *p; p++) *p = toupper(*p);
4294 return (s);
4295 }
4296
4297
4298
4299
kpath_changeunit(double tv[4][4],double tv0[4][4],int Band_Nkpath,double *** Band_kpath)4300 void kpath_changeunit( double tv[4][4], double tv0[4][4], int Band_Nkpath,
4301 double ***Band_kpath )
4302 {
4303 /***********************************************************************
4304 k1 rtv0[0] + k2 rtv0[1] + k3 rtv0[2] = l rtv[0] + m rtv[1] + n rtv[2]
4305 rtv = reciptical vector of tv
4306 rtv0 = reciptical vector of tv0
4307 e.g. l is given by
4308 tv[0] ( k1 rtv0[0] + k2 rtv0[1] + k3 rtv0[2]) = l tv[0] rtv[0]
4309 ************************************************************************/
4310
4311 double tmp[4], CellV;
4312 double rtv[4][4],rtv0[4][4];
4313 int i,j;
4314 double r;
4315 int myid;
4316
4317 MPI_Comm_rank(MPI_COMM_WORLD1,&myid);
4318
4319
4320 Cross_Product(tv[2],tv[3],tmp);
4321 CellV = Dot_Product(tv[1],tmp);
4322
4323 Cross_Product(tv[2],tv[3],tmp);
4324 rtv[1][1] = 2.0*PI*tmp[1]/CellV;
4325 rtv[1][2] = 2.0*PI*tmp[2]/CellV;
4326 rtv[1][3] = 2.0*PI*tmp[3]/CellV;
4327
4328 Cross_Product(tv[3],tv[1],tmp);
4329 rtv[2][1] = 2.0*PI*tmp[1]/CellV;
4330 rtv[2][2] = 2.0*PI*tmp[2]/CellV;
4331 rtv[2][3] = 2.0*PI*tmp[3]/CellV;
4332
4333 Cross_Product(tv[1],tv[2],tmp);
4334 rtv[3][1] = 2.0*PI*tmp[1]/CellV;
4335 rtv[3][2] = 2.0*PI*tmp[2]/CellV;
4336 rtv[3][3] = 2.0*PI*tmp[3]/CellV;
4337
4338 Cross_Product(tv0[2],tv0[3],tmp);
4339 CellV = Dot_Product(tv0[1],tmp);
4340
4341 Cross_Product(tv0[2],tv0[3],tmp);
4342 rtv0[1][1] = 2.0*PI*tmp[1]/CellV;
4343 rtv0[1][2] = 2.0*PI*tmp[2]/CellV;
4344 rtv0[1][3] = 2.0*PI*tmp[3]/CellV;
4345
4346 Cross_Product(tv0[3],tv0[1],tmp);
4347 rtv0[2][1] = 2.0*PI*tmp[1]/CellV;
4348 rtv0[2][2] = 2.0*PI*tmp[2]/CellV;
4349 rtv0[2][3] = 2.0*PI*tmp[3]/CellV;
4350
4351 Cross_Product(tv0[1],tv0[2],tmp);
4352 rtv0[3][1] = 2.0*PI*tmp[1]/CellV;
4353 rtv0[3][2] = 2.0*PI*tmp[2]/CellV;
4354 rtv0[3][3] = 2.0*PI*tmp[3]/CellV;
4355
4356 if (myid==Host_ID){
4357 printf("kpath (converted)\n");
4358 }
4359
4360 for (i=1;i<=Band_Nkpath;i++) {
4361 for (j=1;j<=3;j++) tmp[j]=Band_kpath[i][1][j];
4362 for (j=1;j<=3;j++) {
4363 r = tmp[1]* Dot_Product(tv[j],rtv0[1])
4364 + tmp[2]* Dot_Product(tv[j],rtv0[2])
4365 + tmp[3]* Dot_Product(tv[j],rtv0[3]);
4366 Band_kpath[i][1][j] = r/ Dot_Product(tv[j],rtv[j]);
4367 }
4368 for (j=1;j<=3;j++) tmp[j]=Band_kpath[i][2][j];
4369 for (j=1;j<=3;j++) {
4370 r = tmp[1]* Dot_Product(tv[j],rtv0[1])
4371 + tmp[2]* Dot_Product(tv[j],rtv0[2])
4372 + tmp[3]* Dot_Product(tv[j],rtv0[3]);
4373 Band_kpath[i][2][j] = r/ Dot_Product(tv[j],rtv[j]);
4374 }
4375
4376 if (myid==Host_ID){
4377 printf("(%lf %lf %lf) (%lf %lf %lf)\n",
4378 Band_kpath[i][1][1],Band_kpath[i][1][2],Band_kpath[i][1][3],
4379 Band_kpath[i][2][1],Band_kpath[i][2][2],Band_kpath[i][2][3]);
4380 }
4381 }
4382
4383 }
4384
4385
kpoint_changeunit(double tv[4][4],double tv0[4][4],int MO_Nkpoint,double ** MO_kpoint)4386 void kpoint_changeunit(double tv[4][4],double tv0[4][4],int MO_Nkpoint,
4387 double **MO_kpoint)
4388 {
4389 /***********************************************************************
4390 k1 rtv0[0] + k2 rtv0[1] + k3 rtv0[2] = l rtv[0] + m rtv[1] + n rtv[2]
4391 rtv = reciptical vector of tv
4392 rtv0 = reciptical vector of tv0
4393 e.g. l is given by
4394 tv[0] ( k1 rtv0[0] + k2 rtv0[1] + k3 rtv0[2]) = l tv[0] rtv[0]
4395 ************************************************************************/
4396
4397 double tmp[4], CellV;
4398 double rtv[4][4],rtv0[4][4];
4399 int i,j;
4400 double r;
4401 int myid;
4402
4403 MPI_Comm_rank(MPI_COMM_WORLD1,&myid);
4404
4405 Cross_Product(tv[2],tv[3],tmp);
4406 CellV = Dot_Product(tv[1],tmp);
4407
4408 Cross_Product(tv[2],tv[3],tmp);
4409 rtv[1][1] = 2.0*PI*tmp[1]/CellV;
4410 rtv[1][2] = 2.0*PI*tmp[2]/CellV;
4411 rtv[1][3] = 2.0*PI*tmp[3]/CellV;
4412
4413 Cross_Product(tv[3],tv[1],tmp);
4414 rtv[2][1] = 2.0*PI*tmp[1]/CellV;
4415 rtv[2][2] = 2.0*PI*tmp[2]/CellV;
4416 rtv[2][3] = 2.0*PI*tmp[3]/CellV;
4417
4418 Cross_Product(tv[1],tv[2],tmp);
4419 rtv[3][1] = 2.0*PI*tmp[1]/CellV;
4420 rtv[3][2] = 2.0*PI*tmp[2]/CellV;
4421 rtv[3][3] = 2.0*PI*tmp[3]/CellV;
4422
4423 Cross_Product(tv0[2],tv0[3],tmp);
4424 CellV = Dot_Product(tv0[1],tmp);
4425
4426 Cross_Product(tv0[2],tv0[3],tmp);
4427 rtv0[1][1] = 2.0*PI*tmp[1]/CellV;
4428 rtv0[1][2] = 2.0*PI*tmp[2]/CellV;
4429 rtv0[1][3] = 2.0*PI*tmp[3]/CellV;
4430
4431 Cross_Product(tv0[3],tv0[1],tmp);
4432 rtv0[2][1] = 2.0*PI*tmp[1]/CellV;
4433 rtv0[2][2] = 2.0*PI*tmp[2]/CellV;
4434 rtv0[2][3] = 2.0*PI*tmp[3]/CellV;
4435
4436 Cross_Product(tv0[1],tv0[2],tmp);
4437 rtv0[3][1] = 2.0*PI*tmp[1]/CellV;
4438 rtv0[3][2] = 2.0*PI*tmp[2]/CellV;
4439 rtv0[3][3] = 2.0*PI*tmp[3]/CellV;
4440
4441 if (myid==Host_ID){
4442 printf("kpoint at which wave functions are calculated (converted)\n");
4443 }
4444
4445 for (i=0;i<MO_Nkpoint;i++){
4446 for (j=1;j<=3;j++) tmp[j]= MO_kpoint[i][j];
4447 for (j=1;j<=3;j++) {
4448 r = tmp[1]* Dot_Product(tv[j],rtv0[1])
4449 + tmp[2]* Dot_Product(tv[j],rtv0[2])
4450 + tmp[3]* Dot_Product(tv[j],rtv0[3]);
4451 MO_kpoint[i][j] = r/ Dot_Product(tv[j],rtv[j]);
4452 }
4453 if (myid==Host_ID){
4454 printf("%lf %lf %lf\n",MO_kpoint[i][1],MO_kpoint[i][2],MO_kpoint[i][3]);
4455 }
4456 }
4457 }
4458
4459
4460
4461
4462 /* *** calculate an unit cell of a cluster, ***
4463 * assuming that unit of Gxyz and Rc is A.U.
4464 * cell size = (max[ xyz-Rc ] - min[ xyz+Rc ])*1.01
4465 */
4466
4467
Set_Cluster_UnitCell(double tv[4][4],int unitflag)4468 void Set_Cluster_UnitCell(double tv[4][4], int unitflag)
4469 {
4470 /*
4471 * Species: int SpeciesNum, char SpeName[], char SpeBasis[]
4472 *
4473 * Coordinate: int WhatSpecies[]; double Gxyz[][]
4474 *
4475 * unitflag=0 (Ang.) unitflag=1 (a.u.), used only to print them
4476 *
4477 * tv[][] is always in a.u.
4478 */
4479 FILE *fp;
4480 int i,id,spe,myid;
4481 double *sperc;
4482 double min[4],max[4],gmin[4],gmax[4];
4483 char FN_PAO[YOUSO10];
4484 char ExtPAO[YOUSO10] = ".pao";
4485 char DirPAO[YOUSO10];
4486
4487 double margin=1.10;
4488 double unit;
4489 char *unitstr[2]={"Ang.","a.u."};
4490
4491 MPI_Comm_rank(MPI_COMM_WORLD1,&myid);
4492
4493 /* set DirPAO */
4494
4495 sprintf(DirPAO,"%s/PAO/","/usr/local/share/openmx/DFT_DATA13");
4496
4497 unit=1.0;
4498 if (unitflag==0) unit=BohrR;
4499
4500 sperc=(double*)malloc(sizeof(double)*SpeciesNum);
4501
4502 for (spe=0; spe<SpeciesNum; spe++){
4503 fnjoint2(DirPAO,SpeBasisName[spe],ExtPAO,FN_PAO);
4504 if ((fp = fopen(FN_PAO,"r")) != NULL){
4505 input_open(FN_PAO);
4506 input_double("radial.cutoff.pao",&sperc[spe],(double)0.0);
4507 input_close();
4508 fclose(fp);
4509 }
4510 else {
4511
4512 if (myid==Host_ID){
4513 printf("Set_Cluster_UnitCell: can not open %s\n", FN_PAO);
4514 }
4515
4516 MPI_Finalize();
4517 exit(10);
4518 }
4519 }
4520
4521 if (level_stdout>=2) {
4522 for (spe=0;spe<SpeciesNum;spe++) {
4523 printf("<Set_Cluster_UnitCell> %d %s rc=%lf\n",spe,SpeName[spe],sperc[spe]);
4524 }
4525 }
4526
4527
4528 if (level_stdout>=2) {
4529 printf("<Set_Cluster_UnitCell> x y z Rc\n");
4530 }
4531
4532 for (i=1;i<=3;i++) {
4533 gmax[i]=gmin[i]=Gxyz[1][i];
4534 }
4535
4536 for (i=1;i<=atomnum;i++) {
4537 id=WhatSpecies[i];
4538 /* printf("%d %d %lf\n",i,id,sperc[id]); */
4539 if (level_stdout>=2) {
4540 printf("<Set_Cluster_UnitCell> %lf %lf %lf %lf\n",
4541 Gxyz[i][1],Gxyz[i][2],Gxyz[i][3],sperc[id]);
4542 }
4543 min[1]=Gxyz[i][1]-sperc[id];
4544 min[2]=Gxyz[i][2]-sperc[id];
4545 min[3]=Gxyz[i][3]-sperc[id];
4546 max[1]=Gxyz[i][1]+sperc[id];
4547 max[2]=Gxyz[i][2]+sperc[id];
4548 max[3]=Gxyz[i][3]+sperc[id];
4549 if (min[1]<gmin[1]) gmin[1]=min[1];
4550 if (min[2]<gmin[2]) gmin[2]=min[2];
4551 if (min[3]<gmin[3]) gmin[3]=min[3];
4552 if (max[1]>gmax[1]) gmax[1]=max[1];
4553 if (max[2]>gmax[2]) gmax[2]=max[2];
4554 if (max[3]>gmax[3]) gmax[3]=max[3];
4555 }
4556
4557 /* initialize */
4558 for (id=1;id<=3;id++){
4559 for(i=1;i<=3;i++) {
4560 tv[id][i]=0.0;
4561 }
4562 }
4563
4564 tv[1][1]=(gmax[1]-gmin[1])*margin;
4565 tv[2][2]=(gmax[2]-gmin[2])*margin;
4566 tv[3][3]=(gmax[3]-gmin[3])*margin;
4567
4568 if (level_stdout>=2) {
4569 printf("<Set_Cluster_UnitCell> to determine the unit cell, min and max includes effects of Rc\n");
4570 for (i=1;i<=3;i++)
4571 printf("<Set_Cluster_UnitCell> axis=%d min,max=%lf %lf\n", i, gmin[i]*unit,gmax[i]*unit);
4572 }
4573
4574 if (myid==Host_ID && 0<level_stdout){
4575
4576 printf("<Set_Cluster_UnitCell> automatically determied UnitCell(%s)\n<Set_Cluster_UnitCell> from atomic positions and Rc of PAOs (margin= %5.2lf%%)\n",unitstr[unitflag],(margin-1.0)*100.0);
4577 for(i=1;i<=3;i++) {
4578 printf("<Set_Cluster_UnitCell> %lf %lf %lf\n",
4579 tv[i][1]*unit,tv[i][2]*unit, tv[i][3]*unit);
4580 }
4581 printf("\n");
4582 }
4583
4584 free(sperc);
4585 }
4586
4587
4588
divisible_cheker(int N)4589 int divisible_cheker(int N)
4590 {
4591 /************************
4592 return 0; non-divisible
4593 return 1; divisible
4594 ************************/
4595
4596 int i,po;
4597
4598 if (N!=1){
4599 po = 1;
4600 for (i=0; i<NfundamentalNum; i++){
4601 if ( N!=1 && (N % fundamentalNum[i])==0 ){
4602 po = 0;
4603 N = N/fundamentalNum[i];
4604 }
4605 }
4606 }
4607 else{
4608 po = 0;
4609 }
4610
4611 if (po==0 && N!=1){
4612 divisible_cheker(N);
4613 }
4614
4615 if (po==0) return 1;
4616 else return 0;
4617 }
4618
4619
4620
4621
4622
4623
4624
4625
4626
4627
4628
4629 /* hmweng */
Get_Euler_Rotation_Angle(double zx,double zy,double zz,double xx,double xy,double xz,double * alpha_r,double * beta_r,double * gamma_r)4630 void Get_Euler_Rotation_Angle(
4631 double zx, double zy, double zz,
4632 double xx, double xy, double xz,
4633 double *alpha_r, double *beta_r, double *gamma_r)
4634 {
4635 double norm, coszx, yx,yy,yz, alpha, beta, gamma, tmp1;
4636 int numprocs,myid;
4637
4638 MPI_Comm_size(MPI_COMM_WORLD1,&numprocs);
4639 MPI_Comm_rank(MPI_COMM_WORLD1,&myid);
4640
4641 *alpha_r=0.0;
4642 *beta_r=0.0;
4643 *gamma_r=0.0;
4644 alpha=0.0;
4645 beta=0.0;
4646 gamma=0.0;
4647
4648 /* check orthogonality of z-axis and x-axis and normalise them. */
4649
4650 norm=sqrt(zx*zx+zy*zy+zz*zz);
4651 zx=zx/norm;
4652 zy=zy/norm;
4653 zz=zz/norm;
4654
4655 norm=sqrt(xx*xx+xy*xy+xz*xz);
4656 xx=xx/norm;
4657 xy=xy/norm;
4658 xz=xz/norm;
4659
4660 coszx=(zx*xx+zy*xy+zz*xz);
4661
4662 if(fabs(coszx)>1e-6){ /* not orthogonal */
4663
4664 if (myid==Host_ID){
4665 printf("Error:WF z-axis and x-axis are not orthogonal, please check it and try again.\n");
4666 }
4667
4668 MPI_Finalize();
4669 exit(0);
4670 }
4671 /* orthogonal */
4672 /* then get y-axis = z-axis cross product x-axis to make a right-hand coordinate. */
4673
4674 yx=zy*xz-zz*xy;
4675 yy=zz*xx-zx*xz;
4676 yz=zx*xy-zy*xx;
4677
4678 norm=sqrt(yx*yx+yy*yy+yz*yz);
4679 yx=yx/norm;
4680 yy=yy/norm;
4681 yz=yz/norm;
4682
4683 if(fabs(fabs(zz)-1.0)<1e-6){ /* new z-axis is along old z-axis */
4684 if(zz>0){ /* along the positive direciton */
4685 beta=0.0;
4686 }else{
4687 beta=PI; /* have a 180-degree rotation */
4688 }
4689 /* in this case, alpha and gamma are the same, rotating around the same z-axis.*/
4690 gamma=0.0;
4691 /* alpha is determined by x-axis */
4692 alpha=asin(xy/cos(beta)); /* sin(alpha)*cos(beta)=xy */
4693 if(xx/zz<0){ /* cos(alpha)*cos(beta)=xx */
4694 alpha=PI-alpha;
4695 }
4696 if(alpha<0){
4697 alpha=2.0*PI+alpha;
4698 }
4699 }
4700
4701 else{
4702
4703 beta=acos(zz); /* beta is always in [0, PI] */
4704 tmp1=sqrt(zx*zx+zy*zy); /* sin(beta) */
4705 alpha=asin(zy/tmp1); /* zy=sin(beta)*sin(alpha) asin() gives vale between [-PI/2, PI/2] */
4706 if(zx<0){ /* zx=sin(beta)*cos(alpha) and sin(beta)>=0.0, this means cos(alpha)<0.0 */
4707 alpha=PI-alpha;
4708 }
4709 if(alpha<0.0){
4710 alpha=2.0*PI+alpha; /* make alpha in [0,2PI] */
4711 }
4712
4713 /* determin gamma now */
4714
4715 if(fabs(fabs(-xz/sin(beta))-1.0)<1e-5){
4716 if(-xz/sin(beta)<0.0){
4717 gamma=PI;
4718 }else{
4719 gamma=0.0;
4720 }
4721 }else{
4722 gamma=acos(-xz/sin(beta)); /* xz=-cos(gamma)*sin(beta). acos() gieve a value between [0, PI] */
4723 }
4724 /* we need sin(gamma) to finally determin gamma. xx=cos(alpha)*cos(beta)*cos(gamma)-sin(alpha)*sin(gamma) */
4725 tmp1=-(xx-(-xz/sin(beta))*zz*cos(alpha))/sin(alpha); /* tmp1 is sin(gamma) */
4726 if(tmp1<0){
4727 gamma=2*PI-gamma;
4728 }
4729
4730 } /* cos(beta)!=1.0 or -1.0 */
4731
4732 *alpha_r=alpha;
4733 *beta_r=beta;
4734 *gamma_r=gamma;
4735
4736 if (myid==Host_ID){
4737 printf("z-axis %10.5f %10.5f %10.5f, x-axis %10.5f %10.5f %10.5f\n",zx,zy,zz,xx,xy,xz);
4738 printf("y-axis %10.5f %10.5f %10.5f\n",yx,yy,yz);
4739 printf("Euler Angles are %10.5f, %10.5f, %10.5f.(in degree)\n",
4740 alpha/PI*180.0,beta/PI*180.0,gamma/PI*180.0);
4741 printf("Euler Angles are %10.5f, %10.5f, %10.5f.(in rad)\n",alpha,beta,gamma);
4742 }
4743
4744 } /* end of Get_Euler_Rotation_Angle */
4745
4746
4747 /* hmweng */
4748
Calc_Factorial(int arg)4749 int Calc_Factorial(int arg)
4750 {
4751
4752 /* calculate Factorial of arg. (arg)! */
4753 int result,n;
4754
4755 result = 1;
4756
4757 if(arg<0){
4758 printf("Error. For Calc_Factorial, positive integer is needed!\n");
4759 exit(0);
4760 }
4761
4762 if(arg==0 || arg==1){
4763 result=1;
4764 }
4765
4766 else{
4767 for(n=arg;n>0;n--){
4768 result=n*result;
4769 }
4770 }
4771
4772 return result;
4773 }
4774
4775
4776
4777 /* hmweng */
Get_Rotational_Matrix(double alpha,double beta,double gamma,int L,double tmpRotMat[7][7])4778 void Get_Rotational_Matrix(double alpha, double beta, double gamma, int L, double tmpRotMat[7][7])
4779 {
4780 int j,i,k;
4781 int m, mp;
4782 int jm, jmp, maxk;
4783 double fac1,tmp1,tmp2,fac2,sumk,dj[2*L+1][2*L+1],tmp3, sumr,sumi;
4784 dcomplex Dlm[7][7];
4785 dcomplex RotMat_for_Real_Func[7][7];
4786 dcomplex Umat[7][7], Utmp[7][7], Umat_inv[7][7];
4787 int myid;
4788
4789 /* get MPI ID */
4790 MPI_Comm_rank(MPI_COMM_WORLD1,&myid);
4791
4792 j=L;
4793
4794 for(jmp=0;jmp<7;jmp++){
4795 for(jm=0;jm<7;jm++){
4796 Dlm[jmp][jm].r=0.0;
4797 Dlm[jmp][jm].i=0.0;
4798 RotMat_for_Real_Func[jmp][jm].r=0.0;
4799 RotMat_for_Real_Func[jmp][jm].i=0.0;
4800 Umat[jmp][jm].r=0.0;
4801 Umat[jmp][jm].i=0.0;
4802 Umat_inv[jmp][jm].r=0.0;
4803 Umat_inv[jmp][jm].i=0.0;
4804 tmpRotMat[jmp][jm]=0.0;
4805 if(jmp==jm){
4806 Dlm[jmp][jm].r=1.0;
4807 RotMat_for_Real_Func[jmp][jm].r=1.0;
4808 }
4809 }
4810 }
4811
4812
4813 for(jmp=0;jmp<2*j+1;jmp++){
4814 for(jm=jmp;jm<2*j+1;jm++){
4815 mp=j-jmp; /* mp goes from j to -j */
4816 m=j-jm; /* m goes from mp to -j to satisfy condition mp>=m */
4817
4818 tmp1=(double)(Calc_Factorial(j+m)*Calc_Factorial(j-m)*Calc_Factorial(j+mp)*Calc_Factorial(j-mp));
4819
4820 fac1=sqrt(tmp1);
4821 if((j-mp)>(j+m)){
4822 maxk=j+m;
4823 }else{
4824 maxk=j-mp;
4825 }
4826
4827 sumk=0.0;
4828 for(k=0;k<=maxk;k++){
4829 tmp1=(double)(Calc_Factorial(j-mp-k)*Calc_Factorial(j+m-k)*Calc_Factorial(k+mp-m)*Calc_Factorial(k));
4830
4831 if(fabs(cos(beta/2.0))<1e-8 && 2*j+m-mp-2*k==0){
4832 tmp2=1.0;
4833 }else{
4834 tmp2=exp(log(fabs(cos(beta/2.0)))*(double)(2*j+m-mp-2*k));
4835 }
4836 if(fabs(sin(beta/2.0))<1e-8 && 2*k+mp-m==0){
4837 tmp3=1.0;
4838 }else{
4839 tmp3=exp(log(fabs(sin(beta/2.0)))*(double)(2*k+mp-m));
4840 }
4841 fac2=tmp2*tmp3/tmp1;
4842 if((k+mp-m)%2==0){
4843 sumk=sumk+fac2;
4844 }else{
4845 sumk=sumk-fac2;
4846 }
4847 }/* sum over all k */
4848 dj[jmp][jm]=sumk*fac1;
4849 tmp1= cos((double)mp*alpha+(double)m*gamma);
4850 tmp2=-sin((double)m*gamma+(double)mp*alpha);
4851 Dlm[jmp][jm].r=dj[jmp][jm]*tmp1;
4852 Dlm[jmp][jm].i=dj[jmp][jm]*tmp2;
4853 } /* m */
4854 }/* m' */
4855
4856 /* for those when m'<m */
4857 for(jmp=0;jmp<2*j+1;jmp++){
4858 for(jm=0;jm<jmp;jm++){
4859 mp=j-jmp; /* mp runs from j-1 to -j+1 */
4860 m=j-jm; /* m runs from mp+1 to j to satisfy mp<m */
4861 if((mp-m)%2==0){
4862 dj[jmp][jm]=dj[jm][jmp];
4863 }else{
4864 dj[jmp][jm]=-dj[jm][jmp];
4865 }
4866 tmp1= cos((double)mp*alpha)*cos((double)m*gamma)-sin((double)mp*alpha)*sin((double)m*gamma);
4867 tmp2=-cos((double)mp*alpha)*sin((double)m*gamma)-sin((double)mp*alpha)*cos((double)m*gamma);
4868 Dlm[jmp][jm].r=dj[jmp][jm]*tmp1;
4869 Dlm[jmp][jm].i=dj[jmp][jm]*tmp2;
4870 }
4871 }
4872
4873 if (0){
4874 for(jmp=0;jmp<2*j+1;jmp++){
4875 for(jm=0;jm<2*j+1;jm++){
4876 printf("(%10.6f %10.6f) ",Dlm[jmp][jm].r,Dlm[jmp][jm].i);
4877 }
4878 printf("\n");
4879 }
4880 if(j==1){
4881 printf("Compare with those in Rose's books:\n");
4882 printf("(%10.6f %10.6f) ",(cos(alpha)*cos(gamma)-sin(alpha)*sin(gamma))*(1.0+cos(beta))/2.0,
4883 (-cos(alpha)*sin(gamma)-sin(alpha)*cos(gamma))*(1.0+cos(beta))/2.0);
4884
4885 printf("(%10.6f %10.6f) ",-cos(alpha)*sin(beta)/sqrt(2.0),sin(alpha)*sin(beta)/sqrt(2.0));
4886
4887 printf("(%10.6f %10.6f) ",(cos(alpha)*cos(gamma)+sin(alpha)*sin(gamma))*(1.0-cos(beta))/2.0,
4888 (cos(alpha)*sin(gamma)-sin(alpha)*cos(gamma))*(1.0-cos(beta))/2.0);
4889 printf("\n");
4890
4891 printf("(%10.6f %10.6f) ",cos(gamma)*sin(beta)/sqrt(2.0),-sin(gamma)*sin(beta)/sqrt(2.0));
4892
4893 printf("(%10.6f %10.6f) ",cos(beta),0.0);
4894
4895 printf("(%10.6f %10.6f) ",-cos(gamma)*sin(beta)/sqrt(2.0),-sin(gamma)*sin(beta)/sqrt(2.0));
4896 printf("\n");
4897 printf("(%10.6f %10.6f) ",(cos(alpha)*cos(gamma)+sin(alpha)*sin(gamma))*(1.0-cos(beta))/2.0,
4898 (-cos(alpha)*sin(gamma)+sin(alpha)*cos(gamma))*(1.0-cos(beta))/2.0);
4899
4900 printf("(%10.6f %10.6f) ",cos(alpha)*sin(beta)/sqrt(2.0),sin(alpha)*sin(beta)/sqrt(2.0));
4901
4902 printf("(%10.6f %10.6f) ",(cos(alpha)*cos(gamma)-sin(alpha)*sin(gamma))*(1.0+cos(beta))/2.0,
4903 (cos(alpha)*sin(gamma)+sin(alpha)*cos(gamma))*(1.0+cos(beta))/2.0);
4904 printf("\n");
4905 }else if(j==2){
4906 printf("for Comparision:\n");
4907 /* m'= 2 */
4908 printf("(%10.6f %10.6f) ", (1.0+cos(beta))*(1.0+cos(beta))/4.0*cos(2.0*alpha+2.0*gamma),-(1.0+cos(beta))*(1.0+cos(beta))/4.0*sin(2.0*alpha+2.0*gamma));
4909 printf("(%10.6f %10.6f) ", -0.5*sin(beta)*(1.0+cos(beta))*cos(2.0*alpha+gamma),0.5*sin(beta)*(1.0+cos(beta))*sin(2.0*alpha+gamma));
4910 printf("(%10.6f %10.6f) ", sqrt(6.0)/4.0*sin(beta)*sin(beta)*cos(2.0*alpha),-sqrt(6.0)/4.0*sin(beta)*sin(beta)*sin(2.0*alpha));
4911 printf("(%10.6f %10.6f) ", 0.5*sin(beta)*(cos(beta)-1.0)*cos(2.0*alpha-gamma),-0.5*sin(beta)*(cos(beta)-1.0)*sin(2.0*alpha-gamma));
4912 printf("(%10.6f %10.6f)\n", (1.0-cos(beta))*(1.0-cos(beta))/4.0*cos(2.0*alpha-2.0*gamma),-(1.0-cos(beta))*(1.0-cos(beta))/4.0*sin(2.0*alpha-2.0*gamma));
4913 /* m'= 1 */
4914 printf("(%10.6f %10.6f) ",0.5*sin(beta)*(1.0+cos(beta))*cos(alpha+2.0*gamma),-0.5*sin(beta)*(1.0+cos(beta))*sin(alpha+2.0*gamma));
4915 printf("(%10.6f %10.6f) ", ((1.0+cos(beta))*(1.0+cos(beta))/4.0-3.0/4.0*sin(beta)*sin(beta))*cos(alpha+gamma),-((1.0+cos(beta))*(1.0+cos(beta))/4.0-3.0/4.0*sin(beta)*sin(beta))*sin(alpha+gamma));
4916 printf("(%10.6f %10.6f) ", -sqrt(6.0)/4.0*sin(2.0*beta)*cos(alpha), sqrt(6.0)/4.0*sin(2.0*beta)*sin(alpha));
4917 printf("(%10.6f %10.6f) ",-(cos(2.0*beta)-cos(beta))*0.5*cos(alpha-gamma),(cos(2.0*beta)-cos(beta))*0.5*sin(alpha-gamma));
4918 printf("(%10.6f %10.6f)\n", 0.5*sin(beta)*(cos(beta)-1.0)*cos(alpha-2.0*gamma),-0.5*sin(beta)*(cos(beta)-1.0)*sin(alpha-2.0*gamma));
4919 /* m'= 0 */
4920 printf("(%10.6f %10.6f) ",sqrt(6.0)/4.0*sin(beta)*sin(beta)*cos(2.0*gamma),-sqrt(6.0)/4.0*sin(beta)*sin(beta)*sin(2.0*gamma));
4921 printf("(%10.6f %10.6f) ", sqrt(6.0)/4.0*sin(2.0*beta)*cos(gamma),-sqrt(6.0)/4.0*sin(2.0*beta)*sin(gamma));
4922 printf("(%10.6f %10.6f) ", 0.5*(3.0*cos(beta)*cos(beta)-1.0),0.0);
4923 printf("(%10.6f %10.6f) ", -sqrt(6.0)/4.0*sin(2.0*beta)*cos(-gamma),sqrt(6.0)/4.0*sin(2.0*beta)*sin(-gamma));
4924 printf("(%10.6f %10.6f)\n", sqrt(6.0)/4.0*sin(beta)*sin(beta)*cos(-2.0*gamma),-sqrt(6.0)/4.0*sin(beta)*sin(beta)*sin(-2.0*gamma));
4925 /* m'=-1 */
4926 printf("(%10.6f %10.6f) ", -0.5*sin(beta)*(cos(beta)-1.0)*cos(-1.0*alpha+2.0*gamma), 0.5*sin(beta)*(cos(beta)-1.0)*sin(-1.0*alpha+2.0*gamma));
4927 printf("(%10.6f %10.6f) ",-(cos(2.0*beta)-cos(beta))*0.5*cos(-alpha+gamma),(cos(2.0*beta)-cos(beta))*0.5*sin(-alpha+gamma));
4928 printf("(%10.6f %10.6f) ", sqrt(6.0)/4.0*sin(2.0*beta)*cos(-alpha),-sqrt(6.0)/4.0*sin(2.0*beta)*sin(-alpha));
4929 printf("(%10.6f %10.6f) ", (cos(2.0*beta)+cos(beta))*0.5*cos(-alpha-gamma),-(cos(2.0*beta)+cos(beta))*0.5*sin(-alpha-gamma));
4930 printf("(%10.6f %10.6f)\n", -0.5*sin(beta)*(1.0+cos(beta))*cos(-1.0*alpha-2.0*gamma),0.5*sin(beta)*(1.0+cos(beta))*sin(-1.0*alpha-2.0*gamma));
4931 /* m'=-2 */
4932 printf("(%10.6f %10.6f) ", (1.0-cos(beta))*(1.0-cos(beta))/4.0*cos(-2.0*alpha+2.0*gamma),-(1.0-cos(beta))*(1.0-cos(beta))/4.0*sin(-2.0*alpha+2.0*gamma));
4933 printf("(%10.6f %10.6f) ", -0.5*sin(beta)*(cos(beta)-1.0)*cos(-2.0*alpha+1.0*gamma),0.5*sin(beta)*(cos(beta)-1.0)*sin(-2.0*alpha+1.0*gamma));
4934 printf("(%10.6f %10.6f) ", sqrt(6.0)/4.0*sin(beta)*sin(beta)*cos(-2.0*alpha),-sqrt(6.0)/4.0*sin(beta)*sin(beta)*sin(-2.0*alpha));
4935 printf("(%10.6f %10.6f) ", 0.5*sin(beta)*(1.0+cos(beta))*cos(-2.0*alpha-1.0*gamma),-0.5*sin(beta)*(1.0+cos(beta))*sin(-2.0*alpha-1.0*gamma));
4936 printf("(%10.6f %10.6f)\n", (1.0+cos(beta))*(1.0+cos(beta))/4.0*cos(-2.0*alpha-2.0*gamma),-(1.0+cos(beta))*(1.0+cos(beta))/4.0*sin(-2.0*alpha-2.0*gamma));
4937 }
4938
4939 }
4940 /* The rotation matrix connecting real function is defined as M=U^(-1)*Dlm^(T)*U, U is the transfer matrix
4941 from real to imaginary function for orbitals
4942 p1 px
4943 p0 = U* py
4944 p-1 pz
4945
4946 d2 dz2
4947 d1 dx2-y2
4948 d0 = U* dxy
4949 d-1 dxz
4950 d-2 dyz
4951
4952 px' px
4953 py' = M * py
4954 pz' pz
4955
4956 dz2' dz2
4957 dx2-y2' dx2-y2
4958 dxy' = M * dxy
4959 dxz' dxz
4960 dyz' dyz
4961 Here ' means those in the rotated coordinate, without ' means those in original coordinate.
4962 */
4963 switch(j){
4964 case 1:
4965 Umat[0][0].r=-1.0/sqrt(2.0);
4966 Umat[0][0].i=0.0;
4967
4968 Umat[0][1].r=0.0;
4969 Umat[0][1].i=-1.0/sqrt(2.0);
4970
4971 Umat[0][2].r=0.0;
4972 Umat[0][2].i=0.0;
4973
4974 Umat[1][0].r=0.0;
4975 Umat[1][0].i=0.0;
4976
4977 Umat[1][1].r=0.0;
4978 Umat[1][1].i=0.0;
4979
4980 Umat[1][2].r=1.0;
4981 Umat[1][2].i=0.0;
4982
4983 Umat[2][0].r=1.0/sqrt(2.0);
4984 Umat[2][0].i=0.0;
4985
4986 Umat[2][1].r=0.0;
4987 Umat[2][1].i=-1.0/sqrt(2.0);
4988
4989 Umat[2][2].r=0.0;
4990 Umat[2][2].i=0.0;
4991 break;
4992 case 2:
4993 Umat[0][0].r=0.0;
4994 Umat[0][0].i=0.0;
4995
4996 Umat[0][1].r=1.0/sqrt(2.0);
4997 Umat[0][1].i=0.0;
4998
4999 Umat[0][2].r=0.0;
5000 Umat[0][2].i=1.0/sqrt(2.0);
5001
5002 Umat[0][3].r=0.0;
5003 Umat[0][3].i=0.0;
5004
5005 Umat[0][4].r=0.0;
5006 Umat[0][4].i=0.0;
5007
5008 Umat[1][0].r=0.0;
5009 Umat[1][0].i=0.0;
5010
5011 Umat[1][1].r=0.0;
5012 Umat[1][1].i=0.0;
5013
5014 Umat[1][2].r=0.0;
5015 Umat[1][2].i=0.0;
5016
5017 Umat[1][3].r=-1.0/sqrt(2.0);
5018 Umat[1][3].i=0.0;
5019
5020 Umat[1][4].r=0.0;
5021 Umat[1][4].i=-1.0/sqrt(2.0);
5022
5023 Umat[2][0].r=1.0;
5024 Umat[2][0].i=0.0;
5025
5026 Umat[2][1].r=0.0;
5027 Umat[2][1].i=0.0;
5028
5029 Umat[2][2].r=0.0;
5030 Umat[2][2].i=0.0;
5031
5032 Umat[2][3].r=0.0;
5033 Umat[2][3].i=0.0;
5034
5035 Umat[2][4].r=0.0;
5036 Umat[2][4].i=0.0;
5037
5038 Umat[3][0].r=0.0;
5039 Umat[3][0].i=0.0;
5040
5041 Umat[3][1].r=0.0;
5042 Umat[3][1].i=0.0;
5043
5044 Umat[3][2].r=0.0;
5045 Umat[3][2].i=0.0;
5046
5047 Umat[3][3].r=1.0/sqrt(2.0);
5048 Umat[3][3].i=0.0;
5049
5050 Umat[3][4].r=0.0;
5051 Umat[3][4].i=-1.0/sqrt(2.0);
5052
5053 Umat[4][0].r=0.0;
5054 Umat[4][0].i=0.0;
5055
5056 Umat[4][1].r=1.0/sqrt(2.0);
5057 Umat[4][1].i=0.0;
5058
5059 Umat[4][2].r=0.0;
5060 Umat[4][2].i=-1.0/sqrt(2.0);
5061
5062 Umat[4][3].r=0.0;
5063 Umat[4][3].i=0.0;
5064
5065 Umat[4][4].r=0.0;
5066 Umat[4][4].i=0.0;
5067 break;
5068
5069 case 3:
5070
5071 break;
5072
5073 }
5074
5075 for(jmp=0;jmp<2*j+1;jmp++){
5076 for(jm=0;jm<2*j+1;jm++){
5077 Umat_inv[jmp][jm].r=Umat[jm][jmp].r;
5078 Umat_inv[jmp][jm].i=-Umat[jm][jmp].i;
5079 }
5080 }
5081
5082 /* U^{-1}*Dlm^{T} ==> Utmp */
5083 for(jmp=0;jmp<2*j+1;jmp++){
5084 for(jm=0;jm<2*j+1;jm++){
5085 sumr=0.0; sumi=0.0;
5086 for(i=0;i<2*j+1;i++){
5087 sumr=sumr+Umat_inv[jmp][i].r*Dlm[jm][i].r-Umat_inv[jmp][i].i*Dlm[jm][i].i;
5088 sumi=sumi+Umat_inv[jmp][i].r*Dlm[jm][i].i+Umat_inv[jmp][i].i*Dlm[jm][i].r;
5089 }
5090 Utmp[jmp][jm].r=sumr;
5091 Utmp[jmp][jm].i=sumi;
5092 }
5093 }
5094
5095 /* Utmp * U */
5096 for(jmp=0;jmp<2*j+1;jmp++){
5097 for(jm=0;jm<2*j+1;jm++){
5098 sumr=0.0; sumi=0.0;
5099 for(i=0;i<2*j+1;i++){
5100 sumr=sumr+Utmp[jmp][i].r*Umat[i][jm].r-Utmp[jmp][i].i*Umat[i][jm].i;
5101 sumi=sumi+Utmp[jmp][i].r*Umat[i][jm].i+Utmp[jmp][i].i*Umat[i][jm].r;
5102 }
5103 RotMat_for_Real_Func[jmp][jm].r=sumr;
5104 RotMat_for_Real_Func[jmp][jm].i=sumi;
5105 }
5106 }
5107
5108 for(jmp=0;jmp<2*j+1;jmp++){
5109 for(jm=0;jm<2*j+1;jm++){
5110
5111 if(fabs(RotMat_for_Real_Func[jmp][jm].i)>1e-8){
5112
5113 if (myid==Host_ID){
5114 printf("ERROR:WF The rotational matrix should be real.mp=%i,m=%i\n",j-jmp,j-jm);
5115 }
5116
5117 MPI_Finalize();
5118 exit(0);
5119 }
5120
5121 tmpRotMat[jmp][jm]=RotMat_for_Real_Func[jmp][jm].r;
5122 }
5123 }
5124
5125 if(myid==Host_ID && 0<level_stdout){
5126 printf("Final matrix for real functional rotation:\n");
5127 for(jmp=0;jmp<2*j+1;jmp++){
5128 for(jm=0;jm<2*j+1;jm++){
5129 printf("%10.6f ",RotMat_for_Real_Func[jmp][jm].r);
5130 }
5131 printf("\n");
5132 }
5133 printf("\n");
5134 }
5135
5136 } /* end of Get_Rotational_Matrix */
5137
5138
5139
5140
5141
5142
5143