1 /**********************************************************************
2 Set_Aden_Grid.c:
3
4 Set_Aden_Grid.c is a subroutine to calculate a charge density
5 superposed atomic densities on grid.
6
7 Log of Set_Aden_Grid.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 #include <omp.h>
20
21
Set_Aden_Grid()22 double Set_Aden_Grid()
23 {
24 /****************************************************
25 Densities by the atomic superposition
26 densities on grids
27 ****************************************************/
28
29 static int firsttime=1;
30 int i,j,k,Gc_AN,Mc_AN,NN_S,NN_R;
31 int Cwan,Nc,GNc,GRc,mul,gp,spe;
32 int N3[4],AN,BN,CN,DN,LN;
33 unsigned long long int GN,n2D,N2D;
34 int size_AtomDen_Grid;
35 int size_AtomDen_Snd_Grid_A2B;
36 int size_AtomDen_Rcv_Grid_A2B;
37 double DenA,x,dx,dy,dz,r;
38 double Nele,Nu,Nd,M,ocupcy_u,ocupcy_d;
39 double rho,mag,magx,magy,magz,theta,phi;
40 double time0,TStime,TEtime;
41 double Stime_atom, Etime_atom;
42 double Cxyz[4],Nup,Ndown;
43 double **AtomDen_Snd_Grid_A2B;
44 double **AtomDen_Rcv_Grid_A2B;
45 double **PCCDen_Snd_Grid_A2B;
46 double **PCCDen_Rcv_Grid_A2B;
47 double **AtomDen_Grid;
48 double **PCCDen_Grid;
49 double *Work_Array_Snd_Grid_B2C;
50 double *Work_Array_Rcv_Grid_B2C;
51 double *Work_Array_Snd_Grid_B2D;
52 double *Work_Array_Rcv_Grid_B2D;
53 int numprocs,myid,tag=999,ID,IDS,IDR;
54
55 MPI_Status stat;
56 MPI_Request request;
57 MPI_Status *stat_send;
58 MPI_Status *stat_recv;
59 MPI_Request *request_send;
60 MPI_Request *request_recv;
61
62 /* for OpenMP */
63 int OMPID,Nthrds,Nprocs;
64
65 /* MPI */
66 MPI_Comm_size(mpi_comm_level1,&numprocs);
67 MPI_Comm_rank(mpi_comm_level1,&myid);
68
69 dtime(&TStime);
70
71 /****************************************************
72 allocation of arrays:
73 ****************************************************/
74
75 size_AtomDen_Snd_Grid_A2B = 0;
76 AtomDen_Snd_Grid_A2B = (double**)malloc(sizeof(double*)*numprocs);
77 for (ID=0; ID<numprocs; ID++){
78 AtomDen_Snd_Grid_A2B[ID] = (double*)malloc(sizeof(double)*Num_Snd_Grid_A2B[ID]);
79 size_AtomDen_Snd_Grid_A2B += Num_Snd_Grid_A2B[ID];
80 }
81
82 size_AtomDen_Rcv_Grid_A2B = 0;
83 AtomDen_Rcv_Grid_A2B = (double**)malloc(sizeof(double*)*numprocs);
84 for (ID=0; ID<numprocs; ID++){
85 AtomDen_Rcv_Grid_A2B[ID] = (double*)malloc(sizeof(double)*Num_Rcv_Grid_A2B[ID]);
86 size_AtomDen_Rcv_Grid_A2B += Num_Rcv_Grid_A2B[ID];
87 }
88
89 if (PCC_switch==1){
90
91 PCCDen_Snd_Grid_A2B = (double**)malloc(sizeof(double*)*numprocs);
92 for (ID=0; ID<numprocs; ID++){
93 PCCDen_Snd_Grid_A2B[ID] = (double*)malloc(sizeof(double)*Num_Snd_Grid_A2B[ID]);
94 }
95
96 PCCDen_Rcv_Grid_A2B = (double**)malloc(sizeof(double*)*numprocs);
97 for (ID=0; ID<numprocs; ID++){
98 PCCDen_Rcv_Grid_A2B[ID] = (double*)malloc(sizeof(double)*Num_Rcv_Grid_A2B[ID]);
99 }
100 }
101
102 size_AtomDen_Grid = 0;
103 AtomDen_Grid = (double**)malloc(sizeof(double*)*(Matomnum+1));
104 AtomDen_Grid[0] = (double*)malloc(sizeof(double)*1);
105 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
106 Gc_AN = F_M2G[Mc_AN];
107 AtomDen_Grid[Mc_AN] = (double*)malloc(sizeof(double)*GridN_Atom[Gc_AN]);
108 size_AtomDen_Grid += GridN_Atom[Gc_AN];
109 }
110
111 if (PCC_switch==1){
112 PCCDen_Grid = (double**)malloc(sizeof(double*)*(Matomnum+1));
113 PCCDen_Grid[0] = (double*)malloc(sizeof(double)*1);
114 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
115 Gc_AN = F_M2G[Mc_AN];
116 PCCDen_Grid[Mc_AN] = (double*)malloc(sizeof(double)*GridN_Atom[Gc_AN]);
117 }
118 }
119
120 if (SpinP_switch<=1) mul = 2;
121 else if (SpinP_switch==3) mul = 4;
122
123 Work_Array_Snd_Grid_B2C = (double*)malloc(sizeof(double)*GP_B2C_S[NN_B2C_S]*mul);
124 Work_Array_Rcv_Grid_B2C = (double*)malloc(sizeof(double)*GP_B2C_R[NN_B2C_R]*mul);
125 Work_Array_Snd_Grid_B2D = (double*)malloc(sizeof(double)*GP_B2D_S[NN_B2D_S]*(mul+2));
126 Work_Array_Rcv_Grid_B2D = (double*)malloc(sizeof(double)*GP_B2D_R[NN_B2D_R]*(mul+2));
127
128 /* PrintMemory */
129
130 if (firsttime==1){
131 PrintMemory("Set_Aden_Grid: AtomDen_Grid", sizeof(double)*size_AtomDen_Grid, NULL);
132 if (PCC_switch==1){
133 PrintMemory("Set_Aden_Grid: PCCDen_Grid", sizeof(double)*size_AtomDen_Grid, NULL);
134 }
135 PrintMemory("Set_Aden_Grid: AtomDen_Snd_Grid_A2B", sizeof(double)*size_AtomDen_Snd_Grid_A2B, NULL);
136 PrintMemory("Set_Aden_Grid: AtomDen_Rcv_Grid_A2B", sizeof(double)*size_AtomDen_Rcv_Grid_A2B, NULL);
137 if (PCC_switch==1){
138 PrintMemory("Set_Aden_Grid: PCCDen_Snd_Grid_A2B", sizeof(double)*size_AtomDen_Snd_Grid_A2B, NULL);
139 PrintMemory("Set_Aden_Grid: PCCDen_Rcv_Grid_A2B", sizeof(double)*size_AtomDen_Rcv_Grid_A2B, NULL);
140 }
141
142 PrintMemory("Set_Aden_Grid: Work_Array_Snd_Grid_B2C", sizeof(double)*GP_B2C_S[NN_B2C_S]*mul, NULL);
143 PrintMemory("Set_Aden_Grid: Work_Array_Rcv_Grid_B2C", sizeof(double)*GP_B2C_R[NN_B2C_R]*mul, NULL);
144 PrintMemory("Set_Aden_Grid: Work_Array_Snd_Grid_B2D", sizeof(double)*GP_B2D_S[NN_B2D_S]*(mul+2), NULL);
145 PrintMemory("Set_Aden_Grid: Work_Array_Rcv_Grid_B2D", sizeof(double)*GP_B2D_R[NN_B2D_R]*(mul+2), NULL);
146 firsttime = 0;
147 }
148
149 /******************************************************
150 calculation of AtomDen_Grid and PCCDen_Grid
151 ******************************************************/
152
153 /* calculate AtomDen_Grid and PCCDen_Grid */
154
155 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
156
157 dtime(&Stime_atom);
158
159 Gc_AN = M2G[Mc_AN];
160 Cwan = WhatSpecies[Gc_AN];
161
162 #pragma omp parallel shared(Spe_Atomic_PCC,Spe_VPS_RV,Spe_VPS_XV,Spe_Num_Mesh_VPS,Spe_PAO_RV,Spe_Atomic_Den,Spe_PAO_XV,Spe_Num_Mesh_PAO,PCCDen_Grid,PCC_switch,AtomDen_Grid,Cwan,Gxyz,atv,Mc_AN,CellListAtom,GridListAtom,GridN_Atom,Gc_AN) private(OMPID,Nthrds,Nc,GNc,GRc,Cxyz,dx,dy,dz,r,x)
163 {
164
165 /* get info. on OpenMP */
166
167 OMPID = omp_get_thread_num();
168 Nthrds = omp_get_num_threads();
169
170 for (Nc=OMPID; Nc<GridN_Atom[Gc_AN]; Nc+=Nthrds){
171
172 GNc = GridListAtom[Mc_AN][Nc];
173 GRc = CellListAtom[Mc_AN][Nc];
174
175 Get_Grid_XYZ(GNc,Cxyz);
176 dx = Cxyz[1] + atv[GRc][1] - Gxyz[Gc_AN][1];
177 dy = Cxyz[2] + atv[GRc][2] - Gxyz[Gc_AN][2];
178 dz = Cxyz[3] + atv[GRc][3] - Gxyz[Gc_AN][3];
179
180 x = 0.5*log(dx*dx + dy*dy + dz*dz);
181
182 /* atomic density */
183
184 AtomDen_Grid[Mc_AN][Nc] = KumoF( Spe_Num_Mesh_PAO[Cwan], x,
185 Spe_PAO_XV[Cwan], Spe_PAO_RV[Cwan], Spe_Atomic_Den[Cwan]);
186
187 /* partial core correction */
188
189 if (PCC_switch==1) {
190 PCCDen_Grid[Mc_AN][Nc] = KumoF( Spe_Num_Mesh_VPS[Cwan], x,
191 Spe_VPS_XV[Cwan], Spe_VPS_RV[Cwan], Spe_Atomic_PCC[Cwan]);
192 }
193 }
194
195 } /* #pragma omp parallel */
196
197 dtime(&Etime_atom);
198 time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
199 }
200
201 /******************************************************
202 MPI communication from the partitions A to B
203 ******************************************************/
204
205 /* copy AtomDen_Grid to AtomDen_Snd_Grid_A2B */
206 /* copy PCCDen_Grid to PCCDen_Snd_Grid_A2B */
207
208 for (ID=0; ID<numprocs; ID++) Num_Snd_Grid_A2B[ID] = 0;
209
210 N2D = Ngrid1*Ngrid2;
211
212 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
213
214 Gc_AN = M2G[Mc_AN];
215
216 for (AN=0; AN<GridN_Atom[Gc_AN]; AN++){
217
218 GN = GridListAtom[Mc_AN][AN];
219 GN2N(GN,N3);
220 n2D = N3[1]*Ngrid2 + N3[2];
221 ID = (int)(n2D*(unsigned long long int)numprocs/N2D);
222 AtomDen_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]] = AtomDen_Grid[Mc_AN][AN];
223
224 if (PCC_switch==1){
225 PCCDen_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]] = PCCDen_Grid[Mc_AN][AN];
226 }
227
228 Num_Snd_Grid_A2B[ID]++;
229 }
230 }
231
232 /* MPI: A to B */
233
234 request_send = malloc(sizeof(MPI_Request)*NN_A2B_S);
235 request_recv = malloc(sizeof(MPI_Request)*NN_A2B_R);
236 stat_send = malloc(sizeof(MPI_Status)*NN_A2B_S);
237 stat_recv = malloc(sizeof(MPI_Status)*NN_A2B_R);
238
239 NN_S = 0;
240 NN_R = 0;
241
242 tag = 999;
243 for (ID=1; ID<numprocs; ID++){
244
245 IDS = (myid + ID) % numprocs;
246 IDR = (myid - ID + numprocs) % numprocs;
247
248 if (Num_Snd_Grid_A2B[IDS]!=0){
249 MPI_Isend( &AtomDen_Snd_Grid_A2B[IDS][0], Num_Snd_Grid_A2B[IDS],
250 MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request_send[NN_S]);
251 NN_S++;
252 }
253
254 if (Num_Rcv_Grid_A2B[IDR]!=0){
255 MPI_Irecv( &AtomDen_Rcv_Grid_A2B[IDR][0], Num_Rcv_Grid_A2B[IDR],
256 MPI_DOUBLE, IDR, tag, mpi_comm_level1, &request_recv[NN_R]);
257 NN_R++;
258 }
259 }
260
261 if (NN_S!=0) MPI_Waitall(NN_S,request_send,stat_send);
262 if (NN_R!=0) MPI_Waitall(NN_R,request_recv,stat_recv);
263
264 /* for myid */
265 for (i=0; i<Num_Rcv_Grid_A2B[myid]; i++){
266 AtomDen_Rcv_Grid_A2B[myid][i] = AtomDen_Snd_Grid_A2B[myid][i];
267 }
268
269 /* PCC */
270
271 if (PCC_switch==1){
272
273 NN_S = 0;
274 NN_R = 0;
275
276 tag = 999;
277 for (ID=1; ID<numprocs; ID++){
278
279 IDS = (myid + ID) % numprocs;
280 IDR = (myid - ID + numprocs) % numprocs;
281
282 if (Num_Snd_Grid_A2B[IDS]!=0){
283 MPI_Isend( &PCCDen_Snd_Grid_A2B[IDS][0], Num_Snd_Grid_A2B[IDS],
284 MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request_send[NN_S]);
285 NN_S++;
286 }
287
288 if (Num_Rcv_Grid_A2B[IDR]!=0){
289 MPI_Irecv( &PCCDen_Rcv_Grid_A2B[IDR][0], Num_Rcv_Grid_A2B[IDR],
290 MPI_DOUBLE, IDR, tag, mpi_comm_level1, &request_recv[NN_R]);
291 NN_R++;
292 }
293 }
294
295 if (NN_S!=0) MPI_Waitall(NN_S,request_send,stat_send);
296 if (NN_R!=0) MPI_Waitall(NN_R,request_recv,stat_recv);
297
298 /* for myid */
299 for (i=0; i<Num_Rcv_Grid_A2B[myid]; i++){
300 PCCDen_Rcv_Grid_A2B[myid][i] = PCCDen_Snd_Grid_A2B[myid][i];
301 }
302 }
303
304 free(request_send);
305 free(request_recv);
306 free(stat_send);
307 free(stat_recv);
308
309 /******************************************************
310 superposition of atomic densities and PCC densities
311 in the partition B.
312
313 for spin non-collinear
314 1. set rho, mx, my, mz
315 2. later calculate theta and phi
316 3. n_up = (rho+m)/2
317 n_dn = (rho-m)/2
318 ******************************************************/
319
320 /* initialize arrays */
321
322 if (SpinP_switch==3){
323 for (BN=0; BN<My_NumGridB_AB; BN++){
324 Density_Grid_B[0][BN] = 0.0;
325 Density_Grid_B[1][BN] = 0.0;
326 Density_Grid_B[2][BN] = 0.0;
327 Density_Grid_B[3][BN] = 0.0;
328 }
329 }
330 else{
331 for (BN=0; BN<My_NumGridB_AB; BN++){
332 Density_Grid_B[0][BN] = 0.0;
333 Density_Grid_B[1][BN] = 0.0;
334 }
335 }
336
337 if (PCC_switch==1){
338 for (BN=0; BN<My_NumGridB_AB; BN++){
339 PCCDensity_Grid_B[0][BN] = 0.0;
340 PCCDensity_Grid_B[1][BN] = 0.0;
341 }
342 }
343
344 /* superposition of densities */
345
346 for (ID=0; ID<numprocs; ID++){
347
348 for (LN=0; LN<Num_Rcv_Grid_A2B[ID]; LN++){
349
350 BN = Index_Rcv_Grid_A2B[ID][3*LN+0];
351 Gc_AN = Index_Rcv_Grid_A2B[ID][3*LN+1];
352 GRc = Index_Rcv_Grid_A2B[ID][3*LN+2];
353 spe = WhatSpecies[Gc_AN];
354 Nele = InitN_USpin[Gc_AN] + InitN_DSpin[Gc_AN];
355 Nu = InitN_USpin[Gc_AN];
356 Nd = InitN_DSpin[Gc_AN];
357
358 DenA = AtomDen_Rcv_Grid_A2B[ID][LN];
359
360 if (1.0e-15<Nele){
361 ocupcy_u = Nu/Nele;
362 ocupcy_d = Nd/Nele;
363 }
364 else{
365 ocupcy_u = 0.0;
366 ocupcy_d = 0.0;
367 }
368
369 if (Solver!=4 || (Solver==4 && atv_ijk[GRc][1]==0 )){
370
371 /* spin collinear non-polarization */
372 if ( SpinP_switch==0 ){
373 Density_Grid_B[0][BN] += 0.5*DenA;
374 Density_Grid_B[1][BN] += 0.5*DenA;
375 }
376
377 /* spin collinear polarization */
378 else if ( SpinP_switch==1 ){
379 Density_Grid_B[0][BN] += ocupcy_u*DenA;
380 Density_Grid_B[1][BN] += ocupcy_d*DenA;
381 }
382
383 /* spin non-collinear */
384 else if ( SpinP_switch==3 ){
385
386 theta = Angle0_Spin[Gc_AN];
387 phi = Angle1_Spin[Gc_AN];
388
389 rho = DenA;
390 mag = (ocupcy_u - ocupcy_d)*DenA;
391 magx = mag*sin(theta)*cos(phi);
392 magy = mag*sin(theta)*sin(phi);
393 magz = mag*cos(theta);
394
395 Density_Grid_B[0][BN] += rho;
396 Density_Grid_B[1][BN] += magx;
397 Density_Grid_B[2][BN] += magy;
398 Density_Grid_B[3][BN] += magz;
399 }
400
401 /* partial core correction */
402 /* later add this in Set_XC_Grid */
403
404 if (PCC_switch==1){
405 if ( SpinP_switch==0 ){
406 PCCDensity_Grid_B[0][BN] += 0.5*PCCDen_Rcv_Grid_A2B[ID][LN];
407 PCCDensity_Grid_B[1][BN] += 0.5*PCCDen_Rcv_Grid_A2B[ID][LN];
408 }
409 else {
410
411 if (Spe_OpenCore_flag[spe]==0){
412 PCCDensity_Grid_B[0][BN] += 0.5*PCCDen_Rcv_Grid_A2B[ID][LN];
413 PCCDensity_Grid_B[1][BN] += 0.5*PCCDen_Rcv_Grid_A2B[ID][LN];
414 }
415 else if (Spe_OpenCore_flag[spe]==1){
416 PCCDensity_Grid_B[0][BN] += PCCDen_Rcv_Grid_A2B[ID][LN];
417 }
418 else if (Spe_OpenCore_flag[spe]==-1){
419 PCCDensity_Grid_B[1][BN] += PCCDen_Rcv_Grid_A2B[ID][LN];
420 }
421 }
422 }
423
424 } /* if (Solve!=4.....) */
425
426 } /* AN */
427 } /* ID */
428
429 /****************************************************
430 initialize diagonal and off-diagonal densities
431 in case of spin non-collinear DFT
432 ****************************************************/
433
434 if (SpinP_switch==3){
435 for (BN=0; BN<My_NumGridB_AB; BN++){
436
437 rho = Density_Grid_B[0][BN];
438 magx = Density_Grid_B[1][BN];
439 magy = Density_Grid_B[2][BN];
440 magz = Density_Grid_B[3][BN];
441
442 Density_Grid_B[0][BN] = 0.5*(rho + magz);
443 Density_Grid_B[1][BN] = 0.5*(rho - magz);
444 Density_Grid_B[2][BN] = 0.5*magx;
445 Density_Grid_B[3][BN] =-0.5*magy;
446 }
447 }
448
449 /******************************************************
450 Density_Grid_B to ADensity_Grid_B
451 ******************************************************/
452
453 if ( SpinP_switch==0 ){
454 for (BN=0; BN<My_NumGridB_AB; BN++){
455 ADensity_Grid_B[BN] = Density_Grid_B[0][BN];
456 }
457 }
458 else if ( SpinP_switch==1 || SpinP_switch==3 ){
459 for (BN=0; BN<My_NumGridB_AB; BN++){
460 ADensity_Grid_B[BN] = 0.5*(Density_Grid_B[0][BN] + Density_Grid_B[1][BN]);
461 }
462 }
463
464 /******************************************************
465 MPI: from the partitions B to C
466 ******************************************************/
467
468 request_send = malloc(sizeof(MPI_Request)*NN_B2C_S);
469 request_recv = malloc(sizeof(MPI_Request)*NN_B2C_R);
470 stat_send = malloc(sizeof(MPI_Status)*NN_B2C_S);
471 stat_recv = malloc(sizeof(MPI_Status)*NN_B2C_R);
472
473 NN_S = 0;
474 NN_R = 0;
475
476 if (SpinP_switch<=1) mul = 2;
477 else if (SpinP_switch==3) mul = 4;
478
479 /* MPI_Irecv */
480
481 for (ID=0; ID<NN_B2C_R; ID++){
482
483 IDR = ID_NN_B2C_R[ID];
484 gp = GP_B2C_R[ID];
485
486 if (IDR!=myid){
487 MPI_Irecv( &Work_Array_Rcv_Grid_B2C[mul*gp], Num_Rcv_Grid_B2C[IDR]*mul,
488 MPI_DOUBLE, IDR, tag, mpi_comm_level1, &request_recv[NN_R]);
489 NN_R++;
490 }
491 }
492
493 /* MPI_Isend */
494
495 for (ID=0; ID<NN_B2C_S; ID++){
496
497 IDS = ID_NN_B2C_S[ID];
498 gp = GP_B2C_S[ID];
499
500 /* copy Density_Grid_B to Work_Array_Snd_Grid_B2C */
501
502 for (LN=0; LN<Num_Snd_Grid_B2C[IDS]; LN++){
503
504 BN = Index_Snd_Grid_B2C[IDS][LN];
505
506 if (SpinP_switch<=1){
507 Work_Array_Snd_Grid_B2C[2*gp+2*LN+0] = Density_Grid_B[0][BN];
508 Work_Array_Snd_Grid_B2C[2*gp+2*LN+1] = Density_Grid_B[1][BN];
509 }
510 else if (SpinP_switch==3){
511 Work_Array_Snd_Grid_B2C[4*gp+4*LN+0] = Density_Grid_B[0][BN];
512 Work_Array_Snd_Grid_B2C[4*gp+4*LN+1] = Density_Grid_B[1][BN];
513 Work_Array_Snd_Grid_B2C[4*gp+4*LN+2] = Density_Grid_B[2][BN];
514 Work_Array_Snd_Grid_B2C[4*gp+4*LN+3] = Density_Grid_B[3][BN];
515 }
516 } /* LN */
517
518 if (IDS!=myid){
519 MPI_Isend( &Work_Array_Snd_Grid_B2C[mul*gp], Num_Snd_Grid_B2C[IDS]*mul,
520 MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request_send[NN_S]);
521 NN_S++;
522 }
523 }
524
525 if (NN_S!=0) MPI_Waitall(NN_S,request_send,stat_send);
526 if (NN_R!=0) MPI_Waitall(NN_R,request_recv,stat_recv);
527
528 free(request_send);
529 free(request_recv);
530 free(stat_send);
531 free(stat_recv);
532
533 /* copy Work_Array_Snd_Grid_B2C to Density_Grid */
534
535 for (ID=0; ID<NN_B2C_R; ID++){
536
537 IDR = ID_NN_B2C_R[ID];
538
539 if (IDR==myid){
540
541 gp = GP_B2C_S[ID];
542
543 for (LN=0; LN<Num_Rcv_Grid_B2C[myid]; LN++){
544
545 CN = Index_Rcv_Grid_B2C[myid][LN];
546
547 if (SpinP_switch<=1){
548 Density_Grid[0][CN] = Work_Array_Snd_Grid_B2C[2*gp+2*LN+0];
549 Density_Grid[1][CN] = Work_Array_Snd_Grid_B2C[2*gp+2*LN+1];
550 }
551 else if (SpinP_switch==3){
552 Density_Grid[0][CN] = Work_Array_Snd_Grid_B2C[4*gp+4*LN+0];
553 Density_Grid[1][CN] = Work_Array_Snd_Grid_B2C[4*gp+4*LN+1];
554 Density_Grid[2][CN] = Work_Array_Snd_Grid_B2C[4*gp+4*LN+2];
555 Density_Grid[3][CN] = Work_Array_Snd_Grid_B2C[4*gp+4*LN+3];
556 }
557 } /* LN */
558 }
559
560 else{
561
562 gp = GP_B2C_R[ID];
563
564 for (LN=0; LN<Num_Rcv_Grid_B2C[IDR]; LN++){
565
566 CN = Index_Rcv_Grid_B2C[IDR][LN];
567
568 if (SpinP_switch<=1){
569 Density_Grid[0][CN] = Work_Array_Rcv_Grid_B2C[2*gp+2*LN+0];
570 Density_Grid[1][CN] = Work_Array_Rcv_Grid_B2C[2*gp+2*LN+1];
571 }
572 else if (SpinP_switch==3){
573 Density_Grid[0][CN] = Work_Array_Rcv_Grid_B2C[4*gp+4*LN+0];
574 Density_Grid[1][CN] = Work_Array_Rcv_Grid_B2C[4*gp+4*LN+1];
575 Density_Grid[2][CN] = Work_Array_Rcv_Grid_B2C[4*gp+4*LN+2];
576 Density_Grid[3][CN] = Work_Array_Rcv_Grid_B2C[4*gp+4*LN+3];
577 }
578 }
579 }
580 }
581
582 /******************************************************
583 MPI: from the partitions B to D
584 ******************************************************/
585
586 request_send = malloc(sizeof(MPI_Request)*NN_B2D_S);
587 request_recv = malloc(sizeof(MPI_Request)*NN_B2D_R);
588 stat_send = malloc(sizeof(MPI_Status)*NN_B2D_S);
589 stat_recv = malloc(sizeof(MPI_Status)*NN_B2D_R);
590
591 NN_S = 0;
592 NN_R = 0;
593
594 if (SpinP_switch<=1) mul = 4;
595 else if (SpinP_switch==3) mul = 6;
596
597 /* MPI_Irecv */
598
599 for (ID=0; ID<NN_B2D_R; ID++){
600
601 IDR = ID_NN_B2D_R[ID];
602 gp = GP_B2D_R[ID];
603
604 if (IDR!=myid){
605 MPI_Irecv( &Work_Array_Rcv_Grid_B2D[mul*gp], Num_Rcv_Grid_B2D[IDR]*mul,
606 MPI_DOUBLE, IDR, tag, mpi_comm_level1, &request_recv[NN_R]);
607 NN_R++;
608 }
609 }
610
611 /* MPI_Isend */
612
613 for (ID=0; ID<NN_B2D_S; ID++){
614
615 IDS = ID_NN_B2D_S[ID];
616 gp = GP_B2D_S[ID];
617
618 /* copy Density_Grid_B to Work_Array_Snd_Grid_B2D */
619
620 for (LN=0; LN<Num_Snd_Grid_B2D[IDS]; LN++){
621
622 BN = Index_Snd_Grid_B2D[IDS][LN];
623
624 if (SpinP_switch<=1){
625 Work_Array_Snd_Grid_B2D[4*gp+4*LN+0] = Density_Grid_B[0][BN];
626 Work_Array_Snd_Grid_B2D[4*gp+4*LN+1] = Density_Grid_B[1][BN];
627 Work_Array_Snd_Grid_B2D[4*gp+4*LN+2] = PCCDensity_Grid_B[0][BN];
628 Work_Array_Snd_Grid_B2D[4*gp+4*LN+3] = PCCDensity_Grid_B[1][BN];
629 }
630 else if (SpinP_switch==3){
631 Work_Array_Snd_Grid_B2D[6*gp+6*LN+0] = Density_Grid_B[0][BN];
632 Work_Array_Snd_Grid_B2D[6*gp+6*LN+1] = Density_Grid_B[1][BN];
633 Work_Array_Snd_Grid_B2D[6*gp+6*LN+2] = Density_Grid_B[2][BN];
634 Work_Array_Snd_Grid_B2D[6*gp+6*LN+3] = Density_Grid_B[3][BN];
635 Work_Array_Snd_Grid_B2D[6*gp+6*LN+4] = PCCDensity_Grid_B[0][BN];
636 Work_Array_Snd_Grid_B2D[6*gp+6*LN+5] = PCCDensity_Grid_B[1][BN];
637 }
638 } /* LN */
639
640 if (IDS!=myid){
641 MPI_Isend( &Work_Array_Snd_Grid_B2D[mul*gp], Num_Snd_Grid_B2D[IDS]*mul,
642 MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request_send[NN_S]);
643 NN_S++;
644 }
645 }
646
647 if (NN_S!=0) MPI_Waitall(NN_S,request_send,stat_send);
648 if (NN_R!=0) MPI_Waitall(NN_R,request_recv,stat_recv);
649
650 free(request_send);
651 free(request_recv);
652 free(stat_send);
653 free(stat_recv);
654
655 /* copy Work_Array_Snd_Grid_B2D to Density_Grid */
656
657 for (ID=0; ID<NN_B2D_R; ID++){
658
659 IDR = ID_NN_B2D_R[ID];
660
661 if (IDR==myid){
662
663 gp = GP_B2D_S[ID];
664
665 for (LN=0; LN<Num_Rcv_Grid_B2D[myid]; LN++){
666
667 DN = Index_Rcv_Grid_B2D[myid][LN];
668
669 if (SpinP_switch<=1){
670 Density_Grid_D[0][DN] = Work_Array_Snd_Grid_B2D[4*gp+4*LN+0];
671 Density_Grid_D[1][DN] = Work_Array_Snd_Grid_B2D[4*gp+4*LN+1];
672 PCCDensity_Grid_D[0][DN] = Work_Array_Snd_Grid_B2D[4*gp+4*LN+2];
673 PCCDensity_Grid_D[1][DN] = Work_Array_Snd_Grid_B2D[4*gp+4*LN+3];
674 }
675 else if (SpinP_switch==3){
676 Density_Grid_D[0][DN] = Work_Array_Snd_Grid_B2D[6*gp+6*LN+0];
677 Density_Grid_D[1][DN] = Work_Array_Snd_Grid_B2D[6*gp+6*LN+1];
678 Density_Grid_D[2][DN] = Work_Array_Snd_Grid_B2D[6*gp+6*LN+2];
679 Density_Grid_D[3][DN] = Work_Array_Snd_Grid_B2D[6*gp+6*LN+3];
680 PCCDensity_Grid_D[0][DN] = Work_Array_Snd_Grid_B2D[6*gp+6*LN+4];
681 PCCDensity_Grid_D[1][DN] = Work_Array_Snd_Grid_B2D[6*gp+6*LN+5];
682 }
683 }
684 }
685
686 else {
687
688 gp = GP_B2D_R[ID];
689
690 for (LN=0; LN<Num_Rcv_Grid_B2D[IDR]; LN++){
691
692 DN = Index_Rcv_Grid_B2D[IDR][LN];
693
694 if (SpinP_switch<=1){
695 Density_Grid_D[0][DN] = Work_Array_Rcv_Grid_B2D[4*gp+4*LN+0];
696 Density_Grid_D[1][DN] = Work_Array_Rcv_Grid_B2D[4*gp+4*LN+1];
697 PCCDensity_Grid_D[0][DN] = Work_Array_Rcv_Grid_B2D[4*gp+4*LN+2];
698 PCCDensity_Grid_D[1][DN] = Work_Array_Rcv_Grid_B2D[4*gp+4*LN+3];
699 }
700 else if (SpinP_switch==3){
701 Density_Grid_D[0][DN] = Work_Array_Rcv_Grid_B2D[6*gp+6*LN+0];
702 Density_Grid_D[1][DN] = Work_Array_Rcv_Grid_B2D[6*gp+6*LN+1];
703 Density_Grid_D[2][DN] = Work_Array_Rcv_Grid_B2D[6*gp+6*LN+2];
704 Density_Grid_D[3][DN] = Work_Array_Rcv_Grid_B2D[6*gp+6*LN+3];
705 PCCDensity_Grid_D[0][DN] = Work_Array_Rcv_Grid_B2D[6*gp+6*LN+4];
706 PCCDensity_Grid_D[1][DN] = Work_Array_Rcv_Grid_B2D[6*gp+6*LN+5];
707 }
708 }
709 }
710 }
711
712 /****************************************************
713 freeing of arrays:
714 ****************************************************/
715
716 for (ID=0; ID<numprocs; ID++){
717 free(AtomDen_Snd_Grid_A2B[ID]);
718 }
719 free(AtomDen_Snd_Grid_A2B);
720
721 for (ID=0; ID<numprocs; ID++){
722 free(AtomDen_Rcv_Grid_A2B[ID]);
723 }
724 free(AtomDen_Rcv_Grid_A2B);
725
726 if (PCC_switch==1){
727
728 for (ID=0; ID<numprocs; ID++){
729 free(PCCDen_Snd_Grid_A2B[ID]);
730 }
731 free(PCCDen_Snd_Grid_A2B);
732
733 for (ID=0; ID<numprocs; ID++){
734 free(PCCDen_Rcv_Grid_A2B[ID]);
735 }
736 free(PCCDen_Rcv_Grid_A2B);
737 }
738
739 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
740 free(AtomDen_Grid[Mc_AN]);
741 }
742 free(AtomDen_Grid);
743
744 if (PCC_switch==1){
745 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
746 free(PCCDen_Grid[Mc_AN]);
747 }
748 free(PCCDen_Grid);
749 }
750
751 free(Work_Array_Snd_Grid_B2C);
752 free(Work_Array_Rcv_Grid_B2C);
753 free(Work_Array_Snd_Grid_B2D);
754 free(Work_Array_Rcv_Grid_B2D);
755
756 /* elapsed time */
757 dtime(&TEtime);
758 time0 = TEtime - TStime;
759 return time0;
760 }
761