1 /*
2 Copyright (C) 2006-2007 M.A.L. Marques
3
4 This Source Code Form is subject to the terms of the Mozilla Public
5 License, v. 2.0. If a copy of the MPL was not distributed with this
6 file, You can obtain one at http://mozilla.org/MPL/2.0/.
7 */
8
9
10 #include "util.h"
11 #include "funcs_gga.c"
12 #include "funcs_hyb_gga.c"
13
14 /* Some useful formulas:
15
16 sigma_st = grad rho_s . grad rho_t
17 zk = energy density per unit particle
18
19 vrho_s = d zk / d rho_s
20 vsigma_st = d n*zk / d sigma_st
21
22 v2rho2_st = d^2 n*zk / d rho_s d rho_t
23 v2rhosigma_svx = d^2 n*zk / d rho_s d sigma_tv
24 v2sigma2_stvx = d^2 n*zk / d sigma_st d sigma_vx
25
26 v3rho3_stv = d^3 n*zk / d rho_s d rho_t d rho_v
27 v3rho2sigma_stvx = d^3 n*zk / d rho_s d rho_t d sigma_vx
28 v3rhosigma2_svxyz = d^3 n*zk / d rho_s d sigma_vx d sigma_yz
29 v3sigma3_stvxyz = d^3 n*zk / d sigma_st d sigma_vx d sigma_yz
30
31 if nspin == 2
32 rho(2) = (u, d)
33 sigma(3) = (uu, ud, dd)
34
35 vrho(2) = (u, d)
36 vsigma(3) = (uu, ud, dd)
37
38 v2rho2(3) = (u_u, u_d, d_d)
39 v2rhosigma(6) = (u_uu, u_ud, u_dd, d_uu, d_ud, d_dd)
40 v2sigma2(6) = (uu_uu, uu_ud, uu_dd, ud_ud, ud_dd, dd_dd)
41
42 v3rho3(4) = (u_u_u, u_u_d, u_d_d, d_d_d)
43 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)
44 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)
45 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)
46
47 */
xc_gga(const xc_func_type * func,size_t np,const double * rho,const double * sigma,double * zk,GGA_OUT_PARAMS_NO_EXC (double *))48 void xc_gga(const xc_func_type *func, size_t np, const double *rho, const double *sigma,
49 double *zk, GGA_OUT_PARAMS_NO_EXC(double *))
50 {
51 const xc_dimensions *dim = &(func->dim);
52
53 /* sanity check */
54 if(zk != NULL && !(func->info->flags & XC_FLAGS_HAVE_EXC)){
55 fprintf(stderr, "Functional '%s' does not provide an implementation of Exc\n",
56 func->info->name);
57 exit(1);
58 }
59
60 if(vrho != NULL && !(func->info->flags & XC_FLAGS_HAVE_VXC)){
61 fprintf(stderr, "Functional '%s' does not provide an implementation of vxc\n",
62 func->info->name);
63 exit(1);
64 }
65
66 if(v2rho2 != NULL && !(func->info->flags & XC_FLAGS_HAVE_FXC)){
67 fprintf(stderr, "Functional '%s' does not provide an implementation of fxc\n",
68 func->info->name);
69 exit(1);
70 }
71
72 if(v3rho3 != NULL && !(func->info->flags & XC_FLAGS_HAVE_KXC)){
73 fprintf(stderr, "Functional '%s' does not provide an implementation of kxc\n",
74 func->info->name);
75 exit(1);
76 }
77
78 if(v4rho4 != NULL && !(func->info->flags & XC_FLAGS_HAVE_LXC)){
79 fprintf(stderr, "Functional '%s' does not provide an implementation of lxc\n",
80 func->info->name);
81 exit(1);
82 }
83
84 /* initialize output to zero */
85 if(zk != NULL)
86 libxc_memset(zk, 0, dim->zk*np*sizeof(double));
87
88 if(vrho != NULL){
89 assert(vsigma != NULL);
90
91 libxc_memset(vrho, 0, dim->vrho *np*sizeof(double));
92 libxc_memset(vsigma, 0, dim->vsigma*np*sizeof(double));
93 }
94
95 if(v2rho2 != NULL){
96 assert(v2rhosigma!=NULL && v2sigma2!=NULL);
97
98 libxc_memset(v2rho2, 0, dim->v2rho2 *np*sizeof(double));
99 libxc_memset(v2rhosigma, 0, dim->v2rhosigma*np*sizeof(double));
100 libxc_memset(v2sigma2, 0, dim->v2sigma2 *np*sizeof(double));
101 }
102
103 if(v3rho3 != NULL){
104 assert(v3rho2sigma!=NULL && v3rhosigma2!=NULL && v3sigma3!=NULL);
105
106 libxc_memset(v3rho3, 0, dim->v3rho3 *np*sizeof(double));
107 libxc_memset(v3rho2sigma, 0, dim->v3rho2sigma*np*sizeof(double));
108 libxc_memset(v3rhosigma2, 0, dim->v3rhosigma2*np*sizeof(double));
109 libxc_memset(v3sigma3, 0, dim->v3sigma3 *np*sizeof(double));
110 }
111
112 if(v4rho4 != NULL){
113 assert(v4rho3sigma!=NULL && v4rho2sigma2!=NULL && v4rhosigma3!=NULL && v4sigma4!=NULL);
114
115 libxc_memset(v4rho4, 0, dim->v4rho4 *np*sizeof(double));
116 libxc_memset(v4rho3sigma, 0, dim->v4rho3sigma *np*sizeof(double));
117 libxc_memset(v4rho2sigma2, 0, dim->v4rho2sigma2*np*sizeof(double));
118 libxc_memset(v4rhosigma3, 0, dim->v4rhosigma3 *np*sizeof(double));
119 libxc_memset(v4sigma4, 0, dim->v4sigma4 *np*sizeof(double));
120 }
121
122 /* call functional */
123 if(func->info->gga != NULL)
124 func->info->gga(func, np, rho, sigma, zk, GGA_OUT_PARAMS_NO_EXC(XC_NOARG));
125
126 if(func->mix_coef != NULL)
127 xc_mix_func(func, np, rho, sigma, NULL, NULL, zk, vrho, vsigma, NULL, NULL,
128 v2rho2, v2rhosigma, NULL, NULL, v2sigma2, NULL, NULL, NULL, NULL, NULL,
129 v3rho3, v3rho2sigma, NULL, NULL, v3rhosigma2, NULL, NULL, NULL, NULL, NULL,
130 v3sigma3, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
131 v4rho4, v4rho3sigma, NULL, NULL, v4rho2sigma2, NULL, NULL, NULL, NULL, NULL,
132 v4rhosigma3, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
133 v4sigma4, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
134 NULL, NULL, NULL, NULL, NULL);
135 }
136
137 /* specializations */
138 /* returns only energy */
139 void
xc_gga_exc(const xc_func_type * p,size_t np,const double * rho,const double * sigma,double * zk)140 xc_gga_exc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
141 double *zk)
142 {
143 xc_gga(p, np, rho, sigma, zk, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
144 NULL, NULL, NULL, NULL, NULL);
145 }
146
147 /* returns only potential */
148 void
xc_gga_vxc(const xc_func_type * p,size_t np,const double * rho,const double * sigma,double * vrho,double * vsigma)149 xc_gga_vxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
150 double *vrho, double *vsigma)
151 {
152 xc_gga(p, np, rho, sigma, NULL, vrho, vsigma, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
153 NULL, NULL, NULL, NULL, NULL);
154 }
155
156 /* returns both energy and potential (the most common call usually) */
157 void
xc_gga_exc_vxc(const xc_func_type * p,size_t np,const double * rho,const double * sigma,double * zk,double * vrho,double * vsigma)158 xc_gga_exc_vxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
159 double *zk, double *vrho, double *vsigma)
160 {
161 xc_gga(p, np, rho, sigma, zk, vrho, vsigma, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
162 NULL, NULL, NULL, NULL, NULL);
163 }
164
165 /* returns energy, first and second derivatives */
xc_gga_exc_vxc_fxc(const xc_func_type * p,size_t np,const double * rho,const double * sigma,double * zk,double * vrho,double * vsigma,double * v2rho2,double * v2rhosigma,double * v2sigma2)166 void xc_gga_exc_vxc_fxc (const xc_func_type *p, size_t np, const double *rho, const double *sigma,
167 double *zk, double *vrho, double *vsigma,
168 double *v2rho2, double *v2rhosigma, double *v2sigma2) {
169 xc_gga(p, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, NULL, NULL, NULL, NULL,
170 NULL, NULL, NULL, NULL, NULL);
171 }
172
xc_gga_vxc_fxc(const xc_func_type * p,size_t np,const double * rho,const double * sigma,double * vrho,double * vsigma,double * v2rho2,double * v2rhosigma,double * v2sigma2)173 void xc_gga_vxc_fxc (const xc_func_type *p, size_t np, const double *rho, const double *sigma,
174 double *vrho, double *vsigma,
175 double *v2rho2, double *v2rhosigma, double *v2sigma2) {
176 xc_gga(p, np, rho, sigma, NULL, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, NULL, NULL, NULL, NULL,
177 NULL, NULL, NULL, NULL, NULL);
178 }
179
180 /* returns energy, first, second and third derivatives */
xc_gga_exc_vxc_fxc_kxc(const xc_func_type * p,size_t 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)181 void xc_gga_exc_vxc_fxc_kxc (const xc_func_type *p, size_t np, const double *rho, const double *sigma,
182 double *zk, double *vrho, double *vsigma, double *v2rho2, double *v2rhosigma, double *v2sigma2,
183 double *v3rho3, double *v3rho2sigma, double *v3rhosigma2, double *v3sigma3) {
184 xc_gga(p, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3,
185 NULL, NULL, NULL, NULL, NULL);
186 }
187
xc_gga_vxc_fxc_kxc(const xc_func_type * p,size_t np,const double * rho,const double * sigma,double * vrho,double * vsigma,double * v2rho2,double * v2rhosigma,double * v2sigma2,double * v3rho3,double * v3rho2sigma,double * v3rhosigma2,double * v3sigma3)188 void xc_gga_vxc_fxc_kxc (const xc_func_type *p, size_t np, const double *rho, const double *sigma,
189 double *vrho, double *vsigma, double *v2rho2, double *v2rhosigma, double *v2sigma2,
190 double *v3rho3, double *v3rho2sigma, double *v3rhosigma2, double *v3sigma3) {
191 xc_gga(p, np, rho, sigma, NULL, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3,
192 NULL, NULL, NULL, NULL, NULL);
193 }
194
195
196 /* returns second derivatives */
197 void
xc_gga_fxc(const xc_func_type * p,size_t np,const double * rho,const double * sigma,double * v2rho2,double * v2rhosigma,double * v2sigma2)198 xc_gga_fxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
199 double *v2rho2, double *v2rhosigma, double *v2sigma2)
200 {
201 xc_gga(p, np, rho, sigma, NULL, NULL, NULL, v2rho2, v2rhosigma, v2sigma2, NULL, NULL, NULL, NULL,
202 NULL, NULL, NULL, NULL, NULL);
203 }
204
205 /* returns third derivatives */
206 void
xc_gga_kxc(const xc_func_type * p,size_t np,const double * rho,const double * sigma,double * v3rho3,double * v3rho2sigma,double * v3rhosigma2,double * v3sigma3)207 xc_gga_kxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
208 double *v3rho3, double *v3rho2sigma, double *v3rhosigma2, double *v3sigma3)
209 {
210 xc_gga(p, np, rho, sigma, NULL, NULL, NULL, NULL, NULL, NULL, v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3,
211 NULL, NULL, NULL, NULL, NULL);
212 }
213
214 /* returns fourth derivatives */
215 void
xc_gga_lxc(const xc_func_type * p,size_t np,const double * rho,const double * sigma,double * v4rho4,double * v4rho3sigma,double * v4rho2sigma2,double * v4rhosigma3,double * v4sigma4)216 xc_gga_lxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
217 double *v4rho4, double *v4rho3sigma, double *v4rho2sigma2, double *v4rhosigma3, double *v4sigma4)
218 {
219 xc_gga(p, np, rho, sigma, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
220 v4rho4, v4rho3sigma, v4rho2sigma2, v4rhosigma3, v4sigma4);
221 }
222