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 #include <math.h>
10 #include <stdio.h>
11 #include <stdlib.h>
12
13 #include <xc.h>
14
15
16 typedef struct {
17 int functional;
18 int nspin;
19
20 /* Input */
21 double rho[2]; /* rhoa, rhob */
22 double sigma[3]; /* sigmaaa, sigmaab, sigmabb */
23 double lapl[2]; /* lapla, laplb */
24 double tau[2]; /* taua, taub */
25
26 /* Energy */
27 double zk; /* energy density per unit particle */
28
29 /* First derivatives */
30 double vrho[2]; /* vrhoa, vrhob */
31 double vsigma[3]; /* vsigmaaa, vsigmaab, vsigmabb */
32 double vlapl[2]; /* vlapla, vlaplb */
33 double vtau[2]; /* vtaua, vtaub */
34
35 /* Second derivatives */
36 double v2rho2[3]; /* v2rhoa2, v2rhoab, v2rhob2 */
37 double v2rhosigma[6]; /* v2rhoasigmaaa, v2rhoasigmaab, v2rhoasigmabb
38 v2rhobsigmaaa, v2rhobsigmaab, v2rhobsigmabb */
39 double v2rholapl[3]; /* */
40 double v2rhotau[3]; /* */
41 double v2sigma2[6]; /* v2sigmaaa2, v2sigmaaaab, v2sigmaaabb
42 v2sigmaab2, v2sigmaabbb, v2sigmabb2 */
43 double v2sigmalapl[6]; /* v2sigmaaalapla, v2sigmaaalaplb,
44 v2sigmaablapla, v2sigmaablaplb,
45 v2sigmabblapla, v2sigmabblaplb */
46 double v2sigmatau[6]; /* v2sigmaaataua, v2sigmaaataub,
47 v2sigmaabtaua, v2sigmaabtaub,
48 v2sigmabbtaua, v2sigmabbtaub */
49 double v2lapl2[3]; /* v2lapla2, v2laplab, v2laplb2 */
50 double v2lapltau[3]; /* */
51 double v2tau2[3]; /* v2taua2, v2tauab, v2taub2 */
52
53 /* Third derivatives */
54 double v3rho3[4]; /* v3rhoaaa, v3rhoaab, v3rhoabb, v3rhobbb */
55
56 } xc_values_type;
57
58 /*----------------------------------------------------------*/
init_values(xc_values_type * xc_values,char * argv[])59 void init_values(xc_values_type *xc_values, char *argv[])
60 {
61 int i;
62
63 xc_values->functional = atoi(argv[1]);
64 xc_values->nspin = atoi(argv[2]);
65 xc_values->rho[0] = atof(argv[3]);
66 xc_values->rho[1] = atof(argv[4]);
67 xc_values->sigma[0] = atof(argv[5]);
68 xc_values->sigma[1] = atof(argv[6]);
69 xc_values->sigma[2] = atof(argv[7]);
70 xc_values->lapl[0] = atof(argv[8]);
71 xc_values->lapl[1] = atof(argv[9]);
72 xc_values->tau[0] = atof(argv[10]);
73 xc_values->tau[1] = atof(argv[11]);
74
75 xc_values->zk = 0;
76
77 for(i=0; i<2; i++){
78 xc_values->vrho[i] = 0;
79 xc_values->vlapl[i] = 0;
80 xc_values->vtau[i] = 0;
81 }
82
83 for(i=0; i<3; i++){
84 xc_values->vsigma[i] = 0;
85 xc_values->v2rho2[i] = 0;
86 xc_values->v2lapl2[i] = 0;
87 xc_values->v2tau2[i] = 0;
88 xc_values->v2rholapl[i] = 0;
89 xc_values->v2rhotau[i] = 0;
90 xc_values->v2lapltau[i] = 0;
91 }
92
93 for(i=0; i<4; i++){
94 xc_values->v3rho3[i] = 0;
95 }
96
97 for(i=0; i<6; i++){
98 xc_values->v2rhosigma[i] = 0;
99 xc_values->v2sigma2[i] = 0;
100 xc_values->v2sigmalapl[i] = 0;
101 xc_values->v2sigmatau[i] = 0;
102 }
103 }
104
105
106 /*----------------------------------------------------------*/
print_values(xc_values_type * xc)107 void print_values(xc_values_type *xc)
108 {
109 /*int family = xc_family_from_id(xc->functional, NULL, NULL);*/
110
111 printf(" rhoa= %#0.2E rhob= %#0.2E sigmaaa= %#0.2E sigmaab= %#0.2E sigmabb= %#0.2E lapla= %#0.2E laplb= %#0.2E taua= %#0.2E taub= %#0.2E\n\n",
112 xc->rho[0], xc->rho[1],
113 xc->sigma[0], xc->sigma[1], xc->sigma[2],
114 xc->lapl[0], xc->lapl[1],
115 xc->tau[0], xc->tau[1]);
116 printf(" zk = %#19.12E\n\n",
117 xc->zk);
118 printf(" vrhoa = %#19.12E\n"
119 " vrhob = %#19.12E\n"
120 " vsigmaaa = %#19.12E\n"
121 " vsigmaab = %#19.12E\n"
122 " vsigmabb = %#19.12E\n"
123 " vlapla = %#19.12E\n"
124 " vlaplb = %#19.12E\n"
125 " vtaua = %#19.12E\n"
126 " vtaub = %#19.12E\n\n",
127 xc->vrho[0], xc->vrho[1],
128 xc->vsigma[0], xc->vsigma[1], xc->vsigma[2],
129 xc->vlapl[0], xc->vlapl[1],
130 xc->vtau[0], xc->vtau[1]);
131
132 printf(" v2rhoa2 = %#19.12E\n"
133 " v2rhoab = %#19.12E\n"
134 " v2rhob2 = %#19.12E\n"
135 " v2rhoasigmaaa = %#19.12E\n"
136 " v2rhoasigmaab = %#19.12E\n"
137 " v2rhoasigmabb = %#19.12E\n"
138 " v2rhobsigmaaa = %#19.12E\n"
139 " v2rhobsigmaab = %#19.12E\n"
140 " v2rhobsigmabb = %#19.12E\n"
141 " v2sigmaaa2 = %#19.12E\n"
142 " v2sigmaaaab = %#19.12E\n"
143 " v2sigmaaabb = %#19.12E\n"
144 " v2sigmaab2 = %#19.12E\n"
145 " v2sigmaabbb = %#19.12E\n"
146 " v2sigmabb2 = %#19.12E\n\n",
147 xc->v2rho2[0], xc->v2rho2[1], xc->v2rho2[2],
148 xc->v2rhosigma[0], xc->v2rhosigma[1], xc->v2rhosigma[2],
149 xc->v2rhosigma[3], xc->v2rhosigma[4], xc->v2rhosigma[5],
150 xc->v2sigma2[0], xc->v2sigma2[1], xc->v2sigma2[2],
151 xc->v2sigma2[3], xc->v2sigma2[4], xc->v2sigma2[5]
152 );
153 printf(" v3rhoa3 = %#19.12E\n"
154 " v2rhoaab = %#19.12E\n"
155 " v2rhoabb = %#19.12E\n"
156 " v2rhob3 = %#19.12E\n\n",
157 xc->v3rho3[0], xc->v3rho3[1], xc->v3rho3[2], xc->v3rho3[3]
158 );
159
160 }
161
162
163 /*----------------------------------------------------------*/
main(int argc,char * argv[])164 int main(int argc, char *argv[])
165 {
166 xc_values_type xc;
167 xc_func_type func;
168 const xc_func_info_type *info;
169
170 double *pzk = NULL;
171 double *pvrho = NULL;
172 double *pvsigma = NULL;
173 double *pvlapl = NULL;
174 double *pvtau = NULL;
175 double *pv2rho2 = NULL;
176 double *pv2rhosigma = NULL;
177 double *pv2rholapl = NULL;
178 double *pv2rhotau = NULL;
179 double *pv2sigma2 = NULL;
180 double *pv2sigmalapl = NULL;
181 double *pv2sigmatau = NULL;
182 double *pv2lapl2 = NULL;
183 double *pv2lapltau = NULL;
184 double *pv2tau2 = NULL;
185 double *pv3rho3 = NULL;
186
187 if(argc != 12){
188 printf("Usage:\n%s funct pol rhoa rhob sigmaaa sigmaab sigmabb lapla laplb taua taub\n", argv[0]);
189 return 1;
190 }
191
192 init_values(&xc, argv);
193
194 if(xc.nspin == 1){
195 xc.rho[0] += xc.rho[1];
196 xc.sigma[0] += 2.0*xc.sigma[1] + xc.sigma[2];
197 xc.lapl[0] += xc.lapl[1];
198 xc.tau[0] += xc.tau[1];
199 }
200
201 if(xc_func_init(&func, xc.functional, xc.nspin) != 0){
202 fprintf(stderr, "Functional '%d' not found\n", xc.functional);
203 exit(1);
204 }
205 info = func.info;
206
207 if(info->flags & XC_FLAGS_HAVE_EXC){
208 pzk = &xc.zk;
209 }
210 if(info->flags & XC_FLAGS_HAVE_VXC){
211 pvrho = xc.vrho;
212 pvsigma = xc.vsigma;
213 pvlapl = xc.vlapl;
214 pvtau = xc.vtau;
215 }
216 if(info->flags & XC_FLAGS_HAVE_FXC){
217 pv2rho2 = xc.v2rho2;
218 pv2rhosigma = xc.v2rhosigma;
219 pv2rholapl = xc.v2rholapl;
220 pv2rhotau = xc.v2rhotau;
221 pv2sigma2 = xc.v2sigma2;
222 pv2sigmalapl = xc.v2sigmalapl;
223 pv2sigmatau = xc.v2sigmatau;
224 pv2lapl2 = xc.v2lapl2;
225 pv2lapltau = xc.v2lapltau;
226 pv2tau2 = xc.v2tau2;
227 }
228 if(info->flags & XC_FLAGS_HAVE_KXC){
229 pv3rho3 = xc.v3rho3;
230 }
231
232 switch(func.info->family)
233 {
234 case XC_FAMILY_LDA:
235 case XC_FAMILY_HYB_LDA:
236 xc_lda(&func, 1, xc.rho, pzk, pvrho, pv2rho2, pv3rho3, NULL);
237 break;
238
239 case XC_FAMILY_GGA:
240 case XC_FAMILY_HYB_GGA:
241 xc_gga(&func, 1, xc.rho, xc.sigma, pzk, pvrho, pvsigma,
242 pv2rho2, pv2rhosigma, pv2sigma2, NULL, NULL, NULL, NULL,
243 NULL, NULL, NULL, NULL, NULL);
244 break;
245
246 case XC_FAMILY_MGGA:
247 case XC_FAMILY_HYB_MGGA:
248 xc_mgga(&func, 1, xc.rho, xc.sigma, xc.lapl, xc.tau,
249 pzk, pvrho, pvsigma, pvlapl, pvtau,
250 pv2rho2, pv2rhosigma, pv2rholapl, pv2rhotau, pv2sigma2, pv2sigmalapl, pv2sigmatau, pv2lapl2, pv2lapltau, pv2tau2,
251 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
252 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
253 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
254 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
255 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
256 NULL, NULL, NULL, NULL, NULL
257 );
258 break;
259
260 default:
261 fprintf(stderr,"Family %i not supported by xc-get_data.\nEnding program.\n",func.info->family);
262 exit(1);
263 }
264
265
266 xc_func_end(&func);
267
268 if(xc.nspin == 1){
269 xc.zk *= xc.rho[0];
270 }else{
271 xc.zk *= (xc.rho[0] + xc.rho[1]);
272 }
273
274
275 print_values(&xc);
276
277 return 0;
278 }
279
280