1 /**********************************************************************
2 OutData_Binary.c:
3
4 OutData_Binary.c is a subrutine to output values of electron densities,
5 potetianls, wave functions, and etc. on the grids in the format of
6 Gaussian cube, and atomic cartesian coordinates.
7 However, the data are stored in a binary format, and the text format
8 can be recovered by a post processing code "bin2txt.c".
9
10 Log of OutData_Binary.c:
11
12 27/Dec/2012 Released by T. Ozaki
13 ***********************************************************************/
14
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <math.h>
18 #include <string.h>
19 #include <time.h>
20 #include "openmx_common.h"
21 #include "mpi.h"
22 #include <omp.h>
23
24 #define CUBE_EXTENSION ".cube.bin"
25 #define CAT_FLAG 0
26
27
28 static void out_density();
29 static void out_Veff();
30 static void out_Vhart();
31 static void out_Vna();
32 static void out_Vxc();
33 static void out_grid();
34 static void out_atomxyz();
35 static void out_atomxsf();
36 static void out_coordinates_bulk();
37 static void out_Cluster_MO();
38 static void out_Cluster_NC_MO();
39 static void out_Bulk_MO();
40 static void out_OrbOpt(char *inputfile);
41 static void out_Partial_Charge_Density();
42 static void Set_Partial_Density_Grid(double *****CDM);
43 static void Print_CubeTitle(FILE *fp, int EigenValue_flag, double EigenValue);
44 static void Print_CubeData(FILE *fp, char fext[], double *data, double *data1,char *op);
45 static void Print_VectorData(FILE *fp, char fext[],
46 double *data0, double *data1,
47 double *data2, double *data3);
48
OutData_Binary(char * inputfile)49 double OutData_Binary(char *inputfile)
50 {
51 char operate[YOUSO10];
52 int i,c;
53 int numprocs,myid;
54 char fname1[300];
55 char fname2[300];
56 FILE *fp1,*fp2;
57 char buf[fp_bsize]; /* setvbuf */
58 char buf1[fp_bsize]; /* setvbuf */
59 char buf2[fp_bsize]; /* setvbuf */
60 double Stime,Etime;
61
62 MPI_Comm_size(mpi_comm_level1,&numprocs);
63 MPI_Comm_rank(mpi_comm_level1,&myid);
64
65 dtime(&Stime);
66
67 if (myid==Host_ID && 0<level_stdout) printf("\noutputting data on grids to files...\n\n");
68
69 out_atomxyz();
70
71 if (1<=level_fileout){
72 out_density();
73 out_Veff();
74 out_Vhart();
75 }
76
77 if (2<=level_fileout){
78 out_grid();
79 if (ProExpn_VNA==0) out_Vna();
80 out_Vxc();
81 }
82
83 if (specified_system!=0) out_coordinates_bulk();
84
85 if (MO_fileout==1 && specified_system==0){
86
87 /* spin non-collinear */
88 if (SpinP_switch==3)
89 out_Cluster_NC_MO();
90 /* spin collinear */
91 else
92 out_Cluster_MO();
93 }
94 else if (MO_fileout==1 && specified_system!=0){
95 out_Bulk_MO();
96 }
97
98 if (Cnt_switch==1){
99 out_OrbOpt(inputfile);
100 }
101
102 /*************************************************
103 partial charge density for STM simulations
104 The routine, out_Partial_Charge_Density(),
105 has to be called after calling out_density(),
106 since an array Density_Grid is used in
107 out_Partial_Charge_Density().
108 *************************************************/
109
110 if ( cal_partial_charge ){
111 out_Partial_Charge_Density();
112 }
113
114 /* divide-conquer, gdc, or Krylov */
115
116 if (Dos_fileout && (Solver==5 || Solver==6 || Solver==8) ) {
117
118 if (myid==Host_ID){
119
120 if (CAT_FLAG){
121
122 /* merge files */
123
124 for (i=0; i<numprocs; i++){
125 if (i==0)
126 sprintf(operate,"cat %s%s.Dos.vec%i > tmp1",filepath,filename,i);
127 else
128 sprintf(operate,"cat %s%s.Dos.vec%i >> tmp1",filepath,filename,i);
129
130 system(operate);
131
132 sprintf(operate,"rm %s%s.Dos.vec%i",filepath,filename,i);
133 system(operate);
134 }
135
136 sprintf(operate,"mv tmp1 %s%s.Dos.vec",filepath,filename);
137 system(operate);
138
139 }
140
141 else {
142
143 /* merge files */
144
145 sprintf(fname1,"%s%s.Dos.vec",filepath,filename);
146 fp1 = fopen(fname1,"ab");
147
148 if (fp1!=NULL){
149 remove(fname1);
150 fclose(fp1);
151 }
152
153 for (i=0; i<numprocs; i++){
154
155 sprintf(fname1,"%s%s.Dos.vec",filepath,filename);
156 fp1 = fopen(fname1,"ab");
157 fseek(fp1,0,SEEK_END);
158 sprintf(fname2,"%s%s.Dos.vec%i",filepath,filename,i);
159 fp2 = fopen(fname2,"rb");
160
161 if (fp2!=NULL){
162 for (c=getc(fp2); c!=EOF; c=getc(fp2)) putc(c,fp1);
163 fclose(fp2);
164 }
165 fclose(fp1);
166 }
167
168 for (i=0; i<numprocs; i++){
169 sprintf(fname2,"%s%s.Dos.vec%i",filepath,filename,i);
170 remove(fname2);
171 }
172
173 }
174
175 }
176 }
177
178 /* elapsed time */
179 dtime(&Etime);
180 return Etime - Stime;
181 }
182
183
out_density()184 void out_density()
185 {
186 int ct_AN,spe,i1,i2,i3,GN,i,j,k;
187 int MN;
188 double x,y,z,vx,vy,vz;
189 double phi,theta,sden,oden;
190 double xmin,ymin,zmin,xmax,ymax,zmax;
191 char fname[YOUSO10];
192 char file1[YOUSO10] = ".tden";
193 char file2[YOUSO10] = ".den0";
194 char file3[YOUSO10] = ".den1";
195 char file4[YOUSO10] = ".sden";
196 char file9[YOUSO10] = ".nc.xsf";
197 char file10[YOUSO10] = ".ncsden.xsf.bin";
198 char file12[YOUSO10] = ".nco.xsf";
199 char file11[YOUSO10] = ".dden";
200 FILE *fp;
201 char buf[fp_bsize]; /* setvbuf */
202 int numprocs,myid,ID,digit;
203
204 /* MPI */
205 MPI_Comm_size(mpi_comm_level1,&numprocs);
206 MPI_Comm_rank(mpi_comm_level1,&myid);
207 digit = (int)log10(numprocs) + 1;
208
209 strcat(file1,CUBE_EXTENSION);
210 strcat(file2,CUBE_EXTENSION);
211 strcat(file3,CUBE_EXTENSION);
212 strcat(file4,CUBE_EXTENSION);
213 strcat(file11,CUBE_EXTENSION);
214
215 /****************************************************
216 electron density - atomic charge density
217 ****************************************************/
218
219 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file11,digit,myid);
220
221 if ((fp = fopen(fname,"wb")) != NULL){
222
223 #ifdef xt3
224 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
225 #endif
226
227 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
228
229 for (MN=0; MN<My_NumGridB_AB; MN++){
230 ADensity_Grid_B[MN] = Density_Grid_B[0][MN]
231 +Density_Grid_B[1][MN]
232 -2.0*ADensity_Grid_B[MN];
233 }
234
235 Print_CubeData(fp,file11,ADensity_Grid_B,(void*)NULL,(void*)NULL);
236 }
237 else{
238 printf("Failure of saving the electron density\n");
239 }
240
241 /****************************************************
242 total electron density
243 ****************************************************/
244
245 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file1,digit,myid);
246
247 if ((fp = fopen(fname,"wb")) != NULL){
248
249 #ifdef xt3
250 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
251 #endif
252
253 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
254
255 if (SpinP_switch==0) {
256 for (MN=0; MN<My_NumGridB_AB; MN++){
257 Density_Grid_B[0][MN] = 2.0*Density_Grid_B[0][MN];
258 }
259 Print_CubeData(fp,file1,Density_Grid_B[0],(void*)NULL,(void*)NULL);
260 }
261 else {
262 Print_CubeData(fp,file1,Density_Grid_B[0],Density_Grid_B[1],"add");
263 }
264
265 }
266 else{
267 printf("Failure of saving the electron density\n");
268 }
269
270 /* spin polization */
271
272 if (SpinP_switch==1 || SpinP_switch==3){
273
274 /****************************************************
275 up-spin electron density
276 ****************************************************/
277
278 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file2,digit,myid);
279
280 if ((fp = fopen(fname,"wb")) != NULL){
281
282 #ifdef xt3
283 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
284 #endif
285
286 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
287 Print_CubeData(fp,file2,Density_Grid_B[0],NULL,NULL);
288 }
289 else{
290 printf("Failure of saving the electron density\n");
291 }
292
293 /****************************************************
294 down-spin electron density
295 ****************************************************/
296
297 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file3,digit,myid);
298
299 if ((fp = fopen(fname,"wb")) != NULL){
300
301 #ifdef xt3
302 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
303 #endif
304
305 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
306 Print_CubeData(fp,file3,Density_Grid_B[1],NULL,NULL);
307 }
308 else{
309 printf("Failure of saving the electron density\n");
310 }
311
312 /****************************************************
313 spin electron density
314 ****************************************************/
315
316 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file4,digit,myid);
317
318 if ((fp = fopen(fname,"wb")) != NULL){
319
320 #ifdef xt3
321 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
322 #endif
323
324 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
325 Print_CubeData(fp,file4,Density_Grid_B[0],Density_Grid_B[1],"diff");
326 }
327 else{
328 printf("Failure of saving the electron density\n");
329 }
330 }
331
332 /****************************************************
333 spin electron density with a spin vector
334 ****************************************************/
335
336 if (SpinP_switch==3){
337
338 /* for XCrysDen, *.ncsden.xsf */
339
340 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file10,digit,myid);
341
342 if ((fp = fopen(fname,"wb")) != NULL){
343
344 #ifdef xt3
345 setvbuf(fp,buf,_IOFBF,fp_bsize);
346 #endif
347
348 Print_VectorData(fp,file10,Density_Grid_B[0],Density_Grid_B[1],
349 Density_Grid_B[2],Density_Grid_B[3]);
350
351 }
352 else{
353 printf("Failure of saving the electron spin density with a vector\n");
354 }
355
356 /* non-collinear spin density by Mulliken population */
357
358 if (myid==Host_ID){
359
360 /* for XCrysDen, *.nc.xsf */
361
362 sprintf(fname,"%s%s%s",filepath,filename,file9);
363 if ((fp = fopen(fname,"w")) != NULL){
364
365 #ifdef xt3
366 setvbuf(fp,buf,_IOFBF,fp_bsize);
367 #endif
368
369 fprintf(fp,"CRYSTAL\n");
370 fprintf(fp,"PRIMVEC\n");
371 fprintf(fp," %10.6f %10.6f %10.6f\n",BohrR*tv[1][1], BohrR*tv[1][2], BohrR*tv[1][3]);
372 fprintf(fp," %10.6f %10.6f %10.6f\n",BohrR*tv[2][1], BohrR*tv[2][2], BohrR*tv[2][3]);
373 fprintf(fp," %10.6f %10.6f %10.6f\n",BohrR*tv[3][1], BohrR*tv[3][2], BohrR*tv[3][3]);
374 fprintf(fp,"PRIMCOORD\n");
375 fprintf(fp,"%4d 1\n",atomnum);
376
377 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
378 i = WhatSpecies[ct_AN];
379 j = Spe_WhatAtom[i];
380
381 theta = Angle0_Spin[ct_AN];
382 phi = Angle1_Spin[ct_AN];
383 sden = InitN_USpin[ct_AN] - InitN_DSpin[ct_AN];
384
385 vx = sden*sin(theta)*cos(phi);
386 vy = sden*sin(theta)*sin(phi);
387 vz = sden*cos(theta);
388
389 fprintf(fp,"%s %13.3E %13.3E %13.3E %13.3E %13.3E %13.3E\n",
390 Atom_Symbol[j],
391 BohrR*Gxyz[ct_AN][1],
392 BohrR*Gxyz[ct_AN][2],
393 BohrR*Gxyz[ct_AN][3],
394 vx,vy,vz);
395 }
396
397 fclose(fp);
398 }
399 else{
400 printf("Failure of saving the Mulliken spin vector\n");
401 }
402
403 /* non-collinear obital magnetic moment by 'on-site' approximation */
404 /* for XCrysDen .xsf */
405
406 sprintf(fname,"%s%s%s",filepath,filename,file12);
407 if ((fp = fopen(fname,"w")) != NULL){
408
409 #ifdef xt3
410 setvbuf(fp,buf,_IOFBF,fp_bsize);
411 #endif
412
413 fprintf(fp,"CRYSTAL\n");
414 fprintf(fp,"PRIMVEC\n");
415 fprintf(fp," %10.6f %10.6f %10.6f\n",BohrR*tv[1][1], BohrR*tv[1][2], BohrR*tv[1][3]);
416 fprintf(fp," %10.6f %10.6f %10.6f\n",BohrR*tv[2][1], BohrR*tv[2][2], BohrR*tv[2][3]);
417 fprintf(fp," %10.6f %10.6f %10.6f\n",BohrR*tv[3][1], BohrR*tv[3][2], BohrR*tv[3][3]);
418 fprintf(fp,"PRIMCOORD\n");
419 fprintf(fp,"%4d 1\n",atomnum);
420
421 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
422
423 i = WhatSpecies[ct_AN];
424 j = Spe_WhatAtom[i];
425
426 theta = Angle0_Orbital[ct_AN];
427 phi = Angle1_Orbital[ct_AN];
428 oden = OrbitalMoment[ct_AN];
429
430 vx = oden*sin(theta)*cos(phi);
431 vy = oden*sin(theta)*sin(phi);
432 vz = oden*cos(theta);
433
434 fprintf(fp,"%s %13.3E %13.3E %13.3E %13.3E %13.3E %13.3E\n",
435 Atom_Symbol[j],
436 BohrR*Gxyz[ct_AN][1],
437 BohrR*Gxyz[ct_AN][2],
438 BohrR*Gxyz[ct_AN][3],
439 vx,vy,vz);
440 }
441
442 fclose(fp);
443 }
444 else{
445 printf("Failure of saving orbital magnetic moments\n");
446 }
447
448 } /* if (myid==Host_ID) */
449 } /* if (SpinP_switch==3) */
450
451 }
452
453
454
out_Vhart()455 static void out_Vhart()
456 {
457 int ct_AN,spe,i1,i2,i3,GN,i,j,k,digit;
458 char fname[YOUSO10];
459 char file1[YOUSO10] = ".vhart";
460 FILE *fp;
461 char buf[fp_bsize]; /* setvbuf */
462 int numprocs,myid,ID;
463
464 /* MPI */
465 MPI_Comm_size(mpi_comm_level1,&numprocs);
466 MPI_Comm_rank(mpi_comm_level1,&myid);
467 digit = (int)log10(numprocs) + 1;
468
469 strcat(file1,CUBE_EXTENSION);
470
471 /****************************************************
472 Hartree potential
473 ****************************************************/
474
475 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file1,digit,myid);
476 if ((fp = fopen(fname,"wb")) != NULL){
477
478 #ifdef xt3
479 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
480 #endif
481
482 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
483 Print_CubeData(fp,file1,dVHart_Grid_B,NULL,NULL);
484 }
485 else{
486 printf("Failure of saving the electron density\n");
487 }
488
489 }
490
491
492
493
out_Vna()494 static void out_Vna()
495 {
496 int ct_AN,spe,i1,i2,i3,GN,i,j,k;
497 char fname[YOUSO10];
498 char file1[YOUSO10] = ".vna";
499 FILE *fp;
500 char buf[fp_bsize]; /* setvbuf */
501 int numprocs,myid,ID,digit;
502
503 /* MPI */
504 MPI_Comm_size(mpi_comm_level1,&numprocs);
505 MPI_Comm_rank(mpi_comm_level1,&myid);
506 digit = (int)log10(numprocs) + 1;
507
508 strcat(file1,CUBE_EXTENSION);
509
510 /****************************************************
511 neutral atom potential
512 ****************************************************/
513
514 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file1,digit,myid);
515
516 if ((fp = fopen(fname,"wb")) != NULL){
517
518 #ifdef xt3
519 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
520 #endif
521
522 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
523 Print_CubeData(fp,file1,VNA_Grid_B,NULL,NULL);
524 }
525 else{
526 printf("Failure of saving the electron density\n");
527 }
528 }
529
530
531
532
533
out_Vxc()534 static void out_Vxc()
535 {
536 int ct_AN,spe,i1,i2,i3,GN,i,j,k;
537 int LN,BN,CN,Ng1,Ng2,Ng3,max_spin;
538 double *Work_Array_Snd_Grid_C2B;
539 double *Work_Array_Rcv_Grid_C2B;
540 double **Work_Array_B;
541 char fname[YOUSO10];
542 char file1[YOUSO10] = ".vxc0";
543 char file2[YOUSO10] = ".vxc1";
544 char buf[fp_bsize]; /* setvbuf */
545 FILE *fp;
546 int numprocs,myid,ID,IDS,IDR,tag=999;
547 int digit;
548 MPI_Status stat;
549 MPI_Request request;
550
551 /* MPI */
552 MPI_Comm_size(mpi_comm_level1,&numprocs);
553 MPI_Comm_rank(mpi_comm_level1,&myid);
554 digit = (int)log10(numprocs) + 1;
555
556 strcat(file1,CUBE_EXTENSION);
557 strcat(file2,CUBE_EXTENSION);
558
559 /* allocation of arrays */
560
561 if (SpinP_switch==0) max_spin = 0;
562 else if (SpinP_switch==1) max_spin = 1;
563 else if (SpinP_switch==3) max_spin = 1;
564
565 Work_Array_Snd_Grid_C2B = (double*)malloc(sizeof(double)*Max_Num_Rcv_Grid_B2C*(max_spin+1));
566 Work_Array_Rcv_Grid_C2B = (double*)malloc(sizeof(double)*Max_Num_Snd_Grid_B2C*(max_spin+1));
567 Work_Array_B = (double**)malloc(sizeof(double*)*(max_spin+1));
568 for (i=0; i<(max_spin+1); i++){
569 Work_Array_B[i] = (double*)malloc(sizeof(double)*My_NumGridB_AB);
570 for (j=0; j<My_NumGridB_AB; j++){
571 Work_Array_B[i][j] = 0.0;
572 }
573 }
574
575 /****************************************************
576 MPI communication for Vxc_Grid
577 from the partitions C to B
578 ****************************************************/
579
580 Ng1 = Max_Grid_Index[1] - Min_Grid_Index[1] + 1;
581 Ng2 = Max_Grid_Index[2] - Min_Grid_Index[2] + 1;
582 Ng3 = Max_Grid_Index[3] - Min_Grid_Index[3] + 1;
583
584 tag = 999;
585 for (ID=0; ID<numprocs; ID++){
586
587 IDS = (myid + ID) % numprocs;
588 IDR = (myid - ID + numprocs) % numprocs;
589
590 /* copy Vxc_Grid to Work_Array_Snd_Grid_C2B */
591
592 for (LN=0; LN<Num_Rcv_Grid_B2C[IDS]; LN++){
593
594 CN = Index_Rcv_Grid_B2C[IDS][LN];
595
596 i = CN/(Ng2*Ng3);
597 j = (CN-i*Ng2*Ng3)/Ng3;
598 k = CN - i*Ng2*Ng3 - j*Ng3;
599
600 if ( i<=1 || (Ng1-2)<=i || j<=1 || (Ng2-2)<=j || k<=1 || (Ng3-2)<=k ){
601
602 if (max_spin==0){
603 Work_Array_Snd_Grid_C2B[LN] = (1.0e+10)+0.1;
604 }
605 else {
606 Work_Array_Snd_Grid_C2B[2*LN+0] = (1.0e+10)+0.1;
607 Work_Array_Snd_Grid_C2B[2*LN+1] = (1.0e+10)+0.1;
608 }
609 }
610 else{
611
612 if (max_spin==0){
613 Work_Array_Snd_Grid_C2B[LN] = Vxc_Grid[0][CN];
614 }
615 else {
616 Work_Array_Snd_Grid_C2B[2*LN+0] = Vxc_Grid[0][CN];
617 Work_Array_Snd_Grid_C2B[2*LN+1] = Vxc_Grid[1][CN];
618 }
619 }
620 } /* LN */
621
622 /* MPI_Isend and MPI_Recv */
623
624 if (Num_Rcv_Grid_B2C[IDS]!=0){
625 MPI_Isend( &Work_Array_Snd_Grid_C2B[0], Num_Rcv_Grid_B2C[IDS]*(max_spin+1),
626 MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request);
627 }
628
629 if (Num_Snd_Grid_B2C[IDR]!=0){
630 MPI_Recv( &Work_Array_Rcv_Grid_C2B[0], Num_Snd_Grid_B2C[IDR]*(max_spin+1),
631 MPI_DOUBLE, IDR, tag, mpi_comm_level1, &stat);
632 }
633
634 /* MPI_Wait */
635
636 if (Num_Rcv_Grid_B2C[IDS]!=0) MPI_Wait(&request,&stat);
637
638 /* copy Work_Array_Rcv_Grid_C2B to Work_Array_B */
639
640 for (LN=0; LN<Num_Snd_Grid_B2C[IDR]; LN++){
641
642 BN = Index_Snd_Grid_B2C[IDR][LN];
643
644 if (max_spin==0){
645 if (Work_Array_Rcv_Grid_C2B[LN]<1.0e+10){
646 Work_Array_B[0][BN] = Work_Array_Rcv_Grid_C2B[LN];
647 }
648 }
649 else {
650 if (Work_Array_Rcv_Grid_C2B[2*LN]<1.0e+10){
651 Work_Array_B[0][BN] = Work_Array_Rcv_Grid_C2B[2*LN ];
652 Work_Array_B[1][BN] = Work_Array_Rcv_Grid_C2B[2*LN+1];
653 }
654 }
655
656 } /* LN */
657 } /* ID */
658
659 /****************************************************
660 exchange-correlation potential for up-spin
661 ****************************************************/
662
663 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file1,digit,myid);
664 if ((fp = fopen(fname,"wb")) != NULL){
665
666 #ifdef xt3
667 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
668 #endif
669
670 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
671 Print_CubeData(fp,file1,Work_Array_B[0],NULL,NULL);
672 }
673 else{
674 printf("Failure of saving the electron density\n");
675 }
676
677 /****************************************************
678 exchange-correlation potential for down-spin
679 ****************************************************/
680
681 if (SpinP_switch==1){
682
683 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file2,digit,myid);
684
685 if ((fp = fopen(fname,"wb")) != NULL){
686
687 #ifdef xt3
688 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
689 #endif
690
691 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
692 Print_CubeData(fp,file2,Work_Array_B[1],NULL,NULL);
693 }
694 else{
695 printf("Failure of saving the electron density\n");
696 }
697 }
698
699 /* freeing of arrays */
700 free(Work_Array_Snd_Grid_C2B);
701 free(Work_Array_Rcv_Grid_C2B);
702 for (i=0; i<(max_spin+1); i++){
703 free(Work_Array_B[i]);
704 }
705 free(Work_Array_B);
706 }
707
708
709
710
711
out_Veff()712 static void out_Veff()
713 {
714 int ct_AN,spe,i1,i2,i3,GN,i,j,k,digit;
715 int LN,BN,CN,Ng1,Ng2,Ng3,max_spin;
716 double *Work_Array_Snd_Grid_C2B;
717 double *Work_Array_Rcv_Grid_C2B;
718 double **Work_Array_B;
719 char fname[YOUSO10];
720 char file1[YOUSO10] = ".v0";
721 char file2[YOUSO10] = ".v1";
722 FILE *fp;
723 char buf[fp_bsize]; /* setvbuf */
724 int numprocs,myid,ID,IDS,IDR,tag=999;
725 MPI_Status stat;
726 MPI_Request request;
727
728 /* MPI */
729 MPI_Comm_size(mpi_comm_level1,&numprocs);
730 MPI_Comm_rank(mpi_comm_level1,&myid);
731 digit = (int)log10(numprocs) + 1;
732
733 strcat(file1,CUBE_EXTENSION);
734 strcat(file2,CUBE_EXTENSION);
735
736 /* allocation of arrays */
737
738 if (SpinP_switch==0) max_spin = 0;
739 else if (SpinP_switch==1) max_spin = 1;
740 else if (SpinP_switch==3) max_spin = 1;
741
742 Work_Array_Snd_Grid_C2B = (double*)malloc(sizeof(double)*Max_Num_Rcv_Grid_B2C*(max_spin+1));
743 Work_Array_Rcv_Grid_C2B = (double*)malloc(sizeof(double)*Max_Num_Snd_Grid_B2C*(max_spin+1));
744 Work_Array_B = (double**)malloc(sizeof(double*)*(max_spin+1));
745 for (i=0; i<(max_spin+1); i++){
746 Work_Array_B[i] = (double*)malloc(sizeof(double)*My_NumGridB_AB);
747 for (j=0; j<My_NumGridB_AB; j++){
748 Work_Array_B[i][j] = 0.0;
749 }
750 }
751
752 /****************************************************
753 MPI communication for Vpot_Grid
754 from the partitions C to B
755 ****************************************************/
756
757 Ng1 = Max_Grid_Index[1] - Min_Grid_Index[1] + 1;
758 Ng2 = Max_Grid_Index[2] - Min_Grid_Index[2] + 1;
759 Ng3 = Max_Grid_Index[3] - Min_Grid_Index[3] + 1;
760
761 tag = 999;
762 for (ID=0; ID<numprocs; ID++){
763
764 IDS = (myid + ID) % numprocs;
765 IDR = (myid - ID + numprocs) % numprocs;
766
767 /* copy Vpot_Grid to Work_Array_Snd_Grid_C2B */
768
769 for (LN=0; LN<Num_Rcv_Grid_B2C[IDS]; LN++){
770
771 CN = Index_Rcv_Grid_B2C[IDS][LN];
772
773 i = CN/(Ng2*Ng3);
774 j = (CN-i*Ng2*Ng3)/Ng3;
775 k = CN - i*Ng2*Ng3 - j*Ng3;
776
777 if ( i<=1 || (Ng1-2)<=i || j<=1 || (Ng2-2)<=j || k<=1 || (Ng3-2)<=k ){
778
779 if (max_spin==0){
780 Work_Array_Snd_Grid_C2B[LN] = (1.0e+10)+0.1;
781 }
782 else {
783 Work_Array_Snd_Grid_C2B[2*LN+0] = (1.0e+10)+0.1;
784 Work_Array_Snd_Grid_C2B[2*LN+1] = (1.0e+10)+0.1;
785 }
786 }
787 else{
788
789 if (max_spin==0){
790 Work_Array_Snd_Grid_C2B[LN] = Vpot_Grid[0][CN];
791 }
792 else {
793 Work_Array_Snd_Grid_C2B[2*LN+0] = Vpot_Grid[0][CN];
794 Work_Array_Snd_Grid_C2B[2*LN+1] = Vpot_Grid[1][CN];
795 }
796 }
797 } /* LN */
798
799 /* MPI_Isend and MPI_Recv */
800
801 if (Num_Rcv_Grid_B2C[IDS]!=0){
802 MPI_Isend( &Work_Array_Snd_Grid_C2B[0], Num_Rcv_Grid_B2C[IDS]*(max_spin+1),
803 MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request);
804 }
805
806 if (Num_Snd_Grid_B2C[IDR]!=0){
807 MPI_Recv( &Work_Array_Rcv_Grid_C2B[0], Num_Snd_Grid_B2C[IDR]*(max_spin+1),
808 MPI_DOUBLE, IDR, tag, mpi_comm_level1, &stat);
809 }
810
811 /* MPI_Wait */
812
813 if (Num_Rcv_Grid_B2C[IDS]!=0) MPI_Wait(&request,&stat);
814
815 /* copy Work_Array_Rcv_Grid_C2B to Work_Array_B */
816
817 for (LN=0; LN<Num_Snd_Grid_B2C[IDR]; LN++){
818
819 BN = Index_Snd_Grid_B2C[IDR][LN];
820
821 if (max_spin==0){
822 if (Work_Array_Rcv_Grid_C2B[LN]<1.0e+10){
823 Work_Array_B[0][BN] = Work_Array_Rcv_Grid_C2B[LN];
824 }
825 }
826 else {
827 if (Work_Array_Rcv_Grid_C2B[2*LN]<1.0e+10){
828 Work_Array_B[0][BN] = Work_Array_Rcv_Grid_C2B[2*LN ];
829 Work_Array_B[1][BN] = Work_Array_Rcv_Grid_C2B[2*LN+1];
830 }
831 }
832
833 } /* LN */
834 } /* ID */
835
836 /****************************************************
837 Kohn-Sham potential for up-spin
838 ****************************************************/
839
840 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file1,digit,myid);
841 if ((fp = fopen(fname,"wb")) != NULL){
842
843 #ifdef xt3
844 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
845 #endif
846
847 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
848 Print_CubeData(fp,file1,Work_Array_B[0],NULL,NULL);
849 }
850 else{
851 printf("Failure of saving the electron density\n");
852 }
853
854 /****************************************************
855 Kohn-Sham potential for down-spin
856 ****************************************************/
857
858 if (SpinP_switch==1){
859
860 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file2,digit,myid);
861 if ((fp = fopen(fname,"wb")) != NULL){
862
863 #ifdef xt3
864 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
865 #endif
866
867 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
868 Print_CubeData(fp,file2,Work_Array_B[1],NULL,NULL);
869 }
870 else{
871 printf("Failure of saving the electron density\n");
872 }
873 }
874
875 /* freeing of arrays */
876 free(Work_Array_Snd_Grid_C2B);
877 free(Work_Array_Rcv_Grid_C2B);
878 for (i=0; i<(max_spin+1); i++){
879 free(Work_Array_B[i]);
880 }
881 free(Work_Array_B);
882 }
883
884
885
886
out_grid()887 static void out_grid()
888 {
889 int N,n1,n2,n3,i;
890 char file1[YOUSO10] = ".grid.bin";
891 int numprocs,myid,ID;
892 double x,y,z;
893 double *dlist;
894 char buf[fp_bsize]; /* setvbuf */
895 FILE *fp;
896
897 /* MPI */
898 MPI_Comm_size(mpi_comm_level1,&numprocs);
899 MPI_Comm_rank(mpi_comm_level1,&myid);
900
901 /****************************************************
902 output the real space grids to a file, *.grid
903 ****************************************************/
904
905 if (myid==Host_ID){
906
907 dlist = (double*)malloc(sizeof(double)*Ngrid2*Ngrid3*3);
908
909 fnjoint(filepath,filename,file1);
910
911 if ((fp = fopen(file1,"wb")) != NULL){
912
913 #ifdef xt3
914 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
915 #endif
916
917 fwrite(&Ngrid1, sizeof(int), 1, fp);
918 fwrite(&Ngrid2, sizeof(int), 1, fp);
919 fwrite(&Ngrid3, sizeof(int), 1, fp);
920
921 for (n1=0; n1<Ngrid1; n1++){
922
923 i = 0;
924 for (n2=0; n2<Ngrid2; n2++){
925 for (n3=0; n3<Ngrid3; n3++){
926
927 x = (double)n1*gtv[1][1] + (double)n2*gtv[2][1]
928 + (double)n3*gtv[3][1] + Grid_Origin[1];
929 y = (double)n1*gtv[1][2] + (double)n2*gtv[2][2]
930 + (double)n3*gtv[3][2] + Grid_Origin[2];
931 z = (double)n1*gtv[1][3] + (double)n2*gtv[2][3]
932 + (double)n3*gtv[3][3] + Grid_Origin[3];
933
934 dlist[i] = BohrR*x; i++;
935 dlist[i] = BohrR*y; i++;
936 dlist[i] = BohrR*z; i++;
937 }
938 }
939 fwrite(dlist, sizeof(double), i, fp);
940 }
941
942 fclose(fp);
943 }
944 else{
945 printf("Failure of saving grids\n");
946 }
947
948 free(dlist);
949 }
950 }
951
952
953
954
955
956
957
out_atomxyz()958 static void out_atomxyz()
959 {
960 int ct_AN,spe,i1,i2,i3,GN,i,j,k;
961 char filexyz[YOUSO10] = ".xyz";
962 char buf[fp_bsize]; /* setvbuf */
963 FILE *fp;
964 int numprocs,myid,ID;
965
966 /* MPI */
967 MPI_Comm_size(mpi_comm_level1,&numprocs);
968 MPI_Comm_rank(mpi_comm_level1,&myid);
969
970 /****************************************************
971 cartesian coordinates
972 ****************************************************/
973
974 if (myid==Host_ID){
975
976 fnjoint(filepath,filename,filexyz);
977 if ((fp = fopen(filexyz,"w")) != NULL){
978
979 #ifdef xt3
980 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
981 #endif
982
983 fprintf(fp,"%i \n\n",atomnum);
984 for (k=1; k<=atomnum; k++){
985 i = WhatSpecies[k];
986 j = Spe_WhatAtom[i];
987 fprintf(fp,"%s %8.5f %8.5f %8.5f %18.15f %18.15f %18.15f\n",
988 Atom_Symbol[j],
989 Gxyz[k][1]*BohrR,Gxyz[k][2]*BohrR,Gxyz[k][3]*BohrR,
990 Gxyz[k][17],Gxyz[k][18],Gxyz[k][19]);
991 }
992 fclose(fp);
993 }
994 else{
995 printf("could not save the xyz file\n");
996 }
997 }
998 }
999
1000
1001
out_atomxsf()1002 static void out_atomxsf()
1003 {
1004 int ct_AN,spe,i1,i2,i3,GN,i,j,k;
1005 char filexsf[YOUSO10] = ".coord.xsf";
1006 char buf[fp_bsize]; /* setvbuf */
1007 FILE *fp;
1008 int numprocs,myid,ID;
1009
1010 /* MPI */
1011 MPI_Comm_size(mpi_comm_level1,&numprocs);
1012 MPI_Comm_rank(mpi_comm_level1,&myid);
1013
1014 /****************************************************
1015 cartesian coordinates
1016 ****************************************************/
1017
1018 if (myid==Host_ID){
1019
1020 fnjoint(filepath,filename,filexsf);
1021 if ((fp = fopen(filexsf,"w")) != NULL){
1022
1023 #ifdef xt3
1024 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1025 #endif
1026
1027 fprintf(fp,"CRYSTAL\n");
1028 fprintf(fp,"PRIMVEC\n");
1029 fprintf(fp," %10.6f %10.6f %10.6f\n",BohrR*tv[1][1], BohrR*tv[1][2], BohrR*tv[1][3]);
1030 fprintf(fp," %10.6f %10.6f %10.6f\n",BohrR*tv[2][1], BohrR*tv[2][2], BohrR*tv[2][3]);
1031 fprintf(fp," %10.6f %10.6f %10.6f\n",BohrR*tv[3][1], BohrR*tv[3][2], BohrR*tv[3][3]);
1032 fprintf(fp,"PRIMCOORD 1\n");
1033 fprintf(fp,"%4d %d\n",atomnum, 1);
1034
1035 for (k=1; k<=atomnum; k++){
1036 i = WhatSpecies[k];
1037 /*j = Spe_WhatAtom[i];*/
1038 fprintf(fp,"%4d %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f\n",
1039 Spe_WhatAtom[i],
1040 Gxyz[k][1]*BohrR,Gxyz[k][2]*BohrR,Gxyz[k][3]*BohrR,
1041 Gxyz[k][24],Gxyz[k][25],Gxyz[k][26]);
1042 }
1043 fclose(fp);
1044 }
1045 else{
1046 printf("failure of saving coord.xsf file\n");
1047 }
1048 }
1049 }
1050
1051
1052
out_coordinates_bulk()1053 void out_coordinates_bulk()
1054 {
1055 int n,i1,i2,i3,ct_AN,i,j;
1056 double tx,ty,tz,x,y,z;
1057 char file1[YOUSO10] = ".bulk.xyz";
1058 char buf[fp_bsize]; /* setvbuf */
1059 FILE *fp;
1060 int numprocs,myid,ID;
1061
1062 /* MPI */
1063 MPI_Comm_size(mpi_comm_level1,&numprocs);
1064 MPI_Comm_rank(mpi_comm_level1,&myid);
1065
1066 /****************************************************
1067 atomic coordinates including copied cells
1068 ****************************************************/
1069
1070 if (myid==Host_ID){
1071
1072 n = 1;
1073
1074 fnjoint(filepath,filename,file1);
1075
1076 if ((fp = fopen(file1,"w")) != NULL){
1077
1078 #ifdef xt3
1079 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1080 #endif
1081
1082 fprintf(fp,"%d\n\n",atomnum*(2*n+1)*(2*n+1)*(2*n+1));
1083
1084 for (i1=-n; i1<=n; i1++){
1085 for (i2=-n; i2<=n; i2++){
1086 for (i3=-n; i3<=n; i3++){
1087
1088 tx = (double)i1*tv[1][1] + (double)i2*tv[2][1] + (double)i3*tv[3][1];
1089 ty = (double)i1*tv[1][2] + (double)i2*tv[2][2] + (double)i3*tv[3][2];
1090 tz = (double)i1*tv[1][3] + (double)i2*tv[2][3] + (double)i3*tv[3][3];
1091
1092 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
1093 i = WhatSpecies[ct_AN];
1094 j = Spe_WhatAtom[i];
1095
1096 x = BohrR*(Gxyz[ct_AN][1] + tx);
1097 y = BohrR*(Gxyz[ct_AN][2] + ty);
1098 z = BohrR*(Gxyz[ct_AN][3] + tz);
1099 fprintf(fp,"%s %8.5f %8.5f %8.5f\n",Atom_Symbol[j],x,y,z);
1100 }
1101 }
1102 }
1103 }
1104
1105 fclose(fp);
1106 }
1107 else{
1108 printf("Failure of saving atomic coordinates\n");
1109 }
1110 }
1111
1112 }
1113
1114
1115
out_Cluster_MO()1116 void out_Cluster_MO()
1117 {
1118 int Mc_AN,Gc_AN,Cwan,NO0,spin,Nc;
1119 int orbit,GN,spe,i,i1,i2,i3,so;
1120 double *MO_Grid_tmp;
1121 double *MO_Grid;
1122 char file1[YOUSO10];
1123 char buf[fp_bsize]; /* setvbuf */
1124 FILE *fp;
1125 int numprocs,myid,ID;
1126
1127 /* MPI */
1128 MPI_Comm_size(mpi_comm_level1,&numprocs);
1129 MPI_Comm_rank(mpi_comm_level1,&myid);
1130
1131 /****************************************************
1132 allocation of arrays:
1133
1134 double MO_Grid_tmp[TNumGrid];
1135 double MO_Grid[TNumGrid];
1136 ****************************************************/
1137
1138 MO_Grid_tmp = (double*)malloc(sizeof(double)*TNumGrid);
1139 MO_Grid = (double*)malloc(sizeof(double)*TNumGrid);
1140
1141 /*************
1142 HOMOs
1143 *************/
1144
1145 for (so=0; so<=SO_switch; so++){
1146 for (spin=0; spin<=SpinP_switch; spin++){
1147 for (orbit=0; orbit<num_HOMOs; orbit++){
1148
1149 /* calc. MO on grids */
1150
1151 for (GN=0; GN<TNumGrid; GN++) MO_Grid_tmp[GN] = 0.0;
1152
1153 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1154 Gc_AN = M2G[Mc_AN];
1155 Cwan = WhatSpecies[Gc_AN];
1156 NO0 = Spe_Total_CNO[Cwan];
1157
1158 if (so==0){
1159 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
1160 GN = GridListAtom[Mc_AN][Nc];
1161 for (i=0; i<NO0; i++){
1162 MO_Grid_tmp[GN] += HOMOs_Coef[0][spin][orbit][Gc_AN][i].r*
1163 Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1164 }
1165 }
1166 }
1167
1168 else if (so==1){
1169 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
1170 GN = GridListAtom[Mc_AN][Nc];
1171 for (i=0; i<NO0; i++){
1172 MO_Grid_tmp[GN] += HOMOs_Coef[0][spin][orbit][Gc_AN][i].i*
1173 Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1174 }
1175 }
1176 }
1177
1178 }
1179
1180 MPI_Reduce(&MO_Grid_tmp[0], &MO_Grid[0], TNumGrid, MPI_DOUBLE,
1181 MPI_SUM, Host_ID, mpi_comm_level1);
1182
1183 /* output HOMO on grids */
1184
1185 if (myid==Host_ID){
1186
1187 if (so==0)
1188 sprintf(file1,"%s%s.homo%i_%i_r%s",filepath,filename,spin,orbit,CUBE_EXTENSION);
1189 else if (so==1)
1190 sprintf(file1,"%s%s.homo%i_%i_i%s",filepath,filename,spin,orbit,CUBE_EXTENSION);
1191
1192 if ((fp = fopen(file1,"wb")) != NULL){
1193
1194 #ifdef xt3
1195 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1196 #endif
1197
1198 Print_CubeTitle(fp,1,HOMOs_Coef[0][spin][orbit][0][0].r);
1199 fwrite(MO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1200 fclose(fp);
1201 }
1202 else{
1203 printf("Failure of saving MOs\n");
1204 }
1205
1206 }
1207
1208 } /* orbit */
1209 } /* spin */
1210 } /* so */
1211
1212 /*************
1213 LUMOs
1214 *************/
1215
1216 for (so=0; so<=SO_switch; so++){
1217 for (spin=0; spin<=SpinP_switch; spin++){
1218 for (orbit=0; orbit<num_LUMOs; orbit++){
1219
1220 /* calc. MO on grids */
1221
1222 for (GN=0; GN<TNumGrid; GN++) MO_Grid_tmp[GN] = 0.0;
1223
1224 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1225 Gc_AN = M2G[Mc_AN];
1226 Cwan = WhatSpecies[Gc_AN];
1227 NO0 = Spe_Total_CNO[Cwan];
1228
1229 if (so==0){
1230 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
1231 GN = GridListAtom[Mc_AN][Nc];
1232 for (i=0; i<NO0; i++){
1233 MO_Grid_tmp[GN] += LUMOs_Coef[0][spin][orbit][Gc_AN][i].r*
1234 Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1235 }
1236 }
1237 }
1238
1239 else if (so==1){
1240 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
1241 GN = GridListAtom[Mc_AN][Nc];
1242 for (i=0; i<NO0; i++){
1243 MO_Grid_tmp[GN] += LUMOs_Coef[0][spin][orbit][Gc_AN][i].i*
1244 Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1245 }
1246 }
1247 }
1248
1249 }
1250
1251 /* output LUMO on grids */
1252
1253 MPI_Reduce(&MO_Grid_tmp[0], &MO_Grid[0], TNumGrid, MPI_DOUBLE,
1254 MPI_SUM, Host_ID, mpi_comm_level1);
1255
1256 if (myid==Host_ID){
1257
1258 if (so==0)
1259 sprintf(file1,"%s%s.lumo%i_%i_r%s",filepath,filename,spin,orbit,CUBE_EXTENSION);
1260 else if (so==0)
1261 sprintf(file1,"%s%s.lumo%i_%i_i%s",filepath,filename,spin,orbit,CUBE_EXTENSION);
1262
1263
1264 if ((fp = fopen(file1,"wb")) != NULL){
1265
1266 #ifdef xt3
1267 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1268 #endif
1269
1270 Print_CubeTitle(fp,1,LUMOs_Coef[0][spin][orbit][0][0].r);
1271 fwrite(MO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1272 fclose(fp);
1273 }
1274 else{
1275 printf("Failure of saving MOs\n");
1276 }
1277 }
1278
1279 } /* orbit */
1280 } /* spin */
1281 } /* so */
1282
1283 /****************************************************
1284 freeing of arrays:
1285
1286 double MO_Grid[TNumGrid];
1287 double MO_Grid_tmp[TNumGrid];
1288 ****************************************************/
1289
1290 free(MO_Grid);
1291 free(MO_Grid_tmp);
1292 }
1293
1294
out_Cluster_NC_MO()1295 void out_Cluster_NC_MO()
1296 {
1297 int Mc_AN,Gc_AN,Cwan,NO0,spin,Nc;
1298 int orbit,GN,spe,i,i1,i2,i3;
1299 double *RMO_Grid_tmp,*IMO_Grid_tmp;
1300 double *RMO_Grid,*IMO_Grid;
1301 char file1[YOUSO10];
1302 char buf[fp_bsize]; /* setvbuf */
1303 FILE *fp;
1304 int numprocs,myid,ID;
1305
1306 /* MPI */
1307 MPI_Comm_size(mpi_comm_level1,&numprocs);
1308 MPI_Comm_rank(mpi_comm_level1,&myid);
1309
1310 /****************************************************
1311 allocation of arrays:
1312
1313 double RMO_Grid_tmp[TNumGrid];
1314 double IMO_Grid_tmp[TNumGrid];
1315 double RMO_Grid[TNumGrid];
1316 double IMO_Grid[TNumGrid];
1317 ****************************************************/
1318
1319 RMO_Grid_tmp = (double*)malloc(sizeof(double)*TNumGrid);
1320 IMO_Grid_tmp = (double*)malloc(sizeof(double)*TNumGrid);
1321 RMO_Grid = (double*)malloc(sizeof(double)*TNumGrid);
1322 IMO_Grid = (double*)malloc(sizeof(double)*TNumGrid);
1323
1324 /*************
1325 HOMOs
1326 *************/
1327
1328 for (spin=0; spin<=1; spin++){
1329 for (orbit=0; orbit<num_HOMOs; orbit++){
1330
1331 /* calc. MO on grids */
1332
1333 for (GN=0; GN<TNumGrid; GN++){
1334 RMO_Grid_tmp[GN] = 0.0;
1335 IMO_Grid_tmp[GN] = 0.0;
1336 }
1337
1338 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1339 Gc_AN = M2G[Mc_AN];
1340 Cwan = WhatSpecies[Gc_AN];
1341 NO0 = Spe_Total_CNO[Cwan];
1342 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
1343 GN = GridListAtom[Mc_AN][Nc];
1344 for (i=0; i<NO0; i++){
1345 RMO_Grid_tmp[GN] += HOMOs_Coef[0][spin][orbit][Gc_AN][i].r*
1346 Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1347 IMO_Grid_tmp[GN] += HOMOs_Coef[0][spin][orbit][Gc_AN][i].i*
1348 Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1349 }
1350 }
1351 }
1352
1353 MPI_Reduce(&RMO_Grid_tmp[0], &RMO_Grid[0], TNumGrid, MPI_DOUBLE,
1354 MPI_SUM, Host_ID, mpi_comm_level1);
1355 MPI_Reduce(&IMO_Grid_tmp[0], &IMO_Grid[0], TNumGrid, MPI_DOUBLE,
1356 MPI_SUM, Host_ID, mpi_comm_level1);
1357
1358 if (myid==Host_ID){
1359
1360 /* output the real part of HOMOs on grids */
1361
1362 sprintf(file1,"%s%s.homo%i_%i_r%s",
1363 filepath,filename,spin,orbit,CUBE_EXTENSION);
1364
1365 if ((fp = fopen(file1,"wb")) != NULL){
1366
1367 #ifdef xt3
1368 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1369 #endif
1370 Print_CubeTitle(fp,1,HOMOs_Coef[0][spin][orbit][0][0].r);
1371 fwrite(RMO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1372 fclose(fp);
1373 }
1374 else{
1375 printf("Failure of saving MOs\n");
1376 }
1377
1378 /* output the imaginary part of HOMOs on grids */
1379
1380 sprintf(file1,"%s%s.homo%i_%i_i%s",
1381 filepath,filename,spin,orbit,CUBE_EXTENSION);
1382
1383 if ((fp = fopen(file1,"wb")) != NULL){
1384
1385 #ifdef xt3
1386 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1387 #endif
1388
1389 Print_CubeTitle(fp,1,HOMOs_Coef[0][spin][orbit][0][0].r);
1390 fwrite(IMO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1391 fclose(fp);
1392 }
1393 else{
1394 printf("Failure of saving MOs\n");
1395 }
1396
1397 }
1398 } /* orbit */
1399 } /* spin */
1400
1401 /*************
1402 LUMOs
1403 *************/
1404
1405 for (spin=0; spin<=1; spin++){
1406 for (orbit=0; orbit<num_HOMOs; orbit++){
1407
1408 /* calc. MO on grids */
1409
1410 for (GN=0; GN<TNumGrid; GN++){
1411 RMO_Grid_tmp[GN] = 0.0;
1412 IMO_Grid_tmp[GN] = 0.0;
1413 }
1414
1415 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1416 Gc_AN = M2G[Mc_AN];
1417 Cwan = WhatSpecies[Gc_AN];
1418 NO0 = Spe_Total_CNO[Cwan];
1419 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
1420 GN = GridListAtom[Mc_AN][Nc];
1421 for (i=0; i<NO0; i++){
1422 RMO_Grid_tmp[GN] += LUMOs_Coef[0][spin][orbit][Gc_AN][i].r*
1423 Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1424 IMO_Grid_tmp[GN] += LUMOs_Coef[0][spin][orbit][Gc_AN][i].i*
1425 Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1426 }
1427 }
1428 }
1429
1430 MPI_Reduce(&RMO_Grid_tmp[0], &RMO_Grid[0], TNumGrid, MPI_DOUBLE,
1431 MPI_SUM, Host_ID, mpi_comm_level1);
1432 MPI_Reduce(&IMO_Grid_tmp[0], &IMO_Grid[0], TNumGrid, MPI_DOUBLE,
1433 MPI_SUM, Host_ID, mpi_comm_level1);
1434
1435 if (myid==Host_ID){
1436
1437 /* output the real part of LUMOs on grids */
1438
1439 sprintf(file1,"%s%s.lumo%i_%i_r%s",
1440 filepath,filename,spin,orbit,CUBE_EXTENSION);
1441
1442 if ((fp = fopen(file1,"wb")) != NULL){
1443
1444 #ifdef xt3
1445 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1446 #endif
1447
1448 Print_CubeTitle(fp,1,LUMOs_Coef[0][spin][orbit][0][0].r);
1449 fwrite(RMO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1450 fclose(fp);
1451 }
1452 else{
1453 printf("Failure of saving MOs\n");
1454 }
1455
1456 /* output the imaginary part of HOMOs on grids */
1457
1458 sprintf(file1,"%s%s.lumo%i_%i_i%s",
1459 filepath,filename,spin,orbit,CUBE_EXTENSION);
1460
1461 if ((fp = fopen(file1,"w")) != NULL){
1462
1463 #ifdef xt3
1464 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1465 #endif
1466
1467 Print_CubeTitle(fp,1,LUMOs_Coef[0][spin][orbit][0][0].r);
1468 fwrite(IMO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1469 fclose(fp);
1470 }
1471 else{
1472 printf("Failure of saving MOs\n");
1473 }
1474
1475 }
1476 } /* orbit */
1477 } /* spin */
1478
1479 /****************************************************
1480 freeing of arrays:
1481
1482 double RMO_Grid_tmp[TNumGrid];
1483 double IMO_Grid_tmp[TNumGrid];
1484 double RMO_Grid[TNumGrid];
1485 double IMO_Grid[TNumGrid];
1486 ****************************************************/
1487
1488 free(RMO_Grid_tmp);
1489 free(IMO_Grid_tmp);
1490 free(RMO_Grid);
1491 free(IMO_Grid);
1492 }
1493
1494
1495
1496
out_Bulk_MO()1497 void out_Bulk_MO()
1498 {
1499 int Mc_AN,Gc_AN,Cwan,NO0,spin,Nc;
1500 int kloop,Mh_AN,h_AN,Gh_AN,Rnh,Hwan;
1501 int NO1,l1,l2,l3,Nog,RnG,Nh,Rn;
1502 int orbit,GN,spe,i,j,i1,i2,i3,spinmax;
1503 double co,si,k1,k2,k3,kRn,ReCoef,ImCoef;
1504 double *RMO_Grid;
1505 double *IMO_Grid;
1506 double *RMO_Grid_tmp;
1507 double *IMO_Grid_tmp;
1508 char file1[YOUSO10];
1509 char buf[fp_bsize]; /* setvbuf */
1510 FILE *fp;
1511 int numprocs,myid,ID;
1512
1513 /* MPI */
1514 MPI_Comm_size(mpi_comm_level1,&numprocs);
1515 MPI_Comm_rank(mpi_comm_level1,&myid);
1516
1517 /****************************************************
1518 allocation of arrays:
1519
1520 double RMO_Grid[TNumGrid];
1521 double IMO_Grid[TNumGrid];
1522 double RMO_Grid_tmp[TNumGrid];
1523 double IMO_Grid_tmp[TNumGrid];
1524 ****************************************************/
1525
1526 RMO_Grid = (double*)malloc(sizeof(double)*TNumGrid);
1527 IMO_Grid = (double*)malloc(sizeof(double)*TNumGrid);
1528 RMO_Grid_tmp = (double*)malloc(sizeof(double)*TNumGrid);
1529 IMO_Grid_tmp = (double*)malloc(sizeof(double)*TNumGrid);
1530
1531 if (SpinP_switch==0) spinmax = 0;
1532 else if (SpinP_switch==1) spinmax = 1;
1533 else if (SpinP_switch==3) spinmax = 1;
1534
1535 /*************
1536 HOMOs
1537 *************/
1538
1539 for (kloop=0; kloop<MO_Nkpoint; kloop++){
1540
1541 k1 = MO_kpoint[kloop][1];
1542 k2 = MO_kpoint[kloop][2];
1543 k3 = MO_kpoint[kloop][3];
1544
1545 for (spin=0; spin<=spinmax; spin++){
1546 for (orbit=0; orbit<Bulk_Num_HOMOs[kloop]; orbit++){
1547
1548 /* calc. MO on grids */
1549
1550 for (GN=0; GN<TNumGrid; GN++){
1551 RMO_Grid_tmp[GN] = 0.0;
1552 IMO_Grid_tmp[GN] = 0.0;
1553 }
1554
1555 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1556 Gc_AN = M2G[Mc_AN];
1557 Cwan = WhatSpecies[Gc_AN];
1558 NO0 = Spe_Total_CNO[Cwan];
1559
1560 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
1561 GN = GridListAtom[Mc_AN][Nc];
1562 Rn = CellListAtom[Mc_AN][Nc];
1563
1564 l1 =-atv_ijk[Rn][1];
1565 l2 =-atv_ijk[Rn][2];
1566 l3 =-atv_ijk[Rn][3];
1567
1568 kRn = k1*(double)l1 + k2*(double)l2 + k3*(double)l3;
1569 si = sin(2.0*PI*kRn);
1570 co = cos(2.0*PI*kRn);
1571
1572 for (i=0; i<NO0; i++){
1573 ReCoef = co*HOMOs_Coef[kloop][spin][orbit][Gc_AN][i].r
1574 -si*HOMOs_Coef[kloop][spin][orbit][Gc_AN][i].i;
1575 ImCoef = co*HOMOs_Coef[kloop][spin][orbit][Gc_AN][i].i
1576 +si*HOMOs_Coef[kloop][spin][orbit][Gc_AN][i].r;
1577 RMO_Grid_tmp[GN] += ReCoef*Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1578 IMO_Grid_tmp[GN] += ImCoef*Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1579 }
1580 }
1581 }
1582
1583 MPI_Reduce(&RMO_Grid_tmp[0], &RMO_Grid[0], TNumGrid, MPI_DOUBLE,
1584 MPI_SUM, Host_ID, mpi_comm_level1);
1585 MPI_Reduce(&IMO_Grid_tmp[0], &IMO_Grid[0], TNumGrid, MPI_DOUBLE,
1586 MPI_SUM, Host_ID, mpi_comm_level1);
1587
1588 if (myid==Host_ID){
1589
1590 /* output the real part of HOMOs on grids */
1591
1592 sprintf(file1,"%s%s.homo%i_%i_%i_r%s",
1593 filepath,filename,kloop,spin,orbit,CUBE_EXTENSION);
1594
1595 if ((fp = fopen(file1,"wb")) != NULL){
1596
1597 #ifdef xt3
1598 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1599 #endif
1600
1601 Print_CubeTitle(fp,1,HOMOs_Coef[kloop][spin][orbit][0][0].r);
1602 fwrite(RMO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1603 fclose(fp);
1604 }
1605 else{
1606 printf("Failure of saving MOs\n");
1607 }
1608
1609 /* output the imaginary part of HOMOs on grids */
1610
1611 sprintf(file1,"%s%s.homo%i_%i_%i_i%s",
1612 filepath,filename,kloop,spin,orbit,CUBE_EXTENSION);
1613
1614 if ((fp = fopen(file1,"wb")) != NULL){
1615
1616 #ifdef xt3
1617 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1618 #endif
1619
1620 Print_CubeTitle(fp,1,HOMOs_Coef[kloop][spin][orbit][0][0].r);
1621 fwrite(IMO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1622 fclose(fp);
1623 }
1624 else{
1625 printf("Failure of saving MOs\n");
1626 }
1627 }
1628
1629 /* MPI_Barrier */
1630 MPI_Barrier(mpi_comm_level1);
1631
1632 } /* orbit */
1633 } /* spin */
1634 } /* kloop */
1635
1636 /*************
1637 LUMOs
1638 *************/
1639
1640 for (kloop=0; kloop<MO_Nkpoint; kloop++){
1641
1642 k1 = MO_kpoint[kloop][1];
1643 k2 = MO_kpoint[kloop][2];
1644 k3 = MO_kpoint[kloop][3];
1645
1646 for (spin=0; spin<=spinmax; spin++){
1647 for (orbit=0; orbit<Bulk_Num_LUMOs[kloop]; orbit++){
1648
1649 /* calc. MO on grids */
1650
1651 for (GN=0; GN<TNumGrid; GN++){
1652 RMO_Grid_tmp[GN] = 0.0;
1653 IMO_Grid_tmp[GN] = 0.0;
1654 }
1655
1656 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1657 Gc_AN = M2G[Mc_AN];
1658 Cwan = WhatSpecies[Gc_AN];
1659 NO0 = Spe_Total_CNO[Cwan];
1660
1661 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
1662 GN = GridListAtom[Mc_AN][Nc];
1663 Rn = CellListAtom[Mc_AN][Nc];
1664
1665 l1 =-atv_ijk[Rn][1];
1666 l2 =-atv_ijk[Rn][2];
1667 l3 =-atv_ijk[Rn][3];
1668
1669 kRn = k1*(double)l1 + k2*(double)l2 + k3*(double)l3;
1670 si = sin(2.0*PI*kRn);
1671 co = cos(2.0*PI*kRn);
1672
1673 for (i=0; i<NO0; i++){
1674 ReCoef = co*LUMOs_Coef[kloop][spin][orbit][Gc_AN][i].r
1675 -si*LUMOs_Coef[kloop][spin][orbit][Gc_AN][i].i;
1676 ImCoef = co*LUMOs_Coef[kloop][spin][orbit][Gc_AN][i].i
1677 +si*LUMOs_Coef[kloop][spin][orbit][Gc_AN][i].r;
1678 RMO_Grid_tmp[GN] += ReCoef*Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1679 IMO_Grid_tmp[GN] += ImCoef*Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
1680 }
1681 }
1682 }
1683
1684 MPI_Reduce(&RMO_Grid_tmp[0], &RMO_Grid[0], TNumGrid, MPI_DOUBLE,
1685 MPI_SUM, Host_ID, mpi_comm_level1);
1686 MPI_Reduce(&IMO_Grid_tmp[0], &IMO_Grid[0], TNumGrid, MPI_DOUBLE,
1687 MPI_SUM, Host_ID, mpi_comm_level1);
1688
1689 /* output the real part of LUMOs on grids */
1690
1691 if (myid==Host_ID){
1692
1693 sprintf(file1,"%s%s.lumo%i_%i_%i_r%s",
1694 filepath,filename,kloop,spin,orbit,CUBE_EXTENSION);
1695
1696 if ((fp = fopen(file1,"wb")) != NULL){
1697
1698 #ifdef xt3
1699 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1700 #endif
1701
1702 Print_CubeTitle(fp,1,LUMOs_Coef[kloop][spin][orbit][0][0].r);
1703 fwrite(RMO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1704 fclose(fp);
1705 }
1706 else{
1707 printf("Failure of saving MOs\n");
1708 fclose(fp);
1709 }
1710
1711 /* output the imaginary part of LUMOs on grids */
1712
1713 sprintf(file1,"%s%s.lumo%i_%i_%i_i%s",
1714 filepath,filename,kloop,spin,orbit,CUBE_EXTENSION);
1715
1716 if ((fp = fopen(file1,"wb")) != NULL){
1717
1718 #ifdef xt3
1719 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
1720 #endif
1721
1722 Print_CubeTitle(fp,1,LUMOs_Coef[kloop][spin][orbit][0][0].r);
1723 fwrite(IMO_Grid, sizeof(double), Ngrid1*Ngrid2*Ngrid3, fp);
1724 fclose(fp);
1725 }
1726 else{
1727 printf("Failure of saving MOs\n");
1728 }
1729 }
1730
1731 /* MPI_Barrier */
1732 MPI_Barrier(mpi_comm_level1);
1733
1734 } /* orbit */
1735 } /* spin */
1736 } /* kloop */
1737
1738
1739 /****************************************************
1740 allocation of arrays:
1741
1742 double RMO_Grid[TNumGrid];
1743 double IMO_Grid[TNumGrid];
1744 double RMO_Grid_tmp[TNumGrid];
1745 double IMO_Grid_tmp[TNumGrid];
1746 ****************************************************/
1747
1748 free(RMO_Grid);
1749 free(IMO_Grid);
1750 free(RMO_Grid_tmp);
1751 free(IMO_Grid_tmp);
1752 }
1753
1754
1755
1756
Print_CubeTitle(FILE * fp,int EigenValue_flag,double EigenValue)1757 static void Print_CubeTitle(FILE *fp, int EigenValue_flag, double EigenValue)
1758 {
1759 int ct_AN,n;
1760 int spe;
1761 char clist[300];
1762 double *dlist;
1763
1764 if (EigenValue_flag==0){
1765 sprintf(clist," SYS1\n");
1766 fwrite(clist, sizeof(char), 200, fp);
1767 sprintf(clist," SYS1\n");
1768 fwrite(clist, sizeof(char), 200, fp);
1769 }
1770 else {
1771 sprintf(clist," Absolute eigenvalue=%10.7f (Hartree) Relative eigenvalue=%10.7f (Hartree)\n",
1772 EigenValue,EigenValue-ChemP);
1773 fwrite(clist, sizeof(char), 200, fp);
1774
1775 sprintf(clist, " Chemical Potential=%10.7f (Hartree)\n",ChemP);
1776 fwrite(clist, sizeof(char), 200, fp);
1777 }
1778
1779 fwrite(&atomnum, sizeof(int), 1, fp);
1780 fwrite(&Grid_Origin[1], sizeof(double), 3, fp);
1781 fwrite(&Ngrid1, sizeof(int), 1, fp);
1782 fwrite(>v[1][1], sizeof(double), 3, fp);
1783 fwrite(&Ngrid2, sizeof(int), 1, fp);
1784 fwrite(>v[2][1], sizeof(double), 3, fp);
1785 fwrite(&Ngrid3, sizeof(int), 1, fp);
1786 fwrite(>v[3][1], sizeof(double), 3, fp);
1787
1788 dlist = (double*)malloc(sizeof(double)*atomnum*5);
1789
1790 n = 0;
1791 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
1792 spe = WhatSpecies[ct_AN];
1793 dlist[n] = (double)Spe_WhatAtom[spe]; n++;
1794 dlist[n] = Spe_Core_Charge[spe]-InitN_USpin[ct_AN]-InitN_DSpin[ct_AN]; n++;
1795 dlist[n] = Gxyz[ct_AN][1]; n++;
1796 dlist[n] = Gxyz[ct_AN][2]; n++;
1797 dlist[n] = Gxyz[ct_AN][3]; n++;
1798 }
1799
1800 fwrite(dlist, sizeof(double), atomnum*5, fp);
1801 free(dlist);
1802 }
1803
1804
1805
1806
Print_CubeData(FILE * fp,char fext[],double * data0,double * data1,char * op)1807 static void Print_CubeData(FILE *fp, char fext[], double *data0, double *data1, char *op)
1808 {
1809 int myid,numprocs,ID,BN_AB,n3,cmd,digit,fd,c;
1810 char operate[300],operate1[300],operate2[300];
1811 FILE *fp1,*fp2;
1812
1813 MPI_Comm_size(mpi_comm_level1,&numprocs);
1814 MPI_Comm_rank(mpi_comm_level1,&myid);
1815 digit = (int)log10(numprocs) + 1;
1816
1817 /****************************************************
1818 output data
1819 ****************************************************/
1820
1821 if (op==NULL){
1822 fwrite(data0, sizeof(double), My_NumGridB_AB, fp);
1823 }
1824 else if (strcmp(op,"add")==0){
1825
1826 for (BN_AB=0; BN_AB<My_NumGridB_AB; BN_AB++){
1827 data0[BN_AB] += data1[BN_AB];
1828 }
1829
1830 fwrite(data0, sizeof(double), My_NumGridB_AB, fp);
1831
1832 for (BN_AB=0; BN_AB<My_NumGridB_AB; BN_AB++){
1833 data0[BN_AB] -= data1[BN_AB];
1834 }
1835 }
1836
1837 else if (strcmp(op,"diff")==0){
1838
1839 for (BN_AB=0; BN_AB<My_NumGridB_AB; BN_AB++){
1840 data0[BN_AB] -= data1[BN_AB];
1841 }
1842
1843 fwrite(data0, sizeof(double), My_NumGridB_AB, fp);
1844
1845 for (BN_AB=0; BN_AB<My_NumGridB_AB; BN_AB++){
1846 data0[BN_AB] += data1[BN_AB];
1847 }
1848
1849 }
1850 else {
1851 printf("Print_CubeData: op=%s not supported\n",op);
1852 return;
1853 }
1854
1855 /****************************************************
1856 fclose(fp);
1857 ****************************************************/
1858
1859 fd = fileno(fp);
1860 fsync(fd);
1861 fclose(fp);
1862
1863 /****************************************************
1864 merge files
1865 ****************************************************/
1866
1867 MPI_Barrier(mpi_comm_level1);
1868
1869 if (myid==Host_ID){
1870
1871 if (CAT_FLAG) {
1872
1873 sprintf(operate,"cat %s%s%s_* > %s%s%s",
1874 filepath,filename,fext,filepath,filename,fext);
1875 system(operate);
1876
1877 }
1878
1879 else {
1880
1881 /* check whether the file exists, and if it exists, remove it. */
1882
1883 sprintf(operate1,"%s%s%s",filepath,filename,fext);
1884 fp1 = fopen(operate1, "rb");
1885 if (fp1!=NULL){
1886 fclose(fp1);
1887 remove(operate1);
1888 }
1889
1890 /* merge all the fraction files */
1891
1892 for (ID=0; ID<numprocs; ID++){
1893
1894 sprintf(operate1,"%s%s%s",filepath,filename,fext);
1895 fp1 = fopen(operate1, "ab");
1896 fseek(fp1,0,SEEK_END);
1897
1898 sprintf(operate2,"%s%s%s_%0*i",filepath,filename,fext,digit,ID);
1899 fp2 = fopen(operate2,"rb");
1900
1901 if (fp2!=NULL){
1902 for (c=getc(fp2); c!=EOF; c=getc(fp2)) putc(c,fp1);
1903 fclose(fp2);
1904 }
1905
1906 fclose(fp1);
1907 }
1908 }
1909
1910 /* remove all the fraction files */
1911 sprintf(operate,"rm %s%s%s_*",filepath,filename,fext);
1912 system(operate);
1913
1914 }
1915
1916 }
1917
1918
Print_VectorData(FILE * fp,char fext[],double * data0,double * data1,double * data2,double * data3)1919 static void Print_VectorData(FILE *fp, char fext[],
1920 double *data0, double *data1,
1921 double *data2, double *data3)
1922 {
1923 int i,j,k,i1,i2,i3,n,c,GridNum;
1924 int GN_AB,n1,n2,n3,interval,digit;
1925 int BN_AB,N2D,GNs,GN;
1926 double x,y,z,vx,vy,vz;
1927 double sden,Cxyz[4],theta,phi;
1928 double tv0[4][4];
1929 double *dlist;
1930 char operate[300];
1931 char operate1[300];
1932 char operate2[300];
1933 char clist[300];
1934 int numprocs,myid,ID;
1935 MPI_Status stat;
1936 MPI_Request request;
1937 FILE *fp1,*fp2;
1938
1939 /* MPI */
1940 MPI_Comm_size(mpi_comm_level1,&numprocs);
1941 MPI_Comm_rank(mpi_comm_level1,&myid);
1942 digit = (int)log10(numprocs) + 1;
1943
1944 /****************************************************
1945 output data
1946 ****************************************************/
1947
1948 /* for XCrysDen */
1949
1950 interval = 4;
1951
1952 if (myid==Host_ID){
1953
1954 sprintf(clist,"CRYSTAL\n");
1955 fwrite(clist, sizeof(char), 200, fp);
1956 sprintf(clist,"PRIMVEC\n");
1957 fwrite(clist, sizeof(char), 200, fp);
1958
1959 for (i=1; i<=3; i++){
1960 for (j=1; j<=3; j++){
1961 tv0[i][j] = BohrR*tv[i][j];
1962 }
1963 }
1964
1965 fwrite(&tv0[1][1], sizeof(double), 3, fp);
1966 fwrite(&tv0[2][1], sizeof(double), 3, fp);
1967 fwrite(&tv0[3][1], sizeof(double), 3, fp);
1968
1969 GridNum = 0;
1970 for (n1=0; n1<Ngrid1; n1++){
1971 for (n2=0; n2<Ngrid2; n2++){
1972 for (n3=0; n3<Ngrid3; n3++){
1973 if (n1%interval==0 && n2%interval==0 && n3%interval==0) GridNum++;
1974 }
1975 }
1976 }
1977
1978 fwrite(&atomnum,sizeof(int),1, fp);
1979 fwrite(&GridNum,sizeof(int),1, fp);
1980
1981 for (k=1; k<=atomnum; k++){
1982 i = WhatSpecies[k];
1983 j = Spe_WhatAtom[i];
1984
1985 fwrite(Atom_Symbol[j], sizeof(char), 4, fp);
1986
1987 Cxyz[1] = Gxyz[k][1]*BohrR;
1988 Cxyz[2] = Gxyz[k][2]*BohrR;
1989 Cxyz[3] = Gxyz[k][3]*BohrR;
1990
1991 fwrite(&Cxyz[1], sizeof(double), 3, fp);
1992 }
1993 }
1994
1995 /****************************************************
1996 fwrite vector data
1997 ****************************************************/
1998
1999 dlist = (double*)malloc(sizeof(double)*My_NumGridB_AB*6);
2000
2001 N2D = Ngrid1*Ngrid2;
2002 GNs = ((myid*N2D+numprocs-1)/numprocs)*Ngrid3;
2003
2004 n = 0;
2005 for (BN_AB=0; BN_AB<My_NumGridB_AB; BN_AB++){
2006
2007 GN_AB = BN_AB + GNs;
2008 n1 = GN_AB/(Ngrid2*Ngrid3);
2009 n2 = (GN_AB - n1*(Ngrid2*Ngrid3))/Ngrid3;
2010 n3 = GN_AB - n1*(Ngrid2*Ngrid3) - n2*Ngrid3;
2011
2012 if (n1%interval==0 && n2%interval==0 && n3%interval==0){
2013
2014 GN = n1*Ngrid2*Ngrid3 + n2*Ngrid3 + n3;
2015
2016 Get_Grid_XYZ(GN,Cxyz);
2017 x = BohrR*Cxyz[1];
2018 y = BohrR*Cxyz[2];
2019 z = BohrR*Cxyz[3];
2020
2021 sden = data0[BN_AB] - data1[BN_AB];
2022 theta = data2[BN_AB];
2023 phi = data3[BN_AB];
2024
2025 vx = sden*sin(theta)*cos(phi);
2026 vy = sden*sin(theta)*sin(phi);
2027 vz = sden*cos(theta);
2028
2029 dlist[n ] = x;
2030 dlist[n+1] = y;
2031 dlist[n+2] = z;
2032 dlist[n+3] = vx;
2033 dlist[n+4] = vy;
2034 dlist[n+5] = vz;
2035
2036 n += 6;
2037 }
2038 }
2039
2040 fwrite(dlist, sizeof(double), n, fp);
2041 free(dlist);
2042
2043 /****************************************************
2044 fclose(fp);
2045 ****************************************************/
2046
2047 fclose(fp);
2048
2049 /****************************************************
2050 merge files
2051 ****************************************************/
2052
2053 MPI_Barrier(mpi_comm_level1);
2054
2055 if (myid==Host_ID){
2056
2057 if (CAT_FLAG) {
2058
2059 sprintf(operate,"cat %s%s%s_* > %s%s%s",
2060 filepath,filename,fext,filepath,filename,fext);
2061 system(operate);
2062
2063 }
2064
2065 else {
2066
2067 /* check whether the file exists, and if it exists, remove it. */
2068
2069 sprintf(operate1,"%s%s%s",filepath,filename,fext);
2070 fp1 = fopen(operate1, "rb");
2071 if (fp1!=NULL){
2072 fclose(fp1);
2073 remove(operate1);
2074 }
2075
2076 /* merge all the fraction files */
2077
2078 for (ID=0; ID<numprocs; ID++){
2079
2080 sprintf(operate1,"%s%s%s",filepath,filename,fext);
2081 fp1 = fopen(operate1, "ab");
2082 fseek(fp1,0,SEEK_END);
2083
2084 sprintf(operate2,"%s%s%s_%0*i",filepath,filename,fext,digit,ID);
2085 fp2 = fopen(operate2,"rb");
2086
2087 if (fp2!=NULL){
2088 for (c=getc(fp2); c!=EOF; c=getc(fp2)) putc(c,fp1);
2089 fclose(fp2);
2090 }
2091
2092 fclose(fp1);
2093 }
2094 }
2095
2096 sprintf(operate,"rm %s%s%s_*",filepath,filename,fext);
2097 system(operate);
2098
2099 }
2100 }
2101
2102
2103
2104
2105
out_OrbOpt(char * inputfile)2106 void out_OrbOpt(char *inputfile)
2107 {
2108 int i,j,natom,po,al,p,wan,L0,Mul0,M0;
2109 int num,Mc_AN,Gc_AN,pmax,Mul1;
2110 int al1,al0;
2111 double sum,Max0;
2112 double ***tmp_coes;
2113 char file1[YOUSO10] = ".oopt";
2114 char DirPAO[YOUSO10];
2115 char ExtPAO[YOUSO10] = ".pao";
2116 char FN_PAO[YOUSO10];
2117 char EndLine[YOUSO10] = "<pseudo.atomic.orbitals.L=0";
2118 char fname2[YOUSO10];
2119 char command0[YOUSO10];
2120 double tmp0,tmp1;
2121 double *Tmp_Vec,*Tmp_CntCoes;
2122 int ***tmp_index2;
2123 char *sp0;
2124 FILE *fp,*fp2;
2125 char *buf; /* setvbuf */
2126 int numprocs,myid,ID,tag=999;
2127
2128 MPI_Status stat;
2129 MPI_Request request;
2130
2131 /* MPI */
2132 MPI_Comm_size(mpi_comm_level1,&numprocs);
2133 MPI_Comm_rank(mpi_comm_level1,&myid);
2134
2135 /* set DirPAO */
2136
2137 sprintf(DirPAO,"%s/PAO/","/usr/local/share/openmx/DFT_DATA13");
2138
2139 /* allocate buf */
2140 buf = malloc(fp_bsize); /* setvbuf */
2141
2142 /* allocation of Tmp_Vec */
2143 Tmp_Vec = (double*)malloc(sizeof(double)*List_YOUSO[7]*List_YOUSO[24]);
2144
2145 fnjoint(filepath,filename,file1);
2146
2147 if ((fp = fopen(file1,"w")) != NULL){
2148
2149 #ifdef xt3
2150 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
2151 #endif
2152
2153 Tmp_CntCoes = (double*)malloc(sizeof(double)*List_YOUSO[24]);
2154
2155 tmp_index2 = (int***)malloc(sizeof(int**)*(List_YOUSO[25]+1));
2156 for (i=0; i<(List_YOUSO[25]+1); i++){
2157 tmp_index2[i] = (int**)malloc(sizeof(int*)*List_YOUSO[24]);
2158 for (j=0; j<List_YOUSO[24]; j++){
2159 tmp_index2[i][j] = (int*)malloc(sizeof(int)*(2*(List_YOUSO[25]+1)+1));
2160 }
2161 }
2162
2163 /**********************************************
2164 transformation of optimized orbitals by
2165 an extended Gauss elimination and
2166 the Gram-Schmidt orthogonalization
2167 ***********************************************/
2168
2169 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2170
2171 Gc_AN = M2G[Mc_AN];
2172 wan = WhatSpecies[Gc_AN];
2173
2174 al = -1;
2175 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
2176 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
2177 for (M0=0; M0<=2*L0; M0++){
2178 al++;
2179 tmp_index2[L0][Mul0][M0] = al;
2180 }
2181 }
2182 }
2183
2184 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
2185 for (M0=0; M0<=2*L0; M0++){
2186
2187 /**********************************************
2188 extended Gauss elimination
2189 ***********************************************/
2190
2191 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
2192 al0 = tmp_index2[L0][Mul0][M0];
2193 for (Mul1=0; Mul1<Spe_Num_CBasis[wan][L0]; Mul1++){
2194 al1 = tmp_index2[L0][Mul1][M0];
2195
2196 if (Mul1!=Mul0){
2197
2198 tmp0 = CntCoes[Mc_AN][al0][Mul0];
2199 tmp1 = CntCoes[Mc_AN][al1][Mul0];
2200
2201 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
2202 CntCoes[Mc_AN][al1][p] -= CntCoes[Mc_AN][al0][p]/tmp0*tmp1;
2203 }
2204 }
2205
2206 }
2207 }
2208
2209 /*
2210 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
2211 al0 = tmp_index2[L0][Mul0][M0];
2212 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
2213 printf("Mul0=%2d p=%2d Coes=%15.12f\n",Mul0,p,CntCoes[Mc_AN][al0][p]);
2214 }
2215 }
2216 */
2217
2218 /**********************************************
2219 orthonormalization of optimized orbitals
2220 ***********************************************/
2221
2222 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
2223 al0 = tmp_index2[L0][Mul0][M0];
2224
2225 /* x - sum_i <x|e_i>e_i */
2226
2227 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
2228 Tmp_CntCoes[p] = 0.0;
2229 }
2230
2231 for (Mul1=0; Mul1<Mul0; Mul1++){
2232 al1 = tmp_index2[L0][Mul1][M0];
2233
2234 sum = 0.0;
2235 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
2236 sum = sum + CntCoes[Mc_AN][al0][p]*CntCoes[Mc_AN][al1][p];
2237 }
2238
2239 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
2240 Tmp_CntCoes[p] = Tmp_CntCoes[p] + sum*CntCoes[Mc_AN][al1][p];
2241 }
2242 }
2243
2244 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
2245 CntCoes[Mc_AN][al0][p] = CntCoes[Mc_AN][al0][p] - Tmp_CntCoes[p];
2246 }
2247
2248 /* Normalize */
2249
2250 sum = 0.0;
2251 Max0 = -100.0;
2252 pmax = 0;
2253 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
2254 sum = sum + CntCoes[Mc_AN][al0][p]*CntCoes[Mc_AN][al0][p];
2255 if (Max0<fabs(CntCoes[Mc_AN][al0][p])){
2256 Max0 = fabs(CntCoes[Mc_AN][al0][p]);
2257 pmax = p;
2258 }
2259 }
2260
2261 if (fabs(sum)<1.0e-11)
2262 tmp0 = 0.0;
2263 else
2264 tmp0 = 1.0/sqrt(sum);
2265
2266 tmp1 = sgn(CntCoes[Mc_AN][al0][pmax]);
2267
2268 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
2269 CntCoes[Mc_AN][al0][p] = tmp0*tmp1*CntCoes[Mc_AN][al0][p];
2270 }
2271
2272 }
2273 }
2274 }
2275
2276 } /* Mc_AN */
2277
2278 for (i=0; i<(List_YOUSO[25]+1); i++){
2279 for (j=0; j<List_YOUSO[24]; j++){
2280 free(tmp_index2[i][j]);
2281 }
2282 free(tmp_index2[i]);
2283 }
2284 free(tmp_index2);
2285
2286 free(Tmp_CntCoes);
2287
2288 /**********************************************
2289 set Tmp_Vec
2290 ***********************************************/
2291
2292 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
2293
2294 wan = WhatSpecies[Gc_AN];
2295 ID = G2ID[Gc_AN];
2296
2297 if (myid==ID){
2298
2299 Mc_AN = F_G2M[Gc_AN];
2300
2301 al = -1;
2302 num = 0;
2303 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
2304 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
2305 for (M0=0; M0<=2*L0; M0++){
2306 al++;
2307 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
2308 Tmp_Vec[num] = CntCoes[Mc_AN][al][p];
2309 num++;
2310 }
2311 }
2312 }
2313 }
2314
2315 if (myid!=Host_ID){
2316 MPI_Isend(&num, 1, MPI_INT, Host_ID, tag, mpi_comm_level1, &request);
2317 MPI_Wait(&request,&stat);
2318 MPI_Isend(&Tmp_Vec[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1, &request);
2319 MPI_Wait(&request,&stat);
2320 }
2321
2322 }
2323
2324 else if (ID!=myid && myid==Host_ID){
2325 MPI_Recv(&num, 1, MPI_INT, ID, tag, mpi_comm_level1, &stat);
2326 MPI_Recv(&Tmp_Vec[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
2327 }
2328
2329 /**********************************************
2330 write
2331 ***********************************************/
2332
2333 if (myid==Host_ID){
2334
2335 fprintf(fp,"\nAtom=%2d\n",Gc_AN);
2336 fprintf(fp,"Basis specification %s\n",SpeBasis[wan]);
2337
2338 fprintf(fp,"Contraction coefficients p=");
2339 for (i=0; i<List_YOUSO[24]; i++){
2340 fprintf(fp," %i ",i);
2341 }
2342 fprintf(fp,"\n");
2343
2344 al = -1;
2345 num = 0;
2346
2347 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
2348 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
2349 for (M0=0; M0<=2*L0; M0++){
2350 al++;
2351
2352 fprintf(fp,"Atom=%3d L=%2d Mul=%2d M=%2d ",Gc_AN,L0,Mul0,M0);
2353 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
2354 fprintf(fp,"%9.5f ",Tmp_Vec[num]);
2355 num++;
2356 }
2357 for (p=Spe_Specified_Num[wan][al]; p<List_YOUSO[24]; p++){
2358 fprintf(fp," 0.00000 ");
2359 }
2360 fprintf(fp,"\n");
2361 }
2362 }
2363 }
2364 } /* if (myid==Host_ID) */
2365 } /* Gc_AN */
2366 fclose(fp);
2367 }
2368 else{
2369 printf("Failure of saving contraction coefficients of basis orbitals\n");
2370 }
2371
2372 /****************************************************
2373 outputting of contracted orbitals as *.pao
2374 ****************************************************/
2375
2376 if (CntOrb_fileout==1 && Cnt_switch==1 && RCnt_switch==1){
2377
2378 /****************************************************
2379 allocation of array:
2380
2381 tmp_coes[List_YOUSO[25]+1]
2382 [List_YOUSO[24]]
2383 [List_YOUSO[24]]
2384 ****************************************************/
2385
2386 tmp_coes = (double***)malloc(sizeof(double**)*(List_YOUSO[25]+1));
2387 for (i=0; i<(List_YOUSO[25]+1); i++){
2388 tmp_coes[i] = (double**)malloc(sizeof(double*)*List_YOUSO[24]);
2389 for (j=0; j<List_YOUSO[24]; j++){
2390 tmp_coes[i][j] = (double*)malloc(sizeof(double)*List_YOUSO[24]);
2391 }
2392 }
2393
2394 for (natom=0; natom<Num_CntOrb_Atoms; natom++){
2395
2396 Gc_AN = CntOrb_Atoms[natom];
2397 ID = G2ID[Gc_AN];
2398 wan = WhatSpecies[Gc_AN];
2399
2400 if (myid==ID){
2401
2402 Mc_AN = F_G2M[Gc_AN];
2403
2404 al = -1;
2405 num = 0;
2406 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
2407 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
2408 for (M0=0; M0<=2*L0; M0++){
2409 al++;
2410 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
2411 Tmp_Vec[num] = CntCoes[Mc_AN][al][p];
2412 num++;
2413 }
2414 }
2415 }
2416 }
2417
2418 if (myid!=Host_ID){
2419 MPI_Isend(&num, 1, MPI_INT, Host_ID, tag, mpi_comm_level1, &request);
2420 MPI_Wait(&request,&stat);
2421 MPI_Isend(&Tmp_Vec[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1, &request);
2422 MPI_Wait(&request,&stat);
2423 }
2424
2425 }
2426
2427 else if (ID!=myid && myid==Host_ID){
2428 MPI_Recv(&num, 1, MPI_INT, ID, tag, mpi_comm_level1, &stat);
2429 MPI_Recv(&Tmp_Vec[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
2430 }
2431
2432 /****************************************************
2433 generate a pao file
2434 ****************************************************/
2435
2436 if (myid==Host_ID){
2437
2438 /*
2439 setting of initial vectors using contraction
2440 coefficients and unit vectors (1.0)
2441 */
2442
2443 al = -1;
2444 num = 0;
2445
2446 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
2447 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
2448 for (M0=0; M0<=2*L0; M0++){
2449 al++;
2450
2451 if (M0==0){
2452 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
2453 tmp_coes[L0][Mul0][p] = Tmp_Vec[num];
2454 num++;
2455 }
2456 for (p=Spe_Specified_Num[wan][al]; p<Spe_PAO_Mul[wan]; p++){
2457 tmp_coes[L0][Mul0][p] = 0.0;
2458 }
2459 }
2460 else{
2461 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
2462 num++;
2463 }
2464 }
2465 }
2466 }
2467
2468 for (Mul0=Spe_Num_CBasis[wan][L0]; Mul0<Spe_PAO_Mul[wan]; Mul0++){
2469 for (p=0; p<Spe_PAO_Mul[wan]; p++) tmp_coes[L0][Mul0][p] = 0.0;
2470 tmp_coes[L0][Mul0][Mul0] = 1.0;
2471 }
2472 }
2473
2474 /****************************************************
2475 make *.pao
2476 ****************************************************/
2477
2478 sprintf(fname2,"%s%s_%i%s",filepath,SpeName[wan],Gc_AN,ExtPAO);
2479 if ((fp2 = fopen(fname2,"w")) != NULL){
2480
2481 #ifdef xt3
2482 setvbuf(fp2,buf,_IOFBF,fp_bsize); /* setvbuf */
2483 #endif
2484
2485 fnjoint2(DirPAO,SpeBasisName[wan],ExtPAO,FN_PAO);
2486
2487 /****************************************************
2488 Log of the orbital optimization
2489 and
2490 contraction coefficients
2491 ****************************************************/
2492
2493 fprintf(fp2,"*****************************************************\n");
2494 fprintf(fp2,"*****************************************************\n");
2495 fprintf(fp2," The numerical atomic orbitals were generated\n");
2496 fprintf(fp2," by the variational optimization.\n");
2497 fprintf(fp2," The original file of PAOs was %s.\n",FN_PAO);
2498 fprintf(fp2," Basis specification was %s.\n",SpeBasis[wan]);
2499 fprintf(fp2," The input file was %s.\n",inputfile);
2500 fprintf(fp2,"*****************************************************\n");
2501 fprintf(fp2,"*****************************************************\n\n");
2502
2503 fprintf(fp2,"<Contraction.coefficients\n");
2504
2505 al = -1;
2506 num = 0;
2507
2508 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
2509 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
2510 for (M0=0; M0<=2*L0; M0++){
2511
2512 al++;
2513
2514 if (M0==0){
2515
2516 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
2517
2518 fprintf(fp2," Atom=%3d L=%2d Mul=%2d p=%3d %18.15f\n",
2519 Gc_AN,L0,Mul0,p,Tmp_Vec[num]);
2520
2521 num++;
2522 }
2523 for (p=Spe_Specified_Num[wan][al]; p<List_YOUSO[24]; p++){
2524
2525 fprintf(fp2," Atom=%3d L=%2d Mul=%2d p=%3d %18.15f\n",
2526 Gc_AN,L0,Mul0,p,0.0);
2527 }
2528 }
2529 else{
2530 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
2531 num++;
2532 }
2533 }
2534
2535 }
2536 }
2537 }
2538 fprintf(fp2,"Contraction.coefficients>\n");
2539
2540 fprintf(fp2,"\n*****************************************************\n");
2541 fprintf(fp2,"*****************************************************\n\n");
2542
2543 /****************************************************
2544 from the original file
2545 ****************************************************/
2546
2547 if ((fp = fopen(FN_PAO,"r")) != NULL){
2548
2549 po = 0;
2550 do{
2551 if (fgets(command0,YOUSO10,fp)!=NULL){
2552 command0[strlen(command0)-1] = '\0';
2553 sp0 = strstr(command0,EndLine);
2554 if (sp0!=NULL){
2555 po = 1;
2556 }
2557 else {
2558 fprintf(fp2,"%s\n",command0);
2559 }
2560 }
2561 else{
2562 po = 1;
2563 }
2564 }while(po==0);
2565
2566 fclose(fp);
2567 }
2568 else{
2569 printf("Could not find %s\n",FN_PAO);
2570 exit(1);
2571 }
2572
2573 /****************************************************
2574 contracted orbitals
2575 ****************************************************/
2576
2577 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
2578 fprintf(fp2,"<pseudo.atomic.orbitals.L=%d\n",L0);
2579 for (i=0; i<Spe_Num_Mesh_PAO[wan]; i++){
2580 fprintf(fp2,"%18.15f %18.15f ",Spe_PAO_XV[wan][i],Spe_PAO_RV[wan][i]);
2581 for (Mul0=0; Mul0<Spe_PAO_Mul[wan]; Mul0++){
2582 sum = 0.0;
2583 for (p=0; p<Spe_PAO_Mul[wan]; p++){
2584 sum += tmp_coes[L0][Mul0][p]*Spe_PAO_RWF[wan][L0][p][i];
2585 }
2586 fprintf(fp2,"%18.15f ",sum);
2587 }
2588 fprintf(fp2,"\n");
2589 }
2590 fprintf(fp2,"pseudo.atomic.orbitals.L=%d>\n",L0);
2591 }
2592
2593 for (L0=(Spe_MaxL_Basis[wan]+1); L0<=Spe_PAO_LMAX[wan]; L0++){
2594 fprintf(fp2,"<pseudo.atomic.orbitals.L=%d\n",L0);
2595 for (i=0; i<Spe_Num_Mesh_PAO[wan]; i++){
2596 fprintf(fp2,"%18.15f %18.15f ",Spe_PAO_XV[wan][i],Spe_PAO_RV[wan][i]);
2597 for (Mul0=0; Mul0<Spe_PAO_Mul[wan]; Mul0++){
2598 fprintf(fp2,"%18.15f ",Spe_PAO_RWF[wan][L0][Mul0][i]);
2599 }
2600 fprintf(fp2,"\n");
2601 }
2602 fprintf(fp2,"pseudo.atomic.orbitals.L=%d>\n",L0);
2603 }
2604
2605 fclose(fp2);
2606
2607 }
2608 else{
2609 printf("Could not open %s\n",fname2);
2610 exit(0);
2611 }
2612
2613 } /* if (myid==Host_ID) */
2614 } /* for (natom=0; natom<Num_CntOrb_Atoms; natom++) */
2615
2616 /****************************************************
2617 freeing of arrays:
2618 ****************************************************/
2619
2620 for (i=0; i<(List_YOUSO[25]+1); i++){
2621 for (j=0; j<List_YOUSO[24]; j++){
2622 free(tmp_coes[i][j]);
2623 }
2624 free(tmp_coes[i]);
2625 }
2626 free(tmp_coes);
2627
2628 }
2629
2630 free(Tmp_Vec);
2631 free(buf);
2632 }
2633
2634
2635
2636
Print_CubeData_MO(FILE * fp,double * data,double * data1,char * op)2637 static void Print_CubeData_MO(FILE *fp, double *data, double *data1,char *op)
2638 {
2639 int i1,i2,i3;
2640 int GN;
2641 int cmd;
2642
2643 if (op==NULL) { cmd=0; }
2644 else if (strcmp(op,"add")==0) { cmd=1; }
2645 else if (strcmp(op,"diff")==0) { cmd=2; }
2646 else {
2647 printf("Print_CubeData: op=%s not supported\n",op);
2648 return;
2649 }
2650
2651 for (i1=0; i1<Ngrid1; i1++){
2652 for (i2=0; i2<Ngrid2; i2++){
2653 for (i3=0; i3<Ngrid3; i3++){
2654 GN = i1*Ngrid2*Ngrid3 + i2*Ngrid3 + i3;
2655 switch (cmd) {
2656 case 0:
2657 fprintf(fp,"%13.3E",data[GN]);
2658 break;
2659 case 1:
2660 fprintf(fp,"%13.3E",data[GN]+data1[GN]);
2661 break;
2662 case 2:
2663 fprintf(fp,"%13.3E",data[GN]-data1[GN]);
2664 break;
2665 }
2666 if ((i3+1)%6==0) { fprintf(fp,"\n"); }
2667 }
2668 /* avoid double \n\n when Ngrid3%6 == 0 */
2669 if (Ngrid3%6!=0) fprintf(fp,"\n");
2670 }
2671 }
2672 }
2673
2674
2675
2676
2677
2678
out_Partial_Charge_Density()2679 void out_Partial_Charge_Density()
2680 {
2681 int ct_AN,spe,i1,i2,i3,GN,i,j,k;
2682 int MN,digit;
2683 double x,y,z,vx,vy,vz;
2684 double phi,theta,sden,oden;
2685 double xmin,ymin,zmin,xmax,ymax,zmax;
2686 char fname[YOUSO10];
2687 char file1[YOUSO10] = ".pden";
2688 char file2[YOUSO10] = ".pden0";
2689 char file3[YOUSO10] = ".pden1";
2690 FILE *fp;
2691 char buf[fp_bsize]; /* setvbuf */
2692 int numprocs,myid,ID;
2693
2694 /* MPI */
2695 MPI_Comm_size(mpi_comm_level1,&numprocs);
2696 MPI_Comm_rank(mpi_comm_level1,&myid);
2697 digit = (int)log10(numprocs) + 1;
2698
2699 strcat(file1,CUBE_EXTENSION);
2700 strcat(file2,CUBE_EXTENSION);
2701 strcat(file3,CUBE_EXTENSION);
2702
2703 /****************************************************
2704 calculation of partial charge density
2705 ****************************************************/
2706
2707 Set_Partial_Density_Grid(Partial_DM);
2708
2709 /****************************************************
2710 partial electron density including both up
2711 and down spin contributions
2712 ****************************************************/
2713
2714 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file1,digit,myid);
2715
2716 if ((fp = fopen(fname,"wb")) != NULL){
2717
2718 #ifdef xt3
2719 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
2720 #endif
2721
2722 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
2723
2724 if (SpinP_switch==0) {
2725 for (MN=0; MN<My_NumGridB_AB; MN++){
2726 Density_Grid_B[0][MN] = 2.0*Density_Grid_B[0][MN];
2727 }
2728 Print_CubeData(fp,file1,Density_Grid_B[0],(void*)NULL,(void*)NULL);
2729 }
2730 else {
2731 Print_CubeData(fp,file1,Density_Grid_B[0],Density_Grid_B[1],"add");
2732 }
2733
2734 }
2735 else{
2736 printf("Failure of saving total partial electron density\n");
2737 }
2738
2739 /* spin polization */
2740
2741 if (SpinP_switch==1){
2742
2743 /****************************************************
2744 up-spin electron density
2745 ****************************************************/
2746
2747 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file2,digit,myid);
2748
2749 if ((fp = fopen(fname,"wb")) != NULL){
2750
2751 #ifdef xt3
2752 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
2753 #endif
2754
2755 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
2756 Print_CubeData(fp,file2,Density_Grid_B[0],NULL,NULL);
2757 }
2758 else{
2759 printf("Failure of saving the electron density\n");
2760 }
2761
2762 /****************************************************
2763 down-spin electron density
2764 ****************************************************/
2765
2766 sprintf(fname,"%s%s%s_%0*i",filepath,filename,file3,digit,myid);
2767
2768 if ((fp = fopen(fname,"wb")) != NULL){
2769
2770 #ifdef xt3
2771 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
2772 #endif
2773
2774 if (myid==Host_ID) Print_CubeTitle(fp,0,0.0);
2775 Print_CubeData(fp,file3,Density_Grid_B[1],NULL,NULL);
2776 }
2777 else{
2778 printf("Failure of saving the electron density\n");
2779 }
2780 }
2781
2782 }
2783
2784
2785
2786
2787
Set_Partial_Density_Grid(double ***** CDM)2788 void Set_Partial_Density_Grid(double *****CDM)
2789 {
2790 int al,L0,Mul0,M0,p,size1,size2;
2791 int Gc_AN,Mc_AN,Mh_AN,MN;
2792 int n1,n2,n3,k1,k2,k3,N3[4];
2793 int Cwan,NO0,NO1,Rn,N,Hwan,i,j,k,n;
2794 int Max_Size,My_Max,top_num;
2795 double time0;
2796 int h_AN,Gh_AN,Rnh,spin,Nc,GRc,Nh,Nog;
2797 int Nc_0,Nc_1,Nc_2,Nc_3,Nh_0,Nh_1,Nh_2,Nh_3;
2798 unsigned long long int LN,BN,AN,n2D,N2D;
2799 unsigned long long int GNc,GN;
2800 double tmp0,tmp1,sk1,sk2,sk3,tot_den,sum;
2801 double tmp0_0,tmp0_1,tmp0_2,tmp0_3;
2802 double sum_0,sum_1,sum_2,sum_3;
2803
2804 double d1,d2,d3,cop,sip,sit,cot;
2805 double Re11,Re22,Re12,Im12,phi,theta;
2806 double x,y,z,Cxyz[4];
2807 double TStime,TEtime;
2808 double ***Tmp_Den_Grid;
2809 double **Den_Snd_Grid_A2B;
2810 double **Den_Rcv_Grid_A2B;
2811 double *orbs0,*orbs1;
2812 double *orbs0_0,*orbs0_1,*orbs0_2,*orbs0_3;
2813 double *orbs1_0,*orbs1_1,*orbs1_2,*orbs1_3;
2814
2815 double ***tmp_CDM;
2816 int numprocs,myid,tag=999,ID,IDS,IDR;
2817 double Stime_atom, Etime_atom;
2818
2819 MPI_Status stat;
2820 MPI_Request request;
2821
2822 /* for OpenMP */
2823 int OMPID,Nthrds,Nprocs;
2824
2825 /* MPI */
2826 MPI_Comm_size(mpi_comm_level1,&numprocs);
2827 MPI_Comm_rank(mpi_comm_level1,&myid);
2828
2829 /* allocation of arrays */
2830
2831 Tmp_Den_Grid = (double***)malloc(sizeof(double**)*(SpinP_switch+1));
2832 for (i=0; i<(SpinP_switch+1); i++){
2833 Tmp_Den_Grid[i] = (double**)malloc(sizeof(double*)*(Matomnum+1));
2834 Tmp_Den_Grid[i][0] = (double*)malloc(sizeof(double)*1);
2835 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
2836 Gc_AN = F_M2G[Mc_AN];
2837 Tmp_Den_Grid[i][Mc_AN] = (double*)malloc(sizeof(double)*GridN_Atom[Gc_AN]);
2838 }
2839 }
2840
2841 Den_Snd_Grid_A2B = (double**)malloc(sizeof(double*)*numprocs);
2842 for (ID=0; ID<numprocs; ID++){
2843 Den_Snd_Grid_A2B[ID] = (double*)malloc(sizeof(double)*Num_Snd_Grid_A2B[ID]*(SpinP_switch+1));
2844 }
2845
2846 Den_Rcv_Grid_A2B = (double**)malloc(sizeof(double*)*numprocs);
2847 for (ID=0; ID<numprocs; ID++){
2848 Den_Rcv_Grid_A2B[ID] = (double*)malloc(sizeof(double)*Num_Rcv_Grid_A2B[ID]*(SpinP_switch+1));
2849 }
2850
2851 /**********************************************
2852 calculate Tmp_Den_Grid
2853 ***********************************************/
2854
2855 #pragma omp parallel shared(FNAN,GridN_Atom,Matomnum,myid,G2ID,Orbs_Grid_FNAN,List_YOUSO,time_per_atom,Tmp_Den_Grid,Orbs_Grid,COrbs_Grid,GListTAtoms2,GListTAtoms1,NumOLG,CDM,SpinP_switch,WhatSpecies,ncn,F_G2M,natn,Spe_Total_CNO,M2G) private(OMPID,Nthrds,Nprocs,Mc_AN,h_AN,Stime_atom,Etime_atom,Gc_AN,Cwan,NO0,Gh_AN,Mh_AN,Rnh,Hwan,NO1,spin,i,j,tmp_CDM,Nog,Nc_0,Nc_1,Nc_2,Nc_3,Nh_0,Nh_1,Nh_2,Nh_3,orbs0_0,orbs0_1,orbs0_2,orbs0_3,orbs1_0,orbs1_1,orbs1_2,orbs1_3,sum_0,sum_1,sum_2,sum_3,tmp0_0,tmp0_1,tmp0_2,tmp0_3,Nc,Nh,orbs0,orbs1,sum,tmp0)
2856 {
2857
2858 orbs0 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2859 orbs1 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2860
2861 orbs0_0 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2862 orbs0_1 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2863 orbs0_2 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2864 orbs0_3 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2865 orbs1_0 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2866 orbs1_1 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2867 orbs1_2 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2868 orbs1_3 = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2869
2870 tmp_CDM = (double***)malloc(sizeof(double**)*(SpinP_switch+1));
2871 for (i=0; i<(SpinP_switch+1); i++){
2872 tmp_CDM[i] = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
2873 for (j=0; j<List_YOUSO[7]; j++){
2874 tmp_CDM[i][j] = (double*)malloc(sizeof(double)*List_YOUSO[7]);
2875 }
2876 }
2877
2878 /* get info. on OpenMP */
2879
2880 OMPID = omp_get_thread_num();
2881 Nthrds = omp_get_num_threads();
2882 Nprocs = omp_get_num_procs();
2883
2884 for (Mc_AN=(OMPID*Matomnum/Nthrds+1); Mc_AN<((OMPID+1)*Matomnum/Nthrds+1); Mc_AN++){
2885
2886 dtime(&Stime_atom);
2887
2888 /* set data on Mc_AN */
2889
2890 Gc_AN = M2G[Mc_AN];
2891 Cwan = WhatSpecies[Gc_AN];
2892 NO0 = Spe_Total_CNO[Cwan];
2893
2894 for (spin=0; spin<=SpinP_switch; spin++){
2895 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
2896 Tmp_Den_Grid[spin][Mc_AN][Nc] = 0.0;
2897 }
2898 }
2899
2900 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
2901
2902 /* set data on h_AN */
2903
2904 Gh_AN = natn[Gc_AN][h_AN];
2905 Mh_AN = F_G2M[Gh_AN];
2906 Rnh = ncn[Gc_AN][h_AN];
2907 Hwan = WhatSpecies[Gh_AN];
2908 NO1 = Spe_Total_CNO[Hwan];
2909
2910 /* store CDM into tmp_CDM */
2911
2912 for (spin=0; spin<=SpinP_switch; spin++){
2913 for (i=0; i<NO0; i++){
2914 for (j=0; j<NO1; j++){
2915 tmp_CDM[spin][i][j] = CDM[spin][Mc_AN][h_AN][i][j];
2916 }
2917 }
2918 }
2919
2920 /* summation of non-zero elements */
2921
2922 for (Nog=0; Nog<NumOLG[Mc_AN][h_AN]-3; Nog+=4){
2923
2924 Nc_0 = GListTAtoms1[Mc_AN][h_AN][Nog];
2925 Nc_1 = GListTAtoms1[Mc_AN][h_AN][Nog+1];
2926 Nc_2 = GListTAtoms1[Mc_AN][h_AN][Nog+2];
2927 Nc_3 = GListTAtoms1[Mc_AN][h_AN][Nog+3];
2928
2929 Nh_0 = GListTAtoms2[Mc_AN][h_AN][Nog];
2930 Nh_1 = GListTAtoms2[Mc_AN][h_AN][Nog+1];
2931 Nh_2 = GListTAtoms2[Mc_AN][h_AN][Nog+2];
2932 Nh_3 = GListTAtoms2[Mc_AN][h_AN][Nog+3];
2933
2934 for (i=0; i<NO0; i++){
2935 orbs0_0[i] = Orbs_Grid[Mc_AN][Nc_0][i];/* AITUNE */
2936 orbs0_1[i] = Orbs_Grid[Mc_AN][Nc_1][i];/* AITUNE */
2937 orbs0_2[i] = Orbs_Grid[Mc_AN][Nc_2][i];/* AITUNE */
2938 orbs0_3[i] = Orbs_Grid[Mc_AN][Nc_3][i];/* AITUNE */
2939 }
2940
2941 if (G2ID[Gh_AN]==myid){
2942 for (j=0; j<NO1; j++){
2943 orbs1_0[j] = Orbs_Grid[Mh_AN][Nh_0][j];/* AITUNE */
2944 orbs1_1[j] = Orbs_Grid[Mh_AN][Nh_1][j];/* AITUNE */
2945 orbs1_2[j] = Orbs_Grid[Mh_AN][Nh_2][j];/* AITUNE */
2946 orbs1_3[j] = Orbs_Grid[Mh_AN][Nh_3][j];/* AITUNE */
2947 }
2948 }
2949 else{
2950 for (j=0; j<NO1; j++){
2951 orbs1_0[j] = Orbs_Grid_FNAN[Mc_AN][h_AN][Nog ][j];/* AITUNE */
2952 orbs1_1[j] = Orbs_Grid_FNAN[Mc_AN][h_AN][Nog+1][j];/* AITUNE */
2953 orbs1_2[j] = Orbs_Grid_FNAN[Mc_AN][h_AN][Nog+2][j];/* AITUNE */
2954 orbs1_3[j] = Orbs_Grid_FNAN[Mc_AN][h_AN][Nog+3][j];/* AITUNE */
2955 }
2956 }
2957
2958 for (spin=0; spin<=SpinP_switch; spin++){
2959
2960 /* Tmp_Den_Grid */
2961
2962 sum_0 = 0.0;
2963 sum_1 = 0.0;
2964 sum_2 = 0.0;
2965 sum_3 = 0.0;
2966
2967 for (i=0; i<NO0; i++){
2968
2969 tmp0_0 = 0.0;
2970 tmp0_1 = 0.0;
2971 tmp0_2 = 0.0;
2972 tmp0_3 = 0.0;
2973
2974 for (j=0; j<NO1; j++){
2975 tmp0_0 += orbs1_0[j]*tmp_CDM[spin][i][j];
2976 tmp0_1 += orbs1_1[j]*tmp_CDM[spin][i][j];
2977 tmp0_2 += orbs1_2[j]*tmp_CDM[spin][i][j];
2978 tmp0_3 += orbs1_3[j]*tmp_CDM[spin][i][j];
2979 }
2980
2981 sum_0 += orbs0_0[i]*tmp0_0;
2982 sum_1 += orbs0_1[i]*tmp0_1;
2983 sum_2 += orbs0_2[i]*tmp0_2;
2984 sum_3 += orbs0_3[i]*tmp0_3;
2985 }
2986
2987 Tmp_Den_Grid[spin][Mc_AN][Nc_0] += sum_0;
2988 Tmp_Den_Grid[spin][Mc_AN][Nc_1] += sum_1;
2989 Tmp_Den_Grid[spin][Mc_AN][Nc_2] += sum_2;
2990 Tmp_Den_Grid[spin][Mc_AN][Nc_3] += sum_3;
2991
2992 } /* spin */
2993 } /* Nog */
2994
2995 for (; Nog<NumOLG[Mc_AN][h_AN]; Nog++){
2996
2997 Nc = GListTAtoms1[Mc_AN][h_AN][Nog];
2998 Nh = GListTAtoms2[Mc_AN][h_AN][Nog];
2999
3000 for (i=0; i<NO0; i++){
3001 orbs0[i] = Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
3002 }
3003
3004 if (G2ID[Gh_AN]==myid){
3005 for (j=0; j<NO1; j++){
3006 orbs1[j] = Orbs_Grid[Mh_AN][Nh][j];/* AITUNE */
3007 }
3008 }
3009 else{
3010 for (j=0; j<NO1; j++){
3011 orbs1[j] = Orbs_Grid_FNAN[Mc_AN][h_AN][Nog][j];/* AITUNE */
3012 }
3013 }
3014
3015 for (spin=0; spin<=SpinP_switch; spin++){
3016
3017 /* Tmp_Den_Grid */
3018
3019 sum = 0.0;
3020 for (i=0; i<NO0; i++){
3021 tmp0 = 0.0;
3022 for (j=0; j<NO1; j++){
3023 tmp0 += orbs1[j]*tmp_CDM[spin][i][j];
3024 }
3025 sum += orbs0[i]*tmp0;
3026 }
3027
3028 Tmp_Den_Grid[spin][Mc_AN][Nc] += sum;
3029 }
3030 }
3031
3032 } /* h_AN */
3033
3034 dtime(&Etime_atom);
3035 time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
3036
3037 } /* Mc_AN */
3038
3039 /* freeing of arrays */
3040
3041 free(orbs0);
3042 free(orbs1);
3043
3044 free(orbs0_0);
3045 free(orbs0_1);
3046 free(orbs0_2);
3047 free(orbs0_3);
3048 free(orbs1_0);
3049 free(orbs1_1);
3050 free(orbs1_2);
3051 free(orbs1_3);
3052
3053 for (i=0; i<(SpinP_switch+1); i++){
3054 for (j=0; j<List_YOUSO[7]; j++){
3055 free(tmp_CDM[i][j]);
3056 }
3057 free(tmp_CDM[i]);
3058 }
3059 free(tmp_CDM);
3060
3061 #pragma omp flush(Tmp_Den_Grid)
3062
3063 } /* #pragma omp parallel */
3064
3065 /******************************************************
3066 MPI communication from the partitions A to B
3067 ******************************************************/
3068
3069 /* copy Tmp_Den_Grid to Den_Snd_Grid_A2B */
3070
3071 for (ID=0; ID<numprocs; ID++) Num_Snd_Grid_A2B[ID] = 0;
3072
3073 N2D = Ngrid1*Ngrid2;
3074
3075 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
3076
3077 Gc_AN = M2G[Mc_AN];
3078
3079 for (AN=0; AN<GridN_Atom[Gc_AN]; AN++){
3080
3081 GN = GridListAtom[Mc_AN][AN];
3082 GN2N(GN,N3);
3083 n2D = N3[1]*Ngrid2 + N3[2];
3084 ID = (int)(n2D*(unsigned long long int)numprocs/N2D);
3085
3086 if (SpinP_switch==0){
3087 Den_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]] = Tmp_Den_Grid[0][Mc_AN][AN];
3088 }
3089 else if (SpinP_switch==1){
3090 Den_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]*2+0] = Tmp_Den_Grid[0][Mc_AN][AN];
3091 Den_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]*2+1] = Tmp_Den_Grid[1][Mc_AN][AN];
3092 }
3093 else if (SpinP_switch==3){
3094 Den_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]*4+0] = Tmp_Den_Grid[0][Mc_AN][AN];
3095 Den_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]*4+1] = Tmp_Den_Grid[1][Mc_AN][AN];
3096 Den_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]*4+2] = Tmp_Den_Grid[2][Mc_AN][AN];
3097 Den_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]*4+3] = Tmp_Den_Grid[3][Mc_AN][AN];
3098 }
3099
3100 Num_Snd_Grid_A2B[ID]++;
3101 }
3102 }
3103
3104 /* MPI: A to B */
3105
3106 tag = 999;
3107 for (ID=0; ID<numprocs; ID++){
3108
3109 IDS = (myid + ID) % numprocs;
3110 IDR = (myid - ID + numprocs) % numprocs;
3111
3112 if (Num_Snd_Grid_A2B[IDS]!=0){
3113 MPI_Isend( &Den_Snd_Grid_A2B[IDS][0], Num_Snd_Grid_A2B[IDS]*(SpinP_switch+1),
3114 MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request);
3115 }
3116
3117 if (Num_Rcv_Grid_A2B[IDR]!=0){
3118 MPI_Recv( &Den_Rcv_Grid_A2B[IDR][0], Num_Rcv_Grid_A2B[IDR]*(SpinP_switch+1),
3119 MPI_DOUBLE, IDR, tag, mpi_comm_level1, &stat);
3120 }
3121
3122 if (Num_Snd_Grid_A2B[IDS]!=0) MPI_Wait(&request,&stat);
3123 }
3124
3125 /******************************************************
3126 superposition of rho_i to calculate charge density
3127 in the partition B.
3128 ******************************************************/
3129
3130 /* initialize arrays */
3131
3132 for (spin=0; spin<(SpinP_switch+1); spin++){
3133 for (BN=0; BN<My_NumGridB_AB; BN++){
3134 Density_Grid_B[spin][BN] = 0.0;
3135 }
3136 }
3137
3138 /* superposition of densities rho_i */
3139
3140 for (ID=0; ID<numprocs; ID++){
3141
3142 for (LN=0; LN<Num_Rcv_Grid_A2B[ID]; LN++){
3143
3144 BN = Index_Rcv_Grid_A2B[ID][3*LN+0];
3145 Gc_AN = Index_Rcv_Grid_A2B[ID][3*LN+1];
3146 GRc = Index_Rcv_Grid_A2B[ID][3*LN+2];
3147
3148 if (Solver!=4 || (Solver==4 && atv_ijk[GRc][1]==0 )){
3149
3150 /* spin collinear non-polarization */
3151 if ( SpinP_switch==0 ){
3152 Density_Grid_B[0][BN] += Den_Rcv_Grid_A2B[ID][LN];
3153 }
3154
3155 /* spin collinear polarization */
3156 else if ( SpinP_switch==1 ){
3157 Density_Grid_B[0][BN] += Den_Rcv_Grid_A2B[ID][LN*2 ];
3158 Density_Grid_B[1][BN] += Den_Rcv_Grid_A2B[ID][LN*2+1];
3159 }
3160
3161 /* spin non-collinear */
3162 else if ( SpinP_switch==3 ){
3163 Density_Grid_B[0][BN] += Den_Rcv_Grid_A2B[ID][LN*4 ];
3164 Density_Grid_B[1][BN] += Den_Rcv_Grid_A2B[ID][LN*4+1];
3165 Density_Grid_B[2][BN] += Den_Rcv_Grid_A2B[ID][LN*4+2];
3166 Density_Grid_B[3][BN] += Den_Rcv_Grid_A2B[ID][LN*4+3];
3167 }
3168
3169 } /* if (Solve!=4.....) */
3170
3171 } /* AN */
3172 } /* ID */
3173
3174 /******************************************************
3175 MPI: from the partitions B to C
3176 ******************************************************/
3177
3178 Data_Grid_Copy_B2C_2( Density_Grid_B, Density_Grid );
3179
3180 /******************************************************
3181 diagonalize spinor density
3182 ******************************************************/
3183
3184 if (SpinP_switch==3) diagonalize_nc_density();
3185
3186 /* freeing of arrays */
3187
3188 for (i=0; i<(SpinP_switch+1); i++){
3189 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
3190 free(Tmp_Den_Grid[i][Mc_AN]);
3191 }
3192 free(Tmp_Den_Grid[i]);
3193 }
3194 free(Tmp_Den_Grid);
3195
3196 for (ID=0; ID<numprocs; ID++){
3197 free(Den_Snd_Grid_A2B[ID]);
3198 }
3199 free(Den_Snd_Grid_A2B);
3200
3201 for (ID=0; ID<numprocs; ID++){
3202 free(Den_Rcv_Grid_A2B[ID]);
3203 }
3204 free(Den_Rcv_Grid_A2B);
3205 }
3206