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