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