1 #include <math.h>
2 #include <stdio.h>
3 #include <assert.h>
4
5 #if defined HAVE_CONFIG_H
6 #include "config.h"
7 #endif
8
9 #include "vdw_core.h"
10 #include "vdw_q0.h"
11
12 // 1e-15 corresponds to kF~1e-5 which is the end of the standard q0 grid.
13 // However realistic values of q0 are more like 1e-3 to 5.
14 // libxc uses 1e-10 while QE uses 1e-12. Let us use 1e-12 then for
15 // consistency.
16 #define RHOCUT (1e-12)
17
18
vdwxc_calculate_q0(int N,double Z_ab,double q_cut,double * rho,double * sigma,double * q0,double * rho_dq0drho,double * rho_dq0dsigma)19 void vdwxc_calculate_q0(int N, double Z_ab, double q_cut,
20 double* rho, double* sigma,
21 double* q0, double* rho_dq0drho,
22 double* rho_dq0dsigma)
23 {
24 int i;
25 for (i=0; i<N; i++) {
26 double mrho = *rho;// * rhofactor;
27 double msigma = *sigma; //) * rhofactor;
28
29 if (mrho < RHOCUT) {
30 // TEMP FIX:
31 // for small densities we should in fact use q_cut.
32 // These things will be property of kernel later - probably.
33 mrho = RHOCUT;
34 *q0 = q_cut;
35 *rho_dq0dsigma = 0.0;
36 *rho_dq0drho = 0.0;
37 rho++;
38 sigma++;
39 q0++;
40 rho_dq0drho++;
41 rho_dq0dsigma++;
42 continue;
43 //if (*rho<0)
44 // printf("Warning, negative density %.15f\n", *rho);
45 }
46
47 // kF = (3.0 * np.pi**2 * total_rho)**(1.0 / 3.0)
48 double rho1_root3 = pow(mrho, 1.0 / 3.0);
49 double kF = VDWXC_KF_PREFACTOR * rho1_root3;
50
51 // r_s = (3.0 / (4.0 * np.pi * total_rho))**(1.0 / 3.0)
52 double rs = VDWXC_RS_PREFACTOR / rho1_root3;
53 double rs_sqr = rs*rs;
54 double rs_root2 = sqrt(rs);
55 double rho_sqr = mrho*mrho;
56
57 double gradient_correction = -Z_ab / (36.0 * kF * rho_sqr) * msigma;
58 double LDA_1 = 8.0 * VDWXC_PI / 3.0 * (VDWXC_LDA_A * (1.0 + VDWXC_LDA_a1 * rs));
59 double LDA_2 = 2.0 * VDWXC_LDA_A * (VDWXC_LDA_b1 * rs_root2 + VDWXC_LDA_b2 * rs + \
60 VDWXC_LDA_b3 * rs * rs_root2 + VDWXC_LDA_b4 * rs_sqr);
61
62 double q = kF + LDA_1 * log(1.0 + 1.0 / LDA_2) + gradient_correction;
63 int m;
64 double exponent = 0.0;
65 double dq0_dq = 0.0;
66 double qscaled0 = q / q_cut;
67 double qscaled = 1.0;
68
69 for (m=1; m<=12; m++) {
70 dq0_dq += qscaled;
71 qscaled *= qscaled0;
72 exponent += qscaled / m;
73 }
74
75 //printf("exponent %.15f", exponent);
76 //printf("dq0_dq %.15f", dq0_dq);
77
78 //mask = total_rho >= epsr
79 //q0[mask] = (q_cut * (1.0 - np.exp(-exponent)))[mask]
80 //dq0_dq *= np.exp(-exponent)
81 // Ignore mask for now: DOUBLE CHECK THIS XXX
82 dq0_dq *= exp(-exponent);
83 *q0 = q_cut * (1.0 - exp(-exponent));
84
85 *rho_dq0drho = dq0_dq * (kF / 3.0 - 7.0 / 3.0 * gradient_correction
86 - 8.0 * VDWXC_PI / 9.0 * VDWXC_LDA_A * VDWXC_LDA_a1 * rs
87 * log(1.0 + 1.0 / LDA_2)
88 + LDA_1 / (LDA_2 * (1.0 + LDA_2))
89 * (2.0 * VDWXC_LDA_A * (VDWXC_LDA_b1 / 6.0 * rs_root2
90 + VDWXC_LDA_b2 / 3.0 * rs
91 + VDWXC_LDA_b3 / 2.0 * rs * rs_root2
92 + 2.0 * VDWXC_LDA_b4 / 3.0 * rs_sqr)));
93
94 *rho_dq0dsigma = -dq0_dq * Z_ab / (36.0 * kF * mrho);
95
96 rho++;
97 sigma++;
98 q0++;
99 rho_dq0drho++;
100 rho_dq0dsigma++;
101 }
102 }
103
104
105 // Perdew-wang G[A, alpha1, beta1, ..., beta4](rs) function
vdwxc_G_pw(double A,double a1,double b1,double b2,double b3,double b4,double rootrs,double * e,double * dedrs)106 void vdwxc_G_pw(double A, double a1,
107 double b1, double b2, double b3, double b4,
108 double rootrs, double *e, double *dedrs)
109 {
110 double rs = rootrs * rootrs;
111 double rootrs3 = rootrs * rs;
112 double rs2 = rootrs * rootrs3;
113 double g0 = -2.0 * A * (1.0 + a1 * rs);
114 double g1 = 2.0 * A * (b1 * rootrs + b2 * rs
115 + b3 * rootrs3 + b4 * rs2);
116 *e = g0 * log(1.0 + 1.0 / g1);
117 double dg1drs = A * (b1 / rootrs + 2.0 * b2 +
118 rootrs * (3.0 * b3 + 4.0 * b4 * rootrs));
119 *dedrs = -2.0 * A * a1 * (*e) / g0 - g0 * dg1drs / (g1 * (1.0 + g1));
120 }
121
122
vdwxc_G_zeta0(double rootrs,double * e,double * dedn)123 void vdwxc_G_zeta0(double rootrs, double *e, double *dedn)
124 {
125 vdwxc_G_pw(0.031091, 0.21370, 7.5957, 3.5876, 1.6382, 0.49294,
126 rootrs, e, dedn);
127 }
128
vdwxc_G_zeta1(double rootrs,double * e,double * dedn)129 void vdwxc_G_zeta1(double rootrs, double *e, double *dedn)
130 {
131 vdwxc_G_pw(0.015545, 0.20548, 14.1189, 6.1977, 3.3662, 0.62517,
132 rootrs, e, dedn);
133 }
134
vdwxc_G_alphac(double rootrs,double * e,double * dedn)135 void vdwxc_G_alphac(double rootrs, double *e, double *dedn)
136 {
137 vdwxc_G_pw(0.016887, 0.11125, 10.357, 3.6231, 0.88026, 0.49671,
138 rootrs, e, dedn);
139 (*e) *= -1.0;
140 (*dedn) *= -1.0;
141 }
142
143
vdwxc_compute_lda(double n,double * e,double * v)144 void vdwxc_compute_lda(double n, double *e, double *v)
145 {
146 double rs = VDWXC_RS_PREFACTOR / pow(n, 1.0 / 3.0);
147 double rootrs = sqrt(rs);
148 vdwxc_G_zeta0(rootrs, e, v);
149 // e is now energy density epsilon, and v is derivative wrt. rs.
150 // convert to actual energy/particle and its derivative:
151 *v = (*e) - rs * (*v) / 3.0;
152 //*e *= n;
153 }
154
vdwxc_compute_lda_spin_zeta(double n,double zeta,double * e,double * v_up,double * v_dn)155 void vdwxc_compute_lda_spin_zeta(double n, double zeta,
156 double *e, double *v_up, double *v_dn)
157 {
158 //double C0I = 0.238732414637843;
159 //double C1 = -0.45816529328314287;
160 double CC1 = 1.9236610509315362;
161 double CC2 = 2.5648814012420482;
162 double IF2 = 0.58482236226346462;
163
164 double rs = VDWXC_RS_PREFACTOR / pow(n, 1.0 / 3.0);
165 double rootrs = sqrt(rs);
166
167 double e_zeta0, dedrs_zeta0;
168 vdwxc_G_zeta0(rootrs, &e_zeta0, &dedrs_zeta0);
169 double e_zeta1, dedrs_zeta1;
170 vdwxc_G_zeta1(rootrs, &e_zeta1, &dedrs_zeta1);
171 double alpha, dalphadrs;
172 vdwxc_G_alphac(rootrs, &alpha, &dalphadrs);
173
174 double zp = 1.0 + zeta;
175 double zm = 1.0 - zeta;
176 double xp = pow(zp, 1.0 / 3.0);
177 double xm = pow(zm, 1.0 / 3.0);
178 double f = CC1 * (zp * xp + zm * xm - 2.0);
179 double f1 = CC2 * (xp - xm);
180 double zeta3 = zeta * zeta * zeta;
181 double zeta4 = zeta3 * zeta;
182 double x = 1.0 - zeta4;
183
184 double decdrs = (dedrs_zeta0 * (1.0 - f * zeta4) +
185 dedrs_zeta1 * f * zeta4 +
186 dalphadrs * f * x * IF2);
187 //printf("decdrs %f\n", decdrs);
188 double decdzeta = (4.0 * zeta3 * f * (e_zeta1 - e_zeta0 - alpha * IF2) +
189 f1 * (zeta4 * e_zeta1 - zeta4 * e_zeta0
190 + x * alpha * IF2));
191 //printf("decdzeta %f\n", decdzeta);
192 //printf("f %f alpha %f %f\n", f, alpha, e1 - ec0);
193 *e = e_zeta0 + f * (alpha * IF2 * x + (e_zeta1 - e_zeta0) * zeta4);
194 *v_up = (*e) - rs * decdrs / 3.0 - (zeta - 1.0) * decdzeta;
195 *v_dn = (*e) - rs * decdrs / 3.0 - (zeta + 1.0) * decdzeta;
196 //*e *= n;
197 }
198
vdwxc_compute_lda_spin(double n_up,double n_dn,double * e,double * v_up,double * v_dn)199 void vdwxc_compute_lda_spin(double n_up, double n_dn,
200 double *e, double *v_up, double *v_dn)
201 {
202 double n = n_up + n_dn;
203 double zeta = (n_up - n_dn) / n;
204 vdwxc_compute_lda_spin_zeta(n, zeta, e, v_up, v_dn);
205 }
206
vdwxc_compute_q0x(double Z_ab,double rho,double sigma,double * q0x,double * dq0xdrho,double * dq0xdsigma)207 void vdwxc_compute_q0x(double Z_ab, double rho, double sigma, double *q0x,
208 double *dq0xdrho, double *dq0xdsigma)
209 {
210 double rho_root3 = pow(rho, 1.0 / 3.0);
211 double rho2 = rho * rho;
212 double kF = VDWXC_KF_PREFACTOR * rho_root3;
213 double _36_kF_rho2 = 36.0 * kF * rho2;
214
215 *q0x = kF - Z_ab * sigma / _36_kF_rho2;
216 *dq0xdrho = VDWXC_KF_PREFACTOR / (3.0 * rho_root3 * rho_root3) +
217 7.0 * Z_ab * sigma / (3.0 * _36_kF_rho2 * rho);
218 *dq0xdsigma = -Z_ab / _36_kF_rho2;
219 }
220
vdwxc_compute_q0x_spin(double Z_ab,double rho_up,double rho_dn,double sigma_up,double sigma_dn,double * q0xt,double * dq0xtdrho_up,double * dq0xtdrho_dn,double * dq0xtdsigma_up,double * dq0xtdsigma_dn)221 void vdwxc_compute_q0x_spin(double Z_ab, double rho_up, double rho_dn,
222 double sigma_up, double sigma_dn,
223 double *q0xt,
224 double *dq0xtdrho_up, double *dq0xtdrho_dn,
225 double *dq0xtdsigma_up, double *dq0xtdsigma_dn)
226 {
227
228 double rho = rho_up + rho_dn;
229
230 double q0x_2up, q0x_2dn;
231 double dq0xdrho_2up, dq0xdrho_2dn;
232 double dq0xdsigma_2up, dq0xdsigma_2dn;
233
234
235 if(rho_up > RHOCUT / 2) {
236 vdwxc_compute_q0x(Z_ab, 2.0 * rho_up, 4.0 * sigma_up,
237 &q0x_2up, &dq0xdrho_2up, &dq0xdsigma_2up);
238 } else {
239 q0x_2up = 0.0;
240 dq0xdrho_2up = 0.0;
241 dq0xdsigma_2up = 0.0;
242 }
243 if(rho_dn > RHOCUT / 2) {
244 vdwxc_compute_q0x(Z_ab, 2.0 * rho_dn, 4.0 * sigma_dn,
245 &q0x_2dn, &dq0xdrho_2dn, &dq0xdsigma_2dn);
246 } else {
247 q0x_2dn = 0.0;
248 dq0xdrho_2dn = 0.0;
249 dq0xdsigma_2dn = 0.0;
250 }
251
252 // QE 'saturates' the q0x variables as well for each spin.
253 // This code will enable
254 // the QE behaviour (for testing):
255 /* double qcut = 5.0;
256 double qcutbig = 4.0 * qcut;
257 double hderiv;
258
259 vdwxc_hfilter(q0x_2up, qcutbig, &q0x_2up, &hderiv);
260 dq0xdrho_2up *= hderiv;
261 dq0xdsigma_2up *= hderiv;
262
263 vdwxc_hfilter(q0x_2dn, qcutbig, &q0x_2dn, &hderiv);
264 dq0xdrho_2dn *= hderiv;
265 dq0xdsigma_2dn *= hderiv;*/
266 // done with silly extra cutoff
267
268 *q0xt = (rho_up * q0x_2up + rho_dn * q0x_2dn) / rho;
269 *dq0xtdrho_up = (rho_dn * (q0x_2up - q0x_2dn)
270 + 2.0 * rho_up * rho * dq0xdrho_2up) / (rho * rho);
271 *dq0xtdrho_dn = (rho_up * (q0x_2dn - q0x_2up)
272 + 2.0 * rho_dn * rho * dq0xdrho_2dn) / (rho * rho);
273 *dq0xtdsigma_up = 4.0 * rho_up / rho * dq0xdsigma_2up;
274 *dq0xtdsigma_dn = 4.0 * rho_dn / rho * dq0xdsigma_2dn;
275 }
276
277
vdwxc_hfilter(double q0,double qcut,double * hq0,double * dhq0dq0)278 void vdwxc_hfilter(double q0, double qcut, double *hq0, double *dhq0dq0)
279 {
280 //*dhq0dq0 = 999.;
281 //return; // XXXXXXXXXXXXXXXXXXX
282 int m;
283 double qrel = q0 / qcut;
284 double qrel_power = 1.0;
285 double exponent = 0.0, exponent_deriv = 0.0;
286 for (m=1; m<=12; m++) {
287 exponent_deriv += qrel_power;
288 qrel_power *= qrel;
289 exponent += qrel_power / m;
290 }
291 double exp_exponent = exp(-exponent);
292 *hq0 = qcut * (1.0 - exp_exponent);
293 *dhq0dq0 = exp_exponent * exponent_deriv;
294 }
295
296 #define LDA_VOLUME_FACTOR (-4.0 * VDWXC_PI / 3.0)
297
vdwxc_compute_q0_nospin(double Z_ab,double rho,double sigma,double * q0,double * dq0drho,double * dq0dsigma)298 void vdwxc_compute_q0_nospin(double Z_ab,
299 double rho,
300 double sigma,
301 double* q0,
302 double* dq0drho,
303 double* dq0dsigma)
304 {
305 double qcut = 5.0;
306 if(rho < RHOCUT) {
307 // Temp fix for small rho
308 *q0 = qcut;
309 *dq0drho = 0.0;
310 *dq0dsigma = 0.0;
311 // XXX update
312 return;
313 }
314
315 vdwxc_compute_q0x(Z_ab, rho, sigma, q0, dq0drho, dq0dsigma);
316 *dq0drho *= rho;
317 *dq0dsigma *= rho;
318
319 double e_lda;
320 double v_lda;
321 vdwxc_compute_lda(rho, &e_lda, &v_lda);
322 *q0 += LDA_VOLUME_FACTOR * e_lda;
323 // XXX slightly better to calculate and represent rho * dq0drho
324 *dq0drho += LDA_VOLUME_FACTOR * (v_lda - e_lda);// / rho;
325 double dhq0dq0;
326 vdwxc_hfilter(*q0, qcut, q0, &dhq0dq0);
327 *dq0drho *= dhq0dq0;
328 *dq0dsigma *= dhq0dq0;
329 }
330
vdwxc_compute_q0_spin(double Z_ab,double rho_up,double rho_dn,double sigma_up,double sigma_dn,double * q0,double * rho_dq0drho_up,double * rho_dq0drho_dn,double * rho_dq0dsigma_up,double * rho_dq0dsigma_dn)331 void vdwxc_compute_q0_spin(double Z_ab,
332 double rho_up, double rho_dn,
333 double sigma_up, double sigma_dn,
334 double* q0,
335 double* rho_dq0drho_up, double* rho_dq0drho_dn,
336 double* rho_dq0dsigma_up, double* rho_dq0dsigma_dn)
337 {
338 double qcut = 5.0;
339
340 if(rho_up + rho_dn < RHOCUT) {
341 *q0 = qcut;
342 *rho_dq0drho_up = 0.0;
343 *rho_dq0drho_dn = 0.0;
344 *rho_dq0drho_up = 0.0;
345 *rho_dq0drho_dn = 0.0;
346 return;
347 }
348
349 rho_up = fmax(rho_up, RHOCUT / 2.0);
350 rho_dn = fmax(rho_dn, RHOCUT / 2.0);
351 double rho = rho_up + rho_dn;
352
353 vdwxc_compute_q0x_spin(Z_ab, rho_up, rho_dn,
354 sigma_up, sigma_dn,
355 q0,
356 rho_dq0drho_up, rho_dq0drho_dn,
357 rho_dq0dsigma_up, rho_dq0dsigma_dn);
358 *rho_dq0drho_up *= rho;
359 *rho_dq0drho_dn *= rho;
360 *rho_dq0dsigma_up *= rho;
361 *rho_dq0dsigma_dn *= rho;
362
363 double e_lda;
364 double v_lda_up, v_lda_dn;
365 vdwxc_compute_lda_spin(rho_up, rho_dn,
366 &e_lda, &v_lda_up, &v_lda_dn);
367 *q0 += LDA_VOLUME_FACTOR * e_lda;
368 // XXX slightly better to calculate and represent rho * dq0drho
369 *rho_dq0drho_up += LDA_VOLUME_FACTOR * (v_lda_up - e_lda);
370 *rho_dq0drho_dn += LDA_VOLUME_FACTOR * (v_lda_dn - e_lda);
371
372 double dhq0dq0;
373 vdwxc_hfilter(*q0, qcut, q0, &dhq0dq0);
374 *rho_dq0drho_up *= dhq0dq0;
375 *rho_dq0drho_dn *= dhq0dq0;
376 *rho_dq0dsigma_up *= dhq0dq0;
377 *rho_dq0dsigma_dn *= dhq0dq0;
378 }
379
380
vdwxc_calculate_q0_nospin(int N,double Z_ab,double * rho_g,double * sigma_g,double * q0_g,double * rho_dq0drho_g,double * rho_dq0dsigma_g)381 void vdwxc_calculate_q0_nospin(int N, double Z_ab,
382 double *rho_g,
383 double *sigma_g,
384 double *q0_g,
385 double *rho_dq0drho_g,
386 double *rho_dq0dsigma_g)
387 {
388 int i;
389 for (i=0; i<N; i++) {
390 vdwxc_compute_q0_nospin(Z_ab, rho_g[i],
391 sigma_g[i],
392 &q0_g[i],
393 &rho_dq0drho_g[i],
394 &rho_dq0dsigma_g[i]);
395 }
396 }
397
398
399
vdwxc_calculate_q0_spin(int N,double Z_ab,double * rho_up_g,double * rho_dn_g,double * sigma_up_g,double * sigma_dn_g,double * q0_g,double * rho_dq0drho_sg,double * rho_dq0dsigma_sg)400 void vdwxc_calculate_q0_spin(int N, double Z_ab,
401 double *rho_up_g, double *rho_dn_g,
402 double *sigma_up_g, double *sigma_dn_g,
403 double *q0_g,
404 double *rho_dq0drho_sg,
405 double *rho_dq0dsigma_sg)
406 {
407 int i;
408 for (i=0; i<N; i++) {
409 vdwxc_compute_q0_spin(Z_ab, rho_up_g[i], rho_dn_g[i],
410 sigma_up_g[i], sigma_dn_g[i],
411 &q0_g[i],
412 &rho_dq0drho_sg[i], &rho_dq0drho_sg[i + N],
413 &rho_dq0dsigma_sg[i], &rho_dq0dsigma_sg[i + N]);
414 }
415 }
416