1 /*
2 Copyright (C) 2006-2007 M.A.L. Marques
3 Copyright (C) 2014 Susi Lehtola
4
5 This Source Code Form is subject to the terms of the Mozilla Public
6 License, v. 2.0. If a copy of the MPL was not distributed with this
7 file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 */
9
10 #include <stdio.h>
11 #include <stdlib.h>
12 #include <math.h>
13 #include <string.h>
14
15 #include <xc.h>
16 #include <util.h>
17
18 /* Buffer size (line length) for file reads */
19 #define BUFSIZE 1024
20
21 typedef struct {
22 /* Amount of data points */
23 int n;
24
25 /* Input: density, gradient, laplacian and kinetic energy density */
26 double *rho;
27 double *sigma;
28 double *lapl;
29 double *tau;
30
31 /* Output: energy density */
32 double *zk;
33
34 /* .. and potentials for density, gradient, laplacian and tau */
35 double *vrho;
36 double *vsigma;
37 double *vlapl;
38 double *vtau;
39
40 /* ... and second derivatives */
41 double *v2rho2;
42 double *v2tau2;
43 double *v2lapl2;
44 double *v2rhotau;
45 double *v2rholapl;
46 double *v2lapltau;
47 double *v2sigma2;
48 double *v2rhosigma;
49 double *v2sigmatau;
50 double *v2sigmalapl;
51
52 /* ... and third derivatives */
53 double *v3rho3;
54 } values_t;
55
allocate_memory(values_t * data,int nspin,int order)56 void allocate_memory(values_t *data, int nspin, int order)
57 {
58 data->zk = NULL;
59 data->vrho = NULL;
60 data->vsigma = NULL;
61 data->vlapl = NULL;
62 data->vtau = NULL;
63 data->v2rho2 = NULL;
64 data->v2tau2 = NULL;
65 data->v2lapl2 = NULL;
66 data->v2rhotau = NULL;
67 data->v2rholapl = NULL;
68 data->v2lapltau = NULL;
69 data->v2sigma2 = NULL;
70 data->v2rhosigma = NULL;
71 data->v2sigmatau = NULL;
72 data->v2sigmalapl = NULL;
73 data->v3rho3 = NULL;
74
75 switch(nspin) {
76 case (XC_UNPOLARIZED):
77 data->rho = (double*) libxc_calloc(data->n, sizeof(double));
78 data->sigma = (double*) libxc_calloc(data->n, sizeof(double));
79 data->lapl = (double*) libxc_calloc(data->n, sizeof(double));
80 data->tau = (double*) libxc_calloc(data->n, sizeof(double));
81 switch (order) {
82 case (0):
83 data->zk = (double*) libxc_calloc(data->n, sizeof(double));
84 break;
85 case (1):
86 data->vrho = (double*) libxc_calloc(data->n, sizeof(double));
87 data->vsigma = (double*) libxc_calloc(data->n, sizeof(double));
88 data->vlapl = (double*) libxc_calloc(data->n, sizeof(double));
89 data->vtau = (double*) libxc_calloc(data->n, sizeof(double));
90 break;
91 case (2):
92 data->v2rho2 = (double*) libxc_calloc(data->n, sizeof(double));
93 data->v2tau2 = (double*) libxc_calloc(data->n, sizeof(double));
94 data->v2lapl2 = (double*) libxc_calloc(data->n, sizeof(double));
95 data->v2rhotau = (double*) libxc_calloc(data->n, sizeof(double));
96 data->v2rholapl = (double*) libxc_calloc(data->n, sizeof(double));
97 data->v2lapltau = (double*) libxc_calloc(data->n, sizeof(double));
98 data->v2sigma2 = (double*) libxc_calloc(data->n, sizeof(double));
99 data->v2rhosigma = (double*) libxc_calloc(data->n, sizeof(double));
100 data->v2sigmatau = (double*) libxc_calloc(data->n, sizeof(double));
101 data->v2sigmalapl = (double*) libxc_calloc(data->n, sizeof(double));
102 break;
103 case (3):
104 data->v3rho3 = (double*) libxc_calloc(data->n, sizeof(double));
105 break;
106 default:
107 fprintf(stderr, "order = %i not recognized.\n", order);
108 exit(2);
109 }
110 break;
111
112 case (XC_POLARIZED):
113 data->rho = (double*) libxc_calloc(2*data->n, sizeof(double));
114 data->sigma = (double*) libxc_calloc(3*data->n, sizeof(double));
115 data->lapl = (double*) libxc_calloc(2*data->n, sizeof(double));
116 data->tau = (double*) libxc_calloc(2*data->n, sizeof(double));
117 switch (order) {
118 case (0):
119 data->zk = (double*) libxc_calloc(data->n, sizeof(double));
120 break;
121 case (1):
122 data->vrho = (double*) libxc_calloc(2*data->n, sizeof(double));
123 data->vsigma = (double*) libxc_calloc(3*data->n, sizeof(double));
124 data->vlapl = (double*) libxc_calloc(2*data->n, sizeof(double));
125 data->vtau = (double*) libxc_calloc(2*data->n, sizeof(double));
126 break;
127 case (2):
128 data->v2rho2 = (double*) libxc_calloc(3*data->n, sizeof(double));
129 data->v2tau2 = (double*) libxc_calloc(3*data->n, sizeof(double));
130 data->v2lapl2 = (double*) libxc_calloc(3*data->n, sizeof(double));
131 data->v2rhotau = (double*) libxc_calloc(4*data->n, sizeof(double));
132 data->v2rholapl = (double*) libxc_calloc(4*data->n, sizeof(double));
133 data->v2lapltau = (double*) libxc_calloc(4*data->n, sizeof(double));
134 data->v2sigma2 = (double*) libxc_calloc(6*data->n, sizeof(double));
135 data->v2rhosigma = (double*) libxc_calloc(6*data->n, sizeof(double));
136 data->v2sigmatau = (double*) libxc_calloc(6*data->n, sizeof(double));
137 data->v2sigmalapl = (double*) libxc_calloc(6*data->n, sizeof(double));
138 break;
139 case (3):
140 data->v3rho3 = (double*) libxc_calloc(4*data->n, sizeof(double));
141 break;
142 default:
143 fprintf(stderr, "order = %i not recognized.\n", order);
144 exit(2);
145 }
146 break;
147
148 default:
149 fprintf(stderr, "nspin = %i not recognized.\n", nspin);
150 exit(2);
151 }
152 }
153
154
155 #define FREE_NULL(p) {if(p!=NULL) {libxc_free(p);}; p=NULL;}
156
free_memory(values_t * val)157 void free_memory(values_t *val)
158 {
159 FREE_NULL(val->rho);
160 FREE_NULL(val->sigma);
161 FREE_NULL(val->lapl);
162 FREE_NULL(val->tau);
163 FREE_NULL(val->zk);
164 FREE_NULL(val->vrho);
165 FREE_NULL(val->vsigma);
166 FREE_NULL(val->vlapl);
167 FREE_NULL(val->vtau);
168 FREE_NULL(val->v2rho2);
169 FREE_NULL(val->v2tau2);
170 FREE_NULL(val->v2lapl2);
171 FREE_NULL(val->v2rhotau);
172 FREE_NULL(val->v2rholapl);
173 FREE_NULL(val->v2lapltau);
174 FREE_NULL(val->v2sigma2);
175 FREE_NULL(val->v2rhosigma);
176 FREE_NULL(val->v2sigmatau);
177 FREE_NULL(val->v2sigmalapl);
178 FREE_NULL(val->v3rho3);
179 }
180
drop_laplacian(values_t * val)181 void drop_laplacian(values_t *val)
182 {
183 FREE_NULL(val->lapl);
184 FREE_NULL(val->vlapl);
185 FREE_NULL(val->v2lapl2);
186 FREE_NULL(val->v2rholapl);
187 FREE_NULL(val->v2lapltau);
188 FREE_NULL(val->v2sigmalapl);
189 }
190
read_data(const char * file,int nspin,int order)191 values_t read_data(const char *file, int nspin, int order) {
192 /* Format string */
193 static const char fmt[]="%lf %lf %lf %lf %lf %lf %lf %lf %lf";
194
195 /* Data buffer */
196 char buf[BUFSIZE];
197 char *cp;
198 /* Input data file */
199 FILE *in;
200 /* Loop index */
201 int i;
202 /* Amount of points succesfully read */
203 int nsucc;
204 /* Returned data */
205 values_t data;
206
207 /* Helper variables */
208 double rhoa, rhob;
209 double sigmaaa, sigmaab, sigmabb;
210 double lapla, laplb;
211 double taua, taub;
212
213 /* Open file */
214 in=fopen(file,"r");
215 if(!in) {
216 fprintf(stderr,"Error opening input file %s.\n",file);
217 exit(3);
218 }
219
220 /* Read amount of data points */
221 cp=fgets(buf,BUFSIZE,in);
222 if(cp!=buf) {
223 fprintf(stderr,"Error reading amount of data points.\n");
224 exit(5);
225 }
226 nsucc=sscanf(buf,"%i",&data.n);
227 if(nsucc!=1) {
228 fprintf(stderr,"Error reading amount of input data points.\n");
229 exit(4);
230 }
231
232 /* Allocate memory */
233 allocate_memory(&data, nspin, order);
234
235 for(i=0;i<data.n;i++) {
236 /* Next line of input */
237 cp=fgets(buf,BUFSIZE,in);
238 if(cp!=buf) {
239 fprintf(stderr,"Read error on line %i.\n",i+1);
240 free_memory(&data);
241 exit(5);
242 }
243 /* Read data */
244 nsucc=sscanf(buf, fmt, &rhoa, &rhob, &sigmaaa, &sigmaab, &sigmabb, \
245 &lapla, &laplb, &taua, &taub);
246
247 /* Error control */
248 if(nsucc!=9) {
249 fprintf(stderr,"Read error on line %i: only %i entries read.\n",i+1,nsucc);
250 free_memory(&data);
251 exit(5);
252 }
253
254 /* Store data (if clause suboptimal here but better for code clarity) */
255 if(nspin==XC_POLARIZED) {
256 data.rho[2*i]=rhoa;
257 data.rho[2*i+1]=rhob;
258 data.sigma[3*i]=sigmaaa;
259 data.sigma[3*i+1]=sigmaab;
260 data.sigma[3*i+2]=sigmabb;
261 data.lapl[2*i]=lapla;
262 data.lapl[2*i+1]=laplb;
263 data.tau[2*i]=taua;
264 data.tau[2*i+1]=taub;
265 } else {
266 /* Construct full density data from alpha and beta channels */
267 data.rho[i]=rhoa + rhob;
268 data.sigma[i]=sigmaaa + sigmabb + 2.0*sigmaab;
269 data.lapl[i]=lapla + laplb;
270 data.tau[i]=taua + taub;
271 }
272 }
273
274 /* Close input file */
275 fclose(in);
276
277 return data;
278 }
279
280 /* Print helpers */
printe1(FILE * out,double * x,size_t idx)281 void printe1(FILE *out, double *x, size_t idx) {
282 fprintf(out," % .16e",x[idx]);
283 }
284
printe2(FILE * out,double * x,size_t idx)285 void printe2(FILE *out, double *x, size_t idx) {
286 fprintf(out," % .16e % .16e",x[idx],x[idx+1]);
287 }
288
print03(FILE * out,double * x,size_t idx)289 void print03(FILE *out, double *x, size_t idx) {
290 fprintf(out," % .16e % .16e % .16e",0.0,0.0,0.0);
291 }
292
print01(FILE * out,double * x,size_t idx)293 void print01(FILE *out, double *x, size_t idx) {
294 fprintf(out," % .16e",0.0);
295 }
296
print02(FILE * out,double * x,size_t idx)297 void print02(FILE *out, double *x, size_t idx) {
298 fprintf(out," % .16e % .16e",0.0,0.0);
299 }
300
printe3(FILE * out,double * x,size_t idx)301 void printe3(FILE *out, double *x, size_t idx) {
302 fprintf(out," % .16e % .16e % .16e",x[idx],x[idx+1],x[idx+2]);
303 }
304
305 /*----------------------------------------------------------*/
main(int argc,char * argv[])306 int main(int argc, char *argv[])
307 {
308 int func_id, nspin, order, i;
309 /* Helpers for properties that may not have been implemented */
310 double *zk, *vrho, *v2rho2, *v3rho3;
311
312 static const char sfmt[] =" %23s";
313 static const char sfmt2[]=" %23s %23s";
314 static const char sfmt3[]=" %23s %23s %23s";
315
316 /* Data array */
317 values_t d;
318 /* Functional evaluator */
319 xc_func_type func;
320 /* Flags for functional */
321 int flags;
322 /* Functional family */
323 int family;
324 /* Output file */
325 FILE *out;
326 /* Output file name */
327 char *fname;
328
329 /* Normal print functions */
330 void (*print1)(FILE *out, double *x, size_t idx);
331 void (*print2)(FILE *out, double *x, size_t idx);
332 void (*print3)(FILE *out, double *x, size_t idx);
333 /* Laplacian data print functions */
334 void (*plapl1)(FILE *out, double *x, size_t idx);
335 void (*plapl2)(FILE *out, double *x, size_t idx);
336 void (*plapl3)(FILE *out, double *x, size_t idx);
337
338 if(argc != 6) {
339 fprintf(stderr, "Usage:\n%s funct nspin order input output\n", argv[0]);
340 exit(1);
341 }
342
343 /* Get functional id */
344 func_id = xc_functional_get_number(argv[1]);
345 if(func_id <= 0) {
346 fprintf(stderr, "Functional '%s' not found\n", argv[1]);
347 exit(1);
348 }
349
350 /* Spin-polarized or unpolarized ? */
351 nspin = atoi(argv[2]);
352
353 /* Order of derivatives to compute */
354 order = atoi(argv[3]);
355
356 /* Read in data */
357 d = read_data(argv[4], nspin, order);
358
359 /* Initialize functional */
360 if(xc_func_init(&func, func_id, nspin)) {
361 fprintf(stderr, "Functional '%d' (%s) not found.\nPlease report a bug against functional_get_number.\n", func_id, argv[1]);
362 exit(1);
363 }
364 /* Get flags */
365 flags = func.info->flags;
366 family = func.info->family;
367
368 /* Set helpers */
369 zk = (flags & XC_FLAGS_HAVE_EXC) ? d.zk : NULL;
370 vrho = (flags & XC_FLAGS_HAVE_VXC) ? d.vrho : NULL;
371 v2rho2 = (flags & XC_FLAGS_HAVE_FXC) ? d.v2rho2 : NULL;
372 v3rho3 = (flags & XC_FLAGS_HAVE_KXC) ? d.v3rho3 : NULL;
373 print1 = printe1;
374 print2 = printe2;
375 print3 = printe3;
376 plapl1 = printe1;
377 plapl2 = printe2;
378 plapl3 = printe3;
379
380 /* If functional doesn't need laplacian, drop the values and print
381 out zeros for the functional value */
382 if(!(flags & XC_FLAGS_NEEDS_LAPLACIAN)) {
383 drop_laplacian(&d);
384 plapl1 = print01;
385 plapl2 = print02;
386 plapl3 = print03;
387 }
388
389 /* Evaluate xc functional */
390 switch(family) {
391 case XC_FAMILY_LDA:
392 case XC_FAMILY_HYB_LDA:
393 xc_lda(&func, d.n, d.rho, zk, vrho, v2rho2, v3rho3, NULL);
394 break;
395 case XC_FAMILY_GGA:
396 case XC_FAMILY_HYB_GGA:
397 xc_gga(&func, d.n, d.rho, d.sigma, zk, vrho, d.vsigma,
398 v2rho2, d.v2rhosigma, d.v2sigma2, NULL, NULL, NULL, NULL,
399 NULL, NULL, NULL, NULL, NULL);
400 break;
401 case XC_FAMILY_MGGA:
402 case XC_FAMILY_HYB_MGGA:
403 xc_mgga(&func, d.n, d.rho, d.sigma, d.lapl, d.tau, zk, vrho, d.vsigma, d.vlapl, d.vtau,
404 v2rho2, d.v2rhosigma, d.v2rholapl, d.v2rhotau, d.v2sigma2, d.v2sigmalapl, d.v2sigmatau, d.v2lapl2, d.v2lapltau, d.v2tau2,
405 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
406 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
407 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
408 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
409 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
410 NULL, NULL, NULL, NULL, NULL
411 );
412
413 break;
414
415 default:
416 fprintf(stderr,"Support for family %i not implemented.\n",family);
417 free_memory(&d);
418 exit(1);
419 }
420
421 /* Open output file */
422 fname = argv[5];
423 out = fopen(fname,"w");
424 if(!out) {
425 fprintf(stderr,"Error opening output file %s.\n",fname);
426 free_memory(&d);
427 exit(1);
428 }
429
430 /* Functional id and amount of lines in output */
431 fprintf(out, "%i %i %i\n", func_id, d.n, order);
432
433 switch (order) {
434 case (0): /* energy */
435 fprintf(out, sfmt, "zk");
436 break;
437 case (1): /* first order derivatives */
438 if (nspin == XC_POLARIZED) {
439 fprintf(out, sfmt2, "vrho(a)", "vrho(b)");
440 if (family & (XC_FAMILY_GGA | XC_FAMILY_HYB_GGA | XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA))
441 fprintf(out, sfmt3, "vsigma(aa)", "vsigma(ab)", "vsigma(bb)");
442 if (family & (XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
443 fprintf(out, sfmt2, "vlapl(a)", "vlapl(b)");
444 fprintf(out, sfmt2, "vtau(a)", "vtau(b)");
445 }
446 } else {
447 fprintf(out, sfmt, "vrho");
448 if (family & (XC_FAMILY_GGA | XC_FAMILY_HYB_GGA | XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA))
449 fprintf(out, sfmt, "vsigma");
450 if(family & (XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
451 fprintf(out, sfmt, "vlapl");
452 fprintf(out, sfmt, "vtau");
453 }
454 }
455 break;
456
457 case (2): /* second order derivatives */
458 if (nspin == XC_POLARIZED) {
459 fprintf(out,sfmt3,"v2rho(aa)","v2rho(ab)","v2rho(bb)");
460 if(family & (XC_FAMILY_GGA | XC_FAMILY_HYB_GGA | XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
461 fprintf(out, sfmt3, "v2sigma2(aa-aa)", "v2sigma2(aa-ab)", "v2sigma2(aa-bb)");
462 fprintf(out, sfmt3, "v2sigma2(ab-ab)", "v2sigma2(ab-bb)", "v2sigma2(bb-bb)");
463 fprintf(out, sfmt3, "v2rho(a)sigma(aa)", "v2rho(a)sigma(ab)", "v2rho(a)sigma(bb)");
464 fprintf(out, sfmt3, "v2rho(b)sigma(aa)", "v2rho(b)sigma(ab)", "v2rho(b)sigma(bb)");
465 }
466 if(family & (XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
467 fprintf(out, sfmt3, "v2lapl2(aa)", "v2lapl2(ab)", "v2lapl2(bb)");
468 fprintf(out, sfmt3, "v2tau2(aa)", "v2tau2(ab)", "v2tau2(bb)");
469 fprintf(out, sfmt3, "v2rholapl(aa)", "v2rholapl(ab)", "v2rholapl(bb)");
470 fprintf(out, sfmt3, "v2rhotau(aa)", "v2rhotau(ab)", "v2rhotau(bb)");
471 fprintf(out, sfmt3, "v2lapltau(aa)", "v2lapltau(ab)", "v2lapltau(bb)");
472 fprintf(out, sfmt3, "v2sigma(aa)tau(a)", "v2sigma(aa)tau(b)", "v2sigma(ab)tau(a)");
473 fprintf(out, sfmt3, "v2sigma(ab)tau(b)", "v2sigma(bb)tau(a)", "v2sigma(bb)tau(b)");
474 fprintf(out, sfmt3, "v2sigma(aa)lapl(a)", "v2sigma(aa)lapl(b)", "v2sigma(ab)lapl(a)");
475 fprintf(out, sfmt3, "v2sigma(ab)lapl(b)", "v2sigma(bb)lapl(a)", "v2sigma(bb)lapl(b)");
476 }
477 } else {
478 fprintf(out,sfmt,"v2rho");
479 if(family & (XC_FAMILY_GGA | XC_FAMILY_HYB_GGA | XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
480 fprintf(out, sfmt, "v2sigma2");
481 fprintf(out, sfmt, "v2rhosigma");
482 }
483
484 if(family & (XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
485 fprintf(out, sfmt, "v2lapl2");
486 fprintf(out, sfmt, "v2tau2");
487 fprintf(out, sfmt, "v2rholapl");
488 fprintf(out, sfmt, "v2rhotau");
489 fprintf(out, sfmt, "v2lapltau");
490 fprintf(out, sfmt, "v2sigmatau");
491 fprintf(out, sfmt, "v2sigmalapl");
492 }
493 }
494 break;
495
496 default: /* higher order derivatives ... to be done */
497 fprintf(stderr, "order = %i not recognized.\n", order);
498 exit(2);
499 }
500 fprintf(out,"\n");
501
502 /* Loop over data points */
503 for(i=0;i<d.n;i++) {
504
505 switch (order) {
506 case (0): /* energy */
507 print1(out, d.zk, i);
508 break;
509 case (1): /* first order derivatives */
510 if (nspin == XC_POLARIZED) {
511 print2(out, d.vrho, 2 * i);
512 if (family & (XC_FAMILY_GGA | XC_FAMILY_HYB_GGA | XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA))
513 print3(out, d.vsigma, 3 * i);
514 if (family & (XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
515 plapl2(out, d.vlapl, 2 * i);
516 print2(out, d.vtau, 2 * i);
517 }
518 } else {
519 print1(out, d.vrho, i);
520 if (family & (XC_FAMILY_GGA | XC_FAMILY_HYB_GGA | XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA))
521 print1(out, d.vsigma, i);
522 if (family & (XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
523 plapl1(out, d.vlapl, i);
524 print1(out, d.vtau, i);
525 }
526 }
527 break;
528
529 case (2): /* second order derivatives */
530 if (nspin == XC_POLARIZED) {
531 print3(out, d.v2rho2, 3*i);
532 if(family & (XC_FAMILY_GGA | XC_FAMILY_HYB_GGA | XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
533 print3(out, d.v2sigma2, 6*i);
534 print3(out, d.v2sigma2, 6*i + 3);
535 print3(out, d.v2rhosigma, 6*i);
536 print3(out, d.v2rhosigma, 6*i + 3);
537 }
538 if(family & (XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
539 plapl3(out, d.v2lapl2, 3*i);
540 print3(out, d.v2tau2, 3*i);
541 plapl3(out, d.v2rholapl, 3*i);
542 print3(out, d.v2rhotau, 3*i);
543 plapl3(out, d.v2lapltau, 3*i);
544 print3(out, d.v2sigmatau, 3*i);
545 print3(out, d.v2sigmatau, 3*i + 3);
546 plapl3(out, d.v2sigmalapl, 3*i);
547 plapl3(out, d.v2sigmalapl, 3*i + 3);
548 }
549 } else {
550 print1(out, d.v2rho2, i);
551 if(family & (XC_FAMILY_GGA | XC_FAMILY_HYB_GGA | XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
552 print1(out, d.v2sigma2, i);
553 print1(out, d.v2rhosigma, i);
554 }
555 if(family & (XC_FAMILY_MGGA | XC_FAMILY_HYB_MGGA)) {
556 plapl1(out, d.v2lapl2, i);
557 print1(out, d.v2tau2, i);
558 plapl1(out, d.v2rholapl, i);
559 print1(out, d.v2rhotau, i);
560 plapl1(out, d.v2lapltau, i);
561 print1(out, d.v2sigmatau, i);
562 plapl1(out, d.v2sigmalapl, i);
563 }
564 }
565 break;
566
567 default: /* higher order derivatives ... to be done */
568 fprintf(stderr, "order = %i not recognized.\n", order);
569 exit(2);
570 }
571
572 fprintf(out,"\n");
573 }
574
575 xc_func_end(&func);
576 free_memory(&d);
577 fclose(out);
578
579 return 0;
580 }
581