1 /*
2 Copyright (C) 2006-2018 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 /**
11 * @file work_gga.c
12 * @brief This file is to be included in GGA functionals.
13 */
14
15 #ifdef XC_DEBUG
16 #define __USE_GNU
17 #include <fenv.h>
18 #endif
19
20 /* hack to avoid compiler warnings */
21 #define NOARG
22
23 #ifdef XC_NO_EXC
24 #define OUT_PARAMS GGA_OUT_PARAMS_NO_EXC(NOARG)
25 #else
26 #define OUT_PARAMS zk, GGA_OUT_PARAMS_NO_EXC(NOARG)
27 #endif
28
29 #ifdef HAVE_CUDA
30 __global__ static void
31 work_gga_gpu(const XC(func_type) *p, int order, size_t np, const double *rho, const double *sigma, double *zk, GGA_OUT_PARAMS_NO_EXC(double *));
32 #endif
33
34 /**
35 * @param[in,out] func_type: pointer to functional structure
36 */
37 static void
work_gga(const XC (func_type)* p,size_t np,const double * rho,const double * sigma,double * zk,GGA_OUT_PARAMS_NO_EXC (double *))38 work_gga(const XC(func_type) *p, size_t np,
39 const double *rho, const double *sigma,
40 double *zk, GGA_OUT_PARAMS_NO_EXC(double *))
41 {
42
43 int order = -1;
44 if(zk != NULL) order = 0;
45 if(vrho != NULL) order = 1;
46 if(v2rho2 != NULL) order = 2;
47 if(v3rho3 != NULL) order = 3;
48
49 if(order < 0) return;
50
51 #ifdef XC_DEBUG
52 feenableexcept(FE_DIVBYZERO | FE_INVALID);
53 #endif
54
55 #ifdef HAVE_CUDA
56
57 //make a copy of 'p' since it might be in host-only memory
58 XC(func_type) * pcuda = (XC(func_type) *) libxc_malloc(sizeof(XC(func_type)));
59
60 *pcuda = *p;
61
62 auto nblocks = np/CUDA_BLOCK_SIZE;
63 if(np != nblocks*CUDA_BLOCK_SIZE) nblocks++;
64
65 work_gga_gpu<<<nblocks, CUDA_BLOCK_SIZE>>>(pcuda, order, np, rho, sigma, zk, GGA_OUT_PARAMS_NO_EXC(NOARG));
66
67 libxc_free(pcuda);
68
69 #else
70
71 size_t ip;
72 double my_rho[2] = {0.0, 0.0}, my_sigma[3] = {0.0, 0.0, 0.0};
73 double dens, zeta;
74
75 for(ip = 0; ip < np; ip++){
76 /* sanity check on input parameters */
77 my_rho[0] = max(0.0, rho[0]);
78 my_sigma[0] = max(1e-40, sigma[0]);
79 if(p->nspin == XC_POLARIZED){
80 my_rho[1] = max(0.0, rho[1]);
81 my_sigma[2] = max(1e-40, sigma[2]);
82 /* | grad n |^2 = |grad n_up + grad n_down|^2 > 0 */
83 my_sigma[1] = max(sigma[1], -(sigma[0] + sigma[1])/2.0);
84 }
85 xc_rho2dzeta(p->nspin, my_rho, &dens, &zeta);
86
87 if(dens > p->dens_threshold){
88 if(p->nspin == XC_UNPOLARIZED){ /* unpolarized case */
89 func_unpol(p, order, my_rho, my_sigma, OUT_PARAMS);
90
91 }else if(zeta > 1.0 - 1e-10){ /* ferromagnetic case - spin 0 */
92 func_ferr(p, order, my_rho, my_sigma, OUT_PARAMS);
93
94 }else if(zeta < -1.0 + 1e-10){ /* ferromagnetic case - spin 1 */
95 internal_counters_gga_next(&(p->dim), -1, &rho, &sigma, &zk, GGA_OUT_PARAMS_NO_EXC(&));
96 func_ferr(p, order, &my_rho[1], &my_sigma[2], OUT_PARAMS);
97 internal_counters_gga_prev(&(p->dim), -1, &rho, &sigma, &zk, GGA_OUT_PARAMS_NO_EXC(&));
98
99 }else{ /* polarized (general) case */
100 func_pol(p, order, my_rho, my_sigma, OUT_PARAMS);
101 } /* polarization */
102 }
103
104 /* check for NaNs */
105 #ifdef XC_DEBUG
106 {
107 size_t ii;
108 const xc_dimensions *dim = &(p->dim);
109 int is_OK = 1;
110
111 if(zk != NULL)
112 is_OK = is_OK & isfinite(*zk);
113
114 if(vrho != NULL){
115 for(ii=0; ii < dim->vrho; ii++)
116 is_OK = is_OK && isfinite(vrho[ii]);
117 for(ii=0; ii < dim->vsigma; ii++)
118 is_OK = is_OK && isfinite(vsigma[ii]);
119 }
120
121 if(!is_OK){
122 printf("Problem in the evaluation of the functional\n");
123 if(p->nspin == XC_UNPOLARIZED){
124 printf("./xc-get_data %d 1 %le 0.0 %le 0.0 0.0 0.0 0.0 0.0 0.0\n",
125 p->info->number, *rho, *sigma);
126 }else{
127 printf("./xc-get_data %d 2 %le %le %le %le %le 0.0 0.0 0.0 0.0\n",
128 p->info->number, rho[0], rho[1], sigma[0], sigma[1], sigma[2]);
129 }
130 }
131 }
132 #endif
133
134 internal_counters_gga_next(&(p->dim), 0, &rho, &sigma, &zk, GGA_OUT_PARAMS_NO_EXC(&));
135 } /* for(ip) */
136
137 #endif
138 }
139
140 #ifdef HAVE_CUDA
141
142 __global__ static void
work_gga_gpu(const XC (func_type)* p,int order,size_t np,const double * rho,const double * sigma,double * zk,GGA_OUT_PARAMS_NO_EXC (double *))143 work_gga_gpu(const XC(func_type) *p, int order, size_t np, const double *rho, const double *sigma,
144 double *zk, GGA_OUT_PARAMS_NO_EXC(double *))
145 {
146 double my_rho[2] = {0.0, 0.0}, my_sigma[3] = {0.0, 0.0, 0.0};
147 double dens, zeta;
148
149 size_t ip = blockIdx.x*blockDim.x + threadIdx.x;
150
151 if(ip >= np) return;
152
153 internal_counters_gga_random(&(p->dim), ip, 0, &rho, &sigma, &zk, GGA_OUT_PARAMS_NO_EXC(&));
154
155 /* sanity check on input parameters */
156 my_rho[0] = max(0.0, rho[0]);
157 my_sigma[0] = max(1e-40, sigma[0]);
158 if(p->nspin == XC_POLARIZED){
159 my_rho[1] = max(0.0, rho[1]);
160 my_sigma[2] = max(1e-40, sigma[2]);
161 /* | grad n |^2 = |grad n_up + grad n_down|^2 > 0 */
162 my_sigma[1] = max(sigma[1], -(sigma[0] + sigma[1])/2.0);
163 }
164 xc_rho2dzeta(p->nspin, my_rho, &dens, &zeta);
165
166 if(dens > p->dens_threshold){
167 if(p->nspin == XC_UNPOLARIZED){ /* unpolarized case */
168 func_unpol(p, order, my_rho, my_sigma, OUT_PARAMS);
169
170 } else if(zeta > 1.0 - 1e-10){ /* ferromagnetic case - spin 0 */
171 func_ferr(p, order, my_rho, my_sigma, OUT_PARAMS);
172
173 } else if(zeta < -1.0 + 1e-10){ /* ferromagnetic case - spin 1 */
174 internal_counters_gga_next(&(p->dim), -1, &rho, &sigma, &zk, GGA_OUT_PARAMS_NO_EXC(&));
175 func_ferr(p, order, &my_rho[1], &my_sigma[2], OUT_PARAMS);
176
177 } else { /* polarized (general) case */
178 func_pol(p, order, my_rho, my_sigma, OUT_PARAMS);
179
180 } /* polarization */
181 }
182 }
183
184 #endif
185