1 /**********************************************************************
2 Initial_CntCoes.c:
3
4 Initial_CntCoes.c is a subroutine to prepare initial contraction
5 coefficients for the orbital optimization method.
6
7 Log of Initial_CntCoes.c:
8
9 22/Nov/2001 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <math.h>
16 #include <time.h>
17 #include "openmx_common.h"
18 #include "mpi.h"
19
20
Initial_CntCoes(double ***** nh,double ***** OLP)21 void Initial_CntCoes(double *****nh, double *****OLP)
22 {
23 static int firsttime=1;
24 int i,j,l,n,n2,i1,j1;
25 int wan;
26 int P_min,po;
27 int Second_switch;
28 double time0;
29 int Mc_AN,Gc_AN,wanA;
30 int q,q0,al0,al1,pmax;
31 int Mul0,Mul1,deg_on,deg_num;
32 int al,p,L0,M0,p0;
33 int ig,im,jg,ian,jan,kl,m,Gi;
34 int mu,nu,Anum,Bnum,NUM,NUM1;
35 int h_AN,Gh_AN,Hwan,tno1,tno2,Cwan,spin;
36 double Beta0,scaleF;
37 double *ko,*C0;
38 double **S,**Hks,**D,*abs_sum,*M1,**C,**B;
39 int *jun,*ponu;
40 double tmp0,tmp1,Max0,rc1,fugou,MaxV;
41 double sum,TZ;
42 double Num_State,x,FermiF,Dnum;
43 double LChemP_MAX,LChemP_MIN,LChemP;
44 double TStime,TEtime;
45
46 double *tmp_array;
47 double *tmp_array2;
48 int *MP,*dege;
49 int **tmp_index;
50 int ***tmp_index2;
51 double *Tmp_CntCoes;
52 double **Check_ko;
53 double *Weight_ko;
54 double ***CntCoes_Spe;
55 double ***My_CntCoes_Spe;
56 int *Snd_CntCoes_Size;
57 int *Rcv_CntCoes_Size;
58 int *Snd_H_Size,*Rcv_H_Size;
59 int *Snd_S_Size,*Rcv_S_Size;
60 int size1,size2,num;
61 int numprocs,myid,tag=999,ID,IDS,IDR;
62
63 MPI_Status stat;
64 MPI_Request request;
65
66 /* MPI */
67 MPI_Comm_size(mpi_comm_level1,&numprocs);
68 MPI_Comm_rank(mpi_comm_level1,&myid);
69
70 dtime(&TStime);
71
72 /****************************************************
73 allocation of arrays:
74
75 int MP[List_YOUSO[8]];
76
77 int tmp_index[List_YOUSO[25]+1]
78 [2*(List_YOUSO[25]+1)+1];
79 int tmp_index2[List_YOUSO[25]+1]
80 [List_YOUSO[24]]
81 [2*(List_YOUSO[25]+1)+1];
82
83 double Tmp_CntCoes[List_YOUSO[24]]
84
85 double Check_ko[List_YOUSO[25]+1]
86 [2*(List_YOUSO[25]+1)+1];
87
88 double Weight_ko[List_YOUSO[7]];
89
90 ****************************************************/
91
92 MP = (int*)malloc(sizeof(int)*List_YOUSO[8]);
93
94 tmp_index = (int**)malloc(sizeof(int*)*(List_YOUSO[25]+1));
95 for (i=0; i<(List_YOUSO[25]+1); i++){
96 tmp_index[i] = (int*)malloc(sizeof(int)*(2*(List_YOUSO[25]+1)+1));
97 }
98
99 tmp_index2 = (int***)malloc(sizeof(int**)*(List_YOUSO[25]+1));
100 for (i=0; i<(List_YOUSO[25]+1); i++){
101 tmp_index2[i] = (int**)malloc(sizeof(int*)*List_YOUSO[24]);
102 for (j=0; j<List_YOUSO[24]; j++){
103 tmp_index2[i][j] = (int*)malloc(sizeof(int)*(2*(List_YOUSO[25]+1)+1));
104 }
105 }
106
107 Tmp_CntCoes = (double*)malloc(sizeof(double)*List_YOUSO[24]);
108
109 Check_ko = (double**)malloc(sizeof(double*)*(List_YOUSO[25]+1));
110 for (i=0; i<(List_YOUSO[25]+1); i++){
111 Check_ko[i] = (double*)malloc(sizeof(double)*(2*(List_YOUSO[25]+1)+1));
112 }
113
114 Weight_ko = (double*)malloc(sizeof(double)*List_YOUSO[7]);
115
116 Snd_CntCoes_Size = (int*)malloc(sizeof(int)*Num_Procs);
117 Rcv_CntCoes_Size = (int*)malloc(sizeof(int)*Num_Procs);
118 Snd_H_Size = (int*)malloc(sizeof(int)*Num_Procs);
119 Rcv_H_Size = (int*)malloc(sizeof(int)*Num_Procs);
120 Snd_S_Size = (int*)malloc(sizeof(int)*Num_Procs);
121 Rcv_S_Size = (int*)malloc(sizeof(int)*Num_Procs);
122
123 /* PrintMemory */
124
125 if (firsttime) {
126 PrintMemory("Initial_CntCoes: tmp_index",sizeof(int)*(List_YOUSO[25]+1)*
127 (2*(List_YOUSO[25]+1)+1),NULL);
128 PrintMemory("Initial_CntCoes: tmp_index2",sizeof(int)*(List_YOUSO[25]+1)*
129 List_YOUSO[24]*(2*(List_YOUSO[25]+1)+1) ,NULL);
130 PrintMemory("Initial_CntCoes: Check_ko",sizeof(double)*(List_YOUSO[25]+1)*
131 (2*(List_YOUSO[25]+1)+1),NULL);
132 firsttime=0;
133 }
134
135 /****************************************************
136 MPI:
137
138 nh(=H)
139 ****************************************************/
140
141 /***********************************
142 set data size
143 ************************************/
144
145 for (ID=0; ID<numprocs; ID++){
146
147 IDS = (myid + ID) % numprocs;
148 IDR = (myid - ID + numprocs) % numprocs;
149
150 if (ID!=0){
151 tag = 999;
152
153 /* find data size to send block data */
154 if (F_Snd_Num[IDS]!=0){
155
156 size1 = 0;
157 for (spin=0; spin<=SpinP_switch; spin++){
158 for (n=0; n<F_Snd_Num[IDS]; n++){
159 Mc_AN = Snd_MAN[IDS][n];
160 Gc_AN = Snd_GAN[IDS][n];
161 Cwan = WhatSpecies[Gc_AN];
162 tno1 = Spe_Total_NO[Cwan];
163 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
164 Gh_AN = natn[Gc_AN][h_AN];
165 Hwan = WhatSpecies[Gh_AN];
166 tno2 = Spe_Total_NO[Hwan];
167 size1 += tno1*tno2;
168 }
169 }
170 }
171
172 Snd_H_Size[IDS] = size1;
173 MPI_Isend(&size1, 1, MPI_INT, IDS, tag, mpi_comm_level1, &request);
174 }
175 else{
176 Snd_H_Size[IDS] = 0;
177 }
178
179 /* receiving of size of data */
180
181 if (F_Rcv_Num[IDR]!=0){
182 MPI_Recv(&size2, 1, MPI_INT, IDR, tag, mpi_comm_level1, &stat);
183 Rcv_H_Size[IDR] = size2;
184 }
185 else{
186 Rcv_H_Size[IDR] = 0;
187 }
188
189 if (F_Snd_Num[IDS]!=0) MPI_Wait(&request,&stat);
190
191 }
192 else{
193 Snd_H_Size[IDS] = 0;
194 Rcv_H_Size[IDR] = 0;
195 }
196 }
197
198 /***********************************
199 data transfer
200 ************************************/
201
202 tag = 999;
203 for (ID=0; ID<numprocs; ID++){
204
205 IDS = (myid + ID) % numprocs;
206 IDR = (myid - ID + numprocs) % numprocs;
207
208 if (ID!=0){
209
210 /*****************************
211 sending of data
212 *****************************/
213
214 if (F_Snd_Num[IDS]!=0){
215
216 size1 = Snd_H_Size[IDS];
217
218 /* allocation of array */
219
220 tmp_array = (double*)malloc(sizeof(double)*size1);
221
222 /* multidimentional array to vector array */
223
224 num = 0;
225 for (spin=0; spin<=SpinP_switch; spin++){
226 for (n=0; n<F_Snd_Num[IDS]; n++){
227 Mc_AN = Snd_MAN[IDS][n];
228 Gc_AN = Snd_GAN[IDS][n];
229 Cwan = WhatSpecies[Gc_AN];
230 tno1 = Spe_Total_NO[Cwan];
231 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
232 Gh_AN = natn[Gc_AN][h_AN];
233 Hwan = WhatSpecies[Gh_AN];
234 tno2 = Spe_Total_NO[Hwan];
235 for (i=0; i<tno1; i++){
236 for (j=0; j<tno2; j++){
237 tmp_array[num] = nh[spin][Mc_AN][h_AN][i][j];
238 num++;
239 }
240 }
241 }
242 }
243 }
244
245 MPI_Isend(&tmp_array[0], size1, MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request);
246
247 }
248
249 /*****************************
250 receiving of block data
251 *****************************/
252
253 if (F_Rcv_Num[IDR]!=0){
254
255 size2 = Rcv_H_Size[IDR];
256
257 /* allocation of array */
258 tmp_array2 = (double*)malloc(sizeof(double)*size2);
259
260 MPI_Recv(&tmp_array2[0], size2, MPI_DOUBLE, IDR, tag, mpi_comm_level1, &stat);
261
262 num = 0;
263 for (spin=0; spin<=SpinP_switch; spin++){
264 Mc_AN = S_TopMAN[IDR] - 1; /* S_TopMAN should be used. */
265 for (n=0; n<F_Rcv_Num[IDR]; n++){
266 Mc_AN++;
267 Gc_AN = Rcv_GAN[IDR][n];
268 Cwan = WhatSpecies[Gc_AN];
269 tno1 = Spe_Total_NO[Cwan];
270
271 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
272 Gh_AN = natn[Gc_AN][h_AN];
273 Hwan = WhatSpecies[Gh_AN];
274 tno2 = Spe_Total_NO[Hwan];
275
276 for (i=0; i<tno1; i++){
277 for (j=0; j<tno2; j++){
278 nh[spin][Mc_AN][h_AN][i][j] = tmp_array2[num];
279 num++;
280 }
281 }
282 }
283 }
284 }
285
286 /* freeing of array */
287 free(tmp_array2);
288
289 }
290
291 if (F_Snd_Num[IDS]!=0){
292 MPI_Wait(&request,&stat);
293 free(tmp_array); /* freeing of array */
294 }
295 }
296 }
297
298 /****************************************************
299 MPI:
300
301 OLP[0]
302 ****************************************************/
303
304 /***********************************
305 set data size
306 ************************************/
307
308 for (ID=0; ID<numprocs; ID++){
309
310 IDS = (myid + ID) % numprocs;
311 IDR = (myid - ID + numprocs) % numprocs;
312
313 if (ID!=0){
314 tag = 999;
315
316 /* find data size to send block data */
317 if (F_Snd_Num[IDS]!=0){
318
319 size1 = 0;
320 for (n=0; n<F_Snd_Num[IDS]; n++){
321 Mc_AN = Snd_MAN[IDS][n];
322 Gc_AN = Snd_GAN[IDS][n];
323 Cwan = WhatSpecies[Gc_AN];
324 tno1 = Spe_Total_NO[Cwan];
325 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
326 Gh_AN = natn[Gc_AN][h_AN];
327 Hwan = WhatSpecies[Gh_AN];
328 tno2 = Spe_Total_NO[Hwan];
329 size1 += tno1*tno2;
330 }
331 }
332
333 Snd_S_Size[IDS] = size1;
334 MPI_Isend(&size1, 1, MPI_INT, IDS, tag, mpi_comm_level1, &request);
335 }
336 else{
337 Snd_S_Size[IDS] = 0;
338 }
339
340 /* receiving of size of data */
341
342 if (F_Rcv_Num[IDR]!=0){
343 MPI_Recv(&size2, 1, MPI_INT, IDR, tag, mpi_comm_level1, &stat);
344 Rcv_S_Size[IDR] = size2;
345 }
346 else{
347 Rcv_S_Size[IDR] = 0;
348 }
349
350 if (F_Snd_Num[IDS]!=0) MPI_Wait(&request,&stat);
351
352 }
353 else{
354 Snd_S_Size[IDS] = 0;
355 Rcv_S_Size[IDR] = 0;
356 }
357 }
358
359 /***********************************
360 data transfer
361 ************************************/
362
363 tag = 999;
364 for (ID=0; ID<numprocs; ID++){
365
366 IDS = (myid + ID) % numprocs;
367 IDR = (myid - ID + numprocs) % numprocs;
368
369 if (ID!=0){
370
371 /*****************************
372 sending of data
373 *****************************/
374
375 if (F_Snd_Num[IDS]!=0){
376
377 size1 = Snd_S_Size[IDS];
378
379 /* allocation of array */
380
381 tmp_array = (double*)malloc(sizeof(double)*size1);
382
383 /* multidimentional array to vector array */
384
385 num = 0;
386
387 for (n=0; n<F_Snd_Num[IDS]; n++){
388 Mc_AN = Snd_MAN[IDS][n];
389 Gc_AN = Snd_GAN[IDS][n];
390 Cwan = WhatSpecies[Gc_AN];
391 tno1 = Spe_Total_NO[Cwan];
392 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
393 Gh_AN = natn[Gc_AN][h_AN];
394 Hwan = WhatSpecies[Gh_AN];
395 tno2 = Spe_Total_NO[Hwan];
396 for (i=0; i<tno1; i++){
397 for (j=0; j<tno2; j++){
398 tmp_array[num] = OLP[0][Mc_AN][h_AN][i][j];
399 num++;
400 }
401 }
402 }
403 }
404
405 MPI_Isend(&tmp_array[0], size1, MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request);
406 }
407
408 /*****************************
409 receiving of block data
410 *****************************/
411
412 if (F_Rcv_Num[IDR]!=0){
413
414 size2 = Rcv_S_Size[IDR];
415
416 /* allocation of array */
417 tmp_array2 = (double*)malloc(sizeof(double)*size2);
418
419 MPI_Recv(&tmp_array2[0], size2, MPI_DOUBLE, IDR, tag, mpi_comm_level1, &stat);
420
421 num = 0;
422 Mc_AN = S_TopMAN[IDR] - 1; /* S_TopMAN should be used. */
423 for (n=0; n<F_Rcv_Num[IDR]; n++){
424 Mc_AN++;
425 Gc_AN = Rcv_GAN[IDR][n];
426 Cwan = WhatSpecies[Gc_AN];
427 tno1 = Spe_Total_NO[Cwan];
428
429 for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
430 Gh_AN = natn[Gc_AN][h_AN];
431 Hwan = WhatSpecies[Gh_AN];
432 tno2 = Spe_Total_NO[Hwan];
433 for (i=0; i<tno1; i++){
434 for (j=0; j<tno2; j++){
435 OLP[0][Mc_AN][h_AN][i][j] = tmp_array2[num];
436 num++;
437 }
438 }
439 }
440 }
441
442 /* freeing of array */
443 free(tmp_array2);
444
445 }
446
447 if (F_Snd_Num[IDS]!=0){
448 MPI_Wait(&request,&stat);
449 free(tmp_array); /* freeing of array */
450 }
451 }
452 }
453
454 /****************************************************
455 set of "initial" coefficients of the contraction
456 ****************************************************/
457
458 Second_switch = 1;
459
460 Simple_InitCnt[0] = 0;
461 Simple_InitCnt[1] = 0;
462 Simple_InitCnt[2] = 0;
463 Simple_InitCnt[3] = 0;
464 Simple_InitCnt[4] = 0;
465 Simple_InitCnt[5] = 0;
466
467 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
468 Gc_AN = M2G[Mc_AN];
469 wan = WhatSpecies[Gc_AN];
470
471 for (al=0; al<Spe_Total_CNO[wan]; al++){
472 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
473 CntCoes[Mc_AN][al][p] = 0.0;
474 }
475 }
476
477 al = -1;
478 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
479 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
480 for (M0=0; M0<=2*L0; M0++){
481 al++;
482 CntCoes[Mc_AN][al][Mul0] = 1.0;
483 }
484 }
485 }
486 }
487
488 /****************************************************
489 Setting of Hamiltonian and overlap matrices
490
491 MP indicates the starting position of
492 atom i in arraies H and S
493 ****************************************************/
494
495 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
496 Gc_AN = M2G[Mc_AN];
497 wan = WhatSpecies[Gc_AN];
498
499 if ( FNAN[Gc_AN]<40 ) scaleF = 2.0;
500 else if ( FNAN[Gc_AN]<55 ) scaleF = 1.6;
501 else if ( FNAN[Gc_AN]<82 ) scaleF = 1.2;
502 else scaleF = 1.0;
503
504 rc1 = scaleF*Spe_Atom_Cut1[wan];
505
506 /***********************************************
507 MP indicates the starting position of
508 atom i in arraies H and S
509 ***********************************************/
510
511 Anum = 1;
512 TZ = 0.0;
513 for (i=0; i<=FNAN[Gc_AN]; i++){
514 if (Dis[Gc_AN][i]<=rc1){
515 MP[i] = Anum;
516 Gi = natn[Gc_AN][i];
517 wanA = WhatSpecies[Gi];
518 Anum = Anum + Spe_Total_NO[wanA];
519 TZ = TZ + Spe_Core_Charge[wanA];
520 }
521 }
522 NUM = Anum - 1;
523
524 /****************************************************
525 allocation
526 ****************************************************/
527
528 n2 = NUM + 3;
529
530 ko = (double*)malloc(sizeof(double)*n2);
531 abs_sum = (double*)malloc(sizeof(double)*n2);
532 M1 = (double*)malloc(sizeof(double)*n2);
533 dege = (int*)malloc(sizeof(int)*n2);
534 C0 = (double*)malloc(sizeof(double)*n2);
535
536 S = (double**)malloc(sizeof(double*)*n2);
537 Hks = (double**)malloc(sizeof(double*)*n2);
538 D = (double**)malloc(sizeof(double*)*n2);
539 C = (double**)malloc(sizeof(double*)*n2);
540 B = (double**)malloc(sizeof(double*)*n2);
541
542 for (i=0; i<n2; i++){
543 S[i] = (double*)malloc(sizeof(double)*n2);
544 Hks[i] = (double*)malloc(sizeof(double)*n2);
545 D[i] = (double*)malloc(sizeof(double)*n2);
546 C[i] = (double*)malloc(sizeof(double)*n2);
547 B[i] = (double*)malloc(sizeof(double)*n2);
548 }
549
550 jun = (int*)malloc(sizeof(int)*n2);
551 ponu = (int*)malloc(sizeof(int)*n2);
552
553 /****************************************************
554 calc
555 ****************************************************/
556
557 if (3<=level_stdout){
558 printf("<Initial_CntCoes> Mc_AN=%2d Gc_AN=%2d NUM=%2d\n",Mc_AN,Gc_AN,NUM);
559 }
560
561 for (i=0; i<=FNAN[Gc_AN]; i++){
562
563 if (Dis[Gc_AN][i]<=rc1){
564
565 ig = natn[Gc_AN][i];
566 im = S_G2M[ig]; /* S_G2M must be used. */
567
568 ian = Spe_Total_NO[WhatSpecies[ig]];
569 Anum = MP[i];
570
571 for (j=0; j<=FNAN[Gc_AN]; j++){
572
573 if (Dis[Gc_AN][j]<=rc1){
574
575 kl = RMI1[Mc_AN][i][j];
576 jg = natn[Gc_AN][j];
577 jan = Spe_Total_NO[WhatSpecies[jg]];
578 Bnum = MP[j];
579
580 if (0<=kl){
581 for (m=0; m<ian; m++){
582 for (n=0; n<jan; n++){
583
584 S[Anum+m][Bnum+n] = OLP[0][im][kl][m][n];
585
586 if (SpinP_switch==0)
587 Hks[Anum+m][Bnum+n] = nh[0][im][kl][m][n];
588 else
589 Hks[Anum+m][Bnum+n] = 0.5*(nh[0][im][kl][m][n]
590 + nh[1][im][kl][m][n]);
591 }
592 }
593 }
594
595 else{
596 for (m=0; m<ian; m++){
597 for (n=0; n<jan; n++){
598 S[Anum+m][Bnum+n] = 0.0;
599 Hks[Anum+m][Bnum+n] = 0.0;
600 }
601 }
602 }
603 }
604 }
605 }
606 }
607
608 /***********************************************
609 Solve the generalized eigenvalue problem
610 ***********************************************/
611
612 Eigen_lapack(S,ko,NUM,NUM);
613
614 /***********************************************
615 Searching of negative eigenvalues
616 ************************************************/
617
618 P_min = 1;
619 for (l=1; l<=NUM; l++){
620 if (ko[l]<0.0){
621 P_min = l + 1;
622 if (3<=level_stdout){
623 printf("<Init_CntCoes> Negative EV of OLP %2d %15.12f\n",l,ko[l]);
624 }
625 }
626 }
627 for (l=P_min; l<=NUM; l++){
628 M1[l] = 1.0/sqrt(ko[l]);
629 }
630
631 for (i1=1; i1<=NUM; i1++){
632 for (j1=P_min; j1<=NUM; j1++){
633 sum = 0.0;
634 for (l=1; l<=NUM; l++){
635 sum = sum + Hks[i1][l]*S[l][j1]*M1[j1];
636 }
637 C[i1][j1] = sum;
638 }
639 }
640
641 for (i1=P_min; i1<=NUM; i1++){
642 for (j1=1; j1<=NUM; j1++){
643 sum = 0.0;
644 for (l=1; l<=NUM; l++){
645 sum = sum + M1[i1]*S[l][i1]*C[l][j1];
646 }
647 B[i1][j1] = sum;
648 }
649 }
650
651 for (i1=P_min; i1<=NUM; i1++){
652 for (j1=P_min; j1<=NUM; j1++){
653 D[i1-(P_min-1)][j1-(P_min-1)] = B[i1][j1];
654 }
655 }
656
657 NUM1 = NUM - (P_min - 1);
658 Eigen_lapack(D,ko,NUM1,NUM1);
659
660 /***********************************************
661 transformation to the original eigen vectors.
662 NOTE 244P
663 ***********************************************/
664
665 for (i1=1; i1<=NUM; i1++){
666 for (j1=1; j1<=NUM; j1++){
667 C[i1][j1] = 0.0;
668 }
669 }
670
671 for (i1=1; i1<=NUM; i1++){
672 for (j1=1; j1<=NUM1; j1++){
673 sum = 0.0;
674 for (l=P_min; l<=NUM; l++){
675 sum = sum + S[i1][l]*M1[l]*D[l-(P_min-1)][j1];
676 }
677 C[i1][j1] = sum;
678 }
679 }
680
681 /****************************************************
682 searching of a local chemical potential
683 ****************************************************/
684
685 po = 0;
686 LChemP_MAX = 10.0;
687 LChemP_MIN =-10.0;
688
689 Beta0 = 1.0/(kB*1500.0/eV2Hartree);
690
691 do{
692 LChemP = 0.50*(LChemP_MAX + LChemP_MIN);
693 Num_State = 0.0;
694 for (i1=1; i1<=NUM; i1++){
695 x = (ko[i1] - LChemP)*Beta0;
696 if (x<=-30.0) x = -30.0;
697 if (30.0<=x) x = 30.0;
698 FermiF = 2.0/(1.0 + exp(x));
699 Num_State = Num_State + FermiF;
700 }
701 Dnum = TZ - Num_State;
702 if (0.0<=Dnum) LChemP_MIN = LChemP;
703 else LChemP_MAX = LChemP;
704 if (fabs(Dnum)<0.000000000001) po = 1;
705 }
706 while (po==0);
707
708 if (3<=level_stdout){
709 for (i1=1; i1<=NUM; i1++){
710 x = (ko[i1] - LChemP)*Beta0;
711 if (x<=-30.0) x = -30.0;
712 if (30.0<=x) x = 30.0;
713 FermiF = 1.0/(1.0 + exp(x));
714 printf("<Init_CntCoes> %2d eigenvalue=%15.12f FermiF=%15.12f\n",
715 i1,ko[i1],FermiF);
716 }
717 }
718
719 if (3<=level_stdout){
720 printf("\n");
721 printf("first C Gc_AN=%2d\n",Gc_AN);
722 for (i1=1; i1<=NUM; i1++){
723 for (j1=1; j1<=10; j1++){
724 printf("%6.3f ",C[i1][j1]);
725 }
726 printf("\n");
727 }
728 }
729
730 if (3<=level_stdout){
731 printf("<Init_CntCoes> LChemP=%15.12f\n",LChemP);
732 }
733
734 /************************************************
735 find degenerate eigenvalues
736 ************************************************/
737
738 for (i1=0; i1<=NUM1; i1++) dege[i1] = 0;
739
740 i1 = 0;
741 l = 0;
742 deg_on = 0;
743
744 do {
745 i1++;
746
747 if (fabs(ko[i1]-ko[i1+1])<1.0e-7){
748 if (deg_on==0) deg_on = i1;
749 }
750
751 else {
752
753 if (deg_on!=0){
754 l++;
755 for (j1=deg_on; j1<=i1; j1++){
756 dege[j1] = l;
757 }
758 }
759
760 deg_on = 0;
761 }
762
763 } while(i1<NUM1);
764
765 if (3<=level_stdout){
766 for (i1=1; i1<=NUM1; i1++){
767 printf("i1=%i dege=%i\n",i1,dege[i1]);
768 }
769 }
770
771 /************************************************
772 avaraging of degenerate eigenvectors
773 ************************************************/
774
775 deg_on = 0;
776 deg_num = 0;
777
778 for (i1=1; i1<=NUM1; i1++){
779
780 if (dege[i1]!=0) {
781
782 for (j1=1; j1<=NUM; j1++){
783 if (deg_on==0){
784 C0[j1] = C[j1][i1]*C[j1][i1];
785 }
786 else if (dege[i1-1]!=dege[i1]){
787 C0[j1] = C[j1][i1]*C[j1][i1];
788 }
789 else{
790 C0[j1] += C[j1][i1]*C[j1][i1];
791 }
792 }
793
794 if (dege[i1-1]!=dege[i1]) deg_on = 0;
795
796 if (deg_on==0){
797 deg_on = i1;
798 deg_num = 0;
799 }
800
801 deg_num++;
802 }
803
804 if ( deg_on!=0 && (dege[i1]!=dege[i1+1]) ){
805
806 tmp0 = 1.0/sqrt((double)deg_num);
807 for (j1=1; j1<=NUM; j1++){
808 C0[j1] = sqrt(fabs(C0[j1]))*tmp0;
809 }
810
811 for (l=deg_on; l<(deg_on+deg_num); l++){
812 for (j1=1; j1<=NUM; j1++){
813 C[j1][l] = sgn(C[j1][l])*C0[j1];
814 }
815 }
816
817 deg_on = 0;
818 }
819
820 }
821
822 if (3<=level_stdout){
823 printf("\n");
824 printf("averaged C\n");
825 for (i1=1; i1<=NUM; i1++){
826 for (j1=1; j1<=10; j1++){
827 printf("%6.3f ",C[i1][j1]);
828 }
829 printf("\n");
830 }
831 }
832
833 /************************************************
834 Contraction
835 ************************************************/
836
837 for (i=1; i<=NUM; i++){
838 for (j=1; j<=NUM; j++){
839 B[i][j] = 0.0;
840 }
841 }
842
843 al = -1;
844 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
845 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
846 for (M0=0; M0<=2*L0; M0++){
847 al++;
848 tmp_index2[L0][Mul0][M0] = al;
849 }
850 }
851 }
852
853 Beta0 = 1.0/(kB*1500.0/eV2Hartree);
854
855 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
856
857 if (0<Spe_Num_Basis[wan][L0]){
858
859 for (M0=0; M0<=2*L0; M0++){
860
861 for (mu=1; mu<=NUM1; mu++){
862
863 p = tmp_index2[L0][0][M0];
864 fugou = sgn(C[p+1][mu]);
865 x = (ko[mu] - LChemP)*Beta0;
866 if (x<=-30.0) x = -30.0;
867 if (30.0<=x) x = 30.0;
868 FermiF = 1.0/(1.0 + exp(x));
869
870 sum = 0.0;
871 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
872 al = tmp_index2[L0][Mul0][M0];
873 B[al+1][1] = B[al+1][1] + fugou*FermiF*C[al+1][mu];
874 sum = sum + FermiF*fabs(C[al+1][mu]);
875 }
876
877 abs_sum[mu] = sum;
878 }
879
880 if (Second_switch==0){
881
882 for (mu=1; mu<=NUM1; mu++) ponu[mu] = 0;
883 for (mu=1; mu<=NUM1; mu++){
884 MaxV = -1000.0;
885 for (nu=1; nu<=NUM1; nu++){
886 if (MaxV<abs_sum[nu] && ponu[nu]!=1){
887 MaxV = abs_sum[nu];
888 p = nu;
889 }
890 }
891 jun[mu] = p;
892 ponu[p] = 1;
893 }
894
895
896 if (3<=level_stdout){
897 for (mu=1; mu<=NUM1; mu++){
898 printf("L0=%2d M0= %2d mu=%2d jun=%2d\n",L0,M0,mu,jun[mu]);
899 }
900 }
901
902 for (mu=1; mu<=NUM1; mu++){
903 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
904 al = tmp_index2[L0][Mul0][M0];
905 B[al+1][mu+1] = C[al+1][jun[mu]];
906 }
907 }
908 }
909
910 else if (Second_switch==1){
911
912 for (mu=1; mu<=NUM1; mu++){
913 for (Mul0=0; Mul0<Spe_Num_Basis[wan][L0]; Mul0++){
914 al = tmp_index2[L0][Mul0][M0];
915 if (mu==Mul0) B[al+1][mu+1] = 1.0;
916 }
917 }
918 }
919 }
920 }
921 }
922
923 if (3<=level_stdout){
924 printf("B\n");
925 for (i1=1; i1<=NUM; i1++){
926 for (j1=1; j1<=10; j1++){
927 printf("%6.3f ",B[i1][j1]);
928 }
929 printf("\n");
930 }
931 }
932
933
934 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
935 for (M0=0; M0<=2*L0; M0++){
936 tmp_index[L0][M0] = 1;
937 Check_ko[L0][M0] = 10e+10;
938 }
939 }
940
941 al = -1;
942 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
943 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
944 for (M0=0; M0<=2*L0; M0++){
945
946 al++;
947
948 po = 0;
949 i = tmp_index[L0][M0] - 1;
950 do{
951 i++;
952
953 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
954 p0 = Spe_Trans_Orbital[wan][al][p];
955
956 if (0.001<fabs(B[p0+1][i])){
957 tmp_index[L0][M0] = i;
958 Check_ko[L0][M0] = ko[i];
959 Weight_ko[al] = ko[i];
960 po = 1;
961 }
962 if (po==1) break;
963 }
964
965 }while(po==0);
966
967 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
968 p0 = Spe_Trans_Orbital[wan][al][p];
969 i = tmp_index[L0][M0];
970 if (Simple_InitCnt[L0]==0)
971 CntCoes[Mc_AN][al][p] = B[p0+1][i];
972 }
973
974 tmp_index[L0][M0] = tmp_index[L0][M0] + 1;
975
976 }
977 }
978 }
979
980 if (3<=level_stdout){
981 for (al=0; al<Spe_Total_CNO[wan]; al++){
982 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
983 printf("A Init_CntCoes Mc_AN=%2d Gc_AN=%2d al=%2d p=%2d %15.12f\n",
984 Mc_AN,Gc_AN,al,p,CntCoes[Mc_AN][al][p]);
985 }
986 }
987 }
988
989 /****************************************************
990 Orthonormalization by the Gram-Schmidt method
991 ****************************************************/
992
993 al = -1;
994 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
995 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
996 for (M0=0; M0<=2*L0; M0++){
997 al++;
998 tmp_index2[L0][Mul0][M0] = al;
999 }
1000 }
1001 }
1002
1003 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
1004 for (M0=0; M0<=2*L0; M0++){
1005
1006 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
1007 al0 = tmp_index2[L0][Mul0][M0];
1008
1009 /* x - sum_i <x|e_i>e_i */
1010
1011 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1012 Tmp_CntCoes[p] = 0.0;
1013 }
1014
1015 for (Mul1=0; Mul1<Mul0; Mul1++){
1016 al1 = tmp_index2[L0][Mul1][M0];
1017
1018 sum = 0.0;
1019 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1020 sum = sum + CntCoes[Mc_AN][al0][p]*CntCoes[Mc_AN][al1][p];
1021 }
1022
1023 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1024 Tmp_CntCoes[p] = Tmp_CntCoes[p] + sum*CntCoes[Mc_AN][al1][p];
1025 }
1026 }
1027
1028 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1029 CntCoes[Mc_AN][al0][p] = CntCoes[Mc_AN][al0][p] - Tmp_CntCoes[p];
1030 }
1031
1032 /* Normalize */
1033
1034 sum = 0.0;
1035 Max0 = -100.0;
1036 pmax = 0;
1037 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1038 sum = sum + CntCoes[Mc_AN][al0][p]*CntCoes[Mc_AN][al0][p];
1039 if (Max0<fabs(CntCoes[Mc_AN][al0][p])){
1040 Max0 = fabs(CntCoes[Mc_AN][al0][p]);
1041 pmax = p;
1042 }
1043 }
1044
1045 if (fabs(sum)<1.0e-11)
1046 tmp0 = 0.0;
1047 else
1048 tmp0 = 1.0/sqrt(sum);
1049
1050 tmp1 = sgn(CntCoes[Mc_AN][al0][pmax]);
1051
1052 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1053 CntCoes[Mc_AN][al0][p] = tmp0*tmp1*CntCoes[Mc_AN][al0][p];
1054 }
1055
1056 }
1057 }
1058 }
1059
1060 if (3<=level_stdout){
1061 for (al=0; al<Spe_Total_CNO[wan]; al++){
1062 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1063 printf("B Init_CntCoes Mc_AN=%2d Gc_AN=%2d al=%2d p=%2d %15.12f\n",
1064 Mc_AN,Gc_AN,al,p,CntCoes[Mc_AN][al][p]);
1065 }
1066 }
1067 }
1068
1069 /****************************************************
1070 Restricted contraction
1071 ****************************************************/
1072
1073 if (RCnt_switch==1 || SICnt_switch==1){
1074
1075 al = -1;
1076 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
1077 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
1078
1079 for (p=0; p<Spe_Specified_Num[wan][al+1]; p++){
1080 Tmp_CntCoes[p] = 0.0;
1081 }
1082
1083 al0 = al;
1084 for (M0=0; M0<=2*L0; M0++){
1085 al++;
1086
1087 x = (Weight_ko[al] - LChemP)*Beta;
1088 if (x<=-30.0) x = -30.0;
1089 if (30.0<=x) x = 30.0;
1090 FermiF = 1.0/(1.0 + exp(x));
1091
1092 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1093 Tmp_CntCoes[p] = Tmp_CntCoes[p] + CntCoes[Mc_AN][al][p]*FermiF;
1094 }
1095 }
1096 al = al0;
1097
1098 for (M0=0; M0<=2*L0; M0++){
1099 al++;
1100 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1101 CntCoes[Mc_AN][al][p] = Tmp_CntCoes[p];
1102 }
1103 }
1104
1105 }
1106 }
1107
1108 /****************************************************
1109 Orthonormalization by the Gram-Schmidt method
1110 ****************************************************/
1111
1112 al = -1;
1113 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
1114 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
1115 for (M0=0; M0<=2*L0; M0++){
1116 al++;
1117 tmp_index2[L0][Mul0][M0] = al;
1118 }
1119 }
1120 }
1121
1122 for (L0=0; L0<=Spe_MaxL_Basis[wan]; L0++){
1123 for (M0=0; M0<=2*L0; M0++){
1124
1125 for (Mul0=0; Mul0<Spe_Num_CBasis[wan][L0]; Mul0++){
1126 al0 = tmp_index2[L0][Mul0][M0];
1127
1128 /* x - sum_i <x|e_i>e_i */
1129
1130 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1131 Tmp_CntCoes[p] = 0.0;
1132 }
1133
1134 for (Mul1=0; Mul1<Mul0; Mul1++){
1135 al1 = tmp_index2[L0][Mul1][M0];
1136
1137 sum = 0.0;
1138 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1139 sum = sum + CntCoes[Mc_AN][al0][p]*CntCoes[Mc_AN][al1][p];
1140 }
1141
1142 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1143 Tmp_CntCoes[p] = Tmp_CntCoes[p] + sum*CntCoes[Mc_AN][al1][p];
1144 }
1145 }
1146
1147 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1148 CntCoes[Mc_AN][al0][p] = CntCoes[Mc_AN][al0][p] - Tmp_CntCoes[p];
1149 }
1150
1151 /* Normalize */
1152 sum = 0.0;
1153 Max0 = -100.0;
1154 pmax = 0;
1155 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1156 sum = sum + CntCoes[Mc_AN][al0][p]*CntCoes[Mc_AN][al0][p];
1157 if (Max0<fabs(CntCoes[Mc_AN][al0][p])){
1158 Max0 = fabs(CntCoes[Mc_AN][al0][p]);
1159 pmax = p;
1160 }
1161 }
1162 tmp0 = 1.0/sqrt(sum);
1163 tmp1 = sgn(CntCoes[Mc_AN][al0][pmax]);
1164
1165 for (p=0; p<Spe_Specified_Num[wan][al0]; p++){
1166 CntCoes[Mc_AN][al0][p] = tmp0*tmp1*CntCoes[Mc_AN][al0][p];
1167 }
1168
1169 }
1170 }
1171 }
1172 }
1173
1174 if (3<=level_stdout){
1175 for (al=0; al<Spe_Total_CNO[wan]; al++){
1176 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1177 printf("C Init_CntCoes Mc_AN=%2d Gc_AN=%2d al=%2d p=%2d %15.12f\n",
1178 Mc_AN,Gc_AN,al,p,CntCoes[Mc_AN][al][p]);
1179 }
1180 }
1181 }
1182
1183 /****************************************************
1184 free arrays
1185 ****************************************************/
1186
1187 free(ko);
1188 free(abs_sum);
1189 free(M1);
1190 free(dege);
1191 free(jun);
1192 free(ponu);
1193 free(C0);
1194
1195 for (i=0; i<n2; i++){
1196 free(S[i]);
1197 free(Hks[i]);
1198 free(D[i]);
1199 free(C[i]);
1200 free(B[i]);
1201 }
1202 free(S);
1203 free(Hks);
1204 free(D);
1205 free(C);
1206 free(B);
1207
1208 } /* Mc_AN */
1209
1210 /****************************************************
1211 average contraction coefficients
1212 ****************************************************/
1213
1214 if (ACnt_switch==1){
1215
1216 /* allocation */
1217 My_CntCoes_Spe = (double***)malloc(sizeof(double**)*(SpeciesNum+1));
1218 for (i=0; i<=SpeciesNum; i++){
1219 My_CntCoes_Spe[i] = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
1220 for (j=0; j<List_YOUSO[7]; j++){
1221 My_CntCoes_Spe[i][j] = (double*)malloc(sizeof(double)*List_YOUSO[24]);
1222 }
1223 }
1224
1225 CntCoes_Spe = (double***)malloc(sizeof(double**)*(SpeciesNum+1));
1226 for (i=0; i<=SpeciesNum; i++){
1227 CntCoes_Spe[i] = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
1228 for (j=0; j<List_YOUSO[7]; j++){
1229 CntCoes_Spe[i][j] = (double*)malloc(sizeof(double)*List_YOUSO[24]);
1230 }
1231 }
1232
1233 /* initialize */
1234 for (wan=0; wan<SpeciesNum; wan++){
1235 for (i=0; i<List_YOUSO[7]; i++){
1236 for (j=0; j<List_YOUSO[24]; j++){
1237 My_CntCoes_Spe[wan][i][j] = 0.0;
1238 }
1239 }
1240 }
1241
1242 /* local sum in a proccessor */
1243 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1244 Gc_AN = M2G[Mc_AN];
1245 wan = WhatSpecies[Gc_AN];
1246 for (al=0; al<Spe_Total_CNO[wan]; al++){
1247 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1248 My_CntCoes_Spe[wan][al][p] += CntCoes[Mc_AN][al][p];
1249 }
1250 }
1251 }
1252
1253 /* global sum by MPI */
1254 for (wan=0; wan<SpeciesNum; wan++){
1255 for (al=0; al<Spe_Total_CNO[wan]; al++){
1256 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1257 MPI_Allreduce(&My_CntCoes_Spe[wan][al][p], &CntCoes_Spe[wan][al][p],
1258 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
1259 }
1260 }
1261 }
1262
1263 /* CntCoes_Spe to CntCoes */
1264 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1265 Gc_AN = M2G[Mc_AN];
1266 wan = WhatSpecies[Gc_AN];
1267 for (al=0; al<Spe_Total_CNO[wan]; al++){
1268 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1269 CntCoes[Mc_AN][al][p] = CntCoes_Spe[wan][al][p];
1270 }
1271 }
1272 }
1273
1274 /* free */
1275 for (i=0; i<=SpeciesNum; i++){
1276 for (j=0; j<List_YOUSO[7]; j++){
1277 free(My_CntCoes_Spe[i][j]);
1278 }
1279 free(My_CntCoes_Spe[i]);
1280 }
1281 free(My_CntCoes_Spe);
1282
1283 for (i=0; i<=SpeciesNum; i++){
1284 for (j=0; j<List_YOUSO[7]; j++){
1285 free(CntCoes_Spe[i][j]);
1286 }
1287 free(CntCoes_Spe[i]);
1288 }
1289 free(CntCoes_Spe);
1290
1291 }
1292
1293 /****************************************************
1294 Normalization
1295 ****************************************************/
1296
1297 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
1298 Gc_AN = M2G[Mc_AN];
1299 wan = WhatSpecies[Gc_AN];
1300 for (al=0; al<Spe_Total_CNO[wan]; al++){
1301
1302 sum = 0.0;
1303 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1304 p0 = Spe_Trans_Orbital[wan][al][p];
1305
1306 for (q=0; q<Spe_Specified_Num[wan][al]; q++){
1307 q0 = Spe_Trans_Orbital[wan][al][q];
1308
1309 tmp0 = CntCoes[Mc_AN][al][p]*CntCoes[Mc_AN][al][q];
1310 sum = sum + tmp0*OLP[0][Mc_AN][0][p0][q0];
1311 }
1312 }
1313
1314 tmp0 = 1.0/sqrt(sum);
1315 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1316 CntCoes[Mc_AN][al][p] = CntCoes[Mc_AN][al][p]*tmp0;
1317 }
1318
1319 }
1320
1321 if (3<=level_stdout){
1322 for (al=0; al<Spe_Total_CNO[wan]; al++){
1323 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1324 printf("D Init_CntCoes Mc_AN=%2d Gc_AN=%2d al=%2d p=%2d %15.12f\n",
1325 Mc_AN,Gc_AN,al,p,CntCoes[Mc_AN][al][p]);
1326 }
1327 }
1328 }
1329 }
1330
1331 /****************************************************
1332 MPI:
1333
1334 CntCoes
1335 ****************************************************/
1336
1337 /***********************************
1338 set data size
1339 ************************************/
1340
1341 for (ID=0; ID<numprocs; ID++){
1342
1343 IDS = (myid + ID) % numprocs;
1344 IDR = (myid - ID + numprocs) % numprocs;
1345
1346 if (ID!=0){
1347 tag = 999;
1348
1349 /* find data size to send block data */
1350 if (F_Snd_Num[IDS]!=0){
1351
1352 size1 = 0;
1353 for (n=0; n<F_Snd_Num[IDS]; n++){
1354 Mc_AN = Snd_MAN[IDS][n];
1355 Gc_AN = Snd_GAN[IDS][n];
1356 wan = WhatSpecies[Gc_AN];
1357 for (al=0; al<Spe_Total_CNO[wan]; al++){
1358 size1 += Spe_Specified_Num[wan][al];
1359 }
1360 }
1361
1362 Snd_CntCoes_Size[IDS] = size1;
1363 MPI_Isend(&size1, 1, MPI_INT, IDS, tag, mpi_comm_level1, &request);
1364 }
1365 else{
1366 Snd_CntCoes_Size[IDS] = 0;
1367 }
1368
1369 /* receiving of size of data */
1370
1371 if (F_Rcv_Num[IDR]!=0){
1372 MPI_Recv(&size2, 1, MPI_INT, IDR, tag, mpi_comm_level1, &stat);
1373 Rcv_CntCoes_Size[IDR] = size2;
1374 }
1375 else{
1376 Rcv_CntCoes_Size[IDR] = 0;
1377 }
1378
1379 if (F_Snd_Num[IDS]!=0) MPI_Wait(&request,&stat);
1380
1381 }
1382 else {
1383 Snd_CntCoes_Size[myid] = 0;
1384 Rcv_CntCoes_Size[myid] = 0;
1385 }
1386 }
1387
1388 /***********************************
1389 data transfer
1390 ************************************/
1391
1392 tag = 999;
1393 for (ID=0; ID<numprocs; ID++){
1394
1395 IDS = (myid + ID) % numprocs;
1396 IDR = (myid - ID + numprocs) % numprocs;
1397
1398 if (ID!=0){
1399
1400 /*****************************
1401 sending of data
1402 *****************************/
1403
1404 if (F_Snd_Num[IDS]!=0){
1405
1406 size1 = Snd_CntCoes_Size[IDS];
1407
1408 /* allocation of array */
1409 tmp_array = (double*)malloc(sizeof(double)*size1);
1410
1411 /* multidimentional array to vector array */
1412
1413 num = 0;
1414 for (n=0; n<F_Snd_Num[IDS]; n++){
1415 Mc_AN = Snd_MAN[IDS][n];
1416 Gc_AN = Snd_GAN[IDS][n];
1417 wan = WhatSpecies[Gc_AN];
1418 for (al=0; al<Spe_Total_CNO[wan]; al++){
1419 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1420 tmp_array[num] = CntCoes[Mc_AN][al][p];
1421 num++;
1422 }
1423 }
1424 }
1425
1426 MPI_Isend(&tmp_array[0], size1, MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request);
1427 }
1428
1429 /*****************************
1430 receiving of block data
1431 *****************************/
1432
1433 if (F_Rcv_Num[IDR]!=0){
1434
1435 size2 = Rcv_CntCoes_Size[IDR];
1436
1437 /* allocation of array */
1438 tmp_array2 = (double*)malloc(sizeof(double)*size2);
1439
1440 MPI_Recv(&tmp_array2[0], size2, MPI_DOUBLE, IDR, tag, mpi_comm_level1, &stat);
1441
1442 num = 0;
1443 Mc_AN = F_TopMAN[IDR] - 1;
1444 for (n=0; n<F_Rcv_Num[IDR]; n++){
1445 Mc_AN++;
1446 Gc_AN = Rcv_GAN[IDR][n];
1447 wan = WhatSpecies[Gc_AN];
1448 for (al=0; al<Spe_Total_CNO[wan]; al++){
1449 for (p=0; p<Spe_Specified_Num[wan][al]; p++){
1450 CntCoes[Mc_AN][al][p] = tmp_array2[num];
1451 num++;
1452 }
1453 }
1454 }
1455
1456 /* freeing of array */
1457 free(tmp_array2);
1458 }
1459
1460 if (F_Snd_Num[IDS]!=0){
1461 MPI_Wait(&request,&stat);
1462 free(tmp_array); /* freeing of array */
1463 }
1464 }
1465 }
1466
1467 /****************************************************
1468 freeing of arrays:
1469
1470 int MP[List_YOUSO[8]];
1471
1472 int tmp_index[List_YOUSO[25]+1]
1473 [2*(List_YOUSO[25]+1)+1];
1474 int tmp_index2[List_YOUSO[25]+1]
1475 [List_YOUSO[24]]
1476 [2*(List_YOUSO[25]+1)+1];
1477
1478 double Tmp_CntCoes[List_YOUSO[24]]
1479
1480 double Check_ko[List_YOUSO[25]+1]
1481 [2*(List_YOUSO[25]+1)+1];
1482
1483 double Weight_ko[List_YOUSO[7]];
1484 ****************************************************/
1485
1486 free(MP);
1487
1488 for (i=0; i<(List_YOUSO[25]+1); i++){
1489 free(tmp_index[i]);
1490 }
1491 free(tmp_index);
1492
1493 for (i=0; i<(List_YOUSO[25]+1); i++){
1494 for (j=0; j<List_YOUSO[24]; j++){
1495 free(tmp_index2[i][j]);
1496 }
1497 free(tmp_index2[i]);
1498 }
1499 free(tmp_index2);
1500
1501 free(Tmp_CntCoes);
1502
1503 for (i=0; i<(List_YOUSO[25]+1); i++){
1504 free(Check_ko[i]);
1505 }
1506 free(Check_ko);
1507
1508 free(Weight_ko);
1509
1510 free(Snd_CntCoes_Size);
1511 free(Rcv_CntCoes_Size);
1512
1513 free(Snd_H_Size);
1514 free(Rcv_H_Size);
1515
1516 free(Snd_S_Size);
1517 free(Rcv_S_Size);
1518
1519 /* for elapsed time */
1520 dtime(&TEtime);
1521 time0 = TEtime - TStime;
1522 }
1523
1524