1 /**********************************************************************
2 SCF2File.c:
3
4 SCF2File.c is a subroutine to output connectivity, Hamiltonian,
5 overlap, and etc. to a binary file, filename.scfout.
6
7 Log of SCF2DFT.c:
8
9 1/July/2003 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <math.h>
16 #include <string.h>
17 #include <time.h>
18 #include "openmx_common.h"
19 #include "mpi.h"
20
21 #define MAX_LINE_SIZE 256
22
23 static void Output( FILE *fp, char *inputfile );
24 static void Calc_OLPpo();
25
26 static double ****OLPpox;
27 static double ****OLPpoy;
28 static double ****OLPpoz;
29
30
31
SCF2File(char * mode,char * inputfile)32 void SCF2File(char *mode, char *inputfile)
33 {
34 int Mc_AN,Gc_AN,tno0,tno1;
35 int Cwan,h_AN,Gh_AN,Hwan;
36 int i;
37 char fname[YOUSO10];
38 char operate[300];
39 FILE *fp;
40 int numprocs,myid,ID;
41 int succeed_open;
42 char buf[fp_bsize]; /* setvbuf */
43
44 /* MPI */
45 MPI_Comm_size(mpi_comm_level1,&numprocs);
46 MPI_Comm_rank(mpi_comm_level1,&myid);
47
48 /* allocation of array */
49
50 OLPpox = (double****)malloc(sizeof(double***)*(Matomnum+1));
51 FNAN[0] = 0;
52 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
53
54 if (Mc_AN==0){
55 Gc_AN = 0;
56 tno0 = 1;
57 }
58 else{
59 Gc_AN = S_M2G[Mc_AN];
60 Cwan = WhatSpecies[Gc_AN];
61 tno0 = Spe_Total_CNO[Cwan];
62 }
63
64 OLPpox[Mc_AN] = (double***)malloc(sizeof(double**)*(FNAN[Gc_AN]+1));
65 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
66
67 if (Mc_AN==0){
68 tno1 = 1;
69 }
70 else{
71 Gh_AN = natn[Gc_AN][h_AN];
72 Hwan = WhatSpecies[Gh_AN];
73 tno1 = Spe_Total_CNO[Hwan];
74 }
75
76 OLPpox[Mc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
77 for (i=0; i<tno0; i++){
78 OLPpox[Mc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
79 }
80 }
81 }
82
83 OLPpoy = (double****)malloc(sizeof(double***)*(Matomnum+1));
84 FNAN[0] = 0;
85 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
86
87 if (Mc_AN==0){
88 Gc_AN = 0;
89 tno0 = 1;
90 }
91 else{
92 Gc_AN = S_M2G[Mc_AN];
93 Cwan = WhatSpecies[Gc_AN];
94 tno0 = Spe_Total_CNO[Cwan];
95 }
96
97 OLPpoy[Mc_AN] = (double***)malloc(sizeof(double**)*(FNAN[Gc_AN]+1));
98 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
99
100 if (Mc_AN==0){
101 tno1 = 1;
102 }
103 else{
104 Gh_AN = natn[Gc_AN][h_AN];
105 Hwan = WhatSpecies[Gh_AN];
106 tno1 = Spe_Total_CNO[Hwan];
107 }
108
109 OLPpoy[Mc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
110 for (i=0; i<tno0; i++){
111 OLPpoy[Mc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
112 }
113 }
114 }
115
116 OLPpoz = (double****)malloc(sizeof(double***)*(Matomnum+1));
117 FNAN[0] = 0;
118 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
119
120 if (Mc_AN==0){
121 Gc_AN = 0;
122 tno0 = 1;
123 }
124 else{
125 Gc_AN = S_M2G[Mc_AN];
126 Cwan = WhatSpecies[Gc_AN];
127 tno0 = Spe_Total_CNO[Cwan];
128 }
129
130 OLPpoz[Mc_AN] = (double***)malloc(sizeof(double**)*(FNAN[Gc_AN]+1));
131 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
132
133 if (Mc_AN==0){
134 tno1 = 1;
135 }
136 else{
137 Gh_AN = natn[Gc_AN][h_AN];
138 Hwan = WhatSpecies[Gh_AN];
139 tno1 = Spe_Total_CNO[Hwan];
140 }
141
142 OLPpoz[Mc_AN][h_AN] = (double**)malloc(sizeof(double*)*tno0);
143 for (i=0; i<tno0; i++){
144 OLPpoz[Mc_AN][h_AN][i] = (double*)malloc(sizeof(double)*tno1);
145 }
146 }
147 }
148
149 /* write a file */
150
151 if ( strcasecmp(mode,"write")==0 ) {
152 sprintf(fname,"%s%s.scfout",filepath,filename);
153
154 if (myid==Host_ID){
155
156 if ((fp = fopen(fname,"w")) != NULL){
157
158 #ifdef xt3
159 setvbuf(fp,buf,_IOFBF,fp_bsize); /* setvbuf */
160 #endif
161
162 succeed_open = 1;
163 }
164 else {
165 succeed_open = 0;
166 printf("Failure of making the scfout file (%s).\n",fname);
167 }
168 }
169
170 MPI_Bcast(&succeed_open, 1, MPI_INT, Host_ID, mpi_comm_level1);
171
172 if (succeed_open==1){
173 if (myid==Host_ID) printf("Save the scfout file (%s)\n",fname);
174
175 Calc_OLPpo();
176 Output( fp, inputfile );
177 if (myid==Host_ID) fclose(fp);
178 }
179 }
180
181 /* free array */
182
183 FNAN[0] = 0;
184 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
185
186 if (Mc_AN==0){
187 Gc_AN = 0;
188 tno0 = 1;
189 }
190 else{
191 Gc_AN = S_M2G[Mc_AN];
192 Cwan = WhatSpecies[Gc_AN];
193 tno0 = Spe_Total_CNO[Cwan];
194 }
195
196 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
197
198 if (Mc_AN==0){
199 tno1 = 1;
200 }
201 else{
202 Gh_AN = natn[Gc_AN][h_AN];
203 Hwan = WhatSpecies[Gh_AN];
204 tno1 = Spe_Total_CNO[Hwan];
205 }
206
207 for (i=0; i<tno0; i++){
208 free(OLPpox[Mc_AN][h_AN][i]);
209 }
210 free(OLPpox[Mc_AN][h_AN]);
211 }
212 free(OLPpox[Mc_AN]);
213 }
214 free(OLPpox);
215
216 FNAN[0] = 0;
217 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
218
219 if (Mc_AN==0){
220 Gc_AN = 0;
221 tno0 = 1;
222 }
223 else{
224 Gc_AN = S_M2G[Mc_AN];
225 Cwan = WhatSpecies[Gc_AN];
226 tno0 = Spe_Total_CNO[Cwan];
227 }
228
229 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
230
231 if (Mc_AN==0){
232 tno1 = 1;
233 }
234 else{
235 Gh_AN = natn[Gc_AN][h_AN];
236 Hwan = WhatSpecies[Gh_AN];
237 tno1 = Spe_Total_CNO[Hwan];
238 }
239
240 for (i=0; i<tno0; i++){
241 free(OLPpoy[Mc_AN][h_AN][i]);
242 }
243 free(OLPpoy[Mc_AN][h_AN]);
244 }
245 free(OLPpoy[Mc_AN]);
246 }
247 free(OLPpoy);
248
249 FNAN[0] = 0;
250 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
251
252 if (Mc_AN==0){
253 Gc_AN = 0;
254 tno0 = 1;
255 }
256 else{
257 Gc_AN = S_M2G[Mc_AN];
258 Cwan = WhatSpecies[Gc_AN];
259 tno0 = Spe_Total_CNO[Cwan];
260 }
261
262 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
263
264 if (Mc_AN==0){
265 tno1 = 1;
266 }
267 else{
268 Gh_AN = natn[Gc_AN][h_AN];
269 Hwan = WhatSpecies[Gh_AN];
270 tno1 = Spe_Total_CNO[Hwan];
271 }
272
273 for (i=0; i<tno0; i++){
274 free(OLPpoz[Mc_AN][h_AN][i]);
275 }
276 free(OLPpoz[Mc_AN][h_AN]);
277 }
278 free(OLPpoz[Mc_AN]);
279 }
280 free(OLPpoz);
281 }
282
283
284
285
Output(FILE * fp,char * inputfile)286 void Output( FILE *fp, char *inputfile )
287 {
288 int Gc_AN,Mc_AN,ct_AN,h_AN,i,j,can,Gh_AN;
289 int num,wan1,wan2,TNO1,TNO2,spin,Rn,count;
290 int k,Mh_AN,Rnh,Hwan,NO0,NO1,Qwan;
291 int q_AN,Mq_AN,Gq_AN,Rnq;
292 int i_vec[20],*p_vec;
293 int numprocs,myid,ID,tag=999;
294 double *Tmp_Vec,d_vec[20];
295 FILE *fp_inp;
296 char strg[MAX_LINE_SIZE];
297 char buf[fp_bsize]; /* setvbuf */
298
299 MPI_Status stat;
300 MPI_Request request;
301
302 /* MPI */
303 MPI_Comm_size(mpi_comm_level1,&numprocs);
304 MPI_Comm_rank(mpi_comm_level1,&myid);
305
306 /***************************************************************
307 information of connectivity
308 ****************************************************************/
309
310 if (myid==Host_ID){
311
312 /****************************************************
313 atomnum
314 spinP_switch
315 ****************************************************/
316
317 i = 0;
318 i_vec[0] = atomnum;
319 i_vec[1] = SpinP_switch;
320 i_vec[2] = Catomnum;
321 i_vec[3] = Latomnum;
322 i_vec[4] = Ratomnum;
323 i_vec[5] = TCpyCell;
324
325 fwrite(i_vec,sizeof(int),6,fp);
326
327 /****************************************************
328 atv[Rn][4]
329 ****************************************************/
330
331 for (Rn=0; Rn<=TCpyCell; Rn++){
332 fwrite(atv[Rn],sizeof(double),4,fp);
333 }
334
335 /****************************************************
336 atv_ijk[Rn][4]
337 ****************************************************/
338
339 for (Rn=0; Rn<=TCpyCell; Rn++){
340 fwrite(atv_ijk[Rn],sizeof(int),4,fp);
341 }
342
343 /****************************************************
344 # of orbitals in each atom
345 ****************************************************/
346
347 p_vec = (int*)malloc(sizeof(int)*atomnum);
348 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
349 wan1 = WhatSpecies[ct_AN];
350 TNO1 = Spe_Total_CNO[wan1];
351 p_vec[ct_AN-1] = TNO1;
352 }
353 fwrite(p_vec,sizeof(int),atomnum,fp);
354 free(p_vec);
355
356 /****************************************************
357 FNAN[]:
358 number of first nearest neighbouring atoms
359 ****************************************************/
360
361 p_vec = (int*)malloc(sizeof(int)*atomnum);
362 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
363 p_vec[ct_AN-1] = FNAN[ct_AN];
364 }
365 fwrite(p_vec,sizeof(int),atomnum,fp);
366 free(p_vec);
367
368 /****************************************************
369 natn[][]:
370 grobal index of neighboring atoms of an atom ct_AN
371 ****************************************************/
372
373 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
374 fwrite(natn[ct_AN],sizeof(int),FNAN[ct_AN]+1,fp);
375 }
376
377 /****************************************************
378 ncn[][]:
379 grobal index for cell of neighboring atoms
380 of an atom ct_AN
381 ****************************************************/
382
383 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
384 fwrite(ncn[ct_AN],sizeof(int),FNAN[ct_AN]+1,fp);
385 }
386
387 /****************************************************
388 tv[4][4]:
389 unit cell vectors in Bohr
390 ****************************************************/
391
392 fwrite(tv[1],sizeof(double),4,fp);
393 fwrite(tv[2],sizeof(double),4,fp);
394 fwrite(tv[3],sizeof(double),4,fp);
395
396 /****************************************************
397 rtv[4][4]:
398 reciprocal unit cell vectors in Bohr^{-1}
399 ****************************************************/
400
401 fwrite(rtv[1],sizeof(double),4,fp);
402 fwrite(rtv[2],sizeof(double),4,fp);
403 fwrite(rtv[3],sizeof(double),4,fp);
404
405 /****************************************************
406 Gxyz[][1-3]:
407 atomic coordinates in Bohr
408 ****************************************************/
409
410 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
411 fwrite(Gxyz[ct_AN],sizeof(double),4,fp);
412 }
413
414 } /* if (myid==Host_ID) */
415
416 /***************************************************************
417 Hamiltonian matrix
418 ****************************************************************/
419
420 Tmp_Vec = (double*)malloc(sizeof(double)*List_YOUSO[8]*List_YOUSO[7]*List_YOUSO[7]);
421
422 for (spin=0; spin<=SpinP_switch; spin++){
423
424 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
425 ID = G2ID[Gc_AN];
426
427 if (myid==ID){
428
429 num = 0;
430
431 Mc_AN = F_G2M[Gc_AN];
432 wan1 = WhatSpecies[Gc_AN];
433 TNO1 = Spe_Total_CNO[wan1];
434 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
435 Gh_AN = natn[Gc_AN][h_AN];
436 wan2 = WhatSpecies[Gh_AN];
437 TNO2 = Spe_Total_CNO[wan2];
438
439 if (Cnt_switch==0){
440 for (i=0; i<TNO1; i++){
441 for (j=0; j<TNO2; j++){
442 Tmp_Vec[num] = H[spin][Mc_AN][h_AN][i][j];
443 num++;
444 }
445 }
446 }
447 else{
448 for (i=0; i<TNO1; i++){
449 for (j=0; j<TNO2; j++){
450 Tmp_Vec[num] = CntH[spin][Mc_AN][h_AN][i][j];
451 num++;
452 }
453 }
454 }
455 }
456
457 if (myid!=Host_ID){
458 MPI_Isend(&num, 1, MPI_INT, Host_ID, tag, mpi_comm_level1, &request);
459 MPI_Wait(&request,&stat);
460 MPI_Isend(&Tmp_Vec[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1, &request);
461 MPI_Wait(&request,&stat);
462 }
463 else{
464 fwrite(Tmp_Vec, sizeof(double), num, fp);
465 }
466 }
467
468 else if (ID!=myid && myid==Host_ID){
469 MPI_Recv(&num, 1, MPI_INT, ID, tag, mpi_comm_level1, &stat);
470 MPI_Recv(&Tmp_Vec[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
471 fwrite(Tmp_Vec, sizeof(double), num, fp);
472 }
473
474 }
475 }
476
477 /***************************************************************
478 iHNL
479 ****************************************************************/
480
481 if ( SpinP_switch==3 ){
482
483 for (spin=0; spin<List_YOUSO[5]; spin++){
484
485 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
486 ID = G2ID[Gc_AN];
487
488 if (myid==ID){
489
490 num = 0;
491
492 Mc_AN = F_G2M[Gc_AN];
493 wan1 = WhatSpecies[Gc_AN];
494 TNO1 = Spe_Total_CNO[wan1];
495 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
496 Gh_AN = natn[Gc_AN][h_AN];
497 wan2 = WhatSpecies[Gh_AN];
498 TNO2 = Spe_Total_CNO[wan2];
499
500 for (i=0; i<TNO1; i++){
501 for (j=0; j<TNO2; j++){
502 Tmp_Vec[num] = iHNL[spin][Mc_AN][h_AN][i][j];
503 num++;
504 }
505 }
506 }
507
508 if (myid!=Host_ID){
509 MPI_Isend(&num, 1, MPI_INT, Host_ID, tag, mpi_comm_level1, &request);
510 MPI_Wait(&request,&stat);
511 MPI_Isend(&Tmp_Vec[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1, &request);
512 MPI_Wait(&request,&stat);
513 }
514 else{
515 fwrite(Tmp_Vec, sizeof(double), num, fp);
516 }
517 }
518
519 else if (ID!=myid && myid==Host_ID){
520 MPI_Recv(&num, 1, MPI_INT, ID, tag, mpi_comm_level1, &stat);
521 MPI_Recv(&Tmp_Vec[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
522 fwrite(Tmp_Vec, sizeof(double), num, fp);
523 }
524
525 }
526 }
527 }
528
529 /***************************************************************
530 overlap matrix
531 ****************************************************************/
532
533 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
534 ID = G2ID[Gc_AN];
535
536 if (myid==ID){
537
538 num = 0;
539
540 Mc_AN = F_G2M[Gc_AN];
541 wan1 = WhatSpecies[Gc_AN];
542 TNO1 = Spe_Total_CNO[wan1];
543 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
544 Gh_AN = natn[Gc_AN][h_AN];
545 wan2 = WhatSpecies[Gh_AN];
546 TNO2 = Spe_Total_CNO[wan2];
547
548 if (Cnt_switch==0){
549 for (i=0; i<TNO1; i++){
550 for (j=0; j<TNO2; j++){
551 Tmp_Vec[num] = OLP[0][Mc_AN][h_AN][i][j];
552 num++;
553 }
554 }
555 }
556 else{
557 for (i=0; i<TNO1; i++){
558 for (j=0; j<TNO2; j++){
559 Tmp_Vec[num] = CntOLP[0][Mc_AN][h_AN][i][j];
560 num++;
561 }
562 }
563 }
564 }
565
566 if (myid!=Host_ID){
567 MPI_Isend(&num, 1, MPI_INT, Host_ID, tag, mpi_comm_level1, &request);
568 MPI_Wait(&request,&stat);
569 MPI_Isend(&Tmp_Vec[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1, &request);
570 MPI_Wait(&request,&stat);
571 }
572 else{
573 fwrite(Tmp_Vec, sizeof(double), num, fp);
574 }
575 }
576
577 else if (ID!=myid && myid==Host_ID){
578 MPI_Recv(&num, 1, MPI_INT, ID, tag, mpi_comm_level1, &stat);
579 MPI_Recv(&Tmp_Vec[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
580 fwrite(Tmp_Vec, sizeof(double), num, fp);
581 }
582
583 }
584
585 /***************************************************************
586 overlap matrix with position operator x
587 ****************************************************************/
588
589 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
590 ID = G2ID[Gc_AN];
591
592 if (myid==ID){
593
594 num = 0;
595 Mc_AN = F_G2M[Gc_AN];
596 wan1 = WhatSpecies[Gc_AN];
597 TNO1 = Spe_Total_CNO[wan1];
598
599 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
600 Gh_AN = natn[Gc_AN][h_AN];
601 wan2 = WhatSpecies[Gh_AN];
602 TNO2 = Spe_Total_CNO[wan2];
603
604 for (i=0; i<TNO1; i++){
605 for (j=0; j<TNO2; j++){
606 Tmp_Vec[num] = OLPpox[Mc_AN][h_AN][i][j];
607 num++;
608 }
609 }
610 }
611
612 if (myid!=Host_ID){
613 MPI_Isend(&num, 1, MPI_INT, Host_ID, tag, mpi_comm_level1, &request);
614 MPI_Wait(&request,&stat);
615 MPI_Isend(&Tmp_Vec[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1, &request);
616 MPI_Wait(&request,&stat);
617 }
618 else{
619 fwrite(Tmp_Vec, sizeof(double), num, fp);
620 }
621 }
622
623 else if (ID!=myid && myid==Host_ID){
624 MPI_Recv(&num, 1, MPI_INT, ID, tag, mpi_comm_level1, &stat);
625 MPI_Recv(&Tmp_Vec[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
626 fwrite(Tmp_Vec, sizeof(double), num, fp);
627 }
628 }
629
630 /***************************************************************
631 overlap matrix with position operator y
632 ****************************************************************/
633
634 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
635 ID = G2ID[Gc_AN];
636
637 if (myid==ID){
638
639 num = 0;
640 Mc_AN = F_G2M[Gc_AN];
641 wan1 = WhatSpecies[Gc_AN];
642 TNO1 = Spe_Total_CNO[wan1];
643
644 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
645 Gh_AN = natn[Gc_AN][h_AN];
646 wan2 = WhatSpecies[Gh_AN];
647 TNO2 = Spe_Total_CNO[wan2];
648
649 for (i=0; i<TNO1; i++){
650 for (j=0; j<TNO2; j++){
651 Tmp_Vec[num] = OLPpoy[Mc_AN][h_AN][i][j];
652 num++;
653 }
654 }
655 }
656
657 if (myid!=Host_ID){
658 MPI_Isend(&num, 1, MPI_INT, Host_ID, tag, mpi_comm_level1, &request);
659 MPI_Wait(&request,&stat);
660 MPI_Isend(&Tmp_Vec[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1, &request);
661 MPI_Wait(&request,&stat);
662 }
663 else{
664 fwrite(Tmp_Vec, sizeof(double), num, fp);
665 }
666 }
667
668 else if (ID!=myid && myid==Host_ID){
669 MPI_Recv(&num, 1, MPI_INT, ID, tag, mpi_comm_level1, &stat);
670 MPI_Recv(&Tmp_Vec[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
671 fwrite(Tmp_Vec, sizeof(double), num, fp);
672 }
673 }
674
675 /***************************************************************
676 overlap matrix with position operator z
677 ****************************************************************/
678
679 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
680 ID = G2ID[Gc_AN];
681
682 if (myid==ID){
683
684 num = 0;
685 Mc_AN = F_G2M[Gc_AN];
686 wan1 = WhatSpecies[Gc_AN];
687 TNO1 = Spe_Total_CNO[wan1];
688
689 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
690 Gh_AN = natn[Gc_AN][h_AN];
691 wan2 = WhatSpecies[Gh_AN];
692 TNO2 = Spe_Total_CNO[wan2];
693
694 for (i=0; i<TNO1; i++){
695 for (j=0; j<TNO2; j++){
696 Tmp_Vec[num] = OLPpoz[Mc_AN][h_AN][i][j];
697 num++;
698 }
699 }
700 }
701
702 if (myid!=Host_ID){
703 MPI_Isend(&num, 1, MPI_INT, Host_ID, tag, mpi_comm_level1, &request);
704 MPI_Wait(&request,&stat);
705 MPI_Isend(&Tmp_Vec[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1, &request);
706 MPI_Wait(&request,&stat);
707 }
708 else{
709 fwrite(Tmp_Vec, sizeof(double), num, fp);
710 }
711 }
712
713 else if (ID!=myid && myid==Host_ID){
714 MPI_Recv(&num, 1, MPI_INT, ID, tag, mpi_comm_level1, &stat);
715 MPI_Recv(&Tmp_Vec[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
716 fwrite(Tmp_Vec, sizeof(double), num, fp);
717 }
718 }
719
720 /***************************************************************
721 density matrix
722 ****************************************************************/
723
724 for (spin=0; spin<=SpinP_switch; spin++){
725
726 for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
727 ID = G2ID[Gc_AN];
728
729 if (myid==ID){
730
731 num = 0;
732
733 Mc_AN = F_G2M[Gc_AN];
734 wan1 = WhatSpecies[Gc_AN];
735 TNO1 = Spe_Total_CNO[wan1];
736 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
737 Gh_AN = natn[Gc_AN][h_AN];
738 wan2 = WhatSpecies[Gh_AN];
739 TNO2 = Spe_Total_CNO[wan2];
740
741 for (i=0; i<TNO1; i++){
742 for (j=0; j<TNO2; j++){
743 Tmp_Vec[num] = DM[0][spin][Mc_AN][h_AN][i][j];
744 num++;
745 }
746 }
747 }
748
749 if (myid!=Host_ID){
750 MPI_Isend(&num, 1, MPI_INT, Host_ID, tag, mpi_comm_level1, &request);
751 MPI_Wait(&request,&stat);
752 MPI_Isend(&Tmp_Vec[0], num, MPI_DOUBLE, Host_ID, tag, mpi_comm_level1, &request);
753 MPI_Wait(&request,&stat);
754 }
755 else{
756 fwrite(Tmp_Vec, sizeof(double), num, fp);
757 }
758 }
759
760 else if (ID!=myid && myid==Host_ID){
761 MPI_Recv(&num, 1, MPI_INT, ID, tag, mpi_comm_level1, &stat);
762 MPI_Recv(&Tmp_Vec[0], num, MPI_DOUBLE, ID, tag, mpi_comm_level1, &stat);
763 fwrite(Tmp_Vec, sizeof(double), num, fp);
764 }
765
766 }
767 }
768
769 /* freeing of Tmp_Vec */
770 free(Tmp_Vec);
771
772 if (myid==Host_ID){
773
774 /****************************************************
775 Solver
776 ****************************************************/
777
778 if (PeriodicGamma_flag==1) i_vec[0] = 3;
779 else i_vec[0] = Solver;
780 fwrite(i_vec,sizeof(int),1,fp);
781
782 /****************************************************
783 ChemP
784 Electronic Temp
785 dipole moment (x, y, z) from core charge
786 dipole moment (x, y, z) from back ground charge
787 ****************************************************/
788
789 d_vec[0] = ChemP;
790 d_vec[1] = E_Temp*eV2Hartree;
791 d_vec[2] = dipole_moment[1][1];
792 d_vec[3] = dipole_moment[1][2];
793 d_vec[4] = dipole_moment[1][3];
794 d_vec[5] = dipole_moment[3][1];
795 d_vec[6] = dipole_moment[3][2];
796 d_vec[7] = dipole_moment[3][3];
797 d_vec[8] = Total_Num_Electrons;
798 d_vec[9] = Total_SpinS;
799 fwrite(d_vec,sizeof(double),10,fp);
800 }
801
802 /****************************************************
803 input file
804 ****************************************************/
805
806 if (myid==Host_ID){
807
808 /* find the number of lines in the input file */
809
810 if ((fp_inp = fopen(inputfile,"r")) != NULL){
811
812 #ifdef xt3
813 setvbuf(fp_inp,buf,_IOFBF,fp_bsize); /* setvbuf */
814 #endif
815
816 count = 0;
817 /* fgets gets a carriage return */
818 while( fgets(strg, MAX_LINE_SIZE, fp_inp ) != NULL ){
819 count++;
820 }
821
822 fclose(fp_inp);
823 }
824 else{
825 printf("error #1 in saving *.scfout\n");
826 }
827
828 /* write the input file */
829
830 if ((fp_inp = fopen(inputfile,"r")) != NULL){
831
832 #ifdef xt3
833 setvbuf(fp_inp,buf,_IOFBF,fp_bsize); /* setvbuf */
834 #endif
835
836 i_vec[0] = count;
837 fwrite(i_vec, sizeof(int), 1, fp);
838
839 /* fgets gets a carriage return */
840 while( fgets(strg, MAX_LINE_SIZE, fp_inp ) != NULL ){
841 fwrite(strg, sizeof(char), MAX_LINE_SIZE, fp);
842 }
843
844 fclose(fp_inp);
845 }
846 else{
847 printf("error #1 in saving *.scfout\n");
848 }
849 }
850
851 }
852
853
854
855
856
Calc_OLPpo()857 void Calc_OLPpo()
858 {
859 double time0;
860 int Mc_AN,Gc_AN,Mh_AN,h_AN,Gh_AN;
861 int i,j,Cwan,Hwan,NO0,NO1,spinmax;
862 int Rnh,Rnk,spin,N,NumC[4];
863 int n1,n2,n3,L0,Mul0,M0,L1,Mul1,M1;
864 int Nc,GNc,GRc,Nog,Nh,MN,XC_P_switch;
865 double x,y,z,dx,dy,dz,tmpx,tmpy,tmpz;
866 double bc,dv,r,theta,phi,sum,tmp0,tmp1;
867 double xo,yo,zo,S_coordinate[3];
868 double Cxyz[4];
869 double **ChiVx;
870 double **ChiVy;
871 double **ChiVz;
872 double *tmp_ChiVx;
873 double *tmp_ChiVy;
874 double *tmp_ChiVz;
875 double *tmp_Orbs_Grid;
876 double **tmp_OLPpox;
877 double **tmp_OLPpoy;
878 double **tmp_OLPpoz;
879 double TStime,TEtime;
880 double TStime0,TEtime0;
881 double TStime1,TEtime1;
882 double TStime2,TEtime2;
883 double TStime3,TEtime3;
884 int numprocs,myid,tag=999,ID;
885 double Stime_atom, Etime_atom;
886
887 MPI_Comm_size(mpi_comm_level1,&numprocs);
888 MPI_Comm_rank(mpi_comm_level1,&myid);
889 MPI_Barrier(mpi_comm_level1);
890
891 /****************************************************
892 allocation of arrays:
893 ****************************************************/
894
895 ChiVx = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
896 for (i=0; i<List_YOUSO[7]; i++){
897 ChiVx[i] = (double*)malloc(sizeof(double)*List_YOUSO[11]);
898 }
899
900 ChiVy = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
901 for (i=0; i<List_YOUSO[7]; i++){
902 ChiVy[i] = (double*)malloc(sizeof(double)*List_YOUSO[11]);
903 }
904
905 ChiVz = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
906 for (i=0; i<List_YOUSO[7]; i++){
907 ChiVz[i] = (double*)malloc(sizeof(double)*List_YOUSO[11]);
908 }
909
910 tmp_ChiVx = (double*)malloc(sizeof(double)*List_YOUSO[7]);
911 tmp_ChiVy = (double*)malloc(sizeof(double)*List_YOUSO[7]);
912 tmp_ChiVz = (double*)malloc(sizeof(double)*List_YOUSO[7]);
913
914 tmp_Orbs_Grid = (double*)malloc(sizeof(double)*List_YOUSO[7]);
915
916 tmp_OLPpox = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
917 for (i=0; i<List_YOUSO[7]; i++){
918 tmp_OLPpox[i] = (double*)malloc(sizeof(double)*List_YOUSO[7]);
919 }
920
921 tmp_OLPpoy = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
922 for (i=0; i<List_YOUSO[7]; i++){
923 tmp_OLPpoy[i] = (double*)malloc(sizeof(double)*List_YOUSO[7]);
924 }
925
926 tmp_OLPpoz = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
927 for (i=0; i<List_YOUSO[7]; i++){
928 tmp_OLPpoz[i] = (double*)malloc(sizeof(double)*List_YOUSO[7]);
929 }
930
931 /*****************************************************
932 calculation of matrix elements for OLPpox,y,z
933 *****************************************************/
934
935 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
936
937 dtime(&Stime_atom);
938
939 Gc_AN = M2G[Mc_AN];
940 Cwan = WhatSpecies[Gc_AN];
941 NO0 = Spe_Total_CNO[Cwan];
942
943 for (i=0; i<NO0; i++){
944 for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
945
946 GNc = GridListAtom[Mc_AN][Nc];
947 GRc = CellListAtom[Mc_AN][Nc];
948
949 Get_Grid_XYZ(GNc,Cxyz);
950 x = Cxyz[1] + atv[GRc][1] - Gxyz[Gc_AN][1];
951 y = Cxyz[2] + atv[GRc][2] - Gxyz[Gc_AN][2];
952 z = Cxyz[3] + atv[GRc][3] - Gxyz[Gc_AN][3];
953
954 ChiVx[i][Nc] = x*Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
955 ChiVy[i][Nc] = y*Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
956 ChiVz[i][Nc] = z*Orbs_Grid[Mc_AN][Nc][i];/* AITUNE */
957 }
958 }
959
960 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
961
962 Gh_AN = natn[Gc_AN][h_AN];
963 Mh_AN = F_G2M[Gh_AN];
964
965 Rnh = ncn[Gc_AN][h_AN];
966 Hwan = WhatSpecies[Gh_AN];
967 NO1 = Spe_Total_CNO[Hwan];
968
969 /* initialize */
970
971 for (i=0; i<NO0; i++){
972 for (j=0; j<NO1; j++){
973 tmp_OLPpox[i][j] = 0.0;
974 tmp_OLPpoy[i][j] = 0.0;
975 tmp_OLPpoz[i][j] = 0.0;
976 }
977 }
978
979 /* summation of non-zero elements */
980
981 for (Nog=0; Nog<NumOLG[Mc_AN][h_AN]; Nog++){
982
983 Nc = GListTAtoms1[Mc_AN][h_AN][Nog];
984 Nh = GListTAtoms2[Mc_AN][h_AN][Nog];
985
986 /* store ChiVx,y,z in tmp_ChiVx,y,z */
987
988 for (i=0; i<NO0; i++){
989 tmp_ChiVx[i] = ChiVx[i][Nc];
990 tmp_ChiVy[i] = ChiVy[i][Nc];
991 tmp_ChiVz[i] = ChiVz[i][Nc];
992 }
993
994 /* store Orbs_Grid in tmp_Orbs_Grid */
995
996 if (G2ID[Gh_AN]==myid){
997 for (j=0; j<NO1; j++){
998 tmp_Orbs_Grid[j] = Orbs_Grid[Mh_AN][Nh][j];/* AITUNE */
999 }
1000 }
1001 else{
1002 for (j=0; j<NO1; j++){
1003 tmp_Orbs_Grid[j] = Orbs_Grid_FNAN[Mc_AN][h_AN][Nog][j];/* AITUNE */
1004 }
1005 }
1006
1007 /* integration */
1008
1009 for (i=0; i<NO0; i++){
1010 tmpx = tmp_ChiVx[i];
1011 tmpy = tmp_ChiVy[i];
1012 tmpz = tmp_ChiVz[i];
1013 for (j=0; j<NO1; j++){
1014 tmp_OLPpox[i][j] += tmpx*tmp_Orbs_Grid[j];
1015 tmp_OLPpoy[i][j] += tmpy*tmp_Orbs_Grid[j];
1016 tmp_OLPpoz[i][j] += tmpz*tmp_Orbs_Grid[j];
1017 }
1018 }
1019 }
1020
1021 /* OLPpox,y,z */
1022
1023 for (i=0; i<NO0; i++){
1024 for (j=0; j<NO1; j++){
1025 OLPpox[Mc_AN][h_AN][i][j] = tmp_OLPpox[i][j]*GridVol;
1026 OLPpoy[Mc_AN][h_AN][i][j] = tmp_OLPpoy[i][j]*GridVol;
1027 OLPpoz[Mc_AN][h_AN][i][j] = tmp_OLPpoz[i][j]*GridVol;
1028 }
1029 }
1030
1031 }
1032
1033 dtime(&Etime_atom);
1034 time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
1035 }
1036
1037 /****************************************************
1038 freeing of arrays:
1039 ****************************************************/
1040
1041 for (i=0; i<List_YOUSO[7]; i++){
1042 free(ChiVx[i]);
1043 }
1044 free(ChiVx);
1045
1046 for (i=0; i<List_YOUSO[7]; i++){
1047 free(ChiVy[i]);
1048 }
1049 free(ChiVy);
1050
1051 for (i=0; i<List_YOUSO[7]; i++){
1052 free(ChiVz[i]);
1053 }
1054 free(ChiVz);
1055
1056 free(tmp_ChiVx);
1057 free(tmp_ChiVy);
1058 free(tmp_ChiVz);
1059
1060 free(tmp_Orbs_Grid);
1061
1062 for (i=0; i<List_YOUSO[7]; i++){
1063 free(tmp_OLPpox[i]);
1064 }
1065 free(tmp_OLPpox);
1066
1067 for (i=0; i<List_YOUSO[7]; i++){
1068 free(tmp_OLPpoy[i]);
1069 }
1070 free(tmp_OLPpoy);
1071
1072 for (i=0; i<List_YOUSO[7]; i++){
1073 free(tmp_OLPpoz[i]);
1074 }
1075 free(tmp_OLPpoz);
1076 }
1077