1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
31 /** @file fun-tester.cc
32 Program for testing functional routines in the DFT module.
33 (c) Pawel Salek, pawsa@theochem.kth.se, 2001-10-15
34
35 The test build can be done by:
36 g77 -O fun-tester.c -o fun-tester -L. -ldft -lm
37 or
38 cc -O fun-tester.c -o fun-tester -L. -ldft -lm -lg2c
39
40 NOTES: this file is short but in a separate file to reduce the
41 number of dependences and be able to easily compile the code for
42 the TEST_BUILD.
43 */
44
45 //#define _POSIX_SOURCE 1
46
47 #define __CVERSION__
48 #include <cmath>
49 #include <stdarg.h>
50 #include <stdio.h>
51 #include <stdlib.h>
52 #include <string.h>
53 #include <strings.h>
54
55 #include "dft_common.h"
56 #include "functionals.h"
57
58 #if !defined __inline__
59 /* inline some stuff whenever possible */
60 #define __inline__
61 #endif
62
63 typedef void (*DaltonEnFunc)(real* res, const real* rho,
64 const real* rho13, const real*grad);
65 void edrc_(real* drc, const real* rho, const real* rho13, const real* rhogrd);
66 void evwn_(real* vwn, const real* rho, const real* rho13, const real* rhogrd);
67 void ebck_(real* bck, const real* rho, const real* rho13, const real* rhogrd);
68 void elyp_(real* lyp, const real* rho, const real* rho13, const real* rhogrd);
69
70 void dftpot1(SecondDrv *ds, const real *w, const real* rho, const real* grad,
71 const int* triplet);
72 void condft_(void);
73
74 static __inline__
test_var(real comp,real refer,const char * fun,const char * drv,int * counter)75 void test_var(real comp, real refer, const char* fun, const char* drv,
76 int* counter)
77 {
78 if(template_blas_fabs(comp-refer)>2e-7+5e-5*(template_blas_fabs(comp)+template_blas_fabs(refer))) {
79 if(*counter<95)
80 printf("%s %s: fin.diff: %12g found: %12g, diff=%g\n",
81 fun, drv, (double)refer, (double)comp, (double)template_blas_fabs(comp-refer));
82 ++*counter;
83 }/* else printf("Test '%s:%s' passed (expected: %g found: %g).\n",
84 fun,drv, refer, comp); */
85 }
86
87 const int GRID_STEP = 1;
88
89 /* test_first: test first order derivatives of given functional.
90 Note that really the restricted case is only tested...
91 */
92 static int
test_first(const char * fun_name,EnergyFunc func,FirstOrderFun first_func)93 test_first(const char* fun_name, EnergyFunc func, FirstOrderFun first_func)
94 {
95 int i, j, k, failed = 0;
96 real drho, dgra, resp, resm, num;
97 FunFirstFuncDrv gga;
98 for(i=1; i<=40; i+=GRID_STEP) {
99 for(j=1; j<=40; j+=GRID_STEP) {
100 for(k=-19; k<=19; k+=GRID_STEP) {
101 real rho = i/40.0;
102 real ngrad = j/40.0;
103 real gracos= k/20.0;
104 FunDensProp dt, dp = { (real)0.5*rho,(real)0.2*rho, (real)0.2*ngrad,(real)0.2*ngrad };
105 dp.gradab = dp.grada*dp.gradb*gracos;
106 /* TEST df1000 */
107 drho = rho*1e-4;
108 dt = dp; dt.rhoa -= drho; resm = func(&dt);
109 dt = dp; dt.rhoa += drho; resp = func(&dt);
110 drv1_clear(&gga);
111 first_func(&gga, 1, &dp);
112 num = (resp-resm)/(2*drho);
113 test_var(gga.df1000,num, fun_name, "df1000", &failed);
114
115 /* TEST df0010 */
116 dgra = ngrad*1e-7;
117 dt = dp; dt.grada -= dgra; resm = func(&dt);
118 dt = dp; dt.grada += dgra; resp = func(&dt);
119 num = (resp-resm)/(2*dgra);
120 test_var(gga.df0010,num, fun_name, "df0010", &failed);
121
122 /* TEST df00001 */
123 if(template_blas_fabs(gracos)<1e-5) continue;
124 dgra = gracos*1e-7;
125 dt = dp; dt.gradab -= dgra; resm = func(&dt);
126 dt = dp; dt.gradab += dgra; resp = func(&dt);
127 num = (resp-resm)/(2*dgra);
128 test_var(gga.df00001,num, fun_name, "df00001", &failed);
129 }
130 }
131 }
132 if(failed==0) printf("%-5s (first order derivatives): OK\n", fun_name);
133 return failed;
134 }
135
136 /* test_second:
137 test second order derivatives of given functional.
138 It is assumed that the first order derivatives are OK.
139 */
140 #define COMP_DER(c,delta,eps,field, fun, ord) \
141 eps = c*delta; \
142 drv##ord##_clear(&m); dt = dp; dt.field -= eps; fun ## _fun(&m, 1, &dt); \
143 drv##ord##_clear(&p); dt = dp; dt.field += eps; fun ## _fun(&p, 1, &dt)
144
145 #define T2(der,derdif,eps,label) \
146 num = (p.derdif-m.derdif)/(2*eps);\
147 test_var(d.der, num, fname, label ":" #der, &fail)
148
149 static int
test_second(const char * fname,FirstOrderFun first_fun,SecondOrderFun second_fun)150 test_second(const char* fname,
151 FirstOrderFun first_fun, SecondOrderFun second_fun)
152 {
153 int i, j, k, fail = 0;
154 real drho, dgra, num;
155 FunFirstFuncDrv m, p;
156 FunSecondFuncDrv d;
157 for(i=1; i<=40; i+=GRID_STEP) {
158 for(j=1; j<=40; j+=GRID_STEP) {
159 for(k=-19; k<=19; k+=GRID_STEP) {
160 real rho = i/40.0;
161 real ngrad = j/40.0;
162 real gracos= k/20.0;
163 FunDensProp dt, dp = {(real)0.5*rho, (real)0.2*rho, (real)0.5*ngrad, (real)0.3*ngrad};
164 dp.gradab = dp.grada*dp.gradb*gracos;
165 drv2_clear(&d);
166 second_fun(&d, 1, &dp);
167 drv1_clear(&m);
168 first_fun(&m, 1, &dp);
169 test_var(d.df1000, m.df1000, fname, "df1000X", &fail);
170 test_var(d.df0010, m.df0010, fname, "df0010X", &fail);
171 test_var(d.df00001,m.df00001,fname, "df00001X", &fail);
172 test_var(d.df0100, m.df0100, fname, "df0100X", &fail);
173 test_var(d.df0001, m.df0001, fname, "df0001X", &fail);
174 test_var(d.df00001,m.df00001,fname, "df00001X", &fail);
175
176 /* TEST df2000, df1010, df1001 and df10001 */
177 COMP_DER(rho,1e-5,drho,rhoa, first,1);
178 T2(df2000, df1000, drho, "A");
179 T2(df1010, df0010, drho, "A");
180 T2(df1001, df0001, drho, "A");
181 T2(df10001, df00001, drho, "A");
182
183 /* TEST df0200, df0101, df0110 and df01001 */
184 COMP_DER(rho,1e-7,drho,rhob, first,1);
185 T2(df0200, df0100, drho, "A");
186 T2(df0110, df0010, drho, "A");
187 T2(df0101, df0001, drho, "A");
188 T2(df01001, df00001, drho, "A");
189 T2(df1100, df1000, drho, "A");
190
191 /* TEST df1010, df0110, df0020 */
192 COMP_DER(ngrad,1e-5,dgra, grada, first,1);
193 T2(df1010, df1000, dgra, "B");
194 T2(df0020, df0010, dgra, "B");
195 T2(df0110, df0100, dgra, "B");
196
197 /* TEST df1001, df0101, df0002 */
198 COMP_DER(ngrad,1e-5,dgra, gradb, first,1);
199 T2(df1001, df1000, dgra, "B");
200 T2(df0101, df0100, dgra, "B");
201 T2(df0002, df0001, dgra, "B");
202 T2(df00011, df00001, dgra, "B");
203 }
204 }
205 }
206 if(fail==0) printf("%-5s (second order derivatives): OK\n", fname);
207 return fail;
208 }
209
210 /* test_third:
211 test third order derivatives of given functional.
212 It is assumed that the second order derivatives are OK.
213 */
214 static int
test_third(const char * fname,SecondOrderFun second_fun,ThirdOrderFun third_fun)215 test_third(const char* fname,
216 SecondOrderFun second_fun, ThirdOrderFun third_fun)
217 {
218 int i, j, k, fail = 0;
219 real eps, num;
220 FunSecondFuncDrv m, p;
221 FunThirdFuncDrv d;
222 for(i=1; i<=40; i+=GRID_STEP) {
223 for(j=1; j<=40; j+=GRID_STEP) {
224 for(k=-19; k<=19; k+=GRID_STEP) {
225 real rho = i/40.0;
226 real ngrad = j/40.0;
227 real gracos= k/20.0;
228 FunDensProp dt, dp = {(real)0.5*rho, (real)0.3*rho, (real)0.5*ngrad, (real)1*ngrad};
229 dp.gradab = dp.grada*dp.gradb*gracos;
230 drv3_clear(&d);
231 third_fun(&d, 1, &dp);
232 drv2_clear(&m);
233 second_fun(&m, 1, &dp);
234 test_var(d.df1000, m.df1000, fname, "df1000X", &fail);
235 test_var(d.df0100, m.df0100, fname, "df0100X", &fail);
236 test_var(d.df0010, m.df0010, fname, "df0010X", &fail);
237 test_var(d.df0001, m.df0001, fname, "df0001X", &fail);
238 test_var(d.df00001,m.df00001, fname, "df00001X", &fail);
239 test_var(d.df2000, m.df2000, fname, "df2000X", &fail);
240 test_var(d.df1100, m.df1100, fname, "df1100X", &fail);
241 test_var(d.df1010, m.df1010, fname, "df1010X", &fail);
242 test_var(d.df1001, m.df1001, fname, "df1001X", &fail);
243 test_var(d.df10001,m.df10001, fname, "df10001X", &fail);
244 test_var(d.df0200, m.df0200, fname, "df0200X", &fail);
245 test_var(d.df0110, m.df0110, fname, "df0110X", &fail);
246 test_var(d.df0101, m.df0101, fname, "df0101X", &fail);
247 test_var(d.df01001,m.df01001, fname, "df01001X", &fail);
248 test_var(d.df0020, m.df0020, fname, "df0020X", &fail);
249 test_var(d.df0011, m.df0011, fname, "df0011X", &fail);
250 test_var(d.df00101,m.df00101, fname, "df00101X", &fail);
251 test_var(d.df0002, m.df0002, fname, "df0002X", &fail);
252 test_var(d.df00011,m.df00011, fname, "df00011X", &fail);
253 test_var(d.df00002,m.df00002, fname, "df00002X", &fail);
254
255 /* drhoa: test */
256 COMP_DER(rho,1e-7,eps,rhoa, second,2);
257 T2(df3000, df2000, eps, "A");
258 T2(df2100, df1100, eps, "A");
259 T2(df2010, df1010, eps, "A");
260 T2(df2001, df1001, eps, "A");
261 T2(df1200, df0200, eps, "A");
262 T2(df1110, df0110, eps, "A");
263 T2(df1101, df0101, eps, "A");
264 T2(df11001, df01001, eps, "A");
265 T2(df1020, df0020, eps, "A");
266 T2(df1011, df0011, eps, "A");
267 T2(df1002, df0002, eps, "A");
268
269 /*drhob: test */
270 COMP_DER(rho,1e-7,eps,rhob, second,2);
271 T2(df0300, df0200, eps, "A");
272 T2(df0201, df0101, eps, "A");
273 T2(df0210, df0110, eps, "A");
274 T2(df0102, df0002, eps, "A");
275 T2(df0120, df0020, eps, "A");
276 T2(df1200, df1100, eps, "A");
277 T2(df0111, df0011, eps, "A");
278 T2(df1101, df1001, eps, "A");
279 T2(df1110, df1010, eps, "A");
280 T2(df02001, df01001, eps, "A");
281 T2(df11001, df10001, eps, "A");
282
283 /* dgrada: test */
284 COMP_DER(ngrad,1e-7,eps,grada, second,2);
285 T2(df2010, df2000, eps, "A");
286 T2(df0030, df0020, eps, "A");
287 T2(df0021, df0011, eps, "A");
288
289 /* dgradb: test */
290 COMP_DER(ngrad,1e-7,eps,gradb, second,2);
291 T2(df2001, df2000, eps, "A");
292 T2(df0003, df0002, eps, "A");
293 T2(df0012, df0011, eps, "A");
294 }
295 }
296 }
297 if(fail==0) printf("%-5s (third order derivatives): OK\n", fname);
298 return fail;
299 }
300
301 /* test_fourth:
302 test fourth order derivatives of given functional.
303 It is assumed that the second order derivatives are OK.
304 */
305 static int
test_fourth(const char * fname,ThirdOrderFun third_fun,FourthOrderFun fourth_fun)306 test_fourth(const char* fname,
307 ThirdOrderFun third_fun, FourthOrderFun fourth_fun)
308 {
309 int i, j, k, fail = 0;
310 real eps, num;
311 FunThirdFuncDrv m, p;
312 FunFourthFuncDrv d;
313 for(i=1; i<=40; i+=GRID_STEP) {
314 for(j=1; j<=40; j+=GRID_STEP) {
315 for(k=-19; k<=19; k+=GRID_STEP) {
316 real rho = i/40.0;
317 real ngrad = j/40.0;
318 real gracos= k/20.0;
319 FunDensProp dt, dp = {(real)0.5*rho, (real)0.3*rho, (real)0.5*ngrad, (real)1*ngrad};
320 dp.gradab = dp.grada*dp.gradb*gracos;
321 drv4_clear(&d); fourth_fun(&d, 1, &dp);
322 drv3_clear(&m); third_fun(&m, 1, &dp);
323 test_var(d.df1000, m.df1000, fname, "df1000X", &fail);
324 test_var(d.df0100, m.df0100, fname, "df0100X", &fail);
325 test_var(d.df0010, m.df0010, fname, "df0010X", &fail);
326 test_var(d.df0001, m.df0001, fname, "df0001X", &fail);
327 test_var(d.df00001, m.df00001, fname, "df00001X", &fail);
328 test_var(d.df2000, m.df2000, fname, "df2000X", &fail);
329 test_var(d.df1100, m.df1100, fname, "df1100X", &fail);
330 test_var(d.df1010, m.df1010, fname, "df1010X", &fail);
331 test_var(d.df1001, m.df1001, fname, "df1001X", &fail);
332 test_var(d.df10001, m.df10001, fname, "df10001X", &fail);
333 test_var(d.df0200, m.df0200, fname, "df0200X", &fail);
334 test_var(d.df0110, m.df0110, fname, "df0110X", &fail);
335 test_var(d.df0101, m.df0101, fname, "df0101X", &fail);
336 test_var(d.df01001, m.df01001, fname, "df01001X", &fail);
337 test_var(d.df0020, m.df0020, fname, "df0020X", &fail);
338 test_var(d.df0011, m.df0011, fname, "df0011X", &fail);
339 test_var(d.df00101, m.df00101, fname, "df00101X", &fail);
340 test_var(d.df0002, m.df0002, fname, "df0002X", &fail);
341 test_var(d.df00011, m.df00011, fname, "df00011X", &fail);
342 test_var(d.df00002, m.df00002, fname, "df00002X", &fail);
343 test_var(d.df3000, m.df3000, fname, "df3000X", &fail);
344 test_var(d.df2100, m.df2100, fname, "df2100X", &fail);
345 test_var(d.df2010, m.df2010, fname, "df2010X", &fail);
346 test_var(d.df2001, m.df2001, fname, "df2001X", &fail);
347 test_var(d.df20001, m.df20001, fname, "df20001X", &fail);
348 test_var(d.df1200, m.df1200, fname, "df1200X", &fail);
349 test_var(d.df1110, m.df1110, fname, "df1110X", &fail);
350 test_var(d.df1101, m.df1101, fname, "df1101X", &fail);
351 test_var(d.df11001, m.df11001, fname, "df11001X", &fail);
352 test_var(d.df1020, m.df1020, fname, "df1020X", &fail);
353 test_var(d.df1011, m.df1011, fname, "df1011X", &fail);
354 test_var(d.df10101, m.df10101, fname, "df10101X", &fail);
355 test_var(d.df1002, m.df1002, fname, "df1002X", &fail);
356 test_var(d.df10011, m.df10011, fname, "df10011X", &fail);
357 test_var(d.df10002, m.df10002, fname, "df10002X", &fail);
358 test_var(d.df0300, m.df0300, fname, "df0300X", &fail);
359 test_var(d.df0210, m.df0210, fname, "df0210X", &fail);
360 test_var(d.df0201, m.df0201, fname, "df0201X", &fail);
361 test_var(d.df02001, m.df02001, fname, "df02001X", &fail);
362 test_var(d.df0120, m.df0120, fname, "df0120X", &fail);
363 test_var(d.df0111, m.df0111, fname, "df0111X", &fail);
364 test_var(d.df01101, m.df01101, fname, "df01101X", &fail);
365 test_var(d.df0102, m.df0102, fname, "df0102X", &fail);
366 test_var(d.df01011, m.df01011, fname, "df01011X", &fail);
367 test_var(d.df01002, m.df01002, fname, "df01002X", &fail);
368 test_var(d.df0030, m.df0030, fname, "df0030X", &fail);
369 test_var(d.df0021, m.df0021, fname, "df0021X", &fail);
370 test_var(d.df00201, m.df00201, fname, "df00201X", &fail);
371 test_var(d.df0012, m.df0012, fname, "df0012X", &fail);
372 test_var(d.df00111, m.df00111, fname, "df00111X", &fail);
373 test_var(d.df00102, m.df00102, fname, "df00102X", &fail);
374 test_var(d.df0003, m.df0003, fname, "df0003X", &fail);
375 test_var(d.df00021, m.df00021, fname, "df00021X", &fail);
376 test_var(d.df00012, m.df00012, fname, "df00012X", &fail);
377 test_var(d.df00003, m.df00003, fname, "df00003X", &fail);
378
379 /* drhoa: test */
380 COMP_DER(rho,1e-6,eps,rhoa, third,3);
381 T2(df4000, df3000, eps, "A");
382 T2(df3100, df2100, eps, "A");
383 T2(df3010, df2010, eps, "A");
384 T2(df3001, df2001, eps, "A");
385 T2(df30001, df20001,eps, "A");
386 T2(df2200, df1200, eps, "A");
387 T2(df2110, df1110, eps, "A");
388 T2(df2101, df1101, eps, "A");
389 T2(df21001, df11001,eps, "A");
390 T2(df2020, df1020, eps, "A");
391 T2(df2011, df1011, eps, "A");
392 T2(df20101, df10101,eps, "A");
393 T2(df2002, df1002, eps, "A");
394 T2(df20011, df10011,eps, "A");
395 T2(df20002, df10002,eps, "A");
396 T2(df1300, df0300, eps, "A");
397 T2(df1210, df0210, eps, "A");
398 T2(df1201, df0201, eps, "A");
399 T2(df12001, df02001,eps, "A");
400 T2(df1120, df0120, eps, "A");
401 T2(df1111, df0111, eps, "A");
402 T2(df11101, df01101,eps, "A");
403 T2(df1102, df0102, eps, "A");
404 T2(df11011, df01011,eps, "A");
405 T2(df11002, df01002,eps, "A");
406 T2(df1030, df0030, eps, "A");
407 T2(df1021, df0021, eps, "A");
408 T2(df10201, df00201,eps, "A");
409 T2(df1012, df0012, eps, "A");
410 T2(df10111, df00111,eps, "A");
411 T2(df10102, df00102,eps, "A");
412 T2(df1003, df0003, eps, "A");
413 T2(df10021, df00021,eps, "A");
414 T2(df10012, df00012,eps, "A");
415 T2(df10003, df00003,eps, "A");
416
417 /*drhob: test */
418 COMP_DER(rho,1e-7,eps,rhob, third,3);
419 T2(df3100, df3000, eps, "B");
420 T2(df2200, df2100, eps, "B");
421 T2(df2110, df2010, eps, "B");
422 T2(df2101, df2001, eps, "B");
423 T2(df21001, df20001,eps, "B");
424 T2(df1300, df1200, eps, "B");
425 T2(df1210, df1110, eps, "B");
426 T2(df1201, df1101, eps, "B");
427 T2(df12001, df11001,eps, "B");
428 T2(df1120, df1020, eps, "B");
429 T2(df1111, df1011, eps, "B");
430 T2(df11101, df10101,eps, "B");
431 T2(df1102, df1002, eps, "B");
432 T2(df11011, df10011,eps, "B");
433 T2(df11002, df10002,eps, "B");
434 T2(df0400, df0300, eps, "B");
435 T2(df0310, df0210, eps, "B");
436 T2(df0301, df0201, eps, "B");
437 T2(df03001, df02001,eps, "B");
438 T2(df0220, df0120, eps, "B");
439 T2(df0211, df0111, eps, "B");
440 T2(df02101, df01101,eps, "B");
441 T2(df0202, df0102, eps, "B");
442 T2(df02011, df01011,eps, "B");
443 T2(df02002, df01002,eps, "B");
444 T2(df0130, df0030, eps, "B");
445 T2(df0121, df0021, eps, "B");
446 T2(df01201, df00201,eps, "B");
447 T2(df0112, df0012, eps, "B");
448 T2(df01111, df00111,eps, "B");
449 T2(df01102, df00102,eps, "B");
450 T2(df0103, df0003, eps, "B");
451 T2(df01021, df00021,eps, "B");
452 T2(df01012, df00012,eps, "B");
453 T2(df01003, df00003,eps, "B");
454 /* dgrada: test */
455 COMP_DER(ngrad,1e-7,eps,grada, third,3);
456 T2(df3010, df3000, eps, "C");
457 T2(df2110, df2100, eps, "C");
458 T2(df2020, df2010, eps, "C");
459 T2(df2011, df2001, eps, "C");
460 T2(df20101, df20001,eps, "C");
461 T2(df20011, df20001,eps, "C");
462 T2(df1210, df1200, eps, "C");
463 T2(df1120, df1110, eps, "C");
464 T2(df1111, df1101, eps, "C");
465 T2(df11101, df11001,eps, "C");
466 T2(df1030, df1020, eps, "C");
467 T2(df1021, df1011, eps, "C");
468 T2(df10201, df10101,eps, "C");
469 T2(df1012, df1002, eps, "C");
470 T2(df10111, df10011,eps, "C");
471 T2(df10102, df10002,eps, "C");
472 T2(df1003, df1002, eps, "C");
473 T2(df0310, df0300, eps, "C");
474 T2(df0220, df0210, eps, "C");
475 T2(df0211, df0201, eps, "C");
476 T2(df02101, df02001,eps, "C");
477 T2(df0130, df0120, eps, "C");
478 T2(df0121, df0111, eps, "C");
479 T2(df01201, df01101,eps, "C");
480 T2(df0112, df0102, eps, "C");
481 T2(df01111, df01011,eps, "C");
482 T2(df01102, df01002,eps, "C");
483 T2(df0040, df0030, eps, "C");
484 T2(df0031, df0021, eps, "C");
485 T2(df00301, df00201,eps, "C");
486 T2(df0022, df0012, eps, "C");
487 T2(df00211, df00111,eps, "C");
488 T2(df00202, df00102,eps, "C");
489 T2(df0013, df0003, eps, "C");
490 T2(df00121, df00021,eps, "C");
491 T2(df00112, df00012,eps, "C");
492 T2(df00103, df00003,eps, "C");
493
494 /* dgradb: test */
495 COMP_DER(ngrad,1e-7,eps,gradb, third,3);
496 T2(df3001, df3000, eps, "C");
497 T2(df2101, df2100, eps, "C");
498 T2(df2011, df2010, eps, "C");
499 T2(df2002, df2001, eps, "C");
500 T2(df20011, df20001,eps, "C");
501 T2(df1201, df1200, eps, "C");
502 T2(df1111, df1110, eps, "C");
503 T2(df1102, df1101, eps, "C");
504 T2(df11011, df11001,eps, "C");
505 T2(df1021, df1020, eps, "C");
506 T2(df1012, df1011, eps, "C");
507 T2(df10111, df10101,eps, "C");
508 T2(df1003, df1002, eps, "C");
509 T2(df10021, df10011,eps, "C");
510 T2(df10012, df10002,eps, "C");
511 T2(df0301, df0300, eps, "C");
512 T2(df0211, df0210, eps, "C");
513 T2(df0202, df0201, eps, "C");
514 T2(df02011, df02001,eps, "C");
515 T2(df0121, df0120, eps, "C");
516 T2(df0112, df0111, eps, "C");
517 T2(df01111, df01101,eps, "C");
518 T2(df0103, df0102, eps, "C");
519 T2(df01021, df01011,eps, "C");
520 T2(df01012, df01002,eps, "C");
521 T2(df0031, df0030, eps, "C");
522 T2(df0022, df0021, eps, "C");
523 T2(df00211, df00201,eps, "C");
524 T2(df0013, df0012, eps, "C");
525 T2(df00121, df00101,eps, "C");
526 T2(df00112, df00102,eps, "C");
527 T2(df0004, df0003, eps, "C");
528 T2(df00031, df00021,eps, "C");
529 T2(df00022, df00012,eps, "C");
530 T2(df00013, df00003,eps, "C");
531 }
532 }
533 }
534 if(fail==0) printf("%-5s (fourth order derivatives): OK\n", fname);
535 return fail;
536 }
537
538 static int
test_derivatives(Functional * f,int * orders,DaltonEnFunc dal_fun)539 test_derivatives(Functional* f, int *orders, DaltonEnFunc dal_fun)
540 {
541 int res = 0;
542 /* if(dal_fun) res = test_energy(f->name, f->func, dal_fun); */
543 if(!res && (!orders || orders[0]) )
544 res = test_first(f->name, f->func, f->first);
545 if(!res && (!orders || orders[1]) )
546 res = test_second(f->name, f->first, f->second);
547 if(!res && (!orders || orders[2]) )
548 res = test_third(f->name, f->second, f->third);
549 if(!res && (!orders || orders[3]) )
550 res = test_fourth(f->name, f->third, f->fourth);
551 return res;
552 }
553
554
555 /* main:
556 this is the main test program.
557 */
558 int
main(int argc,char * argv[])559 main(int argc, char* argv[])
560 {
561 int res = 0, i, length, argidx, funcsel = 0;
562 static int funco[] = { 0, 0, 0, 0 };
563 char* arg;
564 Functional* func;
565
566 if(argc<=1) {
567 fprintf(stderr,
568 "Functional derivative tester:\n"
569 "usage: fun-tester [-n] <functional> <options>\n"
570 "-n - test only selected order of derivatives\n"
571 "example: fun-tester GGAKey becke=1 lyp=1\n");
572 return 1;
573 }
574 for(argidx=1; argidx<argc && argv[argidx][0] == '-'; argidx++)
575 switch(argv[argidx][1]) {
576 case '1': funcsel = 1; funco[0]=1; break;
577 case '2': funcsel = 1; funco[1]=1; break;
578 case '3': funcsel = 1; funco[2]=1; break;
579 case '4': funcsel = 1; funco[3]=1; break;
580 default: fprintf(stderr, "option %s is unknown.\n", argv[argidx]);
581 }
582
583 for(i=0; available_functionals[i]; i++)
584 if(strcasecmp(argv[argidx], available_functionals[i]->name)==0)
585 break;
586 if(available_functionals[i]==NULL) {
587 fprintf(stderr, "Functional '%s' not found.\n\n"
588 "Available functionals:\n", argv[1]);
589 for(i=0; available_functionals[i]; i++)
590 fprintf(stderr, " %s\n", available_functionals[i]->name);
591 return 2;
592 }
593 func = available_functionals[i];
594 argidx++;
595 for(length=1, i=argidx; i<argc; i++)
596 length += strlen(argv[i])+1;
597
598 arg = (char*)malloc(length);
599 if(argc>argidx+1)
600 strcpy(arg, argv[argidx]);
601 else
602 *arg = '\0';
603
604 for(i=argidx; i<argc; i++) {
605 strcat(arg, " ");
606 strcat(arg, argv[i]);
607 }
608
609 if(!func->read(arg)) {
610 fprintf(stderr, "Reading configuration for %s from '%s' failed.\n",
611 argv[1], arg);
612 free(arg);
613 return 3;
614 }
615 free(arg);
616
617 if(func->report)
618 func->report();
619 res += test_derivatives(func, funcsel ? funco: NULL, NULL);
620
621 if(res>0)
622 printf("%i tests failed.\n", res);
623 else printf("OK\n");
624 return res;
625 }
626