1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <time.h>
5 #include "openmx_common.h"
6 #include "mpi.h"
7 #include <omp.h>
8
9 /*---------- added by TOYODA 14/JAN/2010 */
10 #include "exx_xc.h"
11 /*---------- until here */
12
13
Set_XC_Grid(int XC_P_switch,int XC_switch,double * Den0,double * Den1,double * Den2,double * Den3,double * Vxc0,double * Vxc1,double * Vxc2,double * Vxc3,double *** dEXC_dGD,double *** dDen_Grid)14 void Set_XC_Grid(int XC_P_switch, int XC_switch,
15 double *Den0, double *Den1,
16 double *Den2, double *Den3,
17 double *Vxc0, double *Vxc1,
18 double *Vxc2, double *Vxc3,
19 double ***dEXC_dGD,
20 double ***dDen_Grid)
21 {
22 /***********************************************************************************************
23 XC_P_switch:
24 0 \epsilon_XC (XC energy density)
25 1 \mu_XC (XC potential)
26 2 \epsilon_XC - \mu_XC
27 3 dfxc/d|\nabra n| *|d|\nabra n|/d\nabra n and nabla n
28
29 Switch 3 is implemetend by yshiihara to calculate the
30 terms for GGA stress. In this case, the input and output parameters are:
31
32 INPUT Den0 : density of up spin n_up
33 INPUT Den1 : density of down spin n_down
34 INPUT Den2 : theta
35 INPUT Den3 : phi
36 OUTPUT Vxc0 : NULL
37 OUTPUT Vxc1 : NULL
38 OUTPUT Vxc2 : NULL
39 OUTPUT Vxc3 : NULL
40 OUTPUT Axc[2(spin)][3(x,y,z)][My_NumGridD] : dfxc/d|\nabra n| *|d|\nabra n|/d\nabra n
41 OUTPUT dDen_Grid[2(spin)][3(x,y,z)][My_NumGridD] : dn/dx, dn/dy, dn/dz
42 ***********************************************************************************************/
43
44 static int firsttime=1;
45 int MN,MN1,MN2,i,j,k,ri,ri1,ri2;
46 int i1,i2,j1,j2,k1,k2,n,nmax;
47 int Ng1,Ng2,Ng3;
48 int dDen_Grid_NULL_flag;
49 int dEXC_dGD_NULL_flag;
50 double den_min=1.0e-14;
51 double Ec_unif[1],Vc_unif[2],Exc[2],Vxc[2];
52 double Ex_unif[1],Vx_unif[2],tot_den;
53 double ED[2],GDENS[3][2];
54 double DEXDD[2],DECDD[2];
55 double DEXDGD[3][2],DECDGD[3][2];
56 double up_x_a,up_x_b,up_x_c;
57 double up_y_a,up_y_b,up_y_c;
58 double up_z_a,up_z_b,up_z_c;
59 double dn_x_a,dn_x_b,dn_x_c;
60 double dn_y_a,dn_y_b,dn_y_c;
61 double dn_z_a,dn_z_b,dn_z_c;
62 double up_a,up_b,up_c;
63 double dn_a,dn_b,dn_c;
64
65 double tmp0,tmp1;
66 double cot,sit,sip,cop,phi,theta;
67 double detA,igtv[4][4];
68 int numprocs,myid;
69
70 /* for OpenMP */
71 int OMPID,Nthrds;
72
73 /* MPI */
74 MPI_Comm_size(mpi_comm_level1,&numprocs);
75 MPI_Comm_rank(mpi_comm_level1,&myid);
76
77 /****************************************************
78 allocation of arrays
79 ****************************************************/
80
81 dDen_Grid_NULL_flag = 0;
82
83 if (XC_switch==4 && dDen_Grid==NULL){
84
85 dDen_Grid_NULL_flag = 1;
86
87 dDen_Grid = (double***)malloc(sizeof(double**)*2);
88 for (k=0; k<=1; k++){
89 dDen_Grid[k] = (double**)malloc(sizeof(double*)*3);
90 for (i=0; i<3; i++){
91 dDen_Grid[k][i] = (double*)malloc(sizeof(double)*My_NumGridD);
92 for (j=0; j<My_NumGridD; j++) dDen_Grid[k][i][j] = 0.0;
93 }
94 }
95 }
96
97 dEXC_dGD_NULL_flag = 0;
98
99 if (XC_switch==4 && dEXC_dGD==NULL){
100
101 dEXC_dGD_NULL_flag = 1;
102
103 dEXC_dGD = (double***)malloc(sizeof(double**)*2);
104 for (k=0; k<=1; k++){
105 dEXC_dGD[k] = (double**)malloc(sizeof(double*)*3);
106 for (i=0; i<3; i++){
107 dEXC_dGD[k][i] = (double*)malloc(sizeof(double)*My_NumGridD);
108 for (j=0; j<My_NumGridD; j++) dEXC_dGD[k][i][j] = 0.0;
109 }
110 }
111 }
112
113 /****************************************************
114 calculate dDen_Grid
115 ****************************************************/
116
117 if (XC_switch==4){
118
119 detA = gtv[1][1]*gtv[2][2]*gtv[3][3]
120 + gtv[1][2]*gtv[2][3]*gtv[3][1]
121 + gtv[1][3]*gtv[2][1]*gtv[3][2]
122 - gtv[1][3]*gtv[2][2]*gtv[3][1]
123 - gtv[1][2]*gtv[2][1]*gtv[3][3]
124 - gtv[1][1]*gtv[2][3]*gtv[3][2];
125
126 igtv[1][1] = (gtv[2][2]*gtv[3][3] - gtv[2][3]*gtv[3][2])/detA;
127 igtv[2][1] = -(gtv[2][1]*gtv[3][3] - gtv[2][3]*gtv[3][1])/detA;
128 igtv[3][1] = (gtv[2][1]*gtv[3][2] - gtv[2][2]*gtv[3][1])/detA;
129
130 igtv[1][2] = -(gtv[1][2]*gtv[3][3] - gtv[1][3]*gtv[3][2])/detA;
131 igtv[2][2] = (gtv[1][1]*gtv[3][3] - gtv[1][3]*gtv[3][1])/detA;
132 igtv[3][2] = -(gtv[1][1]*gtv[3][2] - gtv[1][2]*gtv[3][1])/detA;
133
134 igtv[1][3] = (gtv[1][2]*gtv[2][3] - gtv[1][3]*gtv[2][2])/detA;
135 igtv[2][3] = -(gtv[1][1]*gtv[2][3] - gtv[1][3]*gtv[2][1])/detA;
136 igtv[3][3] = (gtv[1][1]*gtv[2][2] - gtv[1][2]*gtv[2][1])/detA;
137
138 #pragma omp parallel shared(My_NumGridD,Min_Grid_Index_D,Max_Grid_Index_D,igtv,dDen_Grid,PCCDensity_Grid_D,PCC_switch,Den0,Den1,Den2,Den3,den_min) private(OMPID,Nthrds,nmax,n,i,j,k,ri,ri1,ri2,i1,i2,j1,j2,k1,k2,MN,MN1,MN2,up_a,dn_a,up_b,dn_b,up_c,dn_c,Ng1,Ng2,Ng3)
139 {
140
141 OMPID = omp_get_thread_num();
142 Nthrds = omp_get_num_threads();
143
144 Ng1 = Max_Grid_Index_D[1] - Min_Grid_Index_D[1] + 1;
145 Ng2 = Max_Grid_Index_D[2] - Min_Grid_Index_D[2] + 1;
146 Ng3 = Max_Grid_Index_D[3] - Min_Grid_Index_D[3] + 1;
147
148 for (MN=OMPID; MN<My_NumGridD; MN+=Nthrds){
149
150 i = MN/(Ng2*Ng3);
151 j = (MN-i*Ng2*Ng3)/Ng3;
152 k = MN - i*Ng2*Ng3 - j*Ng3;
153
154 if ( i==0 || i==(Ng1-1) || j==0 || j==(Ng2-1) || k==0 || k==(Ng3-1) ){
155
156 dDen_Grid[0][0][MN] = 0.0;
157 dDen_Grid[0][1][MN] = 0.0;
158 dDen_Grid[0][2][MN] = 0.0;
159 dDen_Grid[1][0][MN] = 0.0;
160 dDen_Grid[1][1][MN] = 0.0;
161 dDen_Grid[1][2][MN] = 0.0;
162 }
163
164 else {
165
166 /* set i1, i2, j1, j2, k1, and k2 */
167
168 i1 = i - 1;
169 i2 = i + 1;
170
171 j1 = j - 1;
172 j2 = j + 1;
173
174 k1 = k - 1;
175 k2 = k + 1;
176
177 /* set dDen_Grid */
178
179 if ( den_min<(Den0[MN]+Den1[MN]) ){
180
181 /* a-axis */
182
183 MN1 = i1*Ng2*Ng3 + j*Ng3 + k;
184 MN2 = i2*Ng2*Ng3 + j*Ng3 + k;
185
186 if (PCC_switch==0) {
187 up_a = Den0[MN2] - Den0[MN1];
188 dn_a = Den1[MN2] - Den1[MN1];
189 }
190 else if (PCC_switch==1) {
191 up_a = Den0[MN2] + PCCDensity_Grid_D[0][MN2]
192 - Den0[MN1] - PCCDensity_Grid_D[0][MN1];
193 dn_a = Den1[MN2] + PCCDensity_Grid_D[1][MN2]
194 - Den1[MN1] - PCCDensity_Grid_D[1][MN1];
195 }
196
197 /* b-axis */
198
199 MN1 = i*Ng2*Ng3 + j1*Ng3 + k;
200 MN2 = i*Ng2*Ng3 + j2*Ng3 + k;
201
202 if (PCC_switch==0) {
203 up_b = Den0[MN2] - Den0[MN1];
204 dn_b = Den1[MN2] - Den1[MN1];
205 }
206 else if (PCC_switch==1) {
207 up_b = Den0[MN2] + PCCDensity_Grid_D[0][MN2]
208 - Den0[MN1] - PCCDensity_Grid_D[0][MN1];
209 dn_b = Den1[MN2] + PCCDensity_Grid_D[1][MN2]
210 - Den1[MN1] - PCCDensity_Grid_D[1][MN1];
211 }
212
213 /* c-axis */
214
215 MN1 = i*Ng2*Ng3 + j*Ng3 + k1;
216 MN2 = i*Ng2*Ng3 + j*Ng3 + k2;
217
218 if (PCC_switch==0) {
219 up_c = Den0[MN2] - Den0[MN1];
220 dn_c = Den1[MN2] - Den1[MN1];
221 }
222 else if (PCC_switch==1) {
223 up_c = Den0[MN2] + PCCDensity_Grid_D[0][MN2]
224 - Den0[MN1] - PCCDensity_Grid_D[0][MN1];
225 dn_c = Den1[MN2] + PCCDensity_Grid_D[1][MN2]
226 - Den1[MN1] - PCCDensity_Grid_D[1][MN1];
227 }
228
229 /* up */
230
231 dDen_Grid[0][0][MN] = 0.5*(igtv[1][1]*up_a + igtv[1][2]*up_b + igtv[1][3]*up_c);
232 dDen_Grid[0][1][MN] = 0.5*(igtv[2][1]*up_a + igtv[2][2]*up_b + igtv[2][3]*up_c);
233 dDen_Grid[0][2][MN] = 0.5*(igtv[3][1]*up_a + igtv[3][2]*up_b + igtv[3][3]*up_c);
234
235 /* down */
236
237 dDen_Grid[1][0][MN] = 0.5*(igtv[1][1]*dn_a + igtv[1][2]*dn_b + igtv[1][3]*dn_c);
238 dDen_Grid[1][1][MN] = 0.5*(igtv[2][1]*dn_a + igtv[2][2]*dn_b + igtv[2][3]*dn_c);
239 dDen_Grid[1][2][MN] = 0.5*(igtv[3][1]*dn_a + igtv[3][2]*dn_b + igtv[3][3]*dn_c);
240 }
241
242 else{
243 dDen_Grid[0][0][MN] = 0.0;
244 dDen_Grid[0][1][MN] = 0.0;
245 dDen_Grid[0][2][MN] = 0.0;
246 dDen_Grid[1][0][MN] = 0.0;
247 dDen_Grid[1][1][MN] = 0.0;
248 dDen_Grid[1][2][MN] = 0.0;
249 }
250
251 } /* else */
252
253 } /* n */
254
255 #pragma omp flush(dDen_Grid)
256
257 } /* #pragma omp parallel */
258 } /* if (XC_switch==4) */
259
260 /****************************************************
261 loop MN
262 ****************************************************/
263
264 #pragma omp parallel shared(dDen_Grid,dEXC_dGD,den_min,Vxc0,Vxc1,Vxc2,Vxc3,My_NumGridD,XC_P_switch,XC_switch,Den0,Den1,Den2,Den3,PCC_switch,PCCDensity_Grid_D) private(OMPID,Nthrds,MN,tot_den,tmp0,tmp1,ED,Exc,Ec_unif,Vc_unif,Vxc,Ex_unif,Vx_unif,GDENS,DEXDD,DECDD,DEXDGD,DECDGD)
265 {
266
267 OMPID = omp_get_thread_num();
268 Nthrds = omp_get_num_threads();
269
270 for (MN=OMPID; MN<My_NumGridD; MN+=Nthrds){
271
272 switch(XC_switch){
273
274 /******************************************************************
275 LDA (Ceperly-Alder)
276
277 constructed by Ceperly and Alder,
278 ref.
279 D. M. Ceperley, Phys. Rev. B18, 3126 (1978)
280 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett., 45, 566 (1980)
281
282 and parametrized by Perdew and Zunger.
283 ref.
284 J. Perdew and A. Zunger, Phys. Rev. B23, 5048 (1981)
285 ******************************************************************/
286
287 case 1:
288
289 tot_den = Den0[MN] + Den1[MN];
290
291 /* partial core correction */
292 if (PCC_switch==1) {
293 tot_den += PCCDensity_Grid_D[0][MN] + PCCDensity_Grid_D[1][MN];
294 }
295
296 tmp0 = XC_Ceperly_Alder(tot_den,XC_P_switch);
297
298 Vxc0[MN] = tmp0;
299 Vxc1[MN] = tmp0;
300
301 break;
302
303 /******************************************************************
304 LSDA-CA (Ceperly-Alder)
305
306 constructed by Ceperly and Alder,
307 ref.
308 D. M. Ceperley, Phys. Rev. B18, 3126 (1978)
309 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett., 45, 566 (1980)
310
311 and parametrized by Perdew and Zunger.
312 ref.
313 J. Perdew and A. Zunger, Phys. Rev. B23, 5048 (1981)
314 ******************************************************************/
315
316 case 2:
317
318 ED[0] = Den0[MN];
319 ED[1] = Den1[MN];
320
321 /* partial core correction */
322 if (PCC_switch==1) {
323 ED[0] += PCCDensity_Grid_D[0][MN];
324 ED[1] += PCCDensity_Grid_D[1][MN];
325 }
326
327 XC_CA_LSDA(ED[0], ED[1], Exc, XC_P_switch);
328 Vxc0[MN] = Exc[0];
329 Vxc1[MN] = Exc[1];
330
331 break;
332
333 /******************************************************************
334 LSDA-PW (PW92)
335
336 ref.
337 J.P.Perdew and Yue Wang, Phys. Rev. B45, 13244 (1992)
338 ******************************************************************/
339
340 case 3:
341
342 ED[0] = Den0[MN];
343 ED[1] = Den1[MN];
344
345 /* partial core correction */
346 if (PCC_switch==1) {
347 ED[0] += PCCDensity_Grid_D[0][MN];
348 ED[1] += PCCDensity_Grid_D[1][MN];
349 }
350
351 if ((ED[0]+ED[1])<den_min){
352 Vxc0[MN] = 0.0;
353 Vxc1[MN] = 0.0;
354 }
355 else{
356
357 if (XC_P_switch==0){
358
359 XC_PW92C(ED,Ec_unif,Vc_unif);
360
361 Vxc[0] = Vc_unif[0];
362 Vxc[1] = Vc_unif[1];
363 Exc[0] = Ec_unif[0];
364
365 XC_EX(1,2.0*ED[0],ED,Ex_unif,Vx_unif);
366 Vxc[0] = Vxc[0] + Vx_unif[0];
367 Exc[1] = 2.0*ED[0]*Ex_unif[0];
368
369 XC_EX(1,2.0*ED[1],ED,Ex_unif,Vx_unif);
370 Vxc[1] += Vx_unif[0];
371 Exc[1] += 2.0*ED[1]*Ex_unif[0];
372
373 Exc[1] = 0.5*Exc[1]/(ED[0]+ED[1]);
374
375 Vxc0[MN] = Exc[0] + Exc[1];
376 Vxc1[MN] = Exc[0] + Exc[1];
377 }
378
379 else if (XC_P_switch==1){
380 XC_PW92C(ED,Ec_unif,Vc_unif);
381 Vxc0[MN] = Vc_unif[0];
382 Vxc1[MN] = Vc_unif[1];
383
384 XC_EX(1,2.0*ED[0],ED,Ex_unif,Vx_unif);
385 Vxc0[MN] = Vxc0[MN] + Vx_unif[0];
386
387 XC_EX(1,2.0*ED[1],ED,Ex_unif,Vx_unif);
388 Vxc1[MN] = Vxc1[MN] + Vx_unif[0];
389 }
390
391 else if (XC_P_switch==2){
392
393 XC_PW92C(ED,Ec_unif,Vc_unif);
394
395 Vxc[0] = Vc_unif[0];
396 Vxc[1] = Vc_unif[1];
397 Exc[0] = Ec_unif[0];
398
399 XC_EX(1,2.0*ED[0],ED,Ex_unif,Vx_unif);
400 Vxc[0] = Vxc[0] + Vx_unif[0];
401 Exc[1] = 2.0*ED[0]*Ex_unif[0];
402
403 XC_EX(1,2.0*ED[1],ED,Ex_unif,Vx_unif);
404 Vxc[1] += Vx_unif[0];
405 Exc[1] += 2.0*ED[1]*Ex_unif[0];
406
407 Exc[1] = 0.5*Exc[1]/(ED[0]+ED[1]);
408
409 Vxc0[MN] = Exc[0] + Exc[1] - Vxc[0];
410 Vxc1[MN] = Exc[0] + Exc[1] - Vxc[1];
411 }
412 }
413
414 break;
415
416 /******************************************************************
417 GGA-PBE
418 ref.
419 J. P. Perdew, K. Burke, and M. Ernzerhof,
420 Phys. Rev. Lett. 77, 3865 (1996).
421 ******************************************************************/
422
423 case 4:
424
425 /****************************************************
426 ED[0] density of up spin: n_up
427 ED[1] density of down spin: n_down
428
429 GDENS[0][0] derivative (x) of density of up spin
430 GDENS[1][0] derivative (y) of density of up spin
431 GDENS[2][0] derivative (z) of density of up spin
432 GDENS[0][1] derivative (x) of density of down spin
433 GDENS[1][1] derivative (y) of density of down spin
434 GDENS[2][1] derivative (z) of density of down spin
435
436 DEXDD[0] d(fx)/d(n_up)
437 DEXDD[1] d(fx)/d(n_down)
438 DECDD[0] d(fc)/d(n_up)
439 DECDD[1] d(fc)/d(n_down)
440
441 n'_up_x = d(n_up)/d(x)
442 n'_up_y = d(n_up)/d(y)
443 n'_up_z = d(n_up)/d(z)
444 n'_down_x = d(n_down)/d(x)
445 n'_down_y = d(n_down)/d(y)
446 n'_down_z = d(n_down)/d(z)
447
448 DEXDGD[0][0] d(fx)/d(n'_up_x)
449 DEXDGD[1][0] d(fx)/d(n'_up_y)
450 DEXDGD[2][0] d(fx)/d(n'_up_z)
451 DEXDGD[0][1] d(fx)/d(n'_down_x)
452 DEXDGD[1][1] d(fx)/d(n'_down_y)
453 DEXDGD[2][1] d(fx)/d(n'_down_z)
454
455 DECDGD[0][0] d(fc)/d(n'_up_x)
456 DECDGD[1][0] d(fc)/d(n'_up_y)
457 DECDGD[2][0] d(fc)/d(n'_up_z)
458 DECDGD[0][1] d(fc)/d(n'_down_x)
459 DECDGD[1][1] d(fc)/d(n'_down_y)
460 DECDGD[2][1] d(fc)/d(n'_down_z)
461 ****************************************************/
462
463 ED[0] = Den0[MN];
464 ED[1] = Den1[MN];
465
466 if ((ED[0]+ED[1])<den_min){
467
468 if (XC_P_switch!=3){
469 Vxc0[MN] = 0.0;
470 Vxc1[MN] = 0.0;
471 }
472
473 /* later add its derivatives */
474 if (XC_P_switch!=0){
475 dEXC_dGD[0][0][MN] = 0.0;
476 dEXC_dGD[0][1][MN] = 0.0;
477 dEXC_dGD[0][2][MN] = 0.0;
478
479 dEXC_dGD[1][0][MN] = 0.0;
480 dEXC_dGD[1][1][MN] = 0.0;
481 dEXC_dGD[1][2][MN] = 0.0;
482 }
483 }
484
485 else{
486
487 GDENS[0][0] = dDen_Grid[0][0][MN];
488 GDENS[1][0] = dDen_Grid[0][1][MN];
489 GDENS[2][0] = dDen_Grid[0][2][MN];
490 GDENS[0][1] = dDen_Grid[1][0][MN];
491 GDENS[1][1] = dDen_Grid[1][1][MN];
492 GDENS[2][1] = dDen_Grid[1][2][MN];
493
494 if (PCC_switch==1) {
495 ED[0] += PCCDensity_Grid_D[0][MN];
496 ED[1] += PCCDensity_Grid_D[1][MN];
497 }
498
499 XC_PBE(ED, GDENS, Exc, DEXDD, DECDD, DEXDGD, DECDGD);
500
501 /* XC energy density */
502 if (XC_P_switch==0){
503 Vxc0[MN] = Exc[0] + Exc[1];
504 Vxc1[MN] = Exc[0] + Exc[1];
505 }
506
507 /* XC potential */
508 else if (XC_P_switch==1){
509 Vxc0[MN] = DEXDD[0] + DECDD[0];
510 Vxc1[MN] = DEXDD[1] + DECDD[1];
511 }
512
513 /* XC energy density - XC potential */
514 else if (XC_P_switch==2){
515 Vxc0[MN] = Exc[0] + Exc[1] - DEXDD[0] - DECDD[0];
516 Vxc1[MN] = Exc[0] + Exc[1] - DEXDD[1] - DECDD[1];
517 }
518
519 /* later add its derivatives */
520 if (XC_P_switch!=0){
521 dEXC_dGD[0][0][MN] = DEXDGD[0][0] + DECDGD[0][0];
522 dEXC_dGD[0][1][MN] = DEXDGD[1][0] + DECDGD[1][0];
523 dEXC_dGD[0][2][MN] = DEXDGD[2][0] + DECDGD[2][0];
524
525 dEXC_dGD[1][0][MN] = DEXDGD[0][1] + DECDGD[0][1];
526 dEXC_dGD[1][1][MN] = DEXDGD[1][1] + DECDGD[1][1];
527 dEXC_dGD[1][2][MN] = DEXDGD[2][1] + DECDGD[2][1];
528 }
529
530 }
531
532 break;
533
534 /*---------- added by TOYODA 14/JAN/2010 from here */
535 case 5: /* EXX-TEST */
536 /* only X part of CA-LSDA XC is used */
537
538 ED[0] = Den0[MN];
539 ED[1] = Den1[MN];
540
541 /* partial core correction */
542 if (PCC_switch==1) {
543 ED[0] += PCCDensity_Grid_D[0][MN];
544 ED[1] += PCCDensity_Grid_D[1][MN];
545 }
546
547 EXX_XC_CA_LSDA(ED[0], ED[1], Exc, XC_P_switch);
548 Vxc0[MN] = Exc[0];
549 Vxc1[MN] = Exc[1];
550
551 break;
552 /*---------- added by TOYODA 14/JAN/2010 until here */
553
554 } /* switch(XC_switch) */
555 } /* MN */
556
557 if (XC_switch==4){
558 #pragma omp flush(dEXC_dGD)
559 }
560
561 } /* #pragma omp parallel */
562
563 /****************************************************
564 calculate the second part of XC potential
565 when GGA and XC_P_switch!=0
566 ****************************************************/
567
568 if (XC_switch==4 && (XC_P_switch==1 || XC_P_switch==2)){
569
570 #pragma omp parallel shared(Min_Grid_Index_D,Max_Grid_Index_D,My_NumGridD,XC_P_switch,Vxc0,Vxc1,Vxc2,Vxc3,igtv,dEXC_dGD,Den0,Den1,Den2,Den3,den_min) private(OMPID,Nthrds,nmax,i,j,k,ri,ri1,ri2,i1,i2,j1,j2,k1,k2,MN,MN1,MN2,up_x_a,up_y_a,up_z_a,dn_x_a,dn_y_a,dn_z_a,up_x_b,up_y_b,up_z_b,dn_x_b,dn_y_b,dn_z_b,up_x_c,up_y_c,up_z_c,dn_x_c,dn_y_c,dn_z_c,tmp0,tmp1,Ng1,Ng2,Ng3)
571 {
572
573 OMPID = omp_get_thread_num();
574 Nthrds = omp_get_num_threads();
575
576 Ng1 = Max_Grid_Index_D[1] - Min_Grid_Index_D[1] + 1;
577 Ng2 = Max_Grid_Index_D[2] - Min_Grid_Index_D[2] + 1;
578 Ng3 = Max_Grid_Index_D[3] - Min_Grid_Index_D[3] + 1;
579
580 for (MN=OMPID; MN<My_NumGridD; MN+=Nthrds){
581
582 i = MN/(Ng2*Ng3);
583 j = (MN-i*Ng2*Ng3)/Ng3;
584 k = MN - i*Ng2*Ng3 - j*Ng3;
585
586 if ( i<=1 || (Ng1-2)<=i || j<=1 || (Ng2-2)<=j || k<=1 || (Ng3-2)<=k ){
587
588 Vxc0[MN] = 0.0;
589 Vxc1[MN] = 0.0;
590 }
591
592 else {
593
594 /* set i1, i2, j1, j2, k1, and k2 */
595
596 i1 = i - 1;
597 i2 = i + 1;
598
599 j1 = j - 1;
600 j2 = j + 1;
601
602 k1 = k - 1;
603 k2 = k + 1;
604
605 /* set Vxc_Grid */
606
607 if ( den_min<(Den0[MN]+Den1[MN]) ){
608
609 /* a-axis */
610
611 MN1 = i1*Ng2*Ng3 + j*Ng3 + k;
612 MN2 = i2*Ng2*Ng3 + j*Ng3 + k;
613
614 up_x_a = dEXC_dGD[0][0][MN2] - dEXC_dGD[0][0][MN1];
615 up_y_a = dEXC_dGD[0][1][MN2] - dEXC_dGD[0][1][MN1];
616 up_z_a = dEXC_dGD[0][2][MN2] - dEXC_dGD[0][2][MN1];
617
618 dn_x_a = dEXC_dGD[1][0][MN2] - dEXC_dGD[1][0][MN1];
619 dn_y_a = dEXC_dGD[1][1][MN2] - dEXC_dGD[1][1][MN1];
620 dn_z_a = dEXC_dGD[1][2][MN2] - dEXC_dGD[1][2][MN1];
621
622 /* b-axis */
623
624 MN1 = i*Ng2*Ng3 + j1*Ng3 + k;
625 MN2 = i*Ng2*Ng3 + j2*Ng3 + k;
626
627 up_x_b = dEXC_dGD[0][0][MN2] - dEXC_dGD[0][0][MN1];
628 up_y_b = dEXC_dGD[0][1][MN2] - dEXC_dGD[0][1][MN1];
629 up_z_b = dEXC_dGD[0][2][MN2] - dEXC_dGD[0][2][MN1];
630
631 dn_x_b = dEXC_dGD[1][0][MN2] - dEXC_dGD[1][0][MN1];
632 dn_y_b = dEXC_dGD[1][1][MN2] - dEXC_dGD[1][1][MN1];
633 dn_z_b = dEXC_dGD[1][2][MN2] - dEXC_dGD[1][2][MN1];
634
635 /* c-axis */
636
637 MN1 = i*Ng2*Ng3 + j*Ng3 + k1;
638 MN2 = i*Ng2*Ng3 + j*Ng3 + k2;
639
640 up_x_c = dEXC_dGD[0][0][MN2] - dEXC_dGD[0][0][MN1];
641 up_y_c = dEXC_dGD[0][1][MN2] - dEXC_dGD[0][1][MN1];
642 up_z_c = dEXC_dGD[0][2][MN2] - dEXC_dGD[0][2][MN1];
643
644 dn_x_c = dEXC_dGD[1][0][MN2] - dEXC_dGD[1][0][MN1];
645 dn_y_c = dEXC_dGD[1][1][MN2] - dEXC_dGD[1][1][MN1];
646 dn_z_c = dEXC_dGD[1][2][MN2] - dEXC_dGD[1][2][MN1];
647
648 /* up */
649
650 tmp0 = igtv[1][1]*up_x_a + igtv[1][2]*up_x_b + igtv[1][3]*up_x_c
651 + igtv[2][1]*up_y_a + igtv[2][2]*up_y_b + igtv[2][3]*up_y_c
652 + igtv[3][1]*up_z_a + igtv[3][2]*up_z_b + igtv[3][3]*up_z_c;
653 tmp0 = 0.5*tmp0;
654
655 /* down */
656
657 tmp1 = igtv[1][1]*dn_x_a + igtv[1][2]*dn_x_b + igtv[1][3]*dn_x_c
658 + igtv[2][1]*dn_y_a + igtv[2][2]*dn_y_b + igtv[2][3]*dn_y_c
659 + igtv[3][1]*dn_z_a + igtv[3][2]*dn_z_b + igtv[3][3]*dn_z_c;
660 tmp1 = 0.5*tmp1;
661
662 /* XC potential */
663
664 if (XC_P_switch==1){
665 Vxc0[MN] -= tmp0;
666 Vxc1[MN] -= tmp1;
667 }
668
669 /* XC energy density - XC potential */
670
671 else if (XC_P_switch==2){
672 Vxc0[MN] += tmp0;
673 Vxc1[MN] += tmp1;
674 }
675
676 }
677 }
678 }
679
680 #pragma omp flush(Vxc0,Vxc1,Vxc2,Vxc3)
681
682 } /* #pragma omp parallel */
683 } /* if (XC_switch==4 && XC_P_switch!=0) */
684
685 /****************************************************
686 In case of non-collinear spin DFT
687 ****************************************************/
688
689 if (SpinP_switch==3 && (XC_P_switch==1 || XC_P_switch==2)){
690
691 #pragma omp parallel shared(Den0,Den1,Den2,Den3,Vxc0,Vxc1,Vxc2,Vxc3,My_NumGridD) private(OMPID,Nthrds,MN,tmp0,tmp1,theta,phi,sit,cot,sip,cop)
692 {
693
694 OMPID = omp_get_thread_num();
695 Nthrds = omp_get_num_threads();
696
697 for (MN=OMPID; MN<My_NumGridD; MN+=Nthrds){
698
699 tmp0 = 0.5*(Vxc0[MN] + Vxc1[MN]);
700 tmp1 = 0.5*(Vxc0[MN] - Vxc1[MN]);
701 theta = Den2[MN];
702 phi = Den3[MN];
703
704 sit = sin(theta);
705 cot = cos(theta);
706 sip = sin(phi);
707 cop = cos(phi);
708
709 /*******************************************************************
710 Since Set_XC_Grid is called as
711
712 XC_P_switch = 1;
713 Set_XC_Grid( XC_P_switch, 1,
714 ADensity_Grid_D, ADensity_Grid_D,
715 ADensity_Grid_D, ADensity_Grid_D,
716 RefVxc_Grid_D, RefVxc_Grid_D,
717 RefVxc_Grid_D, RefVxc_Grid_D );
718
719 XC_P_switch = 0;
720 Set_XC_Grid( XC_P_switch, 1,
721 ADensity_Grid_D, ADensity_Grid_D,
722 ADensity_Grid_D, ADensity_Grid_D,
723 RefVxc_Grid_D, RefVxc_Grid_D,
724 RefVxc_Grid_D, RefVxc_Grid_D );
725
726 from Force.c and Total_Energy.c, respectively,
727
728 data have to be stored in order of Vxc3, Vxc2, Vxc1, and Vxc0.
729 Note that only Vxc0 will be used for those calling.
730 ********************************************************************/
731
732 Vxc3[MN] = -tmp1*sit*sip; /* Im Vxc12 */
733 Vxc2[MN] = tmp1*sit*cop; /* Re Vxc12 */
734 Vxc1[MN] = tmp0 - cot*tmp1; /* Re Vxc22 */
735 Vxc0[MN] = tmp0 + cot*tmp1; /* Re Vxc11 */
736 }
737
738 #pragma omp flush(Vxc0,Vxc1,Vxc2,Vxc3)
739
740 } /* #pragma omp parallel */
741 }
742
743 /****************************************************
744 freeing of arrays
745 ****************************************************/
746
747 if (dDen_Grid_NULL_flag==1){
748 for (k=0; k<=1; k++){
749 for (i=0; i<3; i++){
750 free(dDen_Grid[k][i]);
751 }
752 free(dDen_Grid[k]);
753 }
754 free(dDen_Grid);
755 }
756
757 if (dEXC_dGD_NULL_flag==1){
758 for (k=0; k<=1; k++){
759 for (i=0; i<3; i++){
760 free(dEXC_dGD[k][i]);
761 }
762 free(dEXC_dGD[k]);
763 }
764 free(dEXC_dGD);
765 }
766
767 }
768