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