1 /**********************************************************************
2 Set_Vpot.c:
3
4 Set_Vpot.c is a subroutine to calculate the value of local potential
5 on each grid point.
6
7 Log of Set_Vpot.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
22 static void Make_VNA_Grid();
23
24
Set_Vpot(int SCF_iter,int XC_P_switch,double ***** CDM)25 void Set_Vpot(int SCF_iter, int XC_P_switch, double *****CDM)
26 {
27 /****************************************************
28 XC_P_switch:
29 0 \epsilon_XC (XC energy density)
30 1 \mu_XC (XC potential)
31 2 \epsilon_XC - \mu_XC
32 ****************************************************/
33
34 int i,j,k,n,nmax,ri,Mc_AN,Gc_AN,Rn,GNc,GRc;
35 int Nc,Nd,n1,n2,n3,Cwan,spin,MN,BN,DN,ct_AN;
36 int N2D,GNs,BN_AB,GN_AB;
37 int h_AN,Gh_AN,Hwan,Rnh,N3[4];
38 int hNgrid1,hNgrid2,hNgrid3,Ng1,Ng2,Ng3;
39 double Gx,Gy,Gz,sum,x,y,z;
40 double Cxyz[4],dx,dy,dz,r,xc,yc,zc;
41 int numprocs,myid;
42 double Stime_atom,Etime_atom;
43
44 /* for OpenMP */
45 int OMPID,Nthrds,Nprocs;
46
47 /* MPI */
48 MPI_Comm_size(mpi_comm_level1,&numprocs);
49 MPI_Comm_rank(mpi_comm_level1,&myid);
50
51 /****************************************************
52 Vxc on grid
53 ****************************************************/
54
55 Set_XC_Grid(XC_P_switch,XC_switch,
56 Density_Grid_D[0],Density_Grid_D[1],
57 Density_Grid_D[2],Density_Grid_D[3],
58 Vxc_Grid_D[0], Vxc_Grid_D[1],
59 Vxc_Grid_D[2], Vxc_Grid_D[3],
60 NULL,NULL);
61
62 /****************************************************
63 copy Vxc_Grid_D to Vxc_Grid_B
64 ****************************************************/
65
66 Ng1 = Max_Grid_Index_D[1] - Min_Grid_Index_D[1] + 1;
67 Ng2 = Max_Grid_Index_D[2] - Min_Grid_Index_D[2] + 1;
68 Ng3 = Max_Grid_Index_D[3] - Min_Grid_Index_D[3] + 1;
69
70 for (n=0; n<Num_Rcv_Grid_B2D[myid]; n++){
71 DN = Index_Rcv_Grid_B2D[myid][n];
72 BN = Index_Snd_Grid_B2D[myid][n];
73
74 i = DN/(Ng2*Ng3);
75 j = (DN-i*Ng2*Ng3)/Ng3;
76 k = DN - i*Ng2*Ng3 - j*Ng3;
77
78 if ( !(i<=1 || (Ng1-2)<=i || j<=1 || (Ng2-2)<=j || k<=1 || (Ng3-2)<=k)){
79 for (spin=0; spin<=SpinP_switch; spin++){
80 Vxc_Grid_B[spin][BN] = Vxc_Grid_D[spin][DN];
81 }
82 }
83 }
84
85 /****************************************************
86 The neutral atom potential on grids
87 ****************************************************/
88
89 if (SCF_iter==1 && ProExpn_VNA==0) Make_VNA_Grid();
90
91 /****************************************************
92 external electric field
93 ****************************************************/
94
95 if ( SCF_iter==1 && E_Field_switch==1 ){
96
97 /* the center of the system */
98
99 xc = 0.0;
100 yc = 0.0;
101 zc = 0.0;
102
103 for (ct_AN=1; ct_AN<=atomnum; ct_AN++){
104 xc += Gxyz[ct_AN][1];
105 yc += Gxyz[ct_AN][2];
106 zc += Gxyz[ct_AN][3];
107 }
108
109 xc = xc/(double)atomnum;
110 yc = yc/(double)atomnum;
111 zc = zc/(double)atomnum;
112
113 hNgrid1 = Ngrid1/2;
114 hNgrid2 = Ngrid2/2;
115 hNgrid3 = Ngrid3/2;
116
117 /* calculate VEF_Grid in the partition B */
118
119 N2D = Ngrid1*Ngrid2;
120 GNs = ((myid*N2D+numprocs-1)/numprocs)*Ngrid3;
121
122 for (BN_AB=0; BN_AB<My_NumGridB_AB; BN_AB++){
123
124 GN_AB = BN_AB + GNs;
125 i = GN_AB/(Ngrid2*Ngrid3);
126 j = (GN_AB - i*(Ngrid2*Ngrid3))/Ngrid3;
127 k = GN_AB - i*(Ngrid2*Ngrid3) - j*Ngrid3;
128
129 Find_CGrids(1,i,j,k,Cxyz,N3);
130 i = N3[1];
131 j = N3[2];
132 k = N3[3];
133
134 dx = (double)i*length_gtv[1] + Grid_Origin[1] - xc;
135 dy = (double)j*length_gtv[2] + Grid_Origin[2] - yc;
136 dz = (double)k*length_gtv[3] + Grid_Origin[3] - zc;
137 VEF_Grid_B[BN_AB] = dx*E_Field[0] + dy*E_Field[1] + dz*E_Field[2];
138 }
139
140 /* MPI: from the partitions B to C */
141 Data_Grid_Copy_B2C_1( VEF_Grid_B, VEF_Grid );
142
143 } /* if ( SCF_iter<=2 && E_Field_switch==1 ) */
144
145 /****************************************************
146 Sum
147 ****************************************************/
148
149 /* spin non-collinear */
150 if (SpinP_switch==3){
151
152 /******************************
153 diagonal part
154 spin=0: 11
155 spin=1: 22
156 ******************************/
157
158 if ( E_Field_switch==1 ){
159
160 if (ProExpn_VNA==0){
161 for (spin=0; spin<=1; spin++){
162 for (MN=0; MN<My_NumGridB_AB; MN++){
163 Vpot_Grid_B[spin][MN] = F_dVHart_flag*dVHart_Grid_B[MN]
164 + F_Vxc_flag*Vxc_Grid_B[spin][MN]
165 + F_VNA_flag*VNA_Grid_B[MN]
166 + F_VEF_flag*VEF_Grid_B[MN];
167 }
168 }
169 }
170 else{
171 for (spin=0; spin<=1; spin++){
172 for (MN=0; MN<My_NumGridB_AB; MN++){
173 Vpot_Grid_B[spin][MN] = F_dVHart_flag*dVHart_Grid_B[MN]
174 + F_Vxc_flag*Vxc_Grid_B[spin][MN]
175 + F_VEF_flag*VEF_Grid_B[MN];
176 }
177 }
178 }
179
180 }
181
182 else{
183
184 if (ProExpn_VNA==0){
185 for (spin=0; spin<=1; spin++){
186 for (MN=0; MN<My_NumGridB_AB; MN++){
187 Vpot_Grid_B[spin][MN] = F_dVHart_flag*dVHart_Grid_B[MN]
188 + F_Vxc_flag*Vxc_Grid_B[spin][MN]
189 + F_VNA_flag*VNA_Grid_B[MN];
190 }
191 }
192 }
193 else{
194
195 for (spin=0; spin<=1; spin++){
196 for (MN=0; MN<My_NumGridB_AB; MN++){
197 Vpot_Grid_B[spin][MN] = F_dVHart_flag*dVHart_Grid_B[MN]
198 + F_Vxc_flag*Vxc_Grid_B[spin][MN];
199 }
200 }
201 }
202 }
203
204 /******************************
205 off-diagonal part
206 spin=2: real 12
207 spin=3: imaginary 12
208 ******************************/
209
210 for (spin=2; spin<=3; spin++){
211 for (MN=0; MN<My_NumGridB_AB; MN++){
212 Vpot_Grid_B[spin][MN] = F_Vxc_flag*Vxc_Grid_B[spin][MN];
213 }
214 }
215 }
216
217 /* spin collinear */
218 else{
219
220 if ( E_Field_switch==1 ){
221
222 if (ProExpn_VNA==0){
223 for (spin=0; spin<=SpinP_switch; spin++){
224 for (MN=0; MN<My_NumGridB_AB; MN++){
225 Vpot_Grid_B[spin][MN] = F_dVHart_flag*dVHart_Grid_B[MN]
226 + F_Vxc_flag*Vxc_Grid_B[spin][MN]
227 + F_VNA_flag*VNA_Grid_B[MN]
228 + F_VEF_flag*VEF_Grid_B[MN];
229 }
230 }
231 }
232 else{
233 for (spin=0; spin<=SpinP_switch; spin++){
234 for (MN=0; MN<My_NumGridB_AB; MN++){
235 Vpot_Grid_B[spin][MN] = F_dVHart_flag*dVHart_Grid_B[MN]
236 + F_Vxc_flag*Vxc_Grid_B[spin][MN]
237 + F_VEF_flag*VEF_Grid_B[MN];
238 }
239 }
240 }
241
242 }
243 else{
244
245 if (ProExpn_VNA==0){
246 for (spin=0; spin<=SpinP_switch; spin++){
247 for (MN=0; MN<My_NumGridB_AB; MN++){
248 Vpot_Grid_B[spin][MN] = F_dVHart_flag*dVHart_Grid_B[MN]
249 + F_Vxc_flag*Vxc_Grid_B[spin][MN]
250 + F_VNA_flag*VNA_Grid_B[MN];
251 }
252 }
253 }
254 else{
255 for (spin=0; spin<=SpinP_switch; spin++){
256 for (MN=0; MN<My_NumGridB_AB; MN++){
257 Vpot_Grid_B[spin][MN] = F_dVHart_flag*dVHart_Grid_B[MN]
258 + F_Vxc_flag*Vxc_Grid_B[spin][MN];
259 }
260 }
261 }
262 }
263 }
264
265 /******************************************************
266 MPI: from the partitions B to C
267 ******************************************************/
268
269 Data_Grid_Copy_B2C_2( Vpot_Grid_B, Vpot_Grid );
270 }
271
272
273
Make_VNA_Grid()274 void Make_VNA_Grid()
275 {
276 static int firsttime=1;
277 unsigned long long int n2D,N2D,GNc,GN;
278 int i,Mc_AN,Gc_AN,BN,CN,LN,GRc,N3[4];
279 int AN,Nc,MN,Cwan,NN_S,NN_R;
280 int size_AtomVNA_Grid;
281 int size_AtomVNA_Snd_Grid_A2B;
282 int size_AtomVNA_Rcv_Grid_A2B;
283 double Cxyz[4];
284 double r,dx,dy,dz;
285 double **AtomVNA_Grid;
286 double **AtomVNA_Snd_Grid_A2B;
287 double **AtomVNA_Rcv_Grid_A2B;
288 double Stime_atom, Etime_atom;
289 int numprocs,myid,tag=999,ID,IDS,IDR;
290 int OMPID,Nthrds,Nprocs;
291
292 MPI_Status stat;
293 MPI_Request request;
294 MPI_Status *stat_send;
295 MPI_Status *stat_recv;
296 MPI_Request *request_send;
297 MPI_Request *request_recv;
298
299 /* MPI */
300 MPI_Comm_size(mpi_comm_level1,&numprocs);
301 MPI_Comm_rank(mpi_comm_level1,&myid);
302
303 /* allocation of arrays */
304
305 size_AtomVNA_Grid = 1;
306 AtomVNA_Grid = (double**)malloc(sizeof(double*)*(Matomnum+1));
307 AtomVNA_Grid[0] = (double*)malloc(sizeof(double)*1);
308 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
309 Gc_AN = F_M2G[Mc_AN];
310 AtomVNA_Grid[Mc_AN] = (double*)malloc(sizeof(double)*GridN_Atom[Gc_AN]);
311 size_AtomVNA_Grid += GridN_Atom[Gc_AN];
312 }
313
314 size_AtomVNA_Snd_Grid_A2B = 0;
315 AtomVNA_Snd_Grid_A2B = (double**)malloc(sizeof(double*)*numprocs);
316 for (ID=0; ID<numprocs; ID++){
317 AtomVNA_Snd_Grid_A2B[ID] = (double*)malloc(sizeof(double)*Num_Snd_Grid_A2B[ID]);
318 size_AtomVNA_Snd_Grid_A2B += Num_Snd_Grid_A2B[ID];
319 }
320
321 size_AtomVNA_Rcv_Grid_A2B = 0;
322 AtomVNA_Rcv_Grid_A2B = (double**)malloc(sizeof(double*)*numprocs);
323 for (ID=0; ID<numprocs; ID++){
324 AtomVNA_Rcv_Grid_A2B[ID] = (double*)malloc(sizeof(double)*Num_Rcv_Grid_A2B[ID]);
325 size_AtomVNA_Rcv_Grid_A2B += Num_Rcv_Grid_A2B[ID];
326 }
327
328 /* PrintMemory */
329 if (firsttime) {
330 PrintMemory("Set_Vpot: AtomVNA_Grid",sizeof(double)*size_AtomVNA_Grid,NULL);
331 PrintMemory("Set_Vpot: AtomVNA_Snd_Grid_A2B",sizeof(double)*size_AtomVNA_Snd_Grid_A2B,NULL);
332 PrintMemory("Set_Vpot: AtomVNA_Rcv_Grid_A2B",sizeof(double)*size_AtomVNA_Rcv_Grid_A2B,NULL);
333 }
334
335 /* calculation of AtomVNA_Grid */
336
337 for (MN=0; MN<My_NumGridC; MN++) VNA_Grid[MN] = 0.0;
338
339 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
340
341 dtime(&Stime_atom);
342
343 Gc_AN = M2G[Mc_AN];
344 Cwan = WhatSpecies[Gc_AN];
345
346 #pragma omp parallel shared(AtomVNA_Grid,GridN_Atom,atv,Gxyz,Gc_AN,Cwan,Mc_AN,GridListAtom,CellListAtom) private(OMPID,Nthrds,Nprocs,Nc,GNc,GRc,Cxyz,dx,dy,dz,r)
347 {
348
349 OMPID = omp_get_thread_num();
350 Nthrds = omp_get_num_threads();
351 Nprocs = omp_get_num_procs();
352
353 for (Nc=OMPID*GridN_Atom[Gc_AN]/Nthrds; Nc<(OMPID+1)*GridN_Atom[Gc_AN]/Nthrds; Nc++){
354
355 GNc = GridListAtom[Mc_AN][Nc];
356 GRc = CellListAtom[Mc_AN][Nc];
357
358 Get_Grid_XYZ(GNc,Cxyz);
359 dx = Cxyz[1] + atv[GRc][1] - Gxyz[Gc_AN][1];
360 dy = Cxyz[2] + atv[GRc][2] - Gxyz[Gc_AN][2];
361 dz = Cxyz[3] + atv[GRc][3] - Gxyz[Gc_AN][3];
362
363 r = sqrt(dx*dx + dy*dy + dz*dz);
364 AtomVNA_Grid[Mc_AN][Nc] = VNAF(Cwan,r);
365 }
366
367 #pragma omp flush(AtomVNA_Grid)
368
369 } /* #pragma omp parallel */
370
371 dtime(&Etime_atom);
372 time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
373
374 } /* Mc_AN */
375
376 /******************************************************
377 MPI communication from the partitions A to B
378 ******************************************************/
379
380 /* copy AtomVNA_Grid to AtomVNA_Snd_Grid_A2B */
381
382 for (ID=0; ID<numprocs; ID++) Num_Snd_Grid_A2B[ID] = 0;
383
384 N2D = Ngrid1*Ngrid2;
385
386 for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
387
388 Gc_AN = M2G[Mc_AN];
389
390 for (AN=0; AN<GridN_Atom[Gc_AN]; AN++){
391
392 GN = GridListAtom[Mc_AN][AN];
393 GN2N(GN,N3);
394 n2D = N3[1]*Ngrid2 + N3[2];
395 ID = (int)(n2D*(unsigned long long int)numprocs/N2D);
396 AtomVNA_Snd_Grid_A2B[ID][Num_Snd_Grid_A2B[ID]] = AtomVNA_Grid[Mc_AN][AN];
397
398 Num_Snd_Grid_A2B[ID]++;
399 }
400 }
401
402 /* MPI: A to B */
403
404 request_send = malloc(sizeof(MPI_Request)*NN_A2B_S);
405 request_recv = malloc(sizeof(MPI_Request)*NN_A2B_R);
406 stat_send = malloc(sizeof(MPI_Status)*NN_A2B_S);
407 stat_recv = malloc(sizeof(MPI_Status)*NN_A2B_R);
408
409 NN_S = 0;
410 NN_R = 0;
411
412 tag = 999;
413 for (ID=1; ID<numprocs; ID++){
414
415 IDS = (myid + ID) % numprocs;
416 IDR = (myid - ID + numprocs) % numprocs;
417
418 if (Num_Snd_Grid_A2B[IDS]!=0){
419 MPI_Isend(&AtomVNA_Snd_Grid_A2B[IDS][0], Num_Snd_Grid_A2B[IDS], MPI_DOUBLE,
420 IDS, tag, mpi_comm_level1, &request_send[NN_S]);
421 NN_S++;
422 }
423
424 if (Num_Rcv_Grid_A2B[IDR]!=0){
425 MPI_Irecv( &AtomVNA_Rcv_Grid_A2B[IDR][0], Num_Rcv_Grid_A2B[IDR],
426 MPI_DOUBLE, IDR, tag, mpi_comm_level1, &request_recv[NN_R]);
427 NN_R++;
428 }
429 }
430
431 if (NN_S!=0) MPI_Waitall(NN_S,request_send,stat_send);
432 if (NN_R!=0) MPI_Waitall(NN_R,request_recv,stat_recv);
433
434 free(request_send);
435 free(request_recv);
436 free(stat_send);
437 free(stat_recv);
438
439 /* for myid */
440 for (i=0; i<Num_Rcv_Grid_A2B[myid]; i++){
441 AtomVNA_Rcv_Grid_A2B[myid][i] = AtomVNA_Snd_Grid_A2B[myid][i];
442 }
443
444 /******************************************************
445 superposition of VNA in the partition B
446 ******************************************************/
447
448 /* initialize VNA_Grid_B */
449
450 for (BN=0; BN<My_NumGridB_AB; BN++) VNA_Grid_B[BN] = 0.0;
451
452 /* superposition of VNA */
453
454 for (ID=0; ID<numprocs; ID++){
455 for (LN=0; LN<Num_Rcv_Grid_A2B[ID]; LN++){
456
457 BN = Index_Rcv_Grid_A2B[ID][3*LN+0];
458 VNA_Grid_B[BN] += AtomVNA_Rcv_Grid_A2B[ID][LN];
459
460 } /* LN */
461 } /* ID */
462
463 /******************************************************
464 MPI: from the partitions B to C
465 ******************************************************/
466
467 Data_Grid_Copy_B2C_1( VNA_Grid_B, VNA_Grid );
468
469 /* freeing of arrays */
470
471 for (Mc_AN=0; Mc_AN<=Matomnum; Mc_AN++){
472 free(AtomVNA_Grid[Mc_AN]);
473 }
474 free(AtomVNA_Grid);
475
476 for (ID=0; ID<numprocs; ID++){
477 free(AtomVNA_Snd_Grid_A2B[ID]);
478 }
479 free(AtomVNA_Snd_Grid_A2B);
480
481 for (ID=0; ID<numprocs; ID++){
482 free(AtomVNA_Rcv_Grid_A2B[ID]);
483 }
484 free(AtomVNA_Rcv_Grid_A2B);
485 }
486