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