1 /*
2 * Copyright (C) 2015-2020 ABINIT group (MT)
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation.
7 * See http://www.gnu.org/copyleft/gpl.txt.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
12 *
13 */
14
15 /* ===============================================================
16 * Set of C functions interfacing the LibXC library.
17 * (see http://www.tddft.org/programs/Libxc)
18 * ===============================================================
19 */
20
21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24
25 #include <stdlib.h>
26 #include <stdio.h>
27
28 #if defined HAVE_LIBXC
29
30 #include "xc.h"
31 #include "xc_version.h"
32 #include "xc_funcs.h"
33
34 /* if version before 4 get config file */
35 #if ( XC_MAJOR_VERSION < 4 )
36 # include "xc_config.h"
37 #else
38 # define FLOAT double
39 #endif
40
41 /* ===============================================================
42 * Get the SINGLE_PRECISION constant
43 * ===============================================================
44 */
xc_get_singleprecision_constant(int * xc_cst_singleprecision)45 void xc_get_singleprecision_constant(int *xc_cst_singleprecision)
46 {
47 if (sizeof(FLOAT)<sizeof(double))
48 {*xc_cst_singleprecision = 1;}
49 else
50 {*xc_cst_singleprecision = 0;}
51 }
52
53 /* ===============================================================
54 * Get the FAMILY constants
55 * ===============================================================
56 */
xc_get_family_constants(int * xc_cst_family_unknown,int * xc_cst_family_lda,int * xc_cst_family_gga,int * xc_cst_family_mgga,int * xc_cst_family_lca,int * xc_cst_family_oep,int * xc_cst_family_hyb_gga,int * xc_cst_family_hyb_mgga)57 void xc_get_family_constants(int *xc_cst_family_unknown,
58 int *xc_cst_family_lda,
59 int *xc_cst_family_gga,
60 int *xc_cst_family_mgga,
61 int *xc_cst_family_lca,
62 int *xc_cst_family_oep,
63 int *xc_cst_family_hyb_gga,
64 int *xc_cst_family_hyb_mgga)
65 {
66 *xc_cst_family_unknown = XC_FAMILY_UNKNOWN;
67 *xc_cst_family_lda = XC_FAMILY_LDA;
68 *xc_cst_family_gga = XC_FAMILY_GGA;
69 *xc_cst_family_mgga = XC_FAMILY_MGGA;
70 *xc_cst_family_lca = XC_FAMILY_LCA;
71 *xc_cst_family_oep = XC_FAMILY_OEP;
72 #if ( XC_MAJOR_VERSION > 5 )
73 /* ==== libXC v6.0 and later ==== */
74 *xc_cst_family_hyb_gga = -11;
75 *xc_cst_family_hyb_mgga = -11;
76 #else
77 /* ==== Before libXC v6.0 ==== */
78 *xc_cst_family_hyb_gga = XC_FAMILY_HYB_GGA;
79 *xc_cst_family_hyb_mgga = XC_FAMILY_HYB_MGGA;
80 #endif
81 }
82
83 /* ===============================================================
84 * Get the FLAGS constants
85 * ===============================================================
86 */
xc_get_flags_constants(int * xc_cst_flags_have_exc,int * xc_cst_flags_have_vxc,int * xc_cst_flags_have_fxc,int * xc_cst_flags_have_kxc,int * xc_cst_flags_have_lxc,int * xc_cst_flags_needs_laplacian)87 void xc_get_flags_constants(int *xc_cst_flags_have_exc,
88 int *xc_cst_flags_have_vxc,
89 int *xc_cst_flags_have_fxc,
90 int *xc_cst_flags_have_kxc,
91 int *xc_cst_flags_have_lxc,
92 int *xc_cst_flags_needs_laplacian)
93 {
94 *xc_cst_flags_have_exc = XC_FLAGS_HAVE_EXC;
95 *xc_cst_flags_have_vxc = XC_FLAGS_HAVE_VXC;
96 *xc_cst_flags_have_fxc = XC_FLAGS_HAVE_FXC;
97 *xc_cst_flags_have_kxc = XC_FLAGS_HAVE_KXC;
98 *xc_cst_flags_have_lxc = XC_FLAGS_HAVE_LXC;
99 #if ( XC_MAJOR_VERSION > 3 )
100 *xc_cst_flags_needs_laplacian = XC_FLAGS_NEEDS_LAPLACIAN;
101 #else
102 *xc_cst_flags_needs_laplacian = 0;
103 #endif
104 }
105
106 /* ===============================================================
107 * Get the KIND constants
108 * ===============================================================
109 */
xc_get_kind_constants(int * xc_cst_exchange,int * xc_cst_correlation,int * xc_cst_exchange_correlation,int * xc_cst_kinetic)110 void xc_get_kind_constants(int *xc_cst_exchange,
111 int *xc_cst_correlation,
112 int *xc_cst_exchange_correlation,
113 int *xc_cst_kinetic)
114 {
115 *xc_cst_exchange = XC_EXCHANGE;
116 *xc_cst_correlation = XC_CORRELATION;
117 *xc_cst_exchange_correlation = XC_EXCHANGE_CORRELATION;
118 *xc_cst_kinetic = XC_KINETIC;
119 }
120
121 /* ===============================================================
122 * Get the HYBRID constants
123 * ===============================================================
124 */
xc_get_hybrid_constants(int * xc_cst_hyb_none,int * xc_cst_hyb_fock,int * xc_cst_hyb_pt2,int * xc_cst_hyb_erf_sr,int * xc_cst_hyb_yukawa_sr,int * xc_cst_hyb_gaussian_sr,int * xc_cst_hyb_semilocal,int * xc_cst_hyb_hybrid,int * xc_cst_hyb_cam,int * xc_cst_hyb_camy,int * xc_cst_hyb_camg,int * xc_cst_hyb_double_hybrid,int * xc_cst_hyb_mixture)125 void xc_get_hybrid_constants(int *xc_cst_hyb_none,
126 int *xc_cst_hyb_fock,
127 int *xc_cst_hyb_pt2,
128 int *xc_cst_hyb_erf_sr,
129 int *xc_cst_hyb_yukawa_sr,
130 int *xc_cst_hyb_gaussian_sr,
131 int *xc_cst_hyb_semilocal,
132 int *xc_cst_hyb_hybrid,
133 int *xc_cst_hyb_cam,
134 int *xc_cst_hyb_camy,
135 int *xc_cst_hyb_camg,
136 int *xc_cst_hyb_double_hybrid,
137 int *xc_cst_hyb_mixture)
138 {
139 #if ( XC_MAJOR_VERSION > 5 )
140 /* ==== libXC v6.0 and later ==== */
141 *xc_cst_hyb_none = XC_HYB_NONE;
142 *xc_cst_hyb_fock = XC_HYB_FOCK;
143 *xc_cst_hyb_pt2 = XC_HYB_PT2;
144 *xc_cst_hyb_erf_sr = XC_HYB_ERF_SR;
145 *xc_cst_hyb_yukawa_sr = XC_HYB_YUKAWA_SR;
146 *xc_cst_hyb_gaussian_sr = XC_HYB_GAUSSIAN_SR;
147 *xc_cst_hyb_semilocal = XC_HYB_SEMILOCAL;
148 *xc_cst_hyb_hybrid = XC_HYB_HYBRID;
149 *xc_cst_hyb_cam = XC_HYB_CAM;
150 *xc_cst_hyb_camy = XC_HYB_CAMY;
151 *xc_cst_hyb_camg = XC_HYB_CAMG;
152 *xc_cst_hyb_double_hybrid = XC_HYB_DOUBLE_HYBRID;
153 *xc_cst_hyb_mixture = XC_HYB_MIXTURE;
154 #else
155 /* ==== Before libXC v6.0 ==== */
156 *xc_cst_hyb_none = -11; *xc_cst_hyb_fock = -11;
157 *xc_cst_hyb_pt2 = -11; *xc_cst_hyb_erf_sr = -11;
158 *xc_cst_hyb_yukawa_sr = -11; *xc_cst_hyb_gaussian_sr = -11;
159 *xc_cst_hyb_semilocal = -11; *xc_cst_hyb_hybrid = -11;
160 *xc_cst_hyb_cam = -11; *xc_cst_hyb_camy = -11;
161 *xc_cst_hyb_camg = -11; *xc_cst_hyb_double_hybrid = -11;
162 *xc_cst_hyb_mixture = -11;
163 #endif
164 }
165
166 /* ===============================================================
167 * Allocate/free xc_func_type pointer
168 * ===============================================================
169 */
XC(func_type)170 XC(func_type) * xc_func_type_malloc()
171 {return (XC(func_type) *) malloc(sizeof(XC(func_type)));}
172
xc_func_type_free(XC (func_type)** xc_func)173 void xc_func_type_free(XC(func_type) **xc_func)
174 {free(*xc_func);*xc_func = NULL;}
175
176 /* ===============================================================
177 * Wrappers to the LDA/GGA/MGGA functionals
178 * ===============================================================
179 */
180 /* ---------------------------------------------------------------
181 ----- LDA ----- */
xc_get_lda(const XC (func_type)* xc_func,int np,const double * rho,double * zk,double * vrho,double * v2rho2,double * v3rho3)182 void xc_get_lda(const XC(func_type) *xc_func, int np, const double *rho,
183 double *zk, double *vrho, double *v2rho2, double *v3rho3)
184 #if ( XC_MAJOR_VERSION > 4 )
185 /* ==== libXC v5.0 and later ==== */
186 {xc_lda(xc_func, np, rho, zk, vrho, v2rho2, v3rho3, NULL);}
187 #else
188 /* ==== Before libXC v5.0 ==== */
189 {xc_lda(xc_func, np, rho, zk, vrho, v2rho2, v3rho3);}
190 #endif
191 /* ---------------------------------------------------------------
192 ----- GGA ----- */
xc_get_gga(const XC (func_type)* xc_func,int np,const double * rho,const double * sigma,double * zk,double * vrho,double * vsigma,double * v2rho2,double * v2rhosigma,double * v2sigma2,double * v3rho3,double * v3rho2sigma,double * v3rhosigma2,double * v3sigma3)193 void xc_get_gga(const XC(func_type) *xc_func, int np,
194 const double *rho, const double *sigma,
195 double *zk, double *vrho, double *vsigma,
196 double *v2rho2, double *v2rhosigma, double *v2sigma2,
197 double *v3rho3, double *v3rho2sigma, double *v3rhosigma2, double *v3sigma3)
198 #if ( XC_MAJOR_VERSION > 4 )
199 /* ==== libXC v5.0 and later ==== */
200 {xc_gga(xc_func, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2,
201 v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3,
202 NULL, NULL, NULL, NULL, NULL);}
203 #else
204 /* ==== Before libXC v5.0 ==== */
205 {xc_gga(xc_func, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2,
206 v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3);}
207 #endif
208 /* ---------------------------------------------------------------
209 ----- meta-GGA ----- */
xc_get_mgga(const XC (func_type)* xc_func,int np,const double * rho,const double * sigma,const double * lapl,const double * tau,double * zk,double * vrho,double * vsigma,double * vlapl,double * vtau,double * v2rho2,double * v2rhosigma,double * v2rholapl,double * v2rhotau,double * v2sigma2,double * v2sigmalapl,double * v2sigmatau,double * v2lapl2,double * v2lapltau,double * v2tau2)210 void xc_get_mgga(const XC(func_type) *xc_func, int np,
211 const double *rho, const double *sigma, const double *lapl, const double *tau,
212 double *zk, double *vrho, double *vsigma, double *vlapl, double *vtau,
213 double *v2rho2, double *v2rhosigma, double *v2rholapl, double *v2rhotau,
214 double *v2sigma2, double *v2sigmalapl, double *v2sigmatau, double *v2lapl2,
215 double *v2lapltau, double *v2tau2)
216 #if ( XC_MAJOR_VERSION > 4 )
217 /* ==== libXC v5.0 and later ==== */
218 {xc_mgga(xc_func, np, rho, sigma, lapl, tau, zk, vrho, vsigma, vlapl, vtau,
219 v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2,
220 v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2,
221 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
222 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
223 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
224 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
225 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);}
226 #else
227 /* ==== Before libXC v5.0 ==== */
228 {xc_mgga(xc_func, np, rho, sigma, lapl, tau, zk, vrho, vsigma, vlapl, vtau,
229 v2rho2, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl, v2rhotau,
230 v2sigmalapl, v2sigmatau, v2lapltau);}
231 #endif
232
233 /* ===============================================================
234 * Get properties from a xc_func_info_type pointer
235 * These accessors where not provided before libXC v3
236 * ===============================================================
237 */
238 #if ( XC_MAJOR_VERSION > 3 )
239 /* ==== libXC v4.0 and later ==== */
xc_get_info_name(XC (func_type)* xc_func)240 char const *xc_get_info_name(XC(func_type) *xc_func)
241 {return xc_func_info_get_name(xc_func->info);}
xc_get_info_flags(XC (func_type)* xc_func)242 int xc_get_info_flags(XC(func_type) *xc_func)
243 {return xc_func_info_get_flags(xc_func->info);}
xc_get_info_kind(XC (func_type)* xc_func)244 int xc_get_info_kind(XC(func_type) *xc_func)
245 {return xc_func_info_get_kind(xc_func->info);}
xc_get_info_refs(XC (func_type)* xc_func,const int * number)246 char const *xc_get_info_refs(XC(func_type) *xc_func, const int *number)
247 {if (*number>=0&&*number<XC_MAX_REFERENCES)
248 {if (xc_func_info_get_references(xc_func->info,*number) != NULL)
249 {return xc_func_info_get_references(xc_func->info,*number)->ref;}}
250 else {return NULL;}
251 return NULL;}
252
253 #elif ( XC_MAJOR_VERSION > 2 )
254 /* ==== libXC v3.0 ==== */
xc_get_info_name(XC (func_type)* xc_func)255 char const *xc_get_info_name(XC(func_type) *xc_func)
256 {return xc_func_info_get_name(xc_func->info);}
xc_get_info_flags(XC (func_type)* xc_func)257 int xc_get_info_flags(XC(func_type) *xc_func)
258 {return xc_func_info_get_flags(xc_func->info);}
xc_get_info_kind(XC (func_type)* xc_func)259 int xc_get_info_kind(XC(func_type) *xc_func)
260 {return xc_func_info_get_kind(xc_func->info);}
xc_get_info_refs(XC (func_type)* xc_func,const int * number)261 char const *xc_get_info_refs(XC(func_type) *xc_func, const int *number)
262 {if (*number>=0&&*number<=4)
263 {if (xc_func_info_get_ref(xc_func->info,*number) != NULL)
264 {return xc_func_info_get_ref(xc_func->info,*number);}}
265 else {return NULL;}
266 return NULL;}
267
268 #else
269 /* ==== Before libXC v3.0 ==== */
xc_get_info_name(XC (func_type)* xc_func)270 char const *xc_get_info_name(XC(func_type) *xc_func)
271 {return xc_func->info->name;}
xc_get_info_flags(XC (func_type)* xc_func)272 int xc_get_info_flags(XC(func_type) *xc_func)
273 {return xc_func->info->flags;}
xc_get_info_kind(XC (func_type)* xc_func)274 int xc_get_info_kind(XC(func_type) *xc_func)
275 {return xc_func->info->kind;}
xc_get_info_refs(XC (func_type)* xc_func,const int * number)276 char const *xc_get_info_refs(XC(func_type) *xc_func, const int *number)
277 {if (*number==0) {return xc_func->info->refs;} else {return NULL;}
278 return NULL;}
279 #endif
280
281 /* ===============================================================
282 * Wrapper to xc_func_set_ext_params for backward compatibility
283 * Allows to change the parameters of a XC functional
284 * ===============================================================
285 */
xc_func_set_params(XC (func_type)* xc_func,double * ext_params,int n_ext_params)286 void xc_func_set_params(XC(func_type) *xc_func, double *ext_params, int n_ext_params)
287 #if ( XC_MAJOR_VERSION > 4 )
288 /* ==== libXC v5.0 and later ==== */
289 {if (n_ext_params == xc_func->info->ext_params.n)
290 {XC(func_set_ext_params)(xc_func, ext_params);}
291 #elif ( XC_MAJOR_VERSION > 3 )
292 /* ==== libXC v4.0 ==== */
293 {if (xc_func->info->number == XC_HYB_GGA_XC_PBEH && n_ext_params == 1)
294 /* set_ext_params function is missing for PBE0 */
295 {xc_func->cam_alpha=ext_params[0];xc_func->mix_coef[0]=1.0-ext_params[0];}
296 else if (xc_func->info->number == XC_MGGA_X_TB09 && n_ext_params >= 1)
297 /* XC_MGGA_X_TB09 has only one parameter */
298 {XC(func_set_ext_params)(xc_func, ext_params);}
299 else if (n_ext_params == xc_func->info->n_ext_params)
300 {XC(func_set_ext_params)(xc_func, ext_params);}
301
302 #else
303 /* ==== Before libXC v4.0 ==== */
304 {if (xc_func->info->number == XC_LDA_C_XALPHA && n_ext_params == 1)
305 {XC(lda_c_xalpha_set_params)(xc_func, *ext_params);}
306 else if (xc_func->info->number == XC_MGGA_X_TB09 && n_ext_params >= 1)
307 {XC(mgga_x_tb09_set_params)(xc_func, *ext_params);}
308 #if ( XC_MAJOR_VERSION > 2 || ( XC_MAJOR_VERSION > 1 && XC_MINOR_VERSION > 0 ) )
309 else if (xc_func->info->number == XC_HYB_GGA_XC_PBEH && n_ext_params == 1)
310 {XC(hyb_gga_xc_pbeh_set_params)(xc_func, *ext_params);}
311 else if (xc_func->info->number == XC_HYB_GGA_XC_HSE03 && n_ext_params == 3)
312 {XC(hyb_gga_xc_hse_set_params)(xc_func, ext_params[0], ext_params[2]);
313 xc_func->cam_omega=ext_params[1];}
314 else if (xc_func->info->number == XC_HYB_GGA_XC_HSE06 && n_ext_params == 3)
315 {XC(hyb_gga_xc_hse_set_params)(xc_func, ext_params[0], ext_params[2]);
316 xc_func->cam_omega=ext_params[1];}
317 #else
318 else if (xc_func->info->number == XC_HYB_GGA_XC_HSE03 && n_ext_params == 3)
319 {XC(hyb_gga_xc_hse_set_params)(xc_func, ext_params[2]);
320 xc_func->cam_omega=ext_params[1];}
321 else if (xc_func->info->number == XC_HYB_GGA_XC_HSE06 && n_ext_params == 3)
322 {XC(hyb_gga_xc_hse_set_params)(xc_func, ext_params[2]);
323 xc_func->cam_omega=ext_params[1];}
324 #endif
325 #endif
326 else
327 {fprintf(stderr, "BUG: invalid entry in set_params!\n");abort();}
328 }
329
330 /* ===============================================================
331 * Wrapper to xc_func_set_dens_threshold for backward compatibility
332 * Allows to change the zero-density threshold of a XC functional
333 * Only available from libXC v4
334 * ===============================================================
335 */
336 void xc_func_set_density_threshold(XC(func_type) *xc_func, double *dens_threshold)
337 #if ( XC_MAJOR_VERSION > 3 )
338 /* ==== libXC v4.0 and later ==== */
339 {XC(func_set_dens_threshold)(xc_func, *dens_threshold);}
340 #else
341 {fprintf(stdout, "WARNING: setting density threshold not available for libXC<4.0!\n");}
342 #endif
343
344 /* ===============================================================
345 * Missing function:
346 * Return 1 if the functional is hybrid, from its id
347 * ===============================================================
348 */
349 int xc_func_is_hybrid_from_id(int func_id)
350 #if ( XC_MAJOR_VERSION > 5 )
351 /* ==== libXC v6.0 and later ==== */
352 {xc_func_type func; int result=0;
353 if(xc_func_init(&func,func_id,XC_UNPOLARIZED)==0)
354 {if (func.hyb_number_terms>0)
355 {if (func.hyb_type[0] != XC_HYB_NONE){result=1;}}}
356 xc_func_end(&func);
357 return result;
358 }
359 #else
360 /* ==== Before libXC v6.0 ==== */
361 {int family; family=xc_family_from_id(func_id, NULL, NULL);
362 if (family==XC_FAMILY_HYB_GGA || family==XC_FAMILY_HYB_MGGA)
363 {return 1;}
364 else
365 {return 0;}
366 }
367 #endif
368
369 #endif
370