1 #include <math.h>
2
3 #include "vdwxc.h"
4
5 #if defined HAVE_CONFIG_H
6 #include "config.h"
7 #endif
8
9 #include "vdw_core.h"
10
11
12 /* We need to implement:
13
14 * PBE R exchange
15 * GGA-RPW86 exchange
16 * Langreth-Vosko thing with interpolation to GGA-RPW86-X
17 * LDA-PW correlation
18
19 */
20
vdwxc_lda_correlation(int Ng,double * rho_g,double * eps_g,double * dedn_g)21 void vdwxc_lda_correlation(int Ng, double* rho_g, double* eps_g, double* dedn_g)
22 {
23 double
24 C0I = 0.238732414637843,
25 gamma = 0.031091,
26 alpha1 = 0.21370,
27 beta1 = 7.5957,
28 beta2 = 3.5876,
29 beta3 = 1.6382,
30 beta4 = 0.49294; // defined also in vdwxc.h... use them from there.
31
32 int i;
33 for(i=0; i < Ng; i++) {
34 double rtrs = pow(C0I / (*rho_g), 1.0 / 6.0);
35 double rs = rtrs * rtrs;
36 double Q0 = -2.0 * gamma * (1.0 + alpha1 * rtrs * rtrs);
37 double Q1 = 2.0 * gamma * rtrs *
38 (beta1 + rtrs * (beta2 + rtrs * (beta3 + rtrs * beta4)));
39 double G1 = Q0 * log(1.0 + 1.0 / Q1);
40
41 double dQ1drs = gamma * (beta1 / rtrs + 2.0 * beta2 +
42 rtrs * (3.0 * beta3 + 4.0 * beta4 * rtrs));
43 double dGdrs = -2.0 * gamma * alpha1 * G1 / Q0 - Q0
44 * dQ1drs / (Q1 * (Q1 + 1.0));
45
46 //double ec = G1;
47 //double decdrs_0 = dGdrs;
48 *eps_g += (*rho_g) * G1;
49 *dedn_g += G1 - rs * dGdrs / 3.0;
50
51 rho_g++;
52 eps_g++;
53 dedn_g++;
54 }
55 }
56
57
58 // sigma = || grad rho ||^2
vdwxc_cx_semilocal_exchange(int Ng,double * rho_g,double * sigma_g,double * eps_g,double * dedn_g,double * dedsigma_g)59 void vdwxc_cx_semilocal_exchange(int Ng, double* rho_g, double* sigma_g,
60 double* eps_g, double* dedn_g, double* dedsigma_g)
61 {
62 //double tol = 1e-20;
63
64 double alp = 0.021789;
65 double beta = 1.15;
66 double a = 1.851;
67 double b = 17.33;
68 double c = 0.163;
69 double mu_LM = 0.09434;
70 double s_prefactor = 6.18733545256027;
71 double Ax = -0.738558766382022; // = -3./4. * (3./pi)**(1./3)
72
73 int i;
74 for(i=0; i < Ng; i++) {
75 double grad_rho = sqrt(*sigma_g);
76
77 // eventually we need s to power 12. Clip to avoid overflow
78 // (We have individual tolerances on both rho and grho, but
79 // they are not sufficient to guarantee this)
80 double rho_pow_1_3 = pow(*rho_g, 1.0 / 3.0);
81 double rho_pow_4_3 = (*rho_g) * rho_pow_1_3;
82 double s_1 = grad_rho / (s_prefactor * rho_pow_4_3);
83 s_1 = s_1 < 0.0 ? 0.0 : (s_1 > 1e20 ? 1e20 : s_1);
84 double s_2 = s_1 * s_1;
85 double s_3 = s_2 * s_1;
86 double s_4 = s_3 * s_1;
87 double s_5 = s_4 * s_1;
88 double s_6 = s_5 * s_1;
89
90 double fs_rPW86 = pow(1.0 + a * s_2 + b * s_4 + c * s_6, 1. / 15.);
91
92 double fs = (1.0 + mu_LM * s_2) / (1.0 + alp * s_6)
93 + alp * s_6 / (beta + alp * s_6) * fs_rPW86;
94
95 // the energy density for the exchange.
96 *eps_g += Ax * rho_pow_4_3 * fs;
97
98 double df_rPW86_ds = (1. / (15. * pow(fs_rPW86, 14.0))) *
99 (2.0 * a * s_1 + 4.0 * b * s_3 + 6.0 * c * s_5);
100
101 double tmp = 1. + alp * s_6;
102 double df_ds = 1. / (tmp * tmp)
103 * (2.0 * mu_LM * s_1 * (1. + alp * s_6)
104 - 6.0 * alp * s_5 * (1. + mu_LM * s_2))
105 + alp * s_6 / (beta + alp * s_6) * df_rPW86_ds
106 + 6.0 * alp * s_5 * fs_rPW86 / (beta + alp * s_6)
107 * (1. - alp * s_6 / (beta + alp * s_6));
108
109 // de/dn. This is the partial derivative of sx wrt. n, for s constant
110 *dedn_g += Ax * (4.0 / 3.0) *
111 (rho_pow_1_3 * fs - grad_rho * df_ds / (s_prefactor * (*rho_g)));
112 // de/d(nabla n). The other partial derivative
113 *dedsigma_g += 0.5 * Ax * df_ds / (s_prefactor * grad_rho);
114 // (We may or may not understand what that grad_rho is doing here.)
115 rho_g++;
116 sigma_g++;
117 eps_g++;
118 dedn_g++;
119 dedsigma_g++;
120 }
121 }
122
vdwxc_add_cx_exchange(vdwxc_data data,double * rho_g,double * sigma_g,double * eps_g,double * dedn_g,double * dedsigma_g)123 void vdwxc_add_cx_exchange(vdwxc_data data, double* rho_g, double* sigma_g,
124 double* eps_g, double* dedn_g, double* dedsigma_g)
125 {
126 vdwxc_cx_semilocal_exchange(data->Ng, rho_g, sigma_g, eps_g, dedn_g, dedsigma_g);
127 }
128