1 /*
2
3
4 !
5 ! Dalton, a molecular electronic structure program
6 ! Copyright (C) by the authors of Dalton.
7 !
8 ! This program is free software; you can redistribute it and/or
9 ! modify it under the terms of the GNU Lesser General Public
10 ! License version 2.1 as published by the Free Software Foundation.
11 !
12 ! This program is distributed in the hope that it will be useful,
13 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
14 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 ! Lesser General Public License for more details.
16 !
17 ! If a copy of the GNU LGPL v2.1 was not distributed with this
18 ! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
19 !
20
21 !
22
23 */
24 /*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
25 /* fun-b86mx.c:
26
27 Automatically generated code implementing B86MX functional and
28 its derivatives. It is generated by func-codegen.pl being a part of
29 a "Automatic code generation framework for analytical functional
30 derivative evaluation", Pawel Salek, 2005
31
32 This functional is connected by making following changes:
33 1. add "extern Functional b86mxFunctional;" to 'functionals.h'
34 2. add "&b86mxFunctional," to 'functionals.c'
35 3. add "fun-b86mx.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.
36
37 This functional has been generated from following input:
38 ------ cut here -------
39 Cx: -3/4*(6/%PI)^(1/3);
40 bet: 0.00375;
41 lam: 0.007;
42
43 xa: grada/(rhoa^(4/3));
44 xb: gradb/(rhob^(4/3));
45
46 da: (1+lam*xa^2)^(4/5);
47 db: (1+lam*xb^2)^(4/5);
48
49 Exa: (-rhoa^(4/3))*(Cx + bet*xa^2/da);
50 Exb: (-rhob^(4/3))*(Cx + bet*xb^2/db);
51
52 K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb);
53
54
55 ------ cut here -------
56 */
57
58
59 /* strictly conform to XOPEN ANSI C standard */
60 #if !defined(SYS_DEC)
61 /* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
62 * is not specified. */
63 #define _XOPEN_SOURCE 500
64 #define _XOPEN_SOURCE_EXTENDED 1
65 #endif
66 #include <math.h>
67 #include <stddef.h>
68 #include "general.h"
69
70 #define __CVERSION__
71
72 #include "functionals.h"
73
74 /* INTERFACE PART */
b86mx_isgga(void)75 static integer b86mx_isgga(void) { return 1; } /* FIXME: detect! */
76 static integer b86mx_read(const char *conf_line);
77 static real b86mx_energy(const FunDensProp* dp);
78 static void b86mx_first(FunFirstFuncDrv *ds, real factor,
79 const FunDensProp* dp);
80 static void b86mx_second(FunSecondFuncDrv *ds, real factor,
81 const FunDensProp* dp);
82 static void b86mx_third(FunThirdFuncDrv *ds, real factor,
83 const FunDensProp* dp);
84 static void b86mx_fourth(FunFourthFuncDrv *ds, real factor,
85 const FunDensProp* dp);
86
87 Functional B86mxFunctional = {
88 "B86mx", /* name */
89 b86mx_isgga, /* gga-corrected */
90 1,
91 b86mx_read,
92 NULL,
93 b86mx_energy,
94 b86mx_first,
95 b86mx_second,
96 b86mx_third,
97 b86mx_fourth
98 };
99
100 /* IMPLEMENTATION PART */
101 static integer
b86mx_read(const char * conf_line)102 b86mx_read(const char *conf_line)
103 {
104 fun_set_hf_weight(0);
105 return 1;
106 }
107
108 static real
b86mx_energy(const FunDensProp * dp)109 b86mx_energy(const FunDensProp *dp)
110 {
111 real res;
112 real rhoa = dp->rhoa, rhob = dp->rhob;
113 real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
114
115 real t1, t2, t3, t4, t5;
116
117 t1 = -0.75*pow(6.0,0.333333333333333)/pow(3.141592653589793,
118 0.333333333333333);
119 t2 = pow(grada,2.0);
120 t3 = 1/pow(rhoa,2.666666666666667);
121 t4 = pow(gradb,2.0);
122 t5 = 1/pow(rhob,2.666666666666667);
123
124 /* code */
125 res = -1.0*(0.00375*t4*t5/pow(0.007*t4*t5+1.0,0.8)+t1)*
126 pow(rhob,1.333333333333333)-1.0*(0.00375*t2*t3/pow(0.007*t2*
127 t3+1.0,0.8)+t1)*pow(rhoa,1.333333333333333);
128
129 return res;
130 }
131
132 static void
b86mx_first_helper(real rhoa,real grada,real * res)133 b86mx_first_helper(real rhoa, real grada, real *res)
134 { real t1, t2, t3, t4, t5, t6;
135
136 t1 = pow(grada,2.0);
137 t2 = 1/pow(rhoa,2.666666666666667);
138 t3 = 0.007*t1*t2+1.0;
139 t4 = 1/pow(t3,0.8);
140 t5 = 1/pow(t3,1.8);
141 t6 = pow(rhoa,1.333333333333333);
142
143 /* code */
144 res[0] = -1.0*t6*(5.6E-5*t5*pow(grada,4.0)/pow(rhoa,6.333333333333333)-
145 0.01*t1*t4/pow(rhoa,3.666666666666667))-1.333333333333333*
146 (0.00375*t1*t4*t2-0.75*pow(6.0,0.333333333333333)/pow(3.141592653589793,
147 0.333333333333333))*pow(rhoa,0.333333333333333);
148 res[1] = -1.0*t6*(0.0075*grada*t4*t2-4.2E-5*t5*pow(grada,
149 3.0)/pow(rhoa,5.333333333333333));
150 }
151
152 static void
b86mx_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)153 b86mx_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
154 {
155 real res[2];
156
157 b86mx_first_helper(dp->rhoa, dp->grada, res);
158 /* Final assignment */
159 ds->df1000 += factor*res[0];
160 ds->df0010 += factor*res[1];
161
162
163 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
164 fabs(dp->grada-dp->gradb)>1e-13)
165 b86mx_first_helper(dp->rhob, dp->gradb, res);
166 ds->df0100 += factor*res[0];
167 ds->df0001 += factor*res[1];
168
169 }
170
171 static void
b86mx_second_helper(real rhoa,real grada,real * res)172 b86mx_second_helper(real rhoa, real grada, real *res)
173 {
174 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
175 real t11, t12, t13, t14, t15, t16;
176
177 t1 = pow(grada,2.0);
178 t2 = 1/pow(rhoa,2.666666666666667);
179 t3 = 0.007*t1*t2+1.0;
180 t4 = 1/pow(t3,0.8);
181 t5 = 0.00375*t1*t4*t2-0.75*pow(6.0,0.333333333333333)/
182 pow(3.141592653589793,0.333333333333333);
183 t6 = pow(rhoa,0.333333333333333);
184 t7 = pow(grada,4.0);
185 t8 = 1/pow(t3,1.8);
186 t9 = 1/pow(rhoa,6.333333333333333);
187 t10 = 1/pow(rhoa,3.666666666666667);
188 t11 = 5.6E-5*t7*t8*t9-0.01*t1*t4*t10;
189 t12 = pow(rhoa,1.333333333333333);
190 t13 = pow(grada,3.0);
191 t14 = 1/pow(rhoa,5.333333333333333);
192 t15 = 0.0075*grada*t4*t2-4.2E-5*t13*t8*t14;
193 t16 = 1/pow(t3,2.8);
194
195 /* code */
196 res[0] = -1.333333333333333*t5*t6-1.0*t11*t12;
197 res[1] = -1.0*t12*t15;
198 res[2] = -1.0*t12*(1.8816E-6*t16*pow(grada,6.0)/pow(rhoa,
199 10.0)-5.04E-4*t7*t8/pow(rhoa,7.333333333333333)+0.036666666666667*
200 t1*t4/pow(rhoa,4.666666666666667))-0.444444444444444*t5/pow(rhoa,
201 0.666666666666667)-2.666666666666667*t11*t6;
202 res[3] = -1.0*t12*(-1.4112E-6*t16*pow(grada,5.0)/pow(rhoa,
203 9.0)+3.36E-4*t13*t8*t9-0.02*grada*t4*t10)-1.333333333333333*
204 t15*t6;
205 res[4] = -1.0*t12*(1.0584000000000001E-6*t16*t7/pow(rhoa,
206 8.0)+0.0075*t4*t2-2.1000000000000004E-4*t1*t8*t14);
207
208 }
209
210 static void
b86mx_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)211 b86mx_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
212 {
213 real res[5];
214
215 b86mx_second_helper(dp->rhoa, dp->grada, res);
216
217 ds->df1000 += factor*res[0];
218 ds->df0010 += factor*res[1];
219
220 ds->df2000 += factor*res[2];
221 ds->df1010 += factor*res[3];
222 ds->df0020 += factor*res[4];
223
224
225 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
226 fabs(dp->grada-dp->gradb)>1e-13)
227 b86mx_second_helper(dp->rhob, dp->gradb, res);
228 ds->df0100 += factor*res[0];
229 ds->df0001 += factor*res[1];
230
231 ds->df0200 += factor*res[2];
232 ds->df0101 += factor*res[3];
233 ds->df0002 += factor*res[4];
234
235 }
236
237 static void
b86mx_third_helper(real rhoa,real grada,real * res)238 b86mx_third_helper(real rhoa, real grada, real *res)
239 {
240 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
241 real t11, t12, t13, t14, t15, t16, t17, t18;
242 real t19, t20, t21, t22, t23, t24, t25, t26;
243 real t27, t28;
244
245 t1 = pow(grada,2.0);
246 t2 = 1/pow(rhoa,2.666666666666667);
247 t3 = 0.007*t1*t2+1.0;
248 t4 = 1/pow(t3,0.8);
249 t5 = 0.00375*t1*t4*t2-0.75*pow(6.0,0.333333333333333)/
250 pow(3.141592653589793,0.333333333333333);
251 t6 = pow(rhoa,0.333333333333333);
252 t7 = pow(grada,4.0);
253 t8 = 1/pow(t3,1.8);
254 t9 = 1/pow(rhoa,6.333333333333333);
255 t10 = 1/pow(rhoa,3.666666666666667);
256 t11 = 5.6E-5*t7*t8*t9-0.01*t1*t4*t10;
257 t12 = pow(rhoa,1.333333333333333);
258 t13 = pow(grada,3.0);
259 t14 = 1/pow(rhoa,5.333333333333333);
260 t15 = 0.0075*grada*t4*t2-4.2E-5*t13*t8*t14;
261 t16 = 1/pow(rhoa,0.666666666666667);
262 t17 = pow(grada,6.0);
263 t18 = 1/pow(t3,2.8);
264 t19 = 1/pow(rhoa,10.0);
265 t20 = 1/pow(rhoa,7.333333333333333);
266 t21 = 1/pow(rhoa,4.666666666666667);
267 t22 = 0.036666666666667*t1*t4*t21-5.04E-4*t7*t8*t20+1.8816E-6*
268 t17*t18*t19;
269 t23 = pow(grada,5.0);
270 t24 = 1/pow(rhoa,9.0);
271 t25 = -0.02*grada*t4*t10+3.36E-4*t13*t8*t9-1.4112E-6*
272 t23*t18*t24;
273 t26 = 1/pow(rhoa,8.0);
274 t27 = 0.0075*t4*t2-2.1000000000000004E-4*t1*t8*t14+1.0584000000000001E-6*
275 t7*t18*t26;
276 t28 = 1/pow(t3,3.8);
277
278 /* code */
279 res[0] = -1.333333333333333*t5*t6-1.0*t11*t12;
280 res[1] = -1.0*t12*t15;
281 res[2] = -2.666666666666667*t11*t6-0.444444444444444*
282 t16*t5-1.0*t12*t22;
283 res[3] = -1.333333333333333*t15*t6-1.0*t12*t25;
284 res[4] = -1.0*t12*t27;
285 res[5] = -1.0*t12*(9.834495999999997E-8*t28*pow(grada,
286 8.0)/pow(rhoa,13.66666666666667)-3.57504E-5*t17*t18/pow(rhoa,
287 11.0)+0.004243555555556*t7*t8/pow(rhoa,8.333333333333334)-
288 0.171111111111111*t1*t4/pow(rhoa,5.666666666666667))+0.296296296296296*
289 t5/pow(rhoa,1.666666666666667)-4.0*t22*t6-1.333333333333333*
290 t11*t16;
291 res[6] = -1.0*t12*(-7.375872E-8*t28*pow(grada,7.0)/pow(rhoa,
292 12.66666666666667)+0.073333333333333*grada*t4*t21-0.002426666666667*
293 t13*t8*t20+2.39904E-5*t23*t18*t19)-2.666666666666667*t25*t6-
294 0.444444444444444*t15*t16;
295 res[7] = -1.0*t12*(5.531903999999999E-8*t17*t28/pow(rhoa,
296 11.66666666666667)+0.001232*t1*t8*t9-1.55232E-5*t7*t18*t24-
297 0.02*t4*t10)-1.333333333333333*t27*t6;
298 res[8] = -1.0*t12*(-4.148928000000001E-8*t23*t28/pow(rhoa,
299 10.66666666666667)+9.525600000000002E-6*t13*t18*t26-5.040000000000001E-4*
300 grada*t8*t14);
301
302 }
303
304 static void
b86mx_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)305 b86mx_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
306 {
307 real res[9];
308
309 b86mx_third_helper(dp->rhoa, dp->grada, res);
310
311 ds->df1000 += factor*res[0];
312 ds->df0010 += factor*res[1];
313
314 ds->df2000 += factor*res[2];
315 ds->df1010 += factor*res[3];
316 ds->df0020 += factor*res[4];
317
318 ds->df3000 += factor*res[5];
319 ds->df2010 += factor*res[6];
320 ds->df1020 += factor*res[7];
321 ds->df0030 += factor*res[8];
322
323
324 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
325 fabs(dp->grada-dp->gradb)>1e-13)
326 b86mx_third_helper(dp->rhob, dp->gradb, res);
327
328 ds->df0100 += factor*res[0];
329 ds->df0001 += factor*res[1];
330
331 ds->df0200 += factor*res[2];
332 ds->df0101 += factor*res[3];
333 ds->df0002 += factor*res[4];
334
335 ds->df0300 += factor*res[5];
336 ds->df0201 += factor*res[6];
337 ds->df0102 += factor*res[7];
338 ds->df0003 += factor*res[8];
339
340 }
341
342 static void
b86mx_fourth_helper(real rhoa,real grada,real * res)343 b86mx_fourth_helper(real rhoa, real grada, real *res)
344 {
345 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
346 real t11, t12, t13, t14, t15, t16, t17, t18;
347 real t19, t20, t21, t22, t23, t24, t25, t26;
348 real t27, t28, t29, t30, t31, t32, t33, t34;
349 real t35, t36, t37, t38, t39, t40, t41, t42;
350 real t43;
351
352 t1 = pow(grada,2.0);
353 t2 = 1/pow(rhoa,2.666666666666667);
354 t3 = 0.007*t1*t2+1.0;
355 t4 = 1/pow(t3,0.8);
356 t5 = 0.00375*t1*t4*t2-0.75*pow(6.0,0.333333333333333)/
357 pow(3.141592653589793,0.333333333333333);
358 t6 = pow(rhoa,0.333333333333333);
359 t7 = pow(grada,4.0);
360 t8 = 1/pow(t3,1.8);
361 t9 = 1/pow(rhoa,6.333333333333333);
362 t10 = 1/pow(rhoa,3.666666666666667);
363 t11 = 5.6E-5*t7*t8*t9-0.01*t1*t4*t10;
364 t12 = pow(rhoa,1.333333333333333);
365 t13 = pow(grada,3.0);
366 t14 = 1/pow(rhoa,5.333333333333333);
367 t15 = 0.0075*grada*t4*t2-4.2E-5*t13*t8*t14;
368 t16 = 1/pow(rhoa,0.666666666666667);
369 t17 = pow(grada,6.0);
370 t18 = 1/pow(t3,2.8);
371 t19 = 1/pow(rhoa,10.0);
372 t20 = 1/pow(rhoa,7.333333333333333);
373 t21 = 1/pow(rhoa,4.666666666666667);
374 t22 = 0.036666666666667*t1*t4*t21-5.04E-4*t7*t8*t20+1.8816E-6*
375 t17*t18*t19;
376 t23 = pow(grada,5.0);
377 t24 = 1/pow(rhoa,9.0);
378 t25 = -0.02*grada*t4*t10+3.36E-4*t13*t8*t9-1.4112E-6*
379 t23*t18*t24;
380 t26 = 1/pow(rhoa,8.0);
381 t27 = 0.0075*t4*t2-2.1000000000000004E-4*t1*t8*t14+1.0584000000000001E-6*
382 t7*t18*t26;
383 t28 = 1/pow(rhoa,1.666666666666667);
384 t29 = pow(grada,8.0);
385 t30 = 1/pow(t3,3.8);
386 t31 = 1/pow(rhoa,13.66666666666667);
387 t32 = 1/pow(rhoa,11.0);
388 t33 = 1/pow(rhoa,8.333333333333334);
389 t34 = 1/pow(rhoa,5.666666666666667);
390 t35 = -0.171111111111111*t1*t4*t34+0.004243555555556*
391 t7*t8*t33-3.57504E-5*t17*t18*t32+9.834495999999997E-8*t29*
392 t30*t31;
393 t36 = pow(grada,7.0);
394 t37 = 1/pow(rhoa,12.66666666666667);
395 t38 = 0.073333333333333*grada*t4*t21-0.002426666666667*
396 t13*t8*t20+2.39904E-5*t23*t18*t19-7.375872E-8*t36*t30*t37;
397 t39 = 1/
398 pow(rhoa,11.66666666666667);
399 t40 = -0.02*t4*t10+0.001232*t1*t8*t9-1.55232E-5*t7*t18*
400 t24+5.531903999999999E-8*t17*t30*t39;
401 t41 = 1/pow(rhoa,10.66666666666667);
402 t42 = -5.040000000000001E-4*grada*t8*t14+9.525600000000002E-6*
403 t13*t18*t26-4.148928000000001E-8*t23*t30*t41;
404 t43 = 1/pow(t3,4.8);
405
406 /* code */
407 res[0] = -1.333333333333333*t5*t6-1.0*t11*t12;
408 res[1] = -1.0*t12*t15;
409 res[2] = -2.666666666666667*t11*t6-0.444444444444444*
410 t16*t5-1.0*t12*t22;
411 res[3] = -1.333333333333333*t15*t6-1.0*t12*t25;
412 res[4] = -1.0*t12*t27;
413 res[5] = -4.0*t22*t6+0.296296296296296*t28*t5-1.0*t12*
414 t35-1.333333333333333*t11*t16;
415 res[6] = -2.666666666666667*t25*t6-1.0*t12*t38-0.444444444444444*
416 t15*t16;
417 res[7] = -1.333333333333333*t27*t6-1.0*t12*t40;
418 res[8] = -1.0*t12*t42;
419 res[9] = -1.0*t12*(6.97593582933333E-9*t43*pow(grada,
420 10.0)/pow(rhoa,17.33333333333333)-3.212602026666666E-6*t29*
421 t30/pow(rhoa,14.66666666666667)+5.358378666666667E-4*t17*t18/
422 pow(rhoa,12.0)-0.037918222222222*t7*t8/pow(rhoa,9.333333333333334)+
423 0.96962962962963*t1*t4/pow(rhoa,6.666666666666667))-5.333333333333333*
424 t35*t6-0.493827160493827*t2*t5+1.185185185185185*t11*t28-2.666666666666667*
425 t16*t22;
426 res[10] = -1.0*t12*(-5.231951871999998E-9*t43*pow(grada,
427 9.0)/pow(rhoa,16.33333333333333)-0.342222222222222*grada*t4*
428 t34+0.018890666666667*t13*t8*t33-3.214399999999999E-4*t23*
429 t18*t32+2.18817536E-6*t36*t30*t31)-4.0*t38*t6+0.296296296296296*
430 t15*t28-1.333333333333333*t16*t25;
431 res[11] = -1.0*t12*(3.923963904E-9*t29*t43/pow(rhoa,15.33333333333333)-
432 1.45673472E-6*t17*t30*t37+0.073333333333333*t4*t21-0.008101333333333*
433 t1*t8*t20+1.81104E-4*t7*t18*t19)-2.666666666666667*t40*t6-
434 0.444444444444444*t16*t27;
435 res[12] = -1.0*t12*(-2.942972928E-9*t36*t43/pow(rhoa,
436 14.33333333333333)+0.002688*grada*t8*t9+9.4042368E-7*t23*t30*
437 t39-9.31392E-5*t13*t18*t24)-1.333333333333333*t42*t6;
438 res[13] = -1.0*t12*(2.207229696E-9*t17*t43/pow(rhoa,13.33333333333333)-
439 5.8084992E-7*t7*t30*t41+4.127760000000001E-5*t1*t18*t26-5.040000000000001E-4*
440 t8*t14);
441
442 }
443
444 static void
b86mx_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)445 b86mx_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
446 {
447 real res[14];
448
449 b86mx_fourth_helper(dp->rhoa, dp->grada, res);
450
451 ds->df1000 += factor*res[0];
452 ds->df0010 += factor*res[1];
453
454 ds->df2000 += factor*res[2];
455 ds->df1010 += factor*res[3];
456 ds->df0020 += factor*res[4];
457
458 ds->df3000 += factor*res[5];
459 ds->df2010 += factor*res[6];
460 ds->df1020 += factor*res[7];
461 ds->df0030 += factor*res[8];
462
463 ds->df4000 += factor*res[9];
464 ds->df3010 += factor*res[10];
465 ds->df2020 += factor*res[11];
466 ds->df1030 += factor*res[12];
467 ds->df0040 += factor*res[13];
468
469
470 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
471 fabs(dp->grada-dp->gradb)>1e-13)
472 b86mx_fourth_helper(dp->rhob, dp->gradb, res);
473
474 ds->df0100 += factor*res[0];
475 ds->df0001 += factor*res[1];
476
477 ds->df0200 += factor*res[2];
478 ds->df0101 += factor*res[3];
479 ds->df0002 += factor*res[4];
480
481 ds->df0300 += factor*res[5];
482 ds->df0201 += factor*res[6];
483 ds->df0102 += factor*res[7];
484 ds->df0003 += factor*res[8];
485
486 ds->df0400 += factor*res[9];
487 ds->df0301 += factor*res[10];
488 ds->df0202 += factor*res[11];
489 ds->df0103 += factor*res[12];
490 ds->df0004 += factor*res[13];
491
492 }
493