1 /**********************************************************************
2 Mulliken_Charge.c:
3
4 Mulliken_Charge.c is a subrutine to calculate Mulliken charge.
5
6 Log of Mulliken_Charge.c:
7
8 27/Dec/2002 Released by T.Ozaki
9
10 ***********************************************************************/
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <math.h>
16 #include <time.h>
17 #include "openmx_common.h"
18 #include "mpi.h"
19
20
21 #define stdout_MulP 1
22
23
24
Mulliken_Charge(char * mode)25 double Mulliken_Charge( char *mode )
26 {
27 int Mc_AN,Gc_AN,tno0,Cwan,num,l,m,mul;
28 int wan1,wan2,i,j,k,Gs_AN,s_AN,spin;
29 int tag=999;
30 double MulP[4],My_Total_SpinS;
31 double MulP_LA[4],MulP_LB[4],MagML;
32 double summx,summy,summz;
33 double x,y,z,TZ;
34 double My_Total_SpinSx;
35 double My_Total_SpinSy;
36 double My_Total_SpinSz;
37 double My_Total_OrbitalMoment;
38 double My_Total_OrbitalMomentx;
39 double My_Total_OrbitalMomenty;
40 double My_Total_OrbitalMomentz;
41 double sden,tmp,tmp0,tmp1,tmp2;
42 double Total_Mul_up,Total_Mul_dn,Total_Mul;
43 double tmpA0,tmpA1,tmpA2;
44 double tmpB0,tmpB1,tmpB2;
45 double Re11,Re22,Re12,Im12,d1,d2,d3;
46 double theta[2],phi[2],Nup[2],Ndown[2],sit,cot,sip,cop;
47 double Stime_atom, Etime_atom;
48 double my_data[4],data[4];
49 double ***DecMulP;
50 double *tmp_array0;
51 double *sum_l,*sum_mul;
52 double S_coordinate[3];
53 double TStime,TEtime,time0;
54 char *Name_Angular[Supported_MaxL+1][2*(Supported_MaxL+1)+1];
55 char file_MC[YOUSO10] = ".MC";
56 FILE *fp_MC;
57 int numprocs,myid,ID;
58 char buf[fp_bsize]; /* setvbuf */
59 MPI_Status stat;
60
61 MPI_Barrier(mpi_comm_level1);
62 dtime(&TStime);
63
64 /* MPI */
65 MPI_Comm_size(mpi_comm_level1,&numprocs);
66 MPI_Comm_rank(mpi_comm_level1,&myid);
67
68 /* calculate TZ */
69
70 TZ = 0.0;
71 for (i=1; i<=atomnum; i++){
72 wan1 = WhatSpecies[i];
73 TZ = TZ + Spe_Core_Charge[wan1];
74 }
75
76 /* allocation of arrays */
77
78 DecMulP = (double***)malloc(sizeof(double**)*(SpinP_switch+1));
79 for (spin=0; spin<(SpinP_switch+1); spin++){
80 DecMulP[spin] = (double**)malloc(sizeof(double*)*(atomnum+1));
81 for (Gc_AN=0; Gc_AN<=atomnum; Gc_AN++){
82 DecMulP[spin][Gc_AN] = (double*)malloc(sizeof(double)*List_YOUSO[7]);
83 for (i=0; i<List_YOUSO[7]; i++) DecMulP[spin][Gc_AN][i] = 0.0;
84 }
85 }
86
87 tmp_array0 = (double*)malloc(sizeof(double)*(List_YOUSO[7]*(SpinP_switch+1)));
88 sum_l = (double*)malloc(sizeof(double)*(SpinP_switch+1));
89 sum_mul = (double*)malloc(sizeof(double)*(SpinP_switch+1));
90
91 My_Total_SpinS = 0.0;
92 My_Total_SpinSx = 0.0;
93 My_Total_SpinSy = 0.0;
94 My_Total_SpinSz = 0.0;
95
96 My_Total_OrbitalMoment = 0.0;
97 My_Total_OrbitalMomentx = 0.0;
98 My_Total_OrbitalMomenty = 0.0;
99 My_Total_OrbitalMomentz = 0.0;
100
101 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
102
103 dtime(&Stime_atom);
104
105 Gc_AN = M2G[Mc_AN];
106 wan1 = WhatSpecies[Gc_AN];
107
108 MulP[0] = 0.0;
109 MulP[1] = 0.0;
110 MulP[2] = 0.0;
111 MulP[3] = 0.0;
112
113 /****************************************************
114 in case of NEGF, add partial decomposed Mulliken
115 population coming from CL or CR
116 ****************************************************/
117
118 if (Solver==4){
119 for (spin=0; spin<=SpinP_switch; spin++){
120 for (i=0; i<Spe_Total_CNO[wan1]; i++){
121 DecMulP[spin][Gc_AN][i] += TRAN_DecMulP[spin][Mc_AN][i];
122 MulP[spin] += TRAN_DecMulP[spin][Mc_AN][i];
123 }
124 }
125 }
126
127 /****************************************************
128 loop for neighbouring atoms
129 ****************************************************/
130
131 for (s_AN=0; s_AN<=FNAN[Gc_AN]; s_AN++){
132 Gs_AN = natn[Gc_AN][s_AN];
133 wan2 = WhatSpecies[Gs_AN];
134
135 if (Cnt_switch==0){
136 for (spin=0; spin<=SpinP_switch; spin++){
137 for (i=0; i<Spe_Total_CNO[wan1]; i++){
138
139 tmp0 = 0.0;
140 for (j=0; j<Spe_Total_CNO[wan2]; j++){
141 tmp0 += DM[0][spin][Mc_AN][s_AN][i][j]*OLP[0][Mc_AN][s_AN][i][j];
142 }
143
144 /* due to difference in the definition between density matrix and density */
145 if (spin==3){
146 DecMulP[spin][Gc_AN][i] -= tmp0;
147 MulP[spin] -= tmp0;
148 }
149 else{
150 DecMulP[spin][Gc_AN][i] += tmp0;
151 MulP[spin] += tmp0;
152 }
153 }
154 }
155 }
156 else{
157 for (spin=0; spin<=SpinP_switch; spin++){
158 for (i=0; i<Spe_Total_CNO[wan1]; i++){
159
160 tmp0 = 0.0;
161 for (j=0; j<Spe_Total_CNO[wan2]; j++){
162 tmp0 += DM[0][spin][Mc_AN][s_AN][i][j]*CntOLP[0][Mc_AN][s_AN][i][j];
163 }
164
165 /* due to difference in the definition between density matrix and density */
166 if (spin==3){
167 DecMulP[spin][Gc_AN][i] -= tmp0;
168 MulP[spin] -= tmp0;
169 }
170 else {
171 DecMulP[spin][Gc_AN][i] += tmp0;
172 MulP[spin] += tmp0;
173 }
174 }
175 }
176 }
177 }
178
179 /****************************************************
180 if (SpinP_switch==3)
181 spin non-collinear
182
183 U \rho U^+ = n
184 ****************************************************/
185
186 if (SpinP_switch==3){
187
188 Re11 = MulP[0];
189 Re22 = MulP[1];
190 Re12 = MulP[2];
191 Im12 = MulP[3];
192
193 EulerAngle_Spin( 1, Re11, Re22, Re12, Im12, Re12, -Im12, Nup, Ndown, theta, phi );
194
195 MulP[0] = Nup[0];
196 MulP[1] = Ndown[0];
197 MulP[2] = theta[0];
198 MulP[3] = phi[0];
199
200 /* decomposed Mulliken populations */
201 for (i=0; i<Spe_Total_CNO[wan1]; i++){
202
203 Re11 = DecMulP[0][Gc_AN][i];
204 Re22 = DecMulP[1][Gc_AN][i];
205 Re12 = DecMulP[2][Gc_AN][i];
206 Im12 = DecMulP[3][Gc_AN][i];
207
208 EulerAngle_Spin( 1, Re11, Re22, Re12, Im12, Re12, -Im12, Nup, Ndown, theta, phi );
209
210 DecMulP[0][Gc_AN][i] = Nup[0];
211 DecMulP[1][Gc_AN][i] = Ndown[0];
212 DecMulP[2][Gc_AN][i] = theta[0];
213 DecMulP[3][Gc_AN][i] = phi[0];
214 }
215
216 }
217
218 /****************************************************
219 set InitN_USpin and InitN_DSpin
220 ****************************************************/
221
222 if (SpinP_switch==0){
223 InitN_USpin[Gc_AN] = MulP[0];
224 InitN_DSpin[Gc_AN] = MulP[0];
225 }
226 else if (SpinP_switch==1){
227 InitN_USpin[Gc_AN] = MulP[0];
228 InitN_DSpin[Gc_AN] = MulP[1];
229 }
230 else if (SpinP_switch==3){
231 InitN_USpin[Gc_AN] = MulP[0];
232 InitN_DSpin[Gc_AN] = MulP[1];
233 Angle0_Spin[Gc_AN] = MulP[2];
234 Angle1_Spin[Gc_AN] = MulP[3];
235 }
236
237 /****************************************************
238 My_Total_SpinS
239 ****************************************************/
240
241 /* spin non-collinear */
242 if (SpinP_switch==3){
243
244 theta[0] = Angle0_Spin[Gc_AN];
245 phi[0] = Angle1_Spin[Gc_AN];
246 sden = 0.5*(InitN_USpin[Gc_AN] - InitN_DSpin[Gc_AN]);
247 My_Total_SpinSx += sden*sin(theta[0])*cos(phi[0]);
248 My_Total_SpinSy += sden*sin(theta[0])*sin(phi[0]);
249 My_Total_SpinSz += sden*cos(theta[0]);
250
251 }
252 /* spin collinear */
253 else {
254 My_Total_SpinS += 0.5*(InitN_USpin[Gc_AN] - InitN_DSpin[Gc_AN]);
255 }
256
257 dtime(&Etime_atom);
258 time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
259
260 } /* Mc_AN */
261
262 /****************************************
263 MPI:
264
265 My_Total_SpinS
266 ****************************************/
267
268 /* spin non-collinear */
269
270 if (SpinP_switch==3){
271
272 /* spin moment */
273
274 my_data[0] = My_Total_SpinSx;
275 my_data[1] = My_Total_SpinSy;
276 my_data[2] = My_Total_SpinSz;
277
278 MPI_Allreduce(&my_data[0], &data[0], 3, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
279
280 Total_SpinSx = data[0];
281 Total_SpinSy = data[1];
282 Total_SpinSz = data[2];
283
284 xyz2spherical( Total_SpinSx,Total_SpinSy,Total_SpinSz, 0.0,0.0,0.0, S_coordinate );
285
286 Total_SpinS = S_coordinate[0];
287 Total_SpinAngle0 = S_coordinate[1];
288 Total_SpinAngle1 = S_coordinate[2];
289 }
290
291 /* spin collinear */
292 else{
293 MPI_Allreduce(&My_Total_SpinS, &Total_SpinS, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
294 }
295
296 /* MPI: InitN_USpin and InitN_DSpin */
297
298 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
299
300 ID = G2ID[Gc_AN];
301
302 if (ID!=Host_ID){
303
304 if (myid==ID){
305
306 /* spin non-collinear */
307 if (SpinP_switch==3){
308 data[0] = InitN_USpin[Gc_AN];
309 data[1] = InitN_DSpin[Gc_AN];
310 data[2] = Angle0_Spin[Gc_AN];
311 data[3] = Angle1_Spin[Gc_AN];
312 MPI_Send(&data[0], 4, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1);
313 }
314 else{
315 data[0] = InitN_USpin[Gc_AN];
316 data[1] = InitN_DSpin[Gc_AN];
317 MPI_Send(&data[0], 2, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1);
318 }
319 }
320
321 else if (myid==Host_ID){
322
323 /* spin non-collinear */
324 if (SpinP_switch==3){
325 MPI_Recv(&data[0], 4, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
326 InitN_USpin[Gc_AN] = data[0];
327 InitN_DSpin[Gc_AN] = data[1];
328 Angle0_Spin[Gc_AN] = data[2];
329 Angle1_Spin[Gc_AN] = data[3];
330 }
331 else{
332 MPI_Recv(&data[0], 2, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
333 InitN_USpin[Gc_AN] = data[0];
334 InitN_DSpin[Gc_AN] = data[1];
335 }
336 }
337 }
338 }
339
340 /* MPI: DecMulP */
341
342 if ( strcasecmp(mode,"write")==0 ){
343
344 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
345
346 ID = G2ID[Gc_AN];
347
348 if (ID!=Host_ID){
349
350 if (myid==ID){
351
352 num = 0;
353 for (spin=0; spin<(SpinP_switch+1); spin++){
354 wan1 = WhatSpecies[Gc_AN];
355 for (i=0; i<Spe_Total_CNO[wan1]; i++){
356 tmp_array0[num] = DecMulP[spin][Gc_AN][i];
357 num++;
358 }
359 }
360
361 MPI_Send(&tmp_array0[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1);
362
363 }
364 else if (myid==Host_ID){
365
366 wan1 = WhatSpecies[Gc_AN];
367 num = (SpinP_switch+1)*Spe_Total_CNO[wan1];
368 MPI_Recv(&tmp_array0[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
369
370 num = 0;
371 for (spin=0; spin<(SpinP_switch+1); spin++){
372 wan1 = WhatSpecies[Gc_AN];
373 for (i=0; i<Spe_Total_CNO[wan1]; i++){
374 DecMulP[spin][Gc_AN][i] = tmp_array0[num];
375 num++;
376 }
377 }
378 }
379 }
380 }
381 }
382
383 /****************************************
384 stdout MulP
385 ****************************************/
386
387 if (myid==Host_ID && stdout_MulP && strcasecmp(mode,"stdout")==0) {
388
389 Total_Mul_up = 0.0;
390 Total_Mul_dn = 0.0;
391
392 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
393
394 wan1 = WhatSpecies[Gc_AN];
395
396 if (SpinP_switch==0){
397
398 Total_Mul_up += InitN_USpin[Gc_AN];
399 Total_Mul_dn += InitN_USpin[Gc_AN];
400
401 if (Gc_AN<=20 && 1<=level_stdout){
402 printf(" %4d %4s MulP %8.4f%8.4f sum %8.4f\n",
403 Gc_AN,SpeName[wan1],
404 InitN_USpin[Gc_AN],InitN_USpin[Gc_AN],
405 InitN_USpin[Gc_AN]+InitN_USpin[Gc_AN]);
406 }
407 }
408 else if (SpinP_switch==1){
409
410 Total_Mul_up += InitN_USpin[Gc_AN];
411 Total_Mul_dn += InitN_DSpin[Gc_AN];
412
413 if (Gc_AN<=20 && level_stdout<=1){
414 printf(" %4d %4s MulP %8.4f %8.4f sum %8.4f diff %8.4f\n",
415 Gc_AN,SpeName[wan1],
416 InitN_USpin[Gc_AN],InitN_DSpin[Gc_AN],
417 InitN_USpin[Gc_AN]+InitN_DSpin[Gc_AN],
418 InitN_USpin[Gc_AN]-InitN_DSpin[Gc_AN]);
419 }
420 }
421 else if (SpinP_switch==3){
422
423 sden = (InitN_USpin[Gc_AN] - InitN_DSpin[Gc_AN]);
424 theta[0] = Angle0_Spin[Gc_AN];
425 phi[0] = Angle1_Spin[Gc_AN];
426 x = sden*sin(theta[0])*cos(phi[0]);
427 y = sden*sin(theta[0])*sin(phi[0]);
428 z = sden*cos(theta[0]);
429
430 MagML = OrbitalMoment[Gc_AN];
431 theta[0] = Angle0_Orbital[Gc_AN];
432 phi[0] = Angle1_Orbital[Gc_AN];
433 x += MagML*sin(theta[0])*cos(phi[0]);
434 y += MagML*sin(theta[0])*sin(phi[0]);
435 z += MagML*cos(theta[0]);
436
437 xyz2spherical( x,y,z, 0.0,0.0,0.0, S_coordinate );
438
439 Total_Mul_up += InitN_USpin[Gc_AN];
440 Total_Mul_dn += InitN_DSpin[Gc_AN];
441
442 if (Gc_AN<=20 && level_stdout<=1){
443 printf(" %4d %4s MulP%5.2f%5.2f sum %5.2f diff %5.2f (%6.2f %6.2f) Ml %4.2f (%6.2f %6.2f) Ml+s %4.2f (%6.2f %6.2f)\n",
444 Gc_AN,SpeName[wan1],
445 InitN_USpin[Gc_AN],InitN_DSpin[Gc_AN],
446 InitN_USpin[Gc_AN]+InitN_DSpin[Gc_AN],
447 InitN_USpin[Gc_AN]-InitN_DSpin[Gc_AN],
448 Angle0_Spin[Gc_AN]/PI*180.0, Angle1_Spin[Gc_AN]/PI*180.0,
449 OrbitalMoment[Gc_AN],
450 Angle0_Orbital[Gc_AN]/PI*180.0, Angle1_Orbital[Gc_AN]/PI*180.0,
451 S_coordinate[0],S_coordinate[1]/PI*180.0,S_coordinate[2]/PI*180.0);
452 }
453
454 }
455
456 } /* Gc_AN */
457
458 if (20<atomnum && level_stdout<=1){
459 printf(" ..........\n");
460 printf(" ......\n\n");
461 }
462
463 if (0<level_stdout){
464
465 printf(" Sum of MulP: up =%12.5f down =%12.5f\n",
466 Total_Mul_up,Total_Mul_dn);
467 printf(" total=%12.5f ideal(neutral)=%12.5f\n",
468 Total_Mul_up+Total_Mul_dn,TZ);
469 }
470
471 Total_Mul = Total_Mul_up + Total_Mul_dn;
472 }
473
474 MPI_Bcast(&Total_Mul, 1, MPI_DOUBLE, Host_ID, mpi_comm_level1);
475
476 /********************************************************
477 check the stability of the eigenvalue solver
478 by looking the total number of electrons.
479 In some cases with degenerate states, the eigenvectors
480 are not properly calculated, resulting in violation of
481 norm of the eigenstates.
482 *********************************************************/
483
484 rediagonalize_flag_overlap_matrix = 0;
485
486 if (Solver!=4 && strcasecmp(mode,"stdout")==0){
487
488 int tmp_flag;
489 double Dnum;
490
491 Dnum = TZ - Total_Mul - system_charge;
492
493 if (1.0e-8<fabs(Dnum)){
494
495 if (dste_flag==2) tmp_flag = 3; /* vx -> qr */
496 else if (dste_flag==3) tmp_flag = 1; /* qr -> dc */
497 else if (dste_flag==1) tmp_flag = 0; /* dc -> gr */
498 else if (dste_flag==0){ /* gr -> ELPA1 */
499 tmp_flag = 2;
500 rediagonalize_flag_overlap_matrix_ELPA1 = 1;
501 }
502
503 dste_flag = tmp_flag;
504
505 /****************************************************
506 In Cluster_DFT, Band_DFT_Col, or Band_DFT_NonCol,
507 the overlap matrix will be rediagonalized.
508 ****************************************************/
509
510 if (myid==Host_ID && 1<level_stdout){
511 printf("Eigensolver changed Dnum=%18.15f dste_flag=%2d rediagonalize_flag_overlap_matrix_ELPA1=%2d\n",
512 Dnum,dste_flag,rediagonalize_flag_overlap_matrix_ELPA1);
513 }
514
515 rediagonalize_flag_overlap_matrix = 1;
516 }
517 }
518
519 /****************************************
520 file, *.MC
521 ****************************************/
522
523 if (strcasecmp(mode,"write")==0 ){
524 /* MPI: InitN_USpin and InitN_DSpin */
525 MPI_Bcast(&InitN_USpin[0], atomnum+1, MPI_DOUBLE, Host_ID, mpi_comm_level1);
526 MPI_Bcast(&InitN_DSpin[0], atomnum+1, MPI_DOUBLE, Host_ID, mpi_comm_level1);
527 }
528
529 if ( myid==Host_ID && strcasecmp(mode,"write")==0 ){
530
531 Total_Mul_up = 0.0;
532 Total_Mul_dn = 0.0;
533
534 sprintf(file_MC,"%s%s.MC",filepath,filename);
535
536 if ((fp_MC = fopen(file_MC,"w")) != NULL){
537
538 #ifdef xt3
539 setvbuf(fp_MC,buf,_IOFBF,fp_bsize); /* setvbuf */
540 #endif
541
542 fprintf(fp_MC,"\n");
543 fprintf(fp_MC,"***********************************************************\n");
544 fprintf(fp_MC,"***********************************************************\n");
545 fprintf(fp_MC," Mulliken populations \n");
546 fprintf(fp_MC,"***********************************************************\n");
547 fprintf(fp_MC,"***********************************************************\n\n");
548
549 /* spin non-collinear */
550 if (SpinP_switch==3){
551
552 /* total Mulliken charge */
553
554 fprintf(fp_MC," Total spin moment (muB) %12.9f Angles (Deg) %12.9f %12.9f\n\n",
555 2.0*Total_SpinS, Total_SpinAngle0/PI*180.0, Total_SpinAngle1/PI*180.0);
556 fprintf(fp_MC," Up Down Sum Diff theta phi\n");
557 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
558 wan1 = WhatSpecies[Gc_AN];
559 fprintf(fp_MC," %4d %4s %8.5f %8.5f %8.5f %8.5f %10.5f %10.5f\n",
560 Gc_AN,SpeName[wan1],InitN_USpin[Gc_AN],InitN_DSpin[Gc_AN],
561 InitN_USpin[Gc_AN]+InitN_DSpin[Gc_AN],
562 InitN_USpin[Gc_AN]-InitN_DSpin[Gc_AN],
563 Angle0_Spin[Gc_AN]/PI*180.0, Angle1_Spin[Gc_AN]/PI*180.0 );
564
565 Total_Mul_up += InitN_USpin[Gc_AN];
566 Total_Mul_dn += InitN_DSpin[Gc_AN];
567
568 }
569
570 fprintf(fp_MC,"\n");
571 fprintf(fp_MC," Sum of MulP: up =%12.5f down =%12.5f\n",
572 Total_Mul_up,Total_Mul_dn);
573 fprintf(fp_MC," total=%12.5f ideal(neutral)=%12.5f\n",
574 Total_Mul_up+Total_Mul_dn,TZ);
575
576 }
577 else{
578
579 /* total Mulliken charge */
580 fprintf(fp_MC," Total spin moment (muB) %12.9f\n\n",2.0*Total_SpinS);
581 fprintf(fp_MC," Up spin Down spin Sum Diff\n");
582 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
583 wan1 = WhatSpecies[Gc_AN];
584 fprintf(fp_MC," %4d %4s %12.9f %12.9f %12.9f %12.9f\n",
585 Gc_AN,SpeName[wan1],InitN_USpin[Gc_AN],InitN_DSpin[Gc_AN],
586 InitN_USpin[Gc_AN]+InitN_DSpin[Gc_AN],
587 InitN_USpin[Gc_AN]-InitN_DSpin[Gc_AN]);
588
589 Total_Mul_up += InitN_USpin[Gc_AN];
590 Total_Mul_dn += InitN_DSpin[Gc_AN];
591 }
592
593 fprintf(fp_MC,"\n");
594 fprintf(fp_MC," Sum of MulP: up =%12.5f down =%12.5f\n",
595 Total_Mul_up,Total_Mul_dn);
596 fprintf(fp_MC," total=%12.5f ideal(neutral)=%12.5f\n",
597 Total_Mul_up+Total_Mul_dn,TZ);
598 }
599
600 /* decomposed Mulliken charge */
601
602 Name_Angular[0][0] = "s ";
603 Name_Angular[1][0] = "px ";
604 Name_Angular[1][1] = "py ";
605 Name_Angular[1][2] = "pz ";
606 Name_Angular[2][0] = "d3z^2-r^2 ";
607 Name_Angular[2][1] = "dx^2-y^2 ";
608 Name_Angular[2][2] = "dxy ";
609 Name_Angular[2][3] = "dxz ";
610 Name_Angular[2][4] = "dyz ";
611 Name_Angular[3][0] = "f5z^2-3r^2 ";
612 Name_Angular[3][1] = "f5xz^2-xr^2";
613 Name_Angular[3][2] = "f5yz^2-yr^2";
614 Name_Angular[3][3] = "fzx^2-zy^2 ";
615 Name_Angular[3][4] = "fxyz ";
616 Name_Angular[3][5] = "fx^3-3*xy^2";
617 Name_Angular[3][6] = "f3yx^2-y^3 ";
618 Name_Angular[4][0] = "g1 ";
619 Name_Angular[4][1] = "g2 ";
620 Name_Angular[4][2] = "g3 ";
621 Name_Angular[4][3] = "g4 ";
622 Name_Angular[4][4] = "g5 ";
623 Name_Angular[4][5] = "g6 ";
624 Name_Angular[4][6] = "g7 ";
625 Name_Angular[4][7] = "g8 ";
626 Name_Angular[4][8] = "g9 ";
627
628 fprintf(fp_MC,"\n\n Decomposed Mulliken populations\n");
629
630 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
631
632 wan1 = WhatSpecies[Gc_AN];
633
634 /* spin collinear */
635 if (SpinP_switch==0 || SpinP_switch==1){
636
637 fprintf(fp_MC,"\n %4d %4s Up spin Down spin Sum Diff\n",
638 Gc_AN,SpeName[wan1]);
639 fprintf(fp_MC," multiple\n");
640 }
641
642 /* spin non-collinear */
643 else{
644 fprintf(fp_MC,"\n %4d %4s Up spin Down spin Sum Diff Angles(Deg)\n",
645 Gc_AN,SpeName[wan1]);
646 fprintf(fp_MC," multiple\n");
647 }
648
649 num = 0;
650 for (l=0; l<=Supported_MaxL; l++){
651
652 if (SpinP_switch==0){
653 sum_l[0] = 0.0;
654 }
655 else if (SpinP_switch==1){
656 sum_l[0] = 0.0;
657 sum_l[1] = 0.0;
658 }
659 else if (SpinP_switch==3){
660 sum_l[0] = 0.0;
661 sum_l[1] = 0.0;
662 }
663
664 for (mul=0; mul<Spe_Num_CBasis[wan1][l]; mul++){
665
666 if (SpinP_switch==0){
667 sum_mul[0] = 0.0;
668 }
669 else if (SpinP_switch==1){
670 sum_mul[0] = 0.0;
671 sum_mul[1] = 0.0;
672 }
673 else if (SpinP_switch==3){
674 sum_mul[0] = 0.0;
675 sum_mul[1] = 0.0;
676 }
677
678 for (m=0; m<(2*l+1); m++){
679
680 if (SpinP_switch==0){
681 fprintf(fp_MC," %s%2d %12.9f %12.9f %12.9f %12.9f\n",
682 Name_Angular[l][m],mul,
683 DecMulP[0][Gc_AN][num],DecMulP[0][Gc_AN][num],
684 DecMulP[0][Gc_AN][num]+DecMulP[0][Gc_AN][num],
685 DecMulP[0][Gc_AN][num]-DecMulP[0][Gc_AN][num]);
686
687 sum_mul[0] += DecMulP[0][Gc_AN][num];
688 }
689 else if (SpinP_switch==1){
690 fprintf(fp_MC," %s%2d %12.9f %12.9f %12.9f %12.9f\n",
691 Name_Angular[l][m],mul,
692 DecMulP[0][Gc_AN][num],DecMulP[1][Gc_AN][num],
693 DecMulP[0][Gc_AN][num]+DecMulP[1][Gc_AN][num],
694 DecMulP[0][Gc_AN][num]-DecMulP[1][Gc_AN][num]);
695
696 sum_mul[0] += DecMulP[0][Gc_AN][num];
697 sum_mul[1] += DecMulP[1][Gc_AN][num];
698 }
699
700 else if (SpinP_switch==3){
701 fprintf(fp_MC," %s%2d %12.9f %12.9f %12.9f %12.9f %9.4f %9.4f\n",
702 Name_Angular[l][m],mul,
703 DecMulP[0][Gc_AN][num],DecMulP[1][Gc_AN][num],
704 DecMulP[0][Gc_AN][num]+DecMulP[1][Gc_AN][num],
705 DecMulP[0][Gc_AN][num]-DecMulP[1][Gc_AN][num],
706 DecMulP[2][Gc_AN][num]/PI*180.0,
707 DecMulP[3][Gc_AN][num]/PI*180.0);
708
709 sum_mul[0] += DecMulP[0][Gc_AN][num];
710 sum_mul[1] += DecMulP[1][Gc_AN][num];
711 }
712
713 num++;
714 }
715
716 if (SpinP_switch==0){
717
718 fprintf(fp_MC," sum over m %12.9f %12.9f %12.9f %12.9f\n",
719 sum_mul[0],sum_mul[0],
720 sum_mul[0]+sum_mul[0],
721 sum_mul[0]-sum_mul[0]);
722
723 sum_l[0] += sum_mul[0];
724 }
725 else if (SpinP_switch==1){
726 fprintf(fp_MC," sum over m %12.9f %12.9f %12.9f %12.9f\n",
727 sum_mul[0],sum_mul[1],
728 sum_mul[0]+sum_mul[1],
729 sum_mul[0]-sum_mul[1]);
730
731 sum_l[0] += sum_mul[0];
732 sum_l[1] += sum_mul[1];
733 }
734
735 else if (SpinP_switch==3){
736 fprintf(fp_MC," sum over m %12.9f %12.9f %12.9f %12.9f\n",
737 sum_mul[0],sum_mul[1],
738 sum_mul[0]+sum_mul[1],
739 sum_mul[0]-sum_mul[1]);
740
741 sum_l[0] += sum_mul[0];
742 sum_l[1] += sum_mul[1];
743 }
744
745
746 }
747
748 if (Spe_Num_CBasis[wan1][l]!=0){
749
750 if (SpinP_switch==0){
751 fprintf(fp_MC," sum over m+mul %12.9f %12.9f %12.9f %12.9f\n",
752 sum_l[0],sum_l[0],
753 sum_l[0]+sum_l[0],
754 sum_l[0]-sum_l[0]);
755 }
756 else if (SpinP_switch==1){
757 fprintf(fp_MC," sum over m+mul %12.9f %12.9f %12.9f %12.9f\n",
758 sum_l[0],sum_l[1],
759 sum_l[0]+sum_l[1],
760 sum_l[0]-sum_l[1]);
761 }
762
763 else if (SpinP_switch==3){
764 fprintf(fp_MC," sum over m+mul %12.9f %12.9f %12.9f %12.9f\n",
765 sum_l[0],sum_l[1],
766 sum_l[0]+sum_l[1],
767 sum_l[0]-sum_l[1]);
768 }
769
770 }
771 }
772 }
773
774 fclose(fp_MC);
775 }
776 else{
777 printf("Failure of saving the MC file.\n");
778 }
779 }
780
781 /* freeing of arrays */
782
783 for (spin=0; spin<(SpinP_switch+1); spin++){
784 for (Gc_AN=0; Gc_AN<=atomnum; Gc_AN++){
785 free(DecMulP[spin][Gc_AN]);
786 }
787 free(DecMulP[spin]);
788 }
789 free(DecMulP);
790
791 free(tmp_array0);
792
793 free(sum_l);
794 free(sum_mul);
795
796 /* for time */
797 MPI_Barrier(mpi_comm_level1);
798 dtime(&TEtime);
799 time0 = TEtime - TStime;
800 return time0;
801 }
802
803
804
805
806
807
808
809