1 #include <stdio.h>
2 #include <assert.h>
3 #include <math.h>
4 #include "vdw_q0.h"
5
main(int argc,char ** argv)6 int main(int argc, char** argv)
7 {
8 double rs = 0.5;
9 double e, dedrs;
10 vdwxc_G_zeta0(rs, &e, &dedrs);
11 double e1, e2, tmp;
12 double drs = 1e-8;
13 vdwxc_G_zeta0(rs - drs, &e1, &tmp);
14 vdwxc_G_zeta0(rs + drs, &e2, &tmp);
15 double deriv = 0.5 * (e2 - e1) / drs;
16 double deriv_err = fabs(deriv - dedrs);
17 printf("G_pw e=%.16f dedrho=%f deriv_err=%e\n\n", e, dedrs, deriv_err);
18 assert(fabs(-0.0949608976455678 - e) < 1e-15);
19 assert(deriv_err < 1e-9);
20
21 double rho_up = 0.85, rho_dn = 0.35;
22 double rho = rho_up + rho_dn;
23 double zeta = (rho_up - rho_dn) / rho;
24
25 double eps, v;
26 vdwxc_compute_lda(rho, &eps, &v);
27 eps *= rho;
28 double ref = -0.0872542585762377;
29 double err = fabs(ref - eps);
30
31 double drho = 1e-8;
32 vdwxc_compute_lda(rho - drho, &e1, &tmp);
33 vdwxc_compute_lda(rho + drho, &e2, &tmp);
34 e1 *= rho - drho;
35 e2 *= rho + drho;
36 double lda_deriv = 0.5 * (e2 - e1) / drho;
37 //printf("lda spinpair deriv
38 double lda_deriv_err = fabs(lda_deriv - v);
39 printf("lda spinpair\neps=%.16f err=%e\nv=%.11f deriv=%.11f deriv_err=%e\n\n",
40 eps, err, v, lda_deriv, lda_deriv_err);
41 assert(err < 1e-15);
42 assert(lda_deriv_err < 1e-8);
43
44 //double zeta = 0.7;
45 double eps_spin;
46 double v_up, v_dn;
47 double eps_spin_ref = -0.0821360927715360;
48 vdwxc_compute_lda_spin_zeta(rho, zeta, &eps_spin, &v_up, &v_dn);
49 eps_spin *= rho;
50 //eps_spin *= rho;
51 double e_err = fabs(eps_spin_ref - eps_spin);
52 printf("lda spinzeta\neps=%.11f err=%e v_up=%.11f v_dn=%.11f\n\n",
53 eps_spin, e_err, v_up, v_dn);
54 assert(e_err < 1e-15);
55
56 vdwxc_compute_lda_spin(rho_up, rho_dn, &eps_spin, &v_up, &v_dn);
57 eps_spin *= rho;
58
59 double espin1, espin2;
60 vdwxc_compute_lda_spin(rho_up - drho, rho_dn, &espin1, &tmp, &tmp);
61 vdwxc_compute_lda_spin(rho_up + drho, rho_dn, &espin2, &tmp, &tmp);
62 espin1 *= rho - drho;
63 espin2 *= rho + drho;
64
65 double deriv_spin_up = 0.5 * (espin2 - espin1) / drho;
66 double deriv_err_up = fabs(deriv_spin_up - v_up);
67 double err_spin2 = fabs(eps_spin_ref - eps_spin);
68 printf("lda spin\n");
69 printf("eps=%.16f err=%e\n", eps_spin, err_spin2);
70 printf("v_up=%.16f v_dn=%.16f\n", v_up, v_dn);
71 assert(err_spin2 < 1e-15);
72 vdwxc_compute_lda_spin(rho_up, rho_dn - drho, &espin1, &tmp, &tmp);
73 vdwxc_compute_lda_spin(rho_up, rho_dn + drho, &espin2, &tmp, &tmp);
74 espin1 *= rho - drho;
75 espin2 *= rho + drho;
76
77 double deriv_spin_dn = 0.5 * (espin2 - espin1) / drho;
78 double deriv_err_dn = fabs(deriv_spin_dn - v_dn);
79 printf("err_up=%e err_dn=%e\n\n", deriv_err_up, deriv_err_dn);
80 assert(deriv_err_up < 1e-8);
81 assert(deriv_err_dn < 1e-8);
82
83 double Z_ab = -75.0;
84 double sigma = 0.3;
85 double dsigma = drho;
86 double q0, q0_1, q0_2, dq0drho, dq0dsigma;
87 double q0x_ref = 3.419535675401194;
88 vdwxc_compute_q0x(Z_ab, rho, sigma, &q0, &dq0drho, &dq0dsigma);
89
90 vdwxc_compute_q0x(Z_ab, rho - drho, sigma, &q0_1, &tmp, &tmp);
91 vdwxc_compute_q0x(Z_ab, rho + drho, sigma, &q0_2, &tmp, &tmp);
92 double dq0drho_numeric = 0.5 * (q0_2 - q0_1) / drho;
93 double dq0drho_err = fabs(dq0drho_numeric - dq0drho);
94
95 vdwxc_compute_q0x(Z_ab, rho, sigma - dsigma, &q0_1, &tmp, &tmp);
96 vdwxc_compute_q0x(Z_ab, rho, sigma + dsigma, &q0_2, &tmp, &tmp);
97 double dq0dsigma_numeric = 0.5 * (q0_2 - q0_1) / dsigma;
98 double dq0dsigma_err = fabs(dq0dsigma_numeric - dq0dsigma);
99
100 double q0x_err = fabs(q0x_ref - q0);
101 printf("q0x %.15f err %e\n", q0, q0x_err);
102 printf("dq0drho=%f err=%e\n", dq0drho, dq0drho_err);
103 printf("dq0dsigma=%f err=%e\n\n", dq0dsigma, dq0dsigma_err);
104 assert(q0x_err < 1e-14);
105 assert(dq0drho_err < 5e-8);
106 assert(dq0dsigma_err < 2e-8);
107 //double ref2 = -0.0650794082170656;
108 //double err2 = fabs(ref2 - eps2);
109 //printf("eps2=%.16f, err=%e\n", eps2, err2);
110 //assert(err2 < 1e-15);
111
112 //double rho_up = 2.2, n_dn = 0.9;
113 double sigma_up = 1.5, sigma_dn = 0.75;
114 double q0xt, dq0xt_drho_up, dq0xt_drho_dn,
115 dq0xt_dsigma_up, dq0xt_dsigma_dn;
116 Z_ab = -25.0;
117 vdwxc_compute_q0x_spin(Z_ab, rho_up, rho_dn, sigma_up, sigma_dn,
118 &q0xt, &dq0xt_drho_up, &dq0xt_drho_dn,
119 &dq0xt_dsigma_up, &dq0xt_dsigma_dn);
120 double q0xt_ref = 4.144551809731131;
121 double q0xt_err = fabs(q0xt - q0xt_ref);
122 printf("q0x_spin %.15f err=%e\n", q0xt, q0xt_err);
123
124 double q0xt_1, q0xt_2;
125 vdwxc_compute_q0x_spin(Z_ab, rho_up - drho, rho_dn, sigma_up, sigma_dn,
126 &q0xt_1, &tmp, &tmp, &tmp, &tmp);
127 vdwxc_compute_q0x_spin(Z_ab, rho_up + drho, rho_dn, sigma_up, sigma_dn,
128 &q0xt_2, &tmp, &tmp, &tmp, &tmp);
129 double dq0xt_drho_up_numeric = 0.5 * (q0xt_2 - q0xt_1) / drho;
130 err = fabs(dq0xt_drho_up_numeric - dq0xt_drho_up);
131 printf("dq0xtdrho_up=%f, err=%e\n", dq0xt_drho_up, err);
132 assert(q0xt_err < 1e-14);
133 assert(err < 1e-7);
134
135 vdwxc_compute_q0x_spin(Z_ab, rho_up, rho_dn - drho, sigma_up, sigma_dn,
136 &q0xt_1, &tmp, &tmp, &tmp, &tmp);
137 vdwxc_compute_q0x_spin(Z_ab, rho_up, rho_dn + drho, sigma_up, sigma_dn,
138 &q0xt_2, &tmp, &tmp, &tmp, &tmp);
139 double dq0xt_drho_dn_numeric = 0.5 * (q0xt_2 - q0xt_1) / drho;
140 err = fabs(dq0xt_drho_dn_numeric - dq0xt_drho_dn);
141 printf("dq0xtdrho_dn=%f, err=%e\n", dq0xt_drho_dn, err);
142 assert(err < 1e-7);
143
144 vdwxc_compute_q0x_spin(Z_ab, rho_up, rho_dn, sigma_up - dsigma, sigma_dn,
145 &q0xt_1, &tmp, &tmp, &tmp, &tmp);
146 vdwxc_compute_q0x_spin(Z_ab, rho_up, rho_dn, sigma_up + dsigma, sigma_dn,
147 &q0xt_2, &tmp, &tmp, &tmp, &tmp);
148 double dq0xt_dsigma_up_numeric = 0.5 * (q0xt_2 - q0xt_1) / dsigma;
149 err = fabs(dq0xt_dsigma_up_numeric - dq0xt_dsigma_up);
150 printf("dq0xt_dsigma_up %f, err=%e\n", dq0xt_dsigma_up, err);
151 assert(err < 5e-8);
152
153 vdwxc_compute_q0x_spin(Z_ab, rho_up, rho_dn, sigma_up, sigma_dn - dsigma,
154 &q0xt_1, &tmp, &tmp, &tmp, &tmp);
155 vdwxc_compute_q0x_spin(Z_ab, rho_up, rho_dn, sigma_up, sigma_dn + dsigma,
156 &q0xt_2, &tmp, &tmp, &tmp, &tmp);
157 double dq0xt_dsigma_dn_numeric = 0.5 * (q0xt_2 - q0xt_1) / dsigma;
158 err = fabs(dq0xt_dsigma_dn_numeric - dq0xt_dsigma_dn);
159 printf("dq0xt_dsigma_dn %f, err=%e\n\n", dq0xt_dsigma_dn, err);
160 assert(err < 1e-7);
161
162 double oldq0, olddq0_drho, olddq0_dsigma;
163 double newq0, newdq0_drho, newdq0_dsigma;
164 vdwxc_calculate_q0(1, Z_ab, 5.0, &rho, &sigma,
165 &oldq0, &olddq0_drho, &olddq0_dsigma);
166 vdwxc_calculate_q0_nospin(1, Z_ab, &rho, &sigma,
167 &newq0, &newdq0_drho, &newdq0_dsigma);
168
169 double dq0drho_up_spinpair, dq0drho_dn_spinpair;
170 double dq0dsigma_up_spinpair, dq0dsigma_dn_spinpair;
171 double q0_spinpair;
172
173 vdwxc_compute_q0_spin(Z_ab, rho / 2., rho / 2.,
174 sigma / 4., sigma / 4.,
175 &q0_spinpair,
176 &dq0drho_up_spinpair, &dq0drho_dn_spinpair,
177 &dq0dsigma_up_spinpair, &dq0dsigma_dn_spinpair);
178
179 printf("q0 spin-paired old/new comparison\n");
180 printf("q0 old %.15e\n", oldq0);
181 printf("q0 new %.15e\n", newq0);
182 printf("q0 spp %.15e\n", q0_spinpair);
183 printf("old dq0dn %.15e\n", olddq0_drho);
184 printf("new dq0dn %.15e\n", newdq0_drho);
185 printf("spp dq0dn %.15e %.15e\n", dq0drho_up_spinpair, dq0drho_dn_spinpair);
186 printf("old dq0dsigma %.15e\n", olddq0_dsigma);
187 printf("new dq0dsigma %.15e\n", newdq0_dsigma);
188 printf("spp dq0dsigma %.15e %.15e\n\n", dq0dsigma_up_spinpair, dq0dsigma_dn_spinpair);
189
190 // And now the big one!
191 double dq0drho_up, dq0drho_dn;
192 double dq0dsigma_up, dq0dsigma_dn;
193 double q0ref = 4.372399873725762;
194 vdwxc_compute_q0_spin(Z_ab,
195 rho_up, rho_dn,
196 sigma_up, sigma_dn,
197 &q0,
198 &dq0drho_up, &dq0drho_dn,
199 &dq0dsigma_up, &dq0dsigma_dn);
200 dq0drho_up /= rho;
201 dq0drho_dn /= rho;
202 dq0dsigma_up /= rho;
203 dq0dsigma_dn /= rho;
204 double q0err = fabs(q0ref - q0);
205
206 vdwxc_compute_q0_spin(Z_ab, rho_up - drho, rho_dn, sigma_up, sigma_dn,
207 &q0_1, &tmp, &tmp, &tmp, &tmp);
208 vdwxc_compute_q0_spin(Z_ab, rho_up + drho, rho_dn, sigma_up, sigma_dn,
209 &q0_2, &tmp, &tmp, &tmp, &tmp);
210 err = fabs(0.5 * (q0_2 - q0_1) / drho - dq0drho_up);
211 printf("q0_spin\nq0 %.15f err %e\n", q0, q0err);
212 printf("dq0drho_up=%f err=%e\n", dq0drho_up, err);
213 assert(q0err < 1e-14);
214 assert(err < 5e-8);
215
216 vdwxc_compute_q0_spin(Z_ab, rho_up, rho_dn - drho, sigma_up, sigma_dn,
217 &q0_1, &tmp, &tmp, &tmp, &tmp);
218 vdwxc_compute_q0_spin(Z_ab, rho_up, rho_dn + drho, sigma_up, sigma_dn,
219 &q0_2, &tmp, &tmp, &tmp, &tmp);
220 err = fabs(0.5 * (q0_2 - q0_1) / drho - dq0drho_dn);
221 printf("dq0drho_dn=%f err=%e\n", dq0drho_dn, err);
222 assert(err < 1e-7);
223
224 vdwxc_compute_q0_spin(Z_ab, rho_up, rho_dn, sigma_up - dsigma, sigma_dn,
225 &q0_1, &tmp, &tmp, &tmp, &tmp);
226 vdwxc_compute_q0_spin(Z_ab, rho_up, rho_dn, sigma_up + dsigma, sigma_dn,
227 &q0_2, &tmp, &tmp, &tmp, &tmp);
228 err = fabs(0.5 * (q0_2 - q0_1) / dsigma - dq0dsigma_up);
229 printf("dq0dsigma_up=%f err=%e\n", dq0dsigma_up, err);
230 assert(err < 1e-7);
231
232 vdwxc_compute_q0_spin(Z_ab, rho_up, rho_dn, sigma_up, sigma_dn - dsigma,
233 &q0_1, &tmp, &tmp, &tmp, &tmp);
234 vdwxc_compute_q0_spin(Z_ab, rho_up, rho_dn, sigma_up, sigma_dn + dsigma,
235 &q0_2, &tmp, &tmp, &tmp, &tmp);
236 err = fabs(0.5 * (q0_2 - q0_1) / dsigma - dq0dsigma_dn);
237 printf("dq0dsigma_dn=%f err=%e\n", dq0dsigma_dn, err);
238 assert(err < 1e-7);
239
240 double test_q0 = 4.4;
241 double test_dq0 = 1e-8;
242 double hq0, hq0_1, hq0_2;
243 double dhq0dq;
244 double qcut = 5.0;
245 vdwxc_hfilter(test_q0, qcut, &hq0, &dhq0dq);
246 vdwxc_hfilter(test_q0 - test_dq0, qcut, &hq0_1, &tmp);
247 vdwxc_hfilter(test_q0 + test_dq0, qcut, &hq0_2, &tmp);
248 double dhq0dq_numerical = 0.5 * (hq0_2 - hq0_1) / test_dq0;
249 double dhq0dq_err = fabs(dhq0dq_numerical - dhq0dq);
250 printf("hq0 %f, dhq0dq %f, numeric %f, err %e\n",
251 hq0, dhq0dq, dhq0dq_numerical, dhq0dq_err);
252 assert(dhq0dq_err < 2e-8);
253 return 0;
254 }
255