1 /* Copyright 2014-2018 The PySCF Developers. All Rights Reserved.
2
3 Licensed under the Apache License, Version 2.0 (the "License");
4 you may not use this file except in compliance with the License.
5 You may obtain a copy of the License at
6
7 http://www.apache.org/licenses/LICENSE-2.0
8
9 Unless required by applicable law or agreed to in writing, software
10 distributed under the License is distributed on an "AS IS" BASIS,
11 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 See the License for the specific language governing permissions and
13 limitations under the License.
14
15 *
16 * Authors: Qiming Sun <osirpt.sun@gmail.com>
17 * Susi Lehtola <susi.lehtola@gmail.com>
18 *
19 * libxc from
20 * http://www.tddft.org/programs/octopus/wiki/index.php/Libxc:manual
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <assert.h>
26 #include <xc.h>
27 #define MAX(X,Y) ((X) > (Y) ? (X) : (Y))
28
29 /* Extracted from comments of libxc:gga.c
30
31 sigma_st = grad rho_s . grad rho_t
32 zk = energy density per unit particle
33
34 vrho_s = d n*zk / d rho_s
35 vsigma_st = d n*zk / d sigma_st
36
37 v2rho2_st = d^2 n*zk / d rho_s d rho_t
38 v2rhosigma_svx = d^2 n*zk / d rho_s d sigma_tv
39 v2sigma2_stvx = d^2 n*zk / d sigma_st d sigma_vx
40
41 v3rho3_stv = d^3 n*zk / d rho_s d rho_t d rho_v
42 v3rho2sigma_stvx = d^3 n*zk / d rho_s d rho_t d sigma_vx
43 v3rhosigma2_svxyz = d^3 n*zk / d rho_s d sigma_vx d sigma_yz
44 v3sigma3_stvxyz = d^3 n*zk / d sigma_st d sigma_vx d sigma_yz
45
46 if nspin == 2
47 rho(2) = (u, d)
48 sigma(3) = (uu, ud, dd)
49
50 * vxc(N*5):
51 vrho(2) = (u, d)
52 vsigma(3) = (uu, ud, dd)
53
54 * fxc(N*45):
55 v2rho2(3) = (u_u, u_d, d_d)
56 v2rhosigma(6) = (u_uu, u_ud, u_dd, d_uu, d_ud, d_dd)
57 v2sigma2(6) = (uu_uu, uu_ud, uu_dd, ud_ud, ud_dd, dd_dd)
58 v2lapl2(3)
59 vtau2(3)
60 v2rholapl(4)
61 v2rhotau(4)
62 v2lapltau(4)
63 v2sigmalapl(6)
64 v2sigmatau(6)
65
66 * kxc(N*35):
67 v3rho3(4) = (u_u_u, u_u_d, u_d_d, d_d_d)
68 v3rho2sigma(9) = (u_u_uu, u_u_ud, u_u_dd, u_d_uu, u_d_ud, u_d_dd, d_d_uu, d_d_ud, d_d_dd)
69 v3rhosigma2(12) = (u_uu_uu, u_uu_ud, u_uu_dd, u_ud_ud, u_ud_dd, u_dd_dd, d_uu_uu, d_uu_ud, d_uu_dd, d_ud_ud, d_ud_dd, d_dd_dd)
70 v3sigma(10) = (uu_uu_uu, uu_uu_ud, uu_uu_dd, uu_ud_ud, uu_ud_dd, uu_dd_dd, ud_ud_ud, ud_ud_dd, ud_dd_dd, dd_dd_dd)
71
72 */
73 /*
74 * rho_u/rho_d = (den,grad_x,grad_y,grad_z,laplacian,tau)
75 * In spin restricted case (spin == 1), rho_u is assumed to be the
76 * spin-free quantities, rho_d is not used.
77 */
_eval_xc(xc_func_type * func_x,int spin,int np,double * rho_u,double * rho_d,double * ex,double * vxc,double * fxc,double * kxc)78 static void _eval_xc(xc_func_type *func_x, int spin, int np,
79 double *rho_u, double *rho_d,
80 double *ex, double *vxc, double *fxc, double *kxc)
81 {
82 int i;
83 double *rho, *sigma, *lapl, *tau;
84 double *gxu, *gyu, *gzu, *gxd, *gyd, *gzd;
85 double *lapl_u, *lapl_d, *tau_u, *tau_d;
86 double *vsigma = NULL;
87 double *vlapl = NULL;
88 double *vtau = NULL;
89 double *v2rhosigma = NULL;
90 double *v2sigma2 = NULL;
91 double *v2lapl2 = NULL;
92 double *v2tau2 = NULL;
93 double *v2rholapl = NULL;
94 double *v2rhotau = NULL;
95 double *v2sigmalapl = NULL;
96 double *v2sigmatau = NULL;
97 double *v2lapltau = NULL;
98 double *v3rho2sigma = NULL;
99 double *v3rhosigma2 = NULL;
100 double *v3sigma3 = NULL;
101
102 switch (func_x->info->family) {
103 case XC_FAMILY_LDA:
104 #ifdef XC_FAMILY_HYB_LDA
105 case XC_FAMILY_HYB_LDA:
106 #endif
107 // ex is the energy density
108 // NOTE libxc library added ex/ec into vrho/vcrho
109 // vrho = rho d ex/d rho + ex, see work_lda.c:L73
110 if (spin == XC_POLARIZED) {
111 rho = malloc(sizeof(double) * np*2);
112 for (i = 0; i < np; i++) {
113 rho[i*2+0] = rho_u[i];
114 rho[i*2+1] = rho_d[i];
115 }
116 #if XC_MAJOR_VERSION >= 5
117 xc_lda_exc_vxc_fxc_kxc(func_x, np, rho, ex, vxc, fxc, kxc);
118 #else
119 xc_lda(func_x, np, rho, ex, vxc, fxc, kxc);
120 #endif
121 free(rho);
122 } else {
123 rho = rho_u;
124 #if XC_MAJOR_VERSION >= 5
125 xc_lda_exc_vxc_fxc_kxc(func_x, np, rho, ex, vxc, fxc, kxc);
126 #else
127 xc_lda(func_x, np, rho, ex, vxc, fxc, kxc);
128 #endif
129 }
130 break;
131 case XC_FAMILY_GGA:
132 #ifdef XC_FAMILY_HYB_GGA
133 case XC_FAMILY_HYB_GGA:
134 #endif
135 if (spin == XC_POLARIZED) {
136 rho = malloc(sizeof(double) * np*2);
137 sigma = malloc(sizeof(double) * np*3);
138 gxu = rho_u + np;
139 gyu = rho_u + np * 2;
140 gzu = rho_u + np * 3;
141 gxd = rho_d + np;
142 gyd = rho_d + np * 2;
143 gzd = rho_d + np * 3;
144 for (i = 0; i < np; i++) {
145 rho[i*2+0] = rho_u[i];
146 rho[i*2+1] = rho_d[i];
147 sigma[i*3+0] = gxu[i]*gxu[i] + gyu[i]*gyu[i] + gzu[i]*gzu[i];
148 sigma[i*3+1] = gxu[i]*gxd[i] + gyu[i]*gyd[i] + gzu[i]*gzd[i];
149 sigma[i*3+2] = gxd[i]*gxd[i] + gyd[i]*gyd[i] + gzd[i]*gzd[i];
150 }
151 if (vxc != NULL) {
152 // vrho = vxc
153 vsigma = vxc + np * 2;
154 }
155 if (fxc != NULL) {
156 // v2rho2 = fxc
157 v2rhosigma = fxc + np * 3;
158 v2sigma2 = v2rhosigma + np * 6; // np*6
159 }
160 if (kxc != NULL) {
161 // v3rho3 = kxc
162 v3rho2sigma = kxc + np * 4;
163 v3rhosigma2 = v3rho2sigma + np * 9;
164 v3sigma3 = v3rhosigma2 + np * 12; // np*10
165 }
166 #if (XC_MAJOR_VERSION == 2 && XC_MINOR_VERSION < 2)
167 xc_gga(func_x, np, rho, sigma, ex,
168 vxc, vsigma, fxc, v2rhosigma, v2sigma2);
169 #elif XC_MAJOR_VERSION < 5
170 xc_gga(func_x, np, rho, sigma, ex,
171 vxc, vsigma, fxc, v2rhosigma, v2sigma2,
172 kxc, v3rho2sigma, v3rhosigma2, v3sigma3);
173 #else
174 xc_gga_exc_vxc_fxc_kxc(func_x, np, rho, sigma, ex,
175 vxc, vsigma, fxc, v2rhosigma, v2sigma2,
176 kxc, v3rho2sigma, v3rhosigma2, v3sigma3);
177 #endif
178 free(rho);
179 free(sigma);
180 } else {
181 rho = rho_u;
182 sigma = malloc(sizeof(double) * np);
183 gxu = rho_u + np;
184 gyu = rho_u + np * 2;
185 gzu = rho_u + np * 3;
186 for (i = 0; i < np; i++) {
187 sigma[i] = gxu[i]*gxu[i] + gyu[i]*gyu[i] + gzu[i]*gzu[i];
188 }
189 if (vxc != NULL) {
190 vsigma = vxc + np;
191 }
192 if (fxc != NULL) {
193 v2rhosigma = fxc + np;
194 v2sigma2 = v2rhosigma + np;
195 }
196 if (kxc != NULL) {
197 v3rho2sigma = kxc + np;
198 v3rhosigma2 = v3rho2sigma + np;
199 v3sigma3 = v3rhosigma2 + np;
200 }
201 #if (XC_MAJOR_VERSION == 2 && XC_MINOR_VERSION < 2)
202 xc_gga(func_x, np, rho, sigma, ex,
203 vxc, vsigma, fxc, v2rhosigma, v2sigma2);
204 #elif XC_MAJOR_VERSION < 5
205 xc_gga(func_x, np, rho, sigma, ex,
206 vxc, vsigma, fxc, v2rhosigma, v2sigma2,
207 kxc, v3rho2sigma, v3rhosigma2, v3sigma3);
208 #else
209 xc_gga_exc_vxc_fxc_kxc(func_x, np, rho, sigma, ex,
210 vxc, vsigma, fxc, v2rhosigma, v2sigma2,
211 kxc, v3rho2sigma, v3rhosigma2, v3sigma3);
212
213 #endif
214 free(sigma);
215 }
216 break;
217 case XC_FAMILY_MGGA:
218 #ifdef XC_FAMILY_HYB_MGGA
219 case XC_FAMILY_HYB_MGGA:
220 #endif
221 if (spin == XC_POLARIZED) {
222 rho = malloc(sizeof(double) * np*2);
223 sigma = malloc(sizeof(double) * np*3);
224 lapl = malloc(sizeof(double) * np*2);
225 tau = malloc(sizeof(double) * np*2);
226 gxu = rho_u + np;
227 gyu = rho_u + np * 2;
228 gzu = rho_u + np * 3;
229 gxd = rho_d + np;
230 gyd = rho_d + np * 2;
231 gzd = rho_d + np * 3;
232 lapl_u = rho_u + np * 4;
233 tau_u = rho_u + np * 5;
234 lapl_d = rho_d + np * 4;
235 tau_d = rho_d + np * 5;
236 for (i = 0; i < np; i++) {
237 rho[i*2+0] = rho_u[i];
238 rho[i*2+1] = rho_d[i];
239 sigma[i*3+0] = gxu[i]*gxu[i] + gyu[i]*gyu[i] + gzu[i]*gzu[i];
240 sigma[i*3+1] = gxu[i]*gxd[i] + gyu[i]*gyd[i] + gzu[i]*gzd[i];
241 sigma[i*3+2] = gxd[i]*gxd[i] + gyd[i]*gyd[i] + gzd[i]*gzd[i];
242 lapl[i*2+0] = lapl_u[i];
243 lapl[i*2+1] = lapl_d[i];
244 tau[i*2+0] = tau_u[i];
245 tau[i*2+1] = tau_d[i];
246 }
247 if (vxc != NULL) {
248 // vrho = vxc
249 vsigma = vxc + np * 2;
250 vlapl = vsigma + np * 3;
251 vtau = vlapl + np * 2; // np*2
252 }
253 if (fxc != NULL) {
254 // v2rho2 = fxc
255 v2rhosigma = fxc + np * 3;
256 v2sigma2 = v2rhosigma + np * 6;
257 v2lapl2 = v2sigma2 + np * 6;
258 v2tau2 = v2lapl2 + np * 3;
259 v2rholapl = v2tau2 + np * 3;
260 v2rhotau = v2rholapl + np * 4;
261 v2lapltau = v2rhotau + np * 4;
262 v2sigmalapl = v2lapltau + np * 4;
263 v2sigmatau = v2sigmalapl + np * 6; // np*6
264 }
265 #if XC_MAJOR_VERSION >=5
266 xc_mgga_exc_vxc_fxc(func_x, np, rho, sigma, lapl, tau, ex,
267 vxc, vsigma, vlapl, vtau,
268 fxc, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl,
269 v2rhotau, v2sigmalapl, v2sigmatau, v2lapltau);
270 #else
271 xc_mgga(func_x, np, rho, sigma, lapl, tau, ex,
272 vxc, vsigma, vlapl, vtau,
273 fxc, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl,
274 v2rhotau, v2sigmalapl, v2sigmatau, v2lapltau);
275 #endif
276 free(rho);
277 free(sigma);
278 free(lapl);
279 free(tau);
280 } else {
281 rho = rho_u;
282 sigma = malloc(sizeof(double) * np);
283 lapl = rho_u + np * 4;
284 tau = rho_u + np * 5;
285 gxu = rho_u + np;
286 gyu = rho_u + np * 2;
287 gzu = rho_u + np * 3;
288 for (i = 0; i < np; i++) {
289 sigma[i] = gxu[i]*gxu[i] + gyu[i]*gyu[i] + gzu[i]*gzu[i];
290 }
291 if (vxc != NULL) {
292 vsigma = vxc + np;
293 vlapl = vsigma + np;
294 vtau = vlapl + np;
295 }
296 if (fxc != NULL) {
297 v2rhosigma = fxc + np;
298 v2sigma2 = v2rhosigma + np;
299 v2lapl2 = v2sigma2 + np;
300 v2tau2 = v2lapl2 + np;
301 v2rholapl = v2tau2 + np;
302 v2rhotau = v2rholapl + np;
303 v2lapltau = v2rhotau + np;
304 v2sigmalapl = v2lapltau + np;
305 v2sigmatau = v2sigmalapl + np;
306 }
307 #if XC_MAJOR_VERSION >= 5
308 xc_mgga_exc_vxc_fxc(func_x, np, rho, sigma, lapl, tau, ex,
309 vxc, vsigma, vlapl, vtau,
310 fxc, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl,
311 v2rhotau, v2sigmalapl, v2sigmatau, v2lapltau);
312 #else
313 xc_mgga(func_x, np, rho, sigma, lapl, tau, ex,
314 vxc, vsigma, vlapl, vtau,
315 fxc, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl,
316 v2rhotau, v2sigmalapl, v2sigmatau, v2lapltau);
317 #endif
318 free(sigma);
319 }
320 break;
321 default:
322 fprintf(stderr, "functional %d '%s' is not implmented\n",
323 func_x->info->number, func_x->info->name);
324 exit(1);
325 }
326 }
327
LIBXC_is_lda(int xc_id)328 int LIBXC_is_lda(int xc_id)
329 {
330 xc_func_type func;
331 int lda;
332 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
333 fprintf(stderr, "XC functional %d not found\n", xc_id);
334 exit(1);
335 }
336 switch(func.info->family)
337 {
338 case XC_FAMILY_LDA:
339 lda = 1;
340 break;
341 default:
342 lda = 0;
343 }
344
345 xc_func_end(&func);
346 return lda;
347 }
348
LIBXC_is_gga(int xc_id)349 int LIBXC_is_gga(int xc_id)
350 {
351 xc_func_type func;
352 int gga;
353 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
354 fprintf(stderr, "XC functional %d not found\n", xc_id);
355 exit(1);
356 }
357 switch(func.info->family)
358 {
359 case XC_FAMILY_GGA:
360 #ifdef XC_FAMILY_HYB_GGA
361 case XC_FAMILY_HYB_GGA:
362 #endif
363 gga = 1;
364 break;
365 default:
366 gga = 0;
367 }
368
369 xc_func_end(&func);
370 return gga;
371 }
372
LIBXC_is_meta_gga(int xc_id)373 int LIBXC_is_meta_gga(int xc_id)
374 {
375 xc_func_type func;
376 int mgga;
377 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
378 fprintf(stderr, "XC functional %d not found\n", xc_id);
379 exit(1);
380 }
381 switch(func.info->family)
382 {
383 case XC_FAMILY_MGGA:
384 #ifdef XC_FAMILY_HYB_MGGA
385 case XC_FAMILY_HYB_MGGA:
386 #endif
387 mgga = 1;
388 break;
389 default:
390 mgga = 0;
391 }
392
393 xc_func_end(&func);
394 return mgga;
395 }
396
LIBXC_needs_laplacian(int xc_id)397 int LIBXC_needs_laplacian(int xc_id)
398 {
399 xc_func_type func;
400 int lapl;
401 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
402 fprintf(stderr, "XC functional %d not found\n", xc_id);
403 exit(1);
404 }
405 lapl = func.info->flags & XC_FLAGS_NEEDS_LAPLACIAN ? 1 : 0;
406 xc_func_end(&func);
407 return lapl;
408 }
409
LIBXC_is_hybrid(int xc_id)410 int LIBXC_is_hybrid(int xc_id)
411 {
412 xc_func_type func;
413 int hyb;
414 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
415 fprintf(stderr, "XC functional %d not found\n", xc_id);
416 exit(1);
417 }
418
419 #if XC_MAJOR_VERSION < 6
420 switch(func.info->family)
421 {
422 #ifdef XC_FAMILY_HYB_LDA
423 case XC_FAMILY_HYB_LDA:
424 #endif
425 case XC_FAMILY_HYB_GGA:
426 case XC_FAMILY_HYB_MGGA:
427 hyb = 1;
428 break;
429 default:
430 hyb = 0;
431 }
432 #else
433 hyb = (xc_hyb_type(&func) == XC_HYB_HYBRID);
434 #endif
435
436 xc_func_end(&func);
437 return hyb;
438 }
439
LIBXC_hybrid_coeff(int xc_id)440 double LIBXC_hybrid_coeff(int xc_id)
441 {
442 xc_func_type func;
443 double factor;
444 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
445 fprintf(stderr, "XC functional %d not found\n", xc_id);
446 exit(1);
447 }
448
449 #if XC_MAJOR_VERSION < 6
450 switch(func.info->family)
451 {
452 #ifdef XC_FAMILY_HYB_LDA
453 case XC_FAMILY_HYB_LDA:
454 #endif
455 case XC_FAMILY_HYB_GGA:
456 case XC_FAMILY_HYB_MGGA:
457 factor = xc_hyb_exx_coef(&func);
458 break;
459 default:
460 factor = 0;
461 }
462
463 #else
464 if(xc_hyb_type(&func) == XC_HYB_HYBRID)
465 factor = xc_hyb_exx_coef(&func);
466 else
467 factor = 0.0;
468 #endif
469
470 xc_func_end(&func);
471 return factor;
472 }
473
LIBXC_nlc_coeff(int xc_id,double * nlc_pars)474 void LIBXC_nlc_coeff(int xc_id, double *nlc_pars) {
475
476 xc_func_type func;
477 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
478 fprintf(stderr, "XC functional %d not found\n", xc_id);
479 exit(1);
480 }
481 XC(nlc_coef)(&func, &nlc_pars[0], &nlc_pars[1]);
482 xc_func_end(&func);
483 }
484
LIBXC_rsh_coeff(int xc_id,double * rsh_pars)485 void LIBXC_rsh_coeff(int xc_id, double *rsh_pars) {
486
487 xc_func_type func;
488 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
489 fprintf(stderr, "XC functional %d not found\n", xc_id);
490 exit(1);
491 }
492 rsh_pars[0] = 0.0;
493 rsh_pars[1] = 0.0;
494 rsh_pars[2] = 0.0;
495
496 #if XC_MAJOR_VERSION < 6
497 XC(hyb_cam_coef)(&func, &rsh_pars[0], &rsh_pars[1], &rsh_pars[2]);
498 #else
499 switch(xc_hyb_type(&func)) {
500 case(XC_HYB_HYBRID):
501 case(XC_HYB_CAM):
502 XC(hyb_cam_coef)(&func, &rsh_pars[0], &rsh_pars[1], &rsh_pars[2]);
503 }
504 #endif
505 xc_func_end(&func);
506 }
507
LIBXC_is_cam_rsh(int xc_id)508 int LIBXC_is_cam_rsh(int xc_id) {
509 xc_func_type func;
510 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
511 fprintf(stderr, "XC functional %d not found\n", xc_id);
512 exit(1);
513 }
514 #if XC_MAJOR_VERSION < 6
515 int is_cam = func.info->flags & XC_FLAGS_HYB_CAM;
516 #else
517 int is_cam = (xc_hyb_type(&func) == XC_HYB_CAM);
518 #endif
519 xc_func_end(&func);
520 return is_cam;
521 }
522
523 /*
524 * XC_FAMILY_LDA 1
525 * XC_FAMILY_GGA 2
526 * XC_FAMILY_MGGA 4
527 * XC_FAMILY_LCA 8
528 * XC_FAMILY_OEP 16
529 * XC_FAMILY_HYB_GGA 32
530 * XC_FAMILY_HYB_MGGA 64
531 * XC_FAMILY_HYB_LDA 128
532 */
LIBXC_xc_type(int fn_id)533 int LIBXC_xc_type(int fn_id)
534 {
535 xc_func_type func;
536 if (xc_func_init(&func, fn_id, 1) != 0) {
537 fprintf(stderr, "XC functional %d not found\n", fn_id);
538 exit(1);
539 }
540 int type = func.info->family;
541 xc_func_end(&func);
542 return type;
543 }
544
xc_output_length(int nvar,int deriv)545 static int xc_output_length(int nvar, int deriv)
546 {
547 int i;
548 int len = 1.;
549 for (i = 1; i <= nvar; i++) {
550 len *= deriv + i;
551 len /= i;
552 }
553 return len;
554 }
555
556 // return value 0 means no functional needs to be evaluated.
LIBXC_input_length(int nfn,int * fn_id,double * fac,int spin)557 int LIBXC_input_length(int nfn, int *fn_id, double *fac, int spin)
558 {
559 int i;
560 int nvar = 0;
561 xc_func_type func;
562 for (i = 0; i < nfn; i++) {
563 if (xc_func_init(&func, fn_id[i], spin) != 0) {
564 fprintf(stderr, "XC functional %d not found\n",
565 fn_id[i]);
566 exit(1);
567 }
568 if (spin == XC_POLARIZED) {
569 switch (func.info->family) {
570 case XC_FAMILY_LDA:
571 #ifdef XC_FAMILY_HYB_LDA
572 case XC_FAMILY_HYB_LDA:
573 #endif
574 nvar = MAX(nvar, 2);
575 break;
576 case XC_FAMILY_GGA:
577 #ifdef XC_FAMILY_HYB_GGA
578 case XC_FAMILY_HYB_GGA:
579 #endif
580 nvar = MAX(nvar, 5);
581 break;
582 case XC_FAMILY_MGGA:
583 #ifdef XC_FAMILY_HYB_MGGA
584 case XC_FAMILY_HYB_MGGA:
585 #endif
586 nvar = MAX(nvar, 9);
587 }
588 } else {
589 switch (func.info->family) {
590 case XC_FAMILY_LDA:
591 #ifdef XC_FAMILY_HYB_LDA
592 case XC_FAMILY_HYB_LDA:
593 #endif
594 nvar = MAX(nvar, 1);
595 break;
596 case XC_FAMILY_GGA:
597 #ifdef XC_FAMILY_HYB_GGA
598 case XC_FAMILY_HYB_GGA:
599 #endif
600 nvar = MAX(nvar, 2);
601 break;
602 case XC_FAMILY_MGGA:
603 #ifdef XC_FAMILY_HYB_MGGA
604 case XC_FAMILY_HYB_MGGA:
605 #endif
606 nvar = MAX(nvar, 4);
607 }
608 }
609 xc_func_end(&func);
610 }
611 return nvar;
612 }
613
axpy(double * dst,double * src,double fac,int np,int ndst,int nsrc)614 static void axpy(double *dst, double *src, double fac,
615 int np, int ndst, int nsrc)
616 {
617 int i, j;
618 for (j = 0; j < nsrc; j++) {
619 for (i = 0; i < np; i++) {
620 dst[j*np+i] += fac * src[i*nsrc+j];
621 }
622 }
623 }
624
merge_xc(double * dst,double * ebuf,double * vbuf,double * fbuf,double * kbuf,double fac,int np,int ndst,int nvar,int spin,int type)625 static void merge_xc(double *dst, double *ebuf, double *vbuf,
626 double *fbuf, double *kbuf, double fac,
627 int np, int ndst, int nvar, int spin, int type)
628 {
629 const int seg0 [] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
630 // LDA | |
631 // GGA | |
632 // MGGA | |
633 const int vseg1[] = {2, 3, 2, 2};
634 // LDA | |
635 // GGA | |
636 // MGGA | |
637 const int fseg1[] = {3, 6, 6, 3, 3, 4, 4, 4, 6, 6};
638 // LDA | |
639 // GGA | |
640 const int kseg1[] = {4, 9,12,10};
641 int vsegtot, fsegtot, ksegtot;
642 const int *vseg, *fseg, *kseg;
643 if (spin == XC_POLARIZED) {
644 vseg = vseg1;
645 fseg = fseg1;
646 kseg = kseg1;
647 } else {
648 vseg = seg0;
649 fseg = seg0;
650 kseg = seg0;
651 }
652
653 switch (type) {
654 case XC_FAMILY_GGA:
655 #ifdef XC_FAMILY_HYB_GGA
656 case XC_FAMILY_HYB_GGA:
657 #endif
658 vsegtot = 2;
659 fsegtot = 3;
660 ksegtot = 4;
661 break;
662 case XC_FAMILY_MGGA:
663 #ifdef XC_FAMILY_HYB_MGGA
664 case XC_FAMILY_HYB_MGGA:
665 #endif
666 vsegtot = 4;
667 fsegtot = 10;
668 ksegtot = 0; // not supported
669 break;
670 default: //case XC_FAMILY_LDA:
671 vsegtot = 1;
672 fsegtot = 1;
673 ksegtot = 1;
674 }
675
676 int i;
677 size_t offset;
678 axpy(dst, ebuf, fac, np, ndst, 1);
679
680 if (vbuf != NULL) {
681 offset = np;
682 for (i = 0; i < vsegtot; i++) {
683 axpy(dst+offset, vbuf, fac, np, ndst, vseg[i]);
684 offset += np * vseg[i];
685 vbuf += np * vseg[i];
686 }
687 }
688
689 if (fbuf != NULL) {
690 offset = np * xc_output_length(nvar, 1);
691 for (i = 0; i < fsegtot; i++) {
692 axpy(dst+offset, fbuf, fac, np, ndst, fseg[i]);
693 offset += np * fseg[i];
694 fbuf += np * fseg[i];
695 }
696 }
697
698 if (kbuf != NULL) {
699 offset = np * xc_output_length(nvar, 2);
700 for (i = 0; i < ksegtot; i++) {
701 axpy(dst+offset, kbuf, fac, np, ndst, kseg[i]);
702 offset += np * kseg[i];
703 kbuf += np * kseg[i];
704 }
705 }
706 }
707
708 // omega is the range separation parameter mu in xcfun
LIBXC_eval_xc(int nfn,int * fn_id,double * fac,double * omega,int spin,int deriv,int np,double * rho_u,double * rho_d,double * output)709 void LIBXC_eval_xc(int nfn, int *fn_id, double *fac, double *omega,
710 int spin, int deriv, int np,
711 double *rho_u, double *rho_d, double *output)
712 {
713 assert(deriv <= 3);
714 int nvar = LIBXC_input_length(nfn, fn_id, fac, spin);
715 if (nvar == 0) { // No functional needs to be evaluated.
716 return;
717 }
718
719 int outlen = xc_output_length(nvar, deriv);
720 // output buffer is zeroed in the Python caller
721 //NPdset0(output, np*outlen);
722
723 double *ebuf = malloc(sizeof(double) * np);
724 double *vbuf = NULL;
725 double *fbuf = NULL;
726 double *kbuf = NULL;
727 if (deriv > 0) {
728 vbuf = malloc(sizeof(double) * np*9);
729 }
730 if (deriv > 1) {
731 fbuf = malloc(sizeof(double) * np*48);
732 }
733 if (deriv > 2) { // *220 if mgga kxc available
734 kbuf = malloc(sizeof(double) * np*35);
735 }
736
737 int i, j;
738 xc_func_type func;
739 for (i = 0; i < nfn; i++) {
740 if (xc_func_init(&func, fn_id[i], spin) != 0) {
741 fprintf(stderr, "XC functional %d not found\n",
742 fn_id[i]);
743 exit(1);
744 }
745
746 // set the range-separated parameter
747 if (omega[i] != 0) {
748 // skip if func is not a RSH functional
749 #if XC_MAJOR_VERSION < 6
750 if (func.cam_omega != 0) {
751 func.cam_omega = omega[i];
752 }
753 #else
754 if (func.hyb_omega[0] != 0) {
755 func.hyb_omega[0] = omega[i];
756 }
757 #endif
758 // Recursively set the sub-functionals if they are RSH
759 // functionals
760 for (j = 0; j < func.n_func_aux; j++) {
761 #if XC_MAJOR_VERSION < 6
762 if (func.func_aux[j]->cam_omega != 0) {
763 func.func_aux[j]->cam_omega = omega[i];
764 }
765 #else
766 if (func.func_aux[j]->hyb_omega[0] != 0) {
767 func.func_aux[j]->hyb_omega[0] = omega[i];
768 }
769 #endif
770 }
771 }
772
773 // alpha and beta are hardcoded in many functionals in the libxc
774 // code, e.g. the coefficients of B88 (=1-alpha) and
775 // ITYH (=-beta) in cam-b3lyp. Overwriting func->cam_alpha and
776 // func->cam_beta does not update the coefficients accordingly.
777 //func->cam_alpha = alpha;
778 //func->cam_beta = beta;
779 // However, the parameters can be set with the libxc function
780 //void xc_func_set_ext_params_name(xc_func_type *p, const char *name, double par);
781 // since libxc 5.1.0
782 #if defined XC_SET_RELATIVITY
783 xc_lda_x_set_params(&func, relativity);
784 #endif
785 _eval_xc(&func, spin, np, rho_u, rho_d, ebuf, vbuf, fbuf, kbuf);
786 merge_xc(output, ebuf, vbuf, fbuf, kbuf, fac[i],
787 np, outlen, nvar, spin, func.info->family);
788 xc_func_end(&func);
789 }
790
791 free(ebuf);
792 if (deriv > 0) {
793 free(vbuf);
794 }
795 if (deriv > 1) {
796 free(fbuf);
797 }
798 if (deriv > 2) {
799 free(kbuf);
800 }
801 }
802
LIBXC_max_deriv_order(int xc_id)803 int LIBXC_max_deriv_order(int xc_id)
804 {
805 xc_func_type func;
806 int ord;
807 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
808 fprintf(stderr, "XC functional %d not found\n", xc_id);
809 exit(1);
810 }
811
812 if (func.info->flags & XC_FLAGS_HAVE_LXC) {
813 ord = 4;
814 } else if(func.info->flags & XC_FLAGS_HAVE_KXC) {
815 ord = 3;
816 } else if(func.info->flags & XC_FLAGS_HAVE_FXC) {
817 ord = 2;
818 } else if(func.info->flags & XC_FLAGS_HAVE_VXC) {
819 ord = 1;
820 } else if(func.info->flags & XC_FLAGS_HAVE_EXC) {
821 ord = 0;
822 } else {
823 ord = -1;
824 }
825
826 xc_func_end(&func);
827 return ord;
828 }
829
LIBXC_number_of_functionals()830 int LIBXC_number_of_functionals()
831 {
832 return xc_number_of_functionals();
833 }
834
LIBXC_functional_numbers(int * list)835 void LIBXC_functional_numbers(int *list)
836 {
837 return xc_available_functional_numbers(list);
838 }
839
LIBXC_functional_name(int ifunc)840 char * LIBXC_functional_name(int ifunc)
841 {
842 return xc_functional_get_name(ifunc);
843 }
844
LIBXC_version()845 const char * LIBXC_version()
846 {
847 return xc_version_string();
848 }
849
LIBXC_reference()850 const char * LIBXC_reference()
851 {
852 return xc_reference();
853 }
854
LIBXC_reference_doi()855 const char * LIBXC_reference_doi()
856 {
857 return xc_reference_doi();
858 }
859
LIBXC_xc_reference(int xc_id,const char ** refs)860 void LIBXC_xc_reference(int xc_id, const char **refs)
861 {
862 xc_func_type func;
863 if(xc_func_init(&func, xc_id, XC_UNPOLARIZED) != 0){
864 fprintf(stderr, "XC functional %d not found\n", xc_id);
865 exit(1);
866 }
867
868 int i;
869 for (i = 0; i < XC_MAX_REFERENCES; i++) {
870 if (func.info->refs[i] == NULL || func.info->refs[i]->ref == NULL) {
871 refs[i] = NULL;
872 break;
873 }
874 refs[i] = func.info->refs[i]->ref;
875 }
876 }
877