1 #include <cassert>
2 #include <cmath>
3 #include <cstdio>
4 #include <cstdlib>
5
6 #include "XCFun/xcfun.h"
7
8 void check(const char * what, int cond);
9 void checknum(const char * what,
10 double x,
11 double xref,
12 double abserr,
13 double relerr);
14 void consistency_test();
15 void gradient_forms_test();
16 void user_setup_test();
17 void xcfun_get_test();
18
19 /*
20 Run all tests for all functionals.
21 This program is in C to test the no stdc++ feature of xcfun.
22 */
23
check(const char * what,int cond)24 void check(const char * what, int cond) {
25 if (!cond) {
26 fprintf(stderr, "Failed check: %s\n", what);
27 exit(-1);
28 }
29 }
30
checknum(const char * what,double x,double xref,double abserr,double relerr)31 void checknum(const char * what,
32 double x,
33 double xref,
34 double abserr,
35 double relerr) {
36 if (fabs(x - xref) < abserr)
37 return;
38 if (fabs(x - xref) < relerr * xref)
39 return;
40 fprintf(stderr,
41 "Numerical check failed: %s\n got %.16e\n expected %.16e\n Absolute "
42 "error is %.4e, relative error is %.4e\n",
43 what,
44 x,
45 xref,
46 x - xref,
47 (x - xref) / xref);
48 exit(-1);
49 }
50
51 /* Test permutation symmetries over variables and modes etc. */
consistency_test()52 void consistency_test() {
53 auto fun = xcfun_new();
54 double d_unpolarized[8] = {1, 1, 2, -3, 4, 2, -3, 4};
55 double d_pol_a[8] = {1, 2.1, 2, -3, 4, 7, -8, 9};
56 double d_pol_b[8] = {2.1, 1, 7, -8, 9, 2, -3, 4};
57 int nout;
58 xcfun_set(fun, "pbe", 1.0);
59 xcfun_eval_setup(fun, XC_A_B_AX_AY_AZ_BX_BY_BZ, XC_PARTIAL_DERIVATIVES, 1);
60 nout = xcfun_output_length(fun);
61 check("correct output length 1", nout == 9); // 1 + 8
62 auto output = new double[nout];
63 auto out2 = new double[nout];
64 xcfun_eval(fun, d_unpolarized, output);
65
66 checknum("unpolarized symmetry 1", output[1] - output[2], 0, 1e-14, 1e-12);
67 checknum("unpolarized symmetry 2", output[3] - output[6], 0, 1e-14, 1e-12);
68 checknum("unpolarized symmetry 3", output[4] - output[7], 0, 1e-14, 1e-12);
69 checknum("unpolarized symmetry 4", output[5] - output[8], 0, 1e-14, 1e-12);
70 xcfun_eval(fun, d_pol_a, output);
71 xcfun_eval(fun, d_pol_b, out2);
72 checknum("polarized symmetry 1", output[1] - out2[2], 0, 1e-14, 1e-12);
73 checknum("polarized symmetry 2", output[3] - out2[6], 0, 1e-14, 1e-12);
74 checknum("polarized symmetry 3", output[4] - out2[7], 0, 1e-14, 1e-12);
75 checknum("polarized symmetry 4", output[5] - out2[8], 0, 1e-14, 1e-12);
76 checknum("polarized symmetry 5", out2[1] - output[2], 0, 1e-14, 1e-12);
77 checknum("polarized symmetry 6", out2[3] - output[6], 0, 1e-14, 1e-12);
78 checknum("polarized symmetry 7", out2[4] - output[7], 0, 1e-14, 1e-12);
79 checknum("polarized symmetry 8", out2[5] - output[8], 0, 1e-14, 1e-12);
80 delete[] output;
81 delete[] out2;
82 xcfun_delete(fun);
83 }
84
85 // Test that gradient square norma and gradient elements modes are consistent
gradient_forms_test()86 void gradient_forms_test() {
87 auto fun = xcfun_new();
88 double d_elements[8] = {1, 2.1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6};
89 double d_sqnorm[5] = {d_elements[0], d_elements[1]};
90 int nout, i;
91 xcfun_set(fun, "blyp", 1.0);
92
93 xcfun_eval_setup(fun, XC_A_B_AX_AY_AZ_BX_BY_BZ, XC_PARTIAL_DERIVATIVES, 1);
94 nout = xcfun_output_length(fun);
95
96 check("correct output length 1", nout == 9); // 1 + 8
97
98 auto output = new double[nout];
99 xcfun_eval(fun, d_elements, output);
100
101 xcfun_eval_setup(fun, XC_A_B_GAA_GAB_GBB, XC_PARTIAL_DERIVATIVES, 1);
102 nout = xcfun_output_length(fun);
103
104 check("correct output length 1", nout == 6); // 1 + 5
105 auto out2 = new double[nout];
106 d_sqnorm[2] = d_elements[2] * d_elements[2] + d_elements[3] * d_elements[3] +
107 d_elements[4] * d_elements[4];
108 d_sqnorm[3] = d_elements[2] * d_elements[5] + d_elements[3] * d_elements[6] +
109 d_elements[4] * d_elements[7];
110 d_sqnorm[4] = d_elements[5] * d_elements[5] + d_elements[6] * d_elements[6] +
111 d_elements[7] * d_elements[7];
112 xcfun_eval(fun, d_sqnorm, out2);
113
114 checknum("Grad modes energy", output[0] - out2[0], 0, 1e-14, 1e-12);
115 checknum("Grad modes density derivs alpha", output[1] - out2[1], 0, 1e-14, 1e-12);
116 checknum("Grad modes density derivs beta", output[2] - out2[2], 0, 1e-14, 1e-12);
117 // d/dg_ax = d/dgaa * g_ax + d/dgab * g_bx
118 for (i = 0; i < 3; i++) {
119 checknum("Grad modes density grad alpha",
120 output[3 + i] -
121 (2 * out2[3] * d_elements[2 + i] + out2[4] * d_elements[5 + i]),
122 0,
123 1e-14,
124 1e-12);
125 checknum("Grad modes density grad beta",
126 output[6 + i] -
127 (2 * out2[5] * d_elements[5 + i] + out2[4] * d_elements[2 + i]),
128 0,
129 1e-14,
130 1e-12);
131 }
132
133 delete[] output;
134 delete[] out2;
135 xcfun_delete(fun);
136 }
137
user_setup_test()138 void user_setup_test() {
139 auto fun1 = xcfun_new();
140 auto fun2 = xcfun_new();
141 auto fun3 = xcfun_new();
142 xcfun_set(fun1, "lda", 1.0);
143 xcfun_set(fun2, "pbe", 1.0);
144 xcfun_set(fun3, "m06l", 1.0);
145 int rval1 = xcfun_user_eval_setup(fun1, 0, 0, 0, 1, 0, 0, 0, 0);
146 int rval2 = xcfun_user_eval_setup(fun2, 1, 1, 1, 1, 0, 0, 0, 1);
147 int rval3 = xcfun_user_eval_setup(fun3, 2, 2, 2, 1, 0, 1, 0, 1);
148 check("Functional 1 correctly set up", rval1 == 0);
149 check("Functional 2 correctly set up", rval2 == 0);
150 check("Functional 3 correctly set up", rval3 == 0);
151 xcfun_delete(fun1);
152 xcfun_delete(fun2);
153 xcfun_delete(fun3);
154 }
155
xcfun_get_test()156 void xcfun_get_test() {
157 auto fun = xcfun_new();
158 xcfun_set(fun, "B3LYP", 1.0);
159
160 double s, b, lyp, vwn, exx, kt, foo;
161 check("SLATERX is a valid functional", xcfun_get(fun, "SLATERX", &s) == 0);
162 check("BECKECORRX is a valid functional", xcfun_get(fun, "BECKECORRX", &b) == 0);
163 check("LYPC is a valid functional", xcfun_get(fun, "LYPC", &lyp) == 0);
164 check("VWN5C is a valid functional", xcfun_get(fun, "VWN5C", &vwn) == 0);
165 check("EXX is a valid functional", xcfun_get(fun, "EXX", &exx) == 0);
166 check("KTX is a valid functional", xcfun_get(fun, "KTX", &kt) == 0);
167 check("FOO is NOT a valid functional", xcfun_get(fun, "FOO", &foo) != 0);
168
169 checknum("B3LYP contains 80% SLATERX", s, 0.80, 1.0e-14, 1.0e-12);
170 checknum("B3LYP contains 72% BECKECORRX", b, 0.72, 1.0e-14, 1.0e-12);
171 checknum("B3LYP contains 81% LYPC", lyp, 0.81, 1.0e-14, 1.0e-12);
172 checknum("B3LYP contains 19% VWN5C", vwn, 0.19, 1.0e-14, 1.0e-12);
173 checknum("B3LYP contains 20% EXX", exx, 0.20, 1.0e-14, 1.0e-12);
174 checknum("B3LYP contains 0% KTX", kt, 0.00, 1.0e-14, 1.0);
175
176 xcfun_delete(fun);
177 }
178
main()179 int main() {
180 int i = 0;
181 const char *n, *s;
182 consistency_test();
183 gradient_forms_test();
184 user_setup_test();
185 xcfun_get_test();
186 printf("%s", xcfun_splash());
187 printf("XCFun version: %s\n", xcfun_version());
188 printf("\nAvailable functionals and other settings:\n");
189 while ((n = xcfun_enumerate_parameters(i++))) {
190 printf("%s \t", n);
191 if ((s = xcfun_describe_short(n)))
192 printf("%s", s);
193 else
194 printf("[No description]");
195 printf("\n");
196 }
197 printf("\nAvailable aliases:\n");
198 i = 0;
199 while ((n = xcfun_enumerate_aliases(i++))) {
200 printf("%s \t", n);
201 if ((s = xcfun_describe_short(n)))
202 printf("%s", s);
203 else
204 printf("[No description]");
205 printf("\n");
206 }
207 printf("\nRunning tests..\n");
208 if (xcfun_test() == 0) {
209 printf("\nAll tests ok\n");
210 return 0;
211 } else {
212 printf("\nSome tests failed\n");
213 return -1;
214 }
215 }
216