1 /*
2 Copyright (C) 2006-2007 M.A.L. Marques
3 Copyright (C) 2019 X. Andrade
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 #ifndef _LDA_H
11 #define _LDA_H
12
13 /* These are generic header files that are needed basically everywhere */
14 #include <math.h>
15 #include <float.h>
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <string.h>
19 #include <assert.h>
20
21 #ifdef HAVE_CUDA
22 #include <cuda.h>
23 #endif
24
25 #include "xc.h"
26 #include "xc_funcs_worker.h"
27
28 /* we include the references also */
29 #include "references.h"
30
31 /* need config to figure out what needs to be defined or not */
32 #include "config.h"
33
34 #ifdef HAVE_CUDA
35 #define GPU_FUNCTION __host__ __device__
36 #define CUDA_BLOCK_SIZE 256
37 #else
38 #define GPU_FUNCTION
39 #endif
40
41 /* This takes care of disabling specific derivatives from the info structures */
42 #define XC_FLAGS_I_HAVE_EXC XC_FLAGS_HAVE_EXC
43
44 #ifdef XC_DONT_COMPILE_VXC
45 # define XC_FLAGS_I_HAVE_VXC 0
46 #else
47 # define XC_FLAGS_I_HAVE_VXC XC_FLAGS_HAVE_VXC
48 #endif
49
50 #ifdef XC_DONT_COMPILE_FXC
51 # define XC_FLAGS_I_HAVE_FXC 0
52 #else
53 # define XC_FLAGS_I_HAVE_FXC XC_FLAGS_HAVE_FXC
54 #endif
55
56 #ifdef XC_DONT_COMPILE_KXC
57 # define XC_FLAGS_I_HAVE_KXC 0
58 #else
59 # define XC_FLAGS_I_HAVE_KXC XC_FLAGS_HAVE_KXC
60 #endif
61
62 #ifdef XC_DONT_COMPILE_LXC
63 # define XC_FLAGS_I_HAVE_LXC 0
64 #else
65 # define XC_FLAGS_I_HAVE_LXC XC_FLAGS_HAVE_LXC
66 #endif
67
68 #define XC_FLAGS_I_HAVE_ALL (XC_FLAGS_HAVE_EXC | XC_FLAGS_I_HAVE_VXC | \
69 XC_FLAGS_I_HAVE_FXC | XC_FLAGS_I_HAVE_KXC | \
70 XC_FLAGS_I_HAVE_LXC)
71
72 /* Useful mathematical expressions */
73 #ifndef M_E
74 # define M_E 2.7182818284590452354 /* e */
75 #endif
76 #ifndef M_PI
77 # define M_PI 3.14159265358979323846 /* pi */
78 #endif
79 #ifndef M_SQRT2
80 # define M_SQRT2 1.41421356237309504880 /* sqrt(2) */
81 #endif
82
83 #define POW_2(x) ((x)*(x))
84 #define POW_3(x) ((x)*(x)*(x))
85
86 #define POW_1_2(x) sqrt(x)
87 #define POW_1_4(x) sqrt(sqrt(x))
88 #define POW_3_2(x) ((x)*sqrt(x))
89
90 #ifdef HAVE_CBRT
91 #define CBRT(x) cbrt(x)
92 #define POW_1_3(x) cbrt(x)
93 #define POW_2_3(x) (cbrt(x)*cbrt(x))
94 #define POW_4_3(x) ((x)*cbrt(x))
95 #define POW_5_3(x) ((x)*cbrt(x)*cbrt(x))
96 #define POW_7_3(x) ((x)*(x)*cbrt(x))
97 #else
98 #define CBRT(x) pow((x), 1.0/3.0)
99 #define POW_1_3(x) pow((x), 1.0/3.0)
100 #define POW_2_3(x) pow((x), 2.0/3.0)
101 #define POW_4_3(x) pow((x), 4.0/3.0)
102 #define POW_5_3(x) pow((x), 5.0/3.0)
103 #define POW_7_3(x) pow((x), 7.0/3.0)
104 #endif
105
106 /* this is the piecewise function used in maple */
107 #define my_piecewise3(c, x1, x2) ((c) ? (x1) : (x2))
108 #define my_piecewise5(c1, x1, c2, x2, x3) ((c) ? (x1) : ((c2) ? (x2) : (x3)))
109
110 #define M_SQRTPI 1.772453850905516027298167483341145182798L
111 #define M_CBRTPI 1.464591887561523263020142527263790391739L
112 #define M_SQRT3 1.732050807568877293527446341505872366943L
113 #define M_CBRT2 1.259921049894873164767210607278228350570L
114 #define M_CBRT3 1.442249570307408382321638310780109588392L
115 #define M_CBRT4 1.587401051968199474751705639272308260391L
116 #define M_CBRT5 1.709975946676696989353108872543860109868L
117 #define M_CBRT6 1.817120592832139658891211756327260502428L
118 #define M_CBRT7 1.912931182772389101199116839548760282862L
119 #define M_CBRT9 2.080083823051904114530056824357885386338L
120
121 /* Very useful macros */
122 #ifndef min
123 #define min(x,y) (((x)<(y)) ? (x) : (y))
124 #endif
125 #ifndef max
126 #define max(x,y) (((x)<(y)) ? (y) : (x))
127 #endif
128
129 /* some useful constants */
130 #define LOG_DBL_MIN (log(DBL_MIN))
131 #define LOG_DBL_MAX (log(DBL_MAX))
132 #define SQRT_DBL_EPSILON (sqrt(DBL_EPSILON))
133
134 /* special functions */
135 #define Heaviside(x) (((x) >= 0) ? 1.0 : 0.0)
136 GPU_FUNCTION double LambertW(double z);
137 GPU_FUNCTION double xc_dilogarithm(const double x);
138 #define xc_E1_scaled(x) xc_expint_e1_impl(x, 1)
139
140 /* we define this function here, so it can be properly inlined by all compilers */
141 GPU_FUNCTION
142 static inline double
xc_cheb_eval(const double x,const double * cs,const int N)143 xc_cheb_eval(const double x, const double *cs, const int N)
144 {
145 int i;
146 double twox, b0, b1, b2;
147
148 b2 = b1 = b0 = 0.0;
149
150 twox = 2.0*x;
151 for(i=N-1; i>=0; i--){
152 b2 = b1;
153 b1 = b0;
154 b0 = twox*b1 - b2 + cs[i];
155 }
156
157 return 0.5*(b0 - b2);
158 }
159
160 GPU_FUNCTION double xc_bessel_I0_scaled(const double x);
161 GPU_FUNCTION double xc_bessel_I0(const double x);
162 GPU_FUNCTION double xc_bessel_I1_scaled(const double x);
163 GPU_FUNCTION double xc_bessel_I1(const double x);
164 GPU_FUNCTION double xc_bessel_K0_scaled(const double x);
165 GPU_FUNCTION double xc_bessel_K0(const double x);
166 GPU_FUNCTION double xc_bessel_K1_scaled(const double x);
167 GPU_FUNCTION double xc_bessel_K1(const double x);
168
169 GPU_FUNCTION double xc_expint_e1_impl(double x, const int scale);
expint_e1(const double x)170 GPU_FUNCTION static inline double expint_e1(const double x) { return xc_expint_e1_impl( x, 0); }
expint_e1_scaled(const double x)171 GPU_FUNCTION static inline double expint_e1_scaled(const double x) { return xc_expint_e1_impl( x, 1); }
expint_Ei(const double x)172 GPU_FUNCTION static inline double expint_Ei(const double x) { return -xc_expint_e1_impl(-x, 0); }
173 #define Ei(x) expint_Ei(x)
expint_Ei_scaled(const double x)174 GPU_FUNCTION static inline double expint_Ei_scaled(const double x) { return -xc_expint_e1_impl(-x, 1); }
175
176 GPU_FUNCTION double xc_erfcx(double x);
177
178 /* integration */
179 typedef void integr_fn(double *x, int n, void *ex);
180
181 GPU_FUNCTION double xc_integrate(integr_fn func, void *ex, double a, double b);
182
183 GPU_FUNCTION
184 void xc_rdqagse(integr_fn f, void *ex, double *a, double *b,
185 double *epsabs, double *epsrel, int *limit, double *result,
186 double *abserr, int *neval, int *ier, double *alist__,
187 double *blist, double *rlist, double *elist, int *iord, int *last);
188
189 /* root finding */
190 typedef double xc_brent_f(double, void *);
191
192 GPU_FUNCTION
193 double xc_math_brent(xc_brent_f f, double lower_bound, double upper_bound,
194 double TOL, double MAX_ITER, void *f_params);
195
196
197 typedef struct xc_functional_key_t {
198 char name[256];
199 int number;
200 } xc_functional_key_t;
201
202
203 #define M_C 137.0359996287515 /* speed of light */
204
205 #define RS_FACTOR 0.6203504908994000166680068120477781673508 /* (3/(4*Pi))^1/3 */
206 #define X_FACTOR_C 0.9305257363491000250020102180716672510262 /* 3/8*cur(3/pi)*4^(2/3) */
207 #define X_FACTOR_2D_C 1.504505556127350098528211870828726895584 /* 8/(3*sqrt(pi)) */
208 #define K_FACTOR_C 4.557799872345597137288163759599305358515 /* 3/10*(6*pi^2)^(2/3) */
209 #define MU_GE 0.1234567901234567901234567901234567901235 /* 10/81 */
210 #define MU_PBE 0.2195149727645171 /* mu = beta*pi^2/3, beta = 0.06672455060314922 */
211 #define X2S 0.1282782438530421943003109254455883701296 /* 1/(2*(6*pi^2)^(1/3)) */
212 #define X2S_2D 0.1410473958869390717370198628901931464610 /* 1/(2*(4*pi)^(1/2)) */
213 #define FZETAFACTOR 0.5198420997897463295344212145564567011405 /* 2^(4/3) - 2 */
214
215 #define RS(x) (RS_FACTOR/CBRT(x))
216 #define FZETA(x) ((pow(1.0 + (x), 4.0/3.0) + pow(1.0 - (x), 4.0/3.0) - 2.0)/FZETAFACTOR)
217 #define DFZETA(x) ((CBRT(1.0 + (x)) - CBRT(1.0 - (x)))*(4.0/3.0)/FZETAFACTOR)
218 #define D2FZETA(x) ((4.0/9.0)/FZETAFACTOR)* \
219 (fabs(x)==1.0 ? (FLT_MAX) : (pow(1.0 + (x), -2.0/3.0) + pow(1.0 - (x), -2.0/3.0)))
220 #define D3FZETA(x) (-(8.0/27.0)/FZETAFACTOR)* \
221 (fabs(x)==1.0 ? (FLT_MAX) : (pow(1.0 + (x), -5.0/3.0) - pow(1.0 - (x), -5.0/3.0)))
222
223
224 /* The following inlines confuse the xlc compiler */
225 GPU_FUNCTION void xc_rho2dzeta(int nspin, const double *rho, double *d, double *zeta);
226
227 /* Functions to handle the internal counters */
228
229 void internal_counters_set_lda (int nspin, xc_dimensions *dim);
230 GPU_FUNCTION void internal_counters_lda_random
231 (const xc_dimensions *dim, int ip, int offset, const double **rho, double **zk, LDA_OUT_PARAMS_NO_EXC(double **));
232 GPU_FUNCTION void internal_counters_lda_next
233 (const xc_dimensions *dim, int offset, const double **rho, double **zk, LDA_OUT_PARAMS_NO_EXC(double **));
234 GPU_FUNCTION void internal_counters_lda_prev
235 (const xc_dimensions *dim, int offset, const double **rho, double **zk, LDA_OUT_PARAMS_NO_EXC(double **));
236
237 void internal_counters_set_gga (int nspin, xc_dimensions *dim);
238 GPU_FUNCTION void internal_counters_gga_random
239 (const xc_dimensions *dim, int pos, int offset, const double **rho, const double **sigma,
240 double **zk, GGA_OUT_PARAMS_NO_EXC(double **));
241 GPU_FUNCTION void internal_counters_gga_next
242 (const xc_dimensions *dim, int offset, const double **rho, const double **sigma,
243 double **zk, GGA_OUT_PARAMS_NO_EXC(double **));
244 GPU_FUNCTION void internal_counters_gga_prev
245 (const xc_dimensions *dim, int offset, const double **rho, const double **sigma,
246 double **zk, GGA_OUT_PARAMS_NO_EXC(double **));
247
248 void internal_counters_set_mgga(int nspin, xc_dimensions *dim);
249
250 GPU_FUNCTION void internal_counters_mgga_random
251 (const xc_dimensions *dim, const int pos, int offset,
252 const double **rho, const double **sigma, const double **lapl, const double **tau,
253 double **zk, MGGA_OUT_PARAMS_NO_EXC(double **));
254
255 GPU_FUNCTION void internal_counters_mgga_next
256 (const xc_dimensions *dim, int offset,
257 const double **rho, const double **sigma, const double **lapl, const double **tau,
258 double **zk, MGGA_OUT_PARAMS_NO_EXC(double **));
259
260 GPU_FUNCTION void internal_counters_mgga_prev
261 (const xc_dimensions *dim, int offset,
262 const double **rho, const double **sigma, const double **lapl, const double **tau,
263 double **zk, MGGA_OUT_PARAMS_NO_EXC(double **));
264
265 /* Functionals that are defined as mixtures of others */
266 void xc_mix_init(xc_func_type *p, int n_funcs, const int *funcs_id, const double *mix_coef);
267 void xc_mix_func
268 (const xc_func_type *func, size_t np,
269 const double *rho, const double *sigma, const double *lapl, const double *tau,
270 double *zk, MGGA_OUT_PARAMS_NO_EXC(double *));
271
272 /* Some useful functions */
273 const char *get_kind(const xc_func_type *func);
274 const char *get_family(const xc_func_type *func);
275 double get_ext_param(const xc_func_type *func, const double *values, int index);
276 void set_ext_params_cpy (xc_func_type *p, const double *ext_params);
277 void set_ext_params_cpy_omega(xc_func_type *p, const double *ext_params);
278 void set_ext_params_cpy_exx(xc_func_type *p, const double *ext_params);
279 void set_ext_params_cpy_cam(xc_func_type *p, const double *ext_params);
280
281 GPU_FUNCTION
282 double xc_mgga_x_br89_get_x(double Q);
283
284 #ifndef HAVE_CUDA
285 #define libxc_malloc malloc
286 #define libxc_calloc calloc
287 #define libxc_free free
288 #define libxc_memset memset
289 #else
290
291 template <class int_type>
libxc_malloc(const int_type size)292 auto libxc_malloc(const int_type size){
293 void * mem;
294 cudaMallocManaged(&mem, size);
295 return mem;
296 }
297
298 template <class int_type1, class int_type2>
libxc_calloc(const int_type1 size1,const int_type2 size2)299 auto libxc_calloc(const int_type1 size1, const int_type2 size2){
300 void * mem;
301 cudaMallocManaged(&mem, size1*size2);
302 cudaMemset(mem, 0, size1*size2);
303 return mem;
304 }
305
306 #define libxc_free cudaFree
307 #define libxc_memset cudaMemset
308
309 #endif
310
311 #endif
312