1 /*-*-mode:
2 
3 !
4 !  Dalton, a molecular electronic structure program
5 !  Copyright (C) by the authors of Dalton.
6 !
7 !  This program is free software; you can redistribute it and/or
8 !  modify it under the terms of the GNU Lesser General Public
9 !  License version 2.1 as published by the Free Software Foundation.
10 !
11 !  This program is distributed in the hope that it will be useful,
12 !  but WITHOUT ANY WARRANTY; without even the implied warranty of
13 !  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 !  Lesser General Public License for more details.
15 !
16 !  If a copy of the GNU LGPL v2.1 was not distributed with this
17 !  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
18 !
19 
20 !
21 
22 */
23 /* fun-pw91x.c:
24 
25    NOTE: This is the equivalent of the pwggaIIx functional.
26 
27    Automatically generated code implementing PW91X 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 PW91xFunctional;" to 'functionals.h'
34     2. add "&PW91xFunctional," to 'functionals.c'
35     3. add "fun-PW91x.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.
36 
37     This functional has been generated from following input:
38     ------ cut here -------
39 rho:  rhoa + rhob;
40 grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
41 zeta: (rhoa-rhob)/(rhoa+rhob);
42 
43 kfa: (6*%PI*%PI*rhoa)^(1/3);
44 kfb: (6*%PI*%PI*rhob)^(1/3);
45 exa: -3*kfa/(4*%PI);
46 exb: -3*kfb/(4*%PI);
47 
48 sa:  grada/(2*kfa*rhoa);
49 sb:  gradb/(2*kfb*rhob);
50 F0a: 1+0.19645*sa*asinh(7.7956*sa);
51 F0b: 1+0.19645*sb*asinh(7.7956*sb);
52 F1a: F0a + (0.2743-0.1508*exp(-100*sa^2))*sa^2;
53 F1b: F0b + (0.2743-0.1508*exp(-100*sb^2))*sb^2;
54 F2a: F0a + 0.004*sa^4;
55 F2b: F0b + 0.004*sb^4;
56 Fa: F1a/F2a;
57 Fb: F1b/F2b;
58 
59 Exa: 2*rhoa*exa*Fa;
60 Exb: 2*rhob*exb*Fb;
61 
62 K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb)/2;
63 
64 
65     ------ cut here -------
66 */
67 
68 
69 /* strictly conform to XOPEN ANSI C standard */
70 #if !defined(SYS_DEC)
71 /* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
72  * is not specified. */
73 #define _XOPEN_SOURCE          500
74 #define _XOPEN_SOURCE_EXTENDED 1
75 #endif
76 #include <math.h>
77 #include <stddef.h>
78 #include "general.h"
79 
80 #define __CVERSION__
81 
82 #include "functionals.h"
83 
84 /* INTERFACE PART */
pw91x_isgga(void)85 static integer pw91x_isgga(void) { return 1; }
86 static integer pw91x_read(const char *conf_line);
87 static real pw91x_energy(const FunDensProp* dp);
88 static void pw91x_first(FunFirstFuncDrv *ds,   real factor,
89                          const FunDensProp* dp);
90 static void pw91x_second(FunSecondFuncDrv *ds, real factor,
91                           const FunDensProp* dp);
92 static void pw91x_third(FunThirdFuncDrv *ds,   real factor,
93                          const FunDensProp* dp);
94 static void pw91x_fourth(FunFourthFuncDrv *ds,   real factor,
95                           const FunDensProp* dp);
96 
97 Functional PW91xFunctional = {
98   "PW91x",       /* name */
99   pw91x_isgga,   /* gga-corrected */
100    1,
101   pw91x_read,
102   NULL,
103   pw91x_energy,
104   pw91x_first,
105   pw91x_second,
106   pw91x_third,
107   pw91x_fourth
108 };
109 
110 /* IMPLEMENTATION PART */
111 static integer
pw91x_read(const char * conf_line)112 pw91x_read(const char *conf_line)
113 {
114     fun_set_hf_weight(0);
115     return 1;
116 }
117 
118 static real
pw91x_energy(const FunDensProp * dp)119 pw91x_energy(const FunDensProp *dp)
120 {
121     real res;
122     real rhoa = dp->rhoa, rhob = dp->rhob;
123     real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
124 
125     real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
126     real t11, t12, t13, t14, t15, t16, t17;
127 
128     t1 = pow(6.0,0.333333333333333);
129     t2 = 1/pow(3.141592653589793,0.333333333333333);
130     t3 = 1/t1;
131     t4 = 1/pow(3.141592653589793,2.666666666666667);
132     t5 = 1/pow(3.141592653589793,0.666666666666667);
133     t6 = pow(rhoa,1.333333333333333);
134     t7 = 1/t6;
135     t8 = 0.098225*t3*t5*grada*asinh(3.8978*t3*t5*grada*t7)*
136         t7;
137     t9 = 1/pow(6.0,0.666666666666667);
138     t10 = 1/pow(3.141592653589793,1.333333333333333);
139     t11 = pow(grada,2.0);
140     t12 = 1/pow(rhoa,2.666666666666667);
141     t13 = pow(rhob,1.333333333333333);
142     t14 = 1/t13;
143     t15 = 0.098225*t3*t5*gradb*asinh(3.8978*t3*t5*gradb*t14)*
144         t14;
145     t16 = pow(gradb,2.0);
146     t17 = 1/pow(rhob,2.666666666666667);
147 
148    /* code */
149     res = 0.5*(-1.5*t1*t13*t2*(0.25*(0.2743-0.1508/pow(2.718281828459045,
150         25.0*t10*t16*t17*t9))*t10*t16*t17*t9+t15+1.0)/(4.166666666666667E-5*
151         t3*t4*pow(gradb,4.0)/pow(rhob,5.333333333333333)+t15+1.0)-
152         1.5*t1*t2*t6*(0.25*(0.2743-0.1508/pow(2.718281828459045,25.0*
153         t10*t11*t12*t9))*t10*t11*t12*t9+t8+1.0)/(4.166666666666667E-5*
154         t3*t4*pow(grada,4.0)/pow(rhoa,5.333333333333333)+t8+1.0));
155 
156     return res;
157 }
158 
159 static void
pw91x_first_helper(real rhoa,real grada,real * res)160 pw91x_first_helper(real rhoa, real grada, real *res)
161 {    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
162     real t11, t12, t13, t14, t15, t16, t17, t18;
163     real t19, t20, t21, t22, t23, t24, t25, t26;
164     real t27, t28, t29;
165 
166     t1 = pow(6.0,0.333333333333333);
167     t2 = 1/pow(3.141592653589793,0.333333333333333);
168     t3 = 1/t1;
169     t4 = 1/pow(3.141592653589793,2.666666666666667);
170     t5 = pow(grada,4.0);
171     t6 = 1/pow(rhoa,5.333333333333333);
172     t7 = 1/pow(3.141592653589793,0.666666666666667);
173     t8 = pow(rhoa,1.333333333333333);
174     t9 = 1/t8;
175     t10 = asinh(3.8978*t3*t7*grada*t9);
176     t11 = 0.098225*t3*t7*grada*t10*t9;
177     t12 = 4.166666666666667E-5*t3*t4*t5*t6+t11+1.0;
178     t13 = 1/t12;
179     t14 = 1/pow(6.0,0.666666666666667);
180     t15 = 1/pow(3.141592653589793,1.333333333333333);
181     t16 = pow(grada,2.0);
182     t17 = 1/pow(rhoa,2.666666666666667);
183     t18 = 1/pow(2.718281828459045,25.0*t14*t15*t16*t17);
184     t19 = 0.2743-
185         0.1508*t18;
186     t20 = 0.25*t14*t15*t16*t17*t19+t11+1.0;
187     t21 = 1/pow(rhoa,6.333333333333333);
188     t22 = 1/sqrt(15.19284484*t14*t15*t16*t17+1.0);
189     t23 = 1/pow(rhoa,3.666666666666667);
190     t24 = -0.510481873333333*t14*t15*t16*t22*t23;
191     t25 = -0.130966666666667*t10*t3*t7*grada/pow(rhoa,2.333333333333333);
192     t26 = 1/
193         pow(t12,2.0);
194     t27 = pow(grada,3.0);
195     t28 = 0.382861405*t14*t15*grada*t22*t17;
196     t29 = 0.098225*t3*t7*t10*t9;
197 
198    /* code */
199     res[0] = 0.5*(-2.0*t1*t13*t2*t20*pow(rhoa,0.333333333333333)+
200         1.5*t1*t2*t20*(t25+t24-2.222222222222222E-4*t3*t4*t5*t21)*
201         t26*t8-1.5*t1*t13*t2*(t25+t24-0.666666666666667*t14*t15*t16*
202         t19*t23-0.418888888888889*t3*t4*t5*t21*t18)*t8);
203     res[1] = 0.5*(1.5*t1*t2*t20*t26*(t29+t28+1.6666666666666666E-4*
204         t3*t4*t27*t6)*t8-1.5*t1*t13*t2*t8*(0.5*t14*t15*t17*t19*grada+
205         t29+t28+0.314166666666667*t3*t4*t27*t6*t18));
206 }
207 
208 static void
pw91x_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)209 pw91x_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
210 {
211     real res[2];
212 
213     pw91x_first_helper(dp->rhoa, dp->grada, res);
214    /* Final assignment */
215     ds->df1000 += factor*res[0];
216     ds->df0010 += factor*res[1];
217 
218 
219     if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
220        fabs(dp->grada-dp->gradb)>1e-13)
221         pw91x_first_helper(dp->rhob, dp->gradb, res);
222     ds->df0100 += factor*res[0];
223     ds->df0001 += factor*res[1];
224 
225 }
226 
227 static void
pw91x_second_helper(real rhoa,real grada,real * res)228 pw91x_second_helper(real rhoa, real grada, real *res)
229 {
230     real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
231     real t11, t12, t13, t14, t15, t16, t17, t18;
232     real t19, t20, t21, t22, t23, t24, t25, t26;
233     real t27, t28, t29, t30, t31, t32, t33, t34;
234     real t35, t36, t37, t38, t39, t40, t41, t42;
235     real t43, t44, t45, t46, t47, t48, t49;
236 
237     t1 = pow(6.0,0.333333333333333);
238     t2 = 1/pow(3.141592653589793,0.333333333333333);
239     t3 = 1/t1;
240     t4 = 1/pow(3.141592653589793,2.666666666666667);
241     t5 = pow(grada,4.0);
242     t6 = 1/pow(rhoa,5.333333333333333);
243     t7 = 1/pow(3.141592653589793,0.666666666666667);
244     t8 = pow(rhoa,1.333333333333333);
245     t9 = 1/t8;
246     t10 = asinh(3.8978*t3*t7*grada*t9);
247     t11 = 0.098225*t3*t7*grada*t10*t9;
248     t12 = 4.166666666666667E-5*t3*t4*t5*t6+t11+1.0;
249     t13 = 1/t12;
250     t14 = pow(rhoa,0.333333333333333);
251     t15 = 1/pow(6.0,0.666666666666667);
252     t16 = 1/pow(3.141592653589793,1.333333333333333);
253     t17 = pow(grada,2.0);
254     t18 = 1/pow(rhoa,2.666666666666667);
255     t19 = 1/pow(2.718281828459045,25.0*t15*t16*t17*t18);
256     t20 = 0.2743-
257         0.1508*t19;
258     t21 = 0.25*t15*t16*t17*t18*t20+t11+1.0;
259     t22 = 1/pow(rhoa,6.333333333333333);
260     t23 = sqrt(15.19284484*t15*t16*t17*t18+1.0);
261     t24 = 1/t23;
262     t25 = 1/pow(rhoa,3.666666666666667);
263     t26 = -0.510481873333333*t15*t16*t17*t24*t25;
264     t27 = 1/pow(rhoa,2.333333333333333);
265     t28 = -0.130966666666667*t3*t7*grada*t10*t27;
266     t29 = t28+t26-2.222222222222222E-4*t3*t4*t5*t22;
267     t30 = 1/pow(t12,2.0);
268     t31 = t28+t26-0.666666666666667*t15*t16*t17*t20*t25-0.418888888888889*
269         t3*t4*t5*t22*t19;
270     t32 = pow(grada,3.0);
271     t33 = 0.382861405*t15*t16*grada*t24*t18;
272     t34 = 0.098225*t3*t7*t10*t9;
273     t35 = t34+t33+1.6666666666666666E-4*t3*t4*t32*t6;
274     t36 = 0.5*t15*t16*t18*t20*grada+t34+t33+0.314166666666667*
275         t3*t4*t32*t6*t19;
276     t37 = 1/pow(t12,3.0);
277     t38 = 1/pow(rhoa,7.333333333333333);
278     t39 = 1/pow(t23,3.0);
279     t40 = -1.723482643374637*t3*t4*t5*t39*t38;
280     t41 = 1/pow(rhoa,4.666666666666667);
281     t42 = 2.552409366666666*t15*t16*t17*t24*t41;
282     t43 = 0.305588888888889*t10*t3*t7*grada/pow(rhoa,3.333333333333333);
283     t44 = 1/
284         pow(3.141592653589793,4.0);
285     t45 = 1.292611982530978*t3*t4*t32*t39*t22;
286     t46 = -1.53144562*t15*t16*grada*t24*t25;
287     t47 = -0.130966666666667*t3*t7*t10*t27;
288     t48 = -0.969458986898234*t3*t4*t17*t39*t6;
289     t49 = 0.76572281*t15*t16*t24*t18;
290 
291    /* code */
292     res[0] = 0.5*(-1.5*t1*t13*t2*t31*t8+1.5*t1*t2*t21*t29*
293         t30*t8-2.0*t1*t13*t14*t2*t21);
294     res[1] = 0.5*(1.5*t1*t2*t21*t30*t35*t8-1.5*t1*t13*t2*
295         t36*t8);
296     res[2] = 0.5*(-1.5*t1*t13*t2*t8*(-4.65432098765432*t19*
297         t44*pow(grada,6.0)/pow(rhoa,10.0)+t43+t42+2.444444444444445*
298         t15*t16*t17*t20*t41+t40+3.77*t3*t4*t5*t38*t19)-0.666666666666667*
299         t1*t13*t2*t21/pow(rhoa,0.666666666666667)+1.5*t1*t2*t21*t30*
300         (t43+t42+t40+0.001407407407407*t3*t4*t5*t38)*t8-3.0*t1*t2*
301         t21*pow(t29,2.0)*t37*t8+3.0*t1*t2*t29*t30*t31*t8-4.0*t1*t13*
302         t14*t2*t31+4.0*t1*t14*t2*t21*t29*t30);
303     res[3] = 0.5*(-1.5*t1*t13*t2*t8*(3.490740740740741*t19*
304         t44*pow(grada,5.0)/pow(rhoa,9.0)-1.333333333333333*t15*t16*
305         t20*t25*grada+t47+t46+t45-2.513333333333333*t3*t4*t32*t22*
306         t19)+1.5*t1*t2*t21*t30*(t47+t46+t45-8.888888888888888E-4*t3*
307         t4*t32*t22)*t8-3.0*t1*t2*t21*t29*t35*t37*t8+1.5*t1*t2*t29*
308         t30*t36*t8+1.5*t1*t2*t30*t31*t35*t8-2.0*t1*t13*t14*t2*t36+
309         2.0*t1*t14*t2*t21*t30*t35);
310     res[4] = 0.5*(-1.5*t1*t13*t2*t8*(-2.618055555555555*t19*
311         t44*t5/pow(rhoa,8.0)+t49+t48+0.5*t15*t16*t18*t20+1.570833333333333*
312         t3*t4*t17*t6*t19)+1.5*t1*t2*t21*t30*(t49+t48+5.0E-4*t3*t4*
313         t17*t6)*t8-3.0*t1*t2*t21*pow(t35,2.0)*t37*t8+3.0*t1*t2*t30*
314         t35*t36*t8);
315 
316 }
317 
318 static void
pw91x_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)319 pw91x_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
320 {
321     real res[5];
322 
323     pw91x_second_helper(dp->rhoa, dp->grada, res);
324 
325     ds->df1000 += factor*res[0];
326     ds->df0010 += factor*res[1];
327 
328     ds->df2000 += factor*res[2];
329     ds->df1010 += factor*res[3];
330     ds->df0020 += factor*res[4];
331 
332 
333     if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
334        fabs(dp->grada-dp->gradb)>1e-13)
335         pw91x_second_helper(dp->rhob, dp->gradb, res);
336     ds->df0100 += factor*res[0];
337     ds->df0001 += factor*res[1];
338 
339     ds->df0200 += factor*res[2];
340     ds->df0101 += factor*res[3];
341     ds->df0002 += factor*res[4];
342 
343 }
344 
345 static void
pw91x_third_helper(real rhoa,real grada,real * res)346 pw91x_third_helper(real rhoa, real grada, real *res)
347 {
348     real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
349     real t11, t12, t13, t14, t15, t16, t17, t18;
350     real t19, t20, t21, t22, t23, t24, t25, t26;
351     real t27, t28, t29, t30, t31, t32, t33, t34;
352     real t35, t36, t37, t38, t39, t40, t41, t42;
353     real t43, t44, t45, t46, t47, t48, t49, t50;
354     real t51, t52, t53, t54, t55, t56, t57, t58;
355     real t59, t60, t61, t62, t63, t64, t65, t66;
356     real t67, t68, t69, t70, t71, t72, t73, t74;
357     real t75, t76, t77, t78, t79, t80, t81, t82;
358     real t83;
359 
360     t1 = pow(6.0,0.333333333333333);
361     t2 = 1/pow(3.141592653589793,0.333333333333333);
362     t3 = 1/t1;
363     t4 = 1/pow(3.141592653589793,2.666666666666667);
364     t5 = pow(grada,4.0);
365     t6 = 1/pow(rhoa,5.333333333333333);
366     t7 = 1/pow(3.141592653589793,0.666666666666667);
367     t8 = pow(rhoa,1.333333333333333);
368     t9 = 1/t8;
369     t10 = asinh(3.8978*t3*t7*grada*t9);
370     t11 = 0.098225*t3*t7*grada*t10*t9;
371     t12 = 4.166666666666667E-5*t3*t4*t5*t6+t11+1.0;
372     t13 = 1/t12;
373     t14 = pow(rhoa,0.333333333333333);
374     t15 = 1/pow(6.0,0.666666666666667);
375     t16 = 1/pow(3.141592653589793,1.333333333333333);
376     t17 = pow(grada,2.0);
377     t18 = 1/pow(rhoa,2.666666666666667);
378     t19 = 1/pow(2.718281828459045,25.0*t15*t16*t17*t18);
379     t20 = 0.2743-
380         0.1508*t19;
381     t21 = 0.25*t15*t16*t17*t18*t20+t11+1.0;
382     t22 = 1/pow(rhoa,6.333333333333333);
383     t23 = sqrt(15.19284484*t15*t16*t17*t18+1.0);
384     t24 = 1/t23;
385     t25 = 1/pow(rhoa,3.666666666666667);
386     t26 = -0.510481873333333*t15*t16*t17*t24*t25;
387     t27 = 1/pow(rhoa,2.333333333333333);
388     t28 = -0.130966666666667*t3*t7*grada*t10*t27;
389     t29 = t28+t26-2.222222222222222E-4*t3*t4*t5*t22;
390     t30 = 1/pow(t12,2.0);
391     t31 = t28+t26-0.666666666666667*t15*t16*t17*t20*t25-0.418888888888889*
392         t3*t4*t5*t22*t19;
393     t32 = pow(grada,3.0);
394     t33 = 0.382861405*t15*t16*grada*t24*t18;
395     t34 = 0.098225*t3*t7*t10*t9;
396     t35 = t34+t33+1.6666666666666666E-4*t3*t4*t32*t6;
397     t36 = 0.5*t15*t16*t18*t20*grada+t34+t33+0.314166666666667*
398         t3*t4*t32*t6*t19;
399     t37 = 1/pow(rhoa,0.666666666666667);
400     t38 = pow(t29,2.0);
401     t39 = 1/pow(t12,3.0);
402     t40 = 1/pow(rhoa,7.333333333333333);
403     t41 = 1/pow(t23,3.0);
404     t42 = -1.723482643374637*t3*t4*t5*t41*t40;
405     t43 = 1/pow(rhoa,4.666666666666667);
406     t44 = 2.552409366666666*t15*t16*t17*t24*t43;
407     t45 = 1/pow(rhoa,3.333333333333333);
408     t46 = 0.305588888888889*t3*t7*grada*t10*t45;
409     t47 = t46+t44+t42+0.001407407407407*t3*t4*t5*t40;
410     t48 = 1/pow(3.141592653589793,4.0);
411     t49 = pow(grada,6.0);
412     t50 = 1/pow(rhoa,10.0);
413     t51 = t46+t44+2.444444444444445*t15*t16*t17*t20*t43+t42-
414         4.65432098765432*t48*t49*t50*t19+3.77*t3*t4*t5*t40*t19;
415     t52 = 1.292611982530978*t3*t4*t32*t41*t22;
416     t53 = -1.53144562*t15*t16*grada*t24*t25;
417     t54 = -0.130966666666667*t3*t7*t10*t27;
418     t55 = t54+t53+t52-8.888888888888888E-4*t3*t4*t32*t22;
419     t56 = pow(grada,
420         5.0);
421     t57 = 1/pow(rhoa,9.0);
422     t58 = -1.333333333333333*t15*t16*t20*t25*grada+t54+t53+
423         t52+3.490740740740741*t48*t56*t57*t19-2.513333333333333*t3*
424         t4*t32*t22*t19;
425     t59 = pow(t35,2.0);
426     t60 = -0.969458986898234*t3*t4*t17*t41*t6;
427     t61 = 0.76572281*t15*t16*t24*t18;
428     t62 = t61+t60+5.0E-4*t3*t4*t17*t6;
429     t63 = 1/pow(rhoa,8.0);
430     t64 = t61+t60+0.5*t15*t16*t18*t20-2.618055555555555*t48*
431         t5*t63*t19+1.570833333333333*t3*t4*t17*t6*t19;
432     t65 = 1/pow(t12,4.0);
433     t66 = 1/pow(t23,5.0);
434     t67 = 1/pow(rhoa,11.0);
435     t68 = -17.45640292348261*t48*t49*t66*t67;
436     t69 = 1/pow(rhoa,8.333333333333334);
437     t70 = 21.25628593495385*t3*t4*t5*t41*t69;
438     t71 = 1/pow(rhoa,5.666666666666667);
439     t72 = -13.49940953925926*t15*t16*t17*t24*t71;
440     t73 = -1.01862962962963*t10*t3*t7*grada/pow(rhoa,4.333333333333333);
441     t74 = 1/
442         pow(3.141592653589793,5.333333333333333);
443     t75 = 13.09230219261196*t48*t56*t66*t50;
444     t76 = -13.35699048615344*t3*t4*t32*t41*t40;
445     t77 = 6.295943104444444*t15*t16*grada*t24*t43;
446     t78 = 0.305588888888889*t3*t7*t10*t45;
447     t79 = -9.81922664445897*t48*t5*t66*t57;
448     t80 = 7.755671895185867*t3*t4*t17*t41*t22;
449     t81 = -2.041927493333334*t15*t16*t24*t25;
450     t82 = 7.364419983344228*t48*t32*t66*t63;
451     t83 = -3.877835947592934*t3*t4*grada*t41*t6;
452 
453    /* code */
454     res[0] = 0.5*(-1.5*t1*t13*t2*t31*t8+1.5*t1*t2*t21*t29*
455         t30*t8-2.0*t1*t13*t14*t2*t21);
456     res[1] = 0.5*(1.5*t1*t2*t21*t30*t35*t8-1.5*t1*t13*t2*
457         t36*t8);
458     res[2] = 0.5*(-1.5*t1*t13*t2*t51*t8+1.5*t1*t2*t21*t30*
459         t47*t8-3.0*t1*t2*t21*t38*t39*t8+3.0*t1*t2*t29*t30*t31*t8-0.666666666666667*
460         t1*t13*t2*t21*t37-4.0*t1*t13*t14*t2*t31+4.0*t1*t14*t2*t21*
461         t29*t30);
462     res[3] = 0.5*(-1.5*t1*t13*t2*t58*t8+1.5*t1*t2*t21*t30*
463         t55*t8-3.0*t1*t2*t21*t29*t35*t39*t8+1.5*t1*t2*t29*t30*t36*
464         t8+1.5*t1*t2*t30*t31*t35*t8-2.0*t1*t13*t14*t2*t36+2.0*t1*t14*
465         t2*t21*t30*t35);
466     res[4] = 0.5*(-1.5*t1*t13*t2*t64*t8+1.5*t1*t2*t21*t30*
467         t62*t8-3.0*t1*t2*t21*t39*t59*t8+3.0*t1*t2*t30*t35*t36*t8);
468     res[5] = 0.5*(-1.5*t1*t13*t2*t8*(-310.2880658436214*t15*
469         t19*t74*pow(grada,8.0)/pow(rhoa,13.66666666666667)+t73+t72-
470         11.40740740740741*t15*t16*t17*t20*t71+t70+t68-31.74246913580246*
471         t3*t4*t5*t69*t19+88.43209876543209*t48*t49*t67*t19)+0.444444444444444*
472         t1*t13*t2*t21/pow(rhoa,1.666666666666667)+1.5*t1*t2*t21*t30*
473         (t73+t72+t70-0.010320987654321*t3*t4*t5*t69+t68)*t8+9.0*t1*
474         t2*t21*pow(t29,3.0)*t65*t8+4.5*t1*t2*t29*t30*t51*t8-9.0*t1*
475         t2*t21*t29*t39*t47*t8+4.5*t1*t2*t30*t31*t47*t8-9.0*t1*t2*t31*
476         t38*t39*t8-6.0*t1*t13*t14*t2*t51+6.0*t1*t14*t2*t21*t30*t47-
477         12.0*t1*t14*t2*t21*t38*t39-2.0*t1*t13*t2*t31*t37+2.0*t1*t2*
478         t21*t29*t30*t37+12.0*t1*t14*t2*t29*t30*t31);
479     res[6] = 0.5*(-1.5*t1*t13*t2*t8*(232.716049382716*t15*
480         t19*t74*pow(grada,7.0)/pow(rhoa,12.66666666666667)+4.88888888888889*
481         t15*t16*t20*t43*grada+t78+t77+t76+t75-59.34259259259259*t48*
482         t56*t50*t19+18.15185185185185*t3*t4*t32*t40*t19)+1.5*t1*t2*
483         t21*t30*(t78+t77+t76+0.00562962962963*t3*t4*t32*t40+t75)*t8+
484         9.0*t1*t2*t21*t35*t38*t65*t8+3.0*t1*t2*t29*t30*t58*t8-6.0*
485         t1*t2*t21*t29*t39*t55*t8+3.0*t1*t2*t30*t31*t55*t8+1.5*t1*t2*
486         t30*t35*t51*t8-3.0*t1*t2*t21*t35*t39*t47*t8+1.5*t1*t2*t30*
487         t36*t47*t8-3.0*t1*t2*t36*t38*t39*t8-6.0*t1*t2*t29*t31*t35*
488         t39*t8-4.0*t1*t13*t14*t2*t58+4.0*t1*t14*t2*t21*t30*t55-8.0*
489         t1*t14*t2*t21*t29*t35*t39-0.666666666666667*t1*t13*t2*t36*
490         t37+0.666666666666667*t1*t2*t21*t30*t35*t37+4.0*t1*t14*t2*
491         t29*t30*t36+4.0*t1*t14*t2*t30*t31*t35);
492     res[7] = 0.5*(-1.5*t1*t13*t2*t8*(-174.537037037037*t15*
493         t19*t49*t74/pow(rhoa,11.66666666666667)+t81+t80+t79-1.333333333333333*
494         t15*t16*t20*t25+38.39814814814815*t48*t5*t57*t19-9.215555555555554*
495         t3*t4*t17*t22*t19)+1.5*t1*t2*t21*t30*t8*(t81+t80-0.002666666666667*
496         t3*t4*t17*t22+t79)+9.0*t1*t2*t21*t29*t59*t65*t8+1.5*t1*t2*
497         t29*t30*t64*t8-3.0*t1*t2*t21*t29*t39*t62*t8+1.5*t1*t2*t30*
498         t31*t62*t8-3.0*t1*t2*t31*t39*t59*t8+3.0*t1*t2*t30*t35*t58*
499         t8-6.0*t1*t2*t21*t35*t39*t55*t8+3.0*t1*t2*t30*t36*t55*t8-6.0*
500         t1*t2*t29*t35*t36*t39*t8-2.0*t1*t13*t14*t2*t64+2.0*t1*t14*
501         t2*t21*t30*t62-4.0*t1*t14*t2*t21*t39*t59+4.0*t1*t14*t2*t30*
502         t35*t36);
503     res[8] = 0.5*(-1.5*t1*t13*t2*t8*(130.9027777777778*t15*
504         t19*t56*t74/pow(rhoa,10.66666666666667)+t83+t82-23.5625*t48*
505         t32*t63*t19+3.769999999999999*t3*t4*grada*t6*t19)+1.5*t1*t2*
506         t21*t30*t8*(t83+0.001*t3*t4*grada*t6+t82)+9.0*t1*t2*t21*pow(t35,
507         3.0)*t65*t8+4.5*t1*t2*t30*t35*t64*t8-9.0*t1*t2*t21*t35*t39*
508         t62*t8+4.5*t1*t2*t30*t36*t62*t8-9.0*t1*t2*t36*t39*t59*t8);
509 
510 }
511 
512 static void
pw91x_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)513 pw91x_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
514 {
515     real res[9];
516 
517     pw91x_third_helper(dp->rhoa, dp->grada, res);
518 
519     ds->df1000 += factor*res[0];
520     ds->df0010 += factor*res[1];
521 
522     ds->df2000 += factor*res[2];
523     ds->df1010 += factor*res[3];
524     ds->df0020 += factor*res[4];
525 
526     ds->df3000 += factor*res[5];
527     ds->df2010 += factor*res[6];
528     ds->df1020 += factor*res[7];
529     ds->df0030 += factor*res[8];
530 
531 
532     if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
533        fabs(dp->grada-dp->gradb)>1e-13)
534         pw91x_third_helper(dp->rhob, dp->gradb, res);
535 
536     ds->df0100 += factor*res[0];
537     ds->df0001 += factor*res[1];
538 
539     ds->df0200 += factor*res[2];
540     ds->df0101 += factor*res[3];
541     ds->df0002 += factor*res[4];
542 
543     ds->df0300 += factor*res[5];
544     ds->df0201 += factor*res[6];
545     ds->df0102 += factor*res[7];
546     ds->df0003 += factor*res[8];
547 
548 }
549 
550 static void
pw91x_fourth_helper(real rhoa,real grada,real * res)551 pw91x_fourth_helper(real rhoa, real grada, real *res)
552 {
553     real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
554     real t11, t12, t13, t14, t15, t16, t17, t18;
555     real t19, t20, t21, t22, t23, t24, t25, t26;
556     real t27, t28, t29, t30, t31, t32, t33, t34;
557     real t35, t36, t37, t38, t39, t40, t41, t42;
558     real t43, t44, t45, t46, t47, t48, t49, t50;
559     real t51, t52, t53, t54, t55, t56, t57, t58;
560     real t59, t60, t61, t62, t63, t64, t65, t66;
561     real t67, t68, t69, t70, t71, t72, t73, t74;
562     real t75, t76, t77, t78, t79, t80, t81, t82;
563     real t83, t84, t85, t86, t87, t88, t89, t90;
564     real t91, t92, t93, t94, t95, t96, t97, t98;
565     real t99, t100, t101, t102, t103, t104, t105;
566     real t106, t107, t108, t109, t110, t111, t112;
567     real t113, t114, t115, t116, t117, t118, t119;
568     real t120, t121, t122, t123, t124, t125, t126;
569     real t127, t128;
570 
571     t1 = pow(6.0,0.333333333333333);
572     t2 = 1/pow(3.141592653589793,0.333333333333333);
573     t3 = 1/t1;
574     t4 = 1/pow(3.141592653589793,2.666666666666667);
575     t5 = pow(grada,4.0);
576     t6 = 1/pow(rhoa,5.333333333333333);
577     t7 = 1/pow(3.141592653589793,0.666666666666667);
578     t8 = pow(rhoa,1.333333333333333);
579     t9 = 1/t8;
580     t10 = asinh(3.8978*t3*t7*grada*t9);
581     t11 = 0.098225*t3*t7*grada*t10*t9;
582     t12 = 4.166666666666667E-5*t3*t4*t5*t6+t11+1.0;
583     t13 = 1/t12;
584     t14 = pow(rhoa,0.333333333333333);
585     t15 = 1/pow(6.0,0.666666666666667);
586     t16 = 1/pow(3.141592653589793,1.333333333333333);
587     t17 = pow(grada,2.0);
588     t18 = 1/pow(rhoa,2.666666666666667);
589     t19 = 1/pow(2.718281828459045,25.0*t15*t16*t17*t18);
590     t20 = 0.2743-
591         0.1508*t19;
592     t21 = 0.25*t15*t16*t17*t18*t20+t11+1.0;
593     t22 = 1/pow(rhoa,6.333333333333333);
594     t23 = sqrt(15.19284484*t15*t16*t17*t18+1.0);
595     t24 = 1/t23;
596     t25 = 1/pow(rhoa,3.666666666666667);
597     t26 = -0.510481873333333*t15*t16*t17*t24*t25;
598     t27 = 1/pow(rhoa,2.333333333333333);
599     t28 = -0.130966666666667*t3*t7*grada*t10*t27;
600     t29 = t28+t26-2.222222222222222E-4*t3*t4*t5*t22;
601     t30 = 1/pow(t12,2.0);
602     t31 = t28+t26-0.666666666666667*t15*t16*t17*t20*t25-0.418888888888889*
603         t3*t4*t5*t22*t19;
604     t32 = pow(grada,3.0);
605     t33 = 0.382861405*t15*t16*grada*t24*t18;
606     t34 = 0.098225*t3*t7*t10*t9;
607     t35 = t34+t33+1.6666666666666666E-4*t3*t4*t32*t6;
608     t36 = 0.5*t15*t16*t18*t20*grada+t34+t33+0.314166666666667*
609         t3*t4*t32*t6*t19;
610     t37 = 1/pow(rhoa,0.666666666666667);
611     t38 = pow(t29,2.0);
612     t39 = 1/pow(t12,3.0);
613     t40 = 1/pow(rhoa,7.333333333333333);
614     t41 = 1/pow(t23,3.0);
615     t42 = -1.723482643374637*t3*t4*t5*t41*t40;
616     t43 = 1/pow(rhoa,4.666666666666667);
617     t44 = 2.552409366666666*t15*t16*t17*t24*t43;
618     t45 = 1/pow(rhoa,3.333333333333333);
619     t46 = 0.305588888888889*t3*t7*grada*t10*t45;
620     t47 = t46+t44+t42+0.001407407407407*t3*t4*t5*t40;
621     t48 = 1/pow(3.141592653589793,4.0);
622     t49 = pow(grada,6.0);
623     t50 = 1/pow(rhoa,10.0);
624     t51 = t46+t44+2.444444444444445*t15*t16*t17*t20*t43+t42-
625         4.65432098765432*t48*t49*t50*t19+3.77*t3*t4*t5*t40*t19;
626     t52 = 1.292611982530978*t3*t4*t32*t41*t22;
627     t53 = -1.53144562*t15*t16*grada*t24*t25;
628     t54 = -0.130966666666667*t3*t7*t10*t27;
629     t55 = t54+t53+t52-8.888888888888888E-4*t3*t4*t32*t22;
630     t56 = pow(grada,
631         5.0);
632     t57 = 1/pow(rhoa,9.0);
633     t58 = -1.333333333333333*t15*t16*t20*t25*grada+t54+t53+
634         t52+3.490740740740741*t48*t56*t57*t19-2.513333333333333*t3*
635         t4*t32*t22*t19;
636     t59 = pow(t35,2.0);
637     t60 = -0.969458986898234*t3*t4*t17*t41*t6;
638     t61 = 0.76572281*t15*t16*t24*t18;
639     t62 = t61+t60+5.0E-4*t3*t4*t17*t6;
640     t63 = 1/pow(rhoa,8.0);
641     t64 = t61+t60+0.5*t15*t16*t18*t20-2.618055555555555*t48*
642         t5*t63*t19+1.570833333333333*t3*t4*t17*t6*t19;
643     t65 = 1/pow(rhoa,1.666666666666667);
644     t66 = pow(t29,3.0);
645     t67 = 1/pow(t12,4.0);
646     t68 = 1/pow(t23,5.0);
647     t69 = 1/pow(rhoa,11.0);
648     t70 = -17.45640292348261*t48*t49*t68*t69;
649     t71 = 1/pow(rhoa,8.333333333333334);
650     t72 = 21.25628593495385*t3*t4*t5*t41*t71;
651     t73 = 1/pow(rhoa,5.666666666666667);
652     t74 = -13.49940953925926*t15*t16*t17*t24*t73;
653     t75 = 1/pow(rhoa,4.333333333333333);
654     t76 = -1.01862962962963*t3*t7*grada*t10*t75;
655     t77 = t76+t74+t72-0.010320987654321*t3*t4*t5*t71+t70;
656     t78 = 1/
657         pow(3.141592653589793,5.333333333333333);
658     t79 = pow(grada,8.0);
659     t80 = 1/pow(rhoa,13.66666666666667);
660     t81 = t76+t74-11.40740740740741*t15*t16*t17*t20*t73+t72+
661         t70-310.2880658436214*t15*t78*t79*t80*t19-31.74246913580246*
662         t3*t4*t5*t71*t19+88.43209876543209*t48*t49*t69*t19;
663     t82 = 13.09230219261196*t48*t56*t68*t50;
664     t83 = -13.35699048615344*t3*t4*t32*t41*t40;
665     t84 = 6.295943104444444*t15*t16*grada*t24*t43;
666     t85 = 0.305588888888889*t3*t7*t10*t45;
667     t86 = t85+t84+t83+0.00562962962963*t3*t4*t32*t40+t82;
668     t87 = pow(grada,
669         7.0);
670     t88 = 1/pow(rhoa,12.66666666666667);
671     t89 = 4.88888888888889*t15*t16*t20*t43*grada+t85+t84+
672         t83+t82+232.716049382716*t15*t78*t87*t88*t19-59.34259259259259*
673         t48*t56*t50*t19+18.15185185185185*t3*t4*t32*t40*t19;
674     t90 = -9.81922664445897*t48*t5*t68*t57;
675     t91 = 7.755671895185867*t3*t4*t17*t41*t22;
676     t92 = -2.041927493333334*t15*t16*t24*t25;
677     t93 = t92+t91-0.002666666666667*t3*t4*t17*t22+t90;
678     t94 = 1/pow(rhoa,11.66666666666667);
679     t95 = t92+t91+t90-1.333333333333333*t15*t16*t20*t25-174.537037037037*
680         t15*t78*t49*t94*t19+38.39814814814815*t48*t5*t57*t19-9.215555555555554*
681         t3*t4*t17*t22*t19;
682     t96 = pow(t35,3.0);
683     t97 = 7.364419983344228*t48*t32*t68*t63;
684     t98 = -3.877835947592934*t3*t4*grada*t41*t6;
685     t99 = t98+0.001*t3*t4*grada*t6+t97;
686     t100 = 1/pow(rhoa,10.66666666666667);
687     t101 = 3.769999999999999*t3*t4*grada*t6*t19-23.5625*t48*
688         t32*t63*t19+130.9027777777778*t15*t78*t56*t100*t19+t98+t97;
689     t102 = 1/
690         pow(t12,5.0);
691     t103 = 1/pow(t23,7.0);
692     t104 = 1/pow(rhoa,14.66666666666667);
693     t105 = -1768.082807206625*t15*t78*t79*t103*t104;
694     t106 = 1/pow(rhoa,12.0);
695     t107 = 407.3160682145942*t48*t49*t68*t106;
696     t108 = 1/pow(rhoa,9.333333333333334);
697     t109 = -222.7122571383003*t3*t4*t5*t41*t108;
698     t110 = 1/pow(rhoa,6.666666666666667);
699     t111 = 81.79054014962962*t15*t16*t17*t24*t110;
700     t112 = 4.414061728395062*t3*t7*grada*t10*t6;
701     t113 = 1/pow(3.141592653589793,6.666666666666667);
702     t114 = 1326.062105404968*t15*t78*t87*t103*t80;
703     t115 = -266.2101445831098*t48*t56*t68*t69;
704     t116 = 119.207549500079*t3*t4*t32*t41*t71;
705     t117 = -30.96923364888889*t15*t16*grada*t24*t73;
706     t118 = -1.01862962962963*t3*t7*t10*t75;
707     t119 = -994.5465790537263*t15*t78*t49*t103*t88;
708     t120 = 166.9268529558024*t48*t5*t68*t50;
709     t121 = -56.0131859096757*t3*t4*t17*t41*t40;
710     t122 = 7.487067475555555*t15*t16*t24*t43;
711     t123 = 745.9099342902949*t15*t78*t56*t103*t94;
712     t124 = -98.1922664445897*t48*t32*t68*t57;
713     t125 = 20.68179172049565*t3*t4*grada*t41*t22;
714     t126 = -559.4324507177213*t15*t78*t5*t103*t100;
715     t127 = 51.5509398834096*t48*t17*t68*t63;
716     t128 = -3.877835947592934*t3*t4*t41*t6;
717 
718    /* code */
719     res[0] = 0.5*(-1.5*t1*t13*t2*t31*t8+1.5*t1*t2*t21*t29*
720         t30*t8-2.0*t1*t13*t14*t2*t21);
721     res[1] = 0.5*(1.5*t1*t2*t21*t30*t35*t8-1.5*t1*t13*t2*
722         t36*t8);
723     res[2] = 0.5*(-1.5*t1*t13*t2*t51*t8+1.5*t1*t2*t21*t30*
724         t47*t8-3.0*t1*t2*t21*t38*t39*t8+3.0*t1*t2*t29*t30*t31*t8-0.666666666666667*
725         t1*t13*t2*t21*t37-4.0*t1*t13*t14*t2*t31+4.0*t1*t14*t2*t21*
726         t29*t30);
727     res[3] = 0.5*(-1.5*t1*t13*t2*t58*t8+1.5*t1*t2*t21*t30*
728         t55*t8-3.0*t1*t2*t21*t29*t35*t39*t8+1.5*t1*t2*t29*t30*t36*
729         t8+1.5*t1*t2*t30*t31*t35*t8-2.0*t1*t13*t14*t2*t36+2.0*t1*t14*
730         t2*t21*t30*t35);
731     res[4] = 0.5*(-1.5*t1*t13*t2*t64*t8+1.5*t1*t2*t21*t30*
732         t62*t8-3.0*t1*t2*t21*t39*t59*t8+3.0*t1*t2*t30*t35*t36*t8);
733     res[5] = 0.5*(-1.5*t1*t13*t2*t8*t81+1.5*t1*t2*t21*t30*
734         t77*t8+9.0*t1*t2*t21*t66*t67*t8+4.5*t1*t2*t29*t30*t51*t8-9.0*
735         t1*t2*t21*t29*t39*t47*t8+4.5*t1*t2*t30*t31*t47*t8-9.0*t1*t2*
736         t31*t38*t39*t8+0.444444444444444*t1*t13*t2*t21*t65-6.0*t1*
737         t13*t14*t2*t51+6.0*t1*t14*t2*t21*t30*t47-12.0*t1*t14*t2*t21*
738         t38*t39-2.0*t1*t13*t2*t31*t37+2.0*t1*t2*t21*t29*t30*t37+12.0*
739         t1*t14*t2*t29*t30*t31);
740     res[6] = 0.5*(-1.5*t1*t13*t2*t8*t89+1.5*t1*t2*t21*t30*
741         t8*t86+9.0*t1*t2*t21*t35*t38*t67*t8+3.0*t1*t2*t29*t30*t58*
742         t8-6.0*t1*t2*t21*t29*t39*t55*t8+3.0*t1*t2*t30*t31*t55*t8+1.5*
743         t1*t2*t30*t35*t51*t8-3.0*t1*t2*t21*t35*t39*t47*t8+1.5*t1*t2*
744         t30*t36*t47*t8-3.0*t1*t2*t36*t38*t39*t8-6.0*t1*t2*t29*t31*
745         t35*t39*t8-4.0*t1*t13*t14*t2*t58+4.0*t1*t14*t2*t21*t30*t55-
746         8.0*t1*t14*t2*t21*t29*t35*t39-0.666666666666667*t1*t13*t2*
747         t36*t37+0.666666666666667*t1*t2*t21*t30*t35*t37+4.0*t1*t14*
748         t2*t29*t30*t36+4.0*t1*t14*t2*t30*t31*t35);
749     res[7] = 0.5*(-1.5*t1*t13*t2*t8*t95+1.5*t1*t2*t21*t30*
750         t8*t93+9.0*t1*t2*t21*t29*t59*t67*t8+1.5*t1*t2*t29*t30*t64*
751         t8-3.0*t1*t2*t21*t29*t39*t62*t8+1.5*t1*t2*t30*t31*t62*t8-3.0*
752         t1*t2*t31*t39*t59*t8+3.0*t1*t2*t30*t35*t58*t8-6.0*t1*t2*t21*
753         t35*t39*t55*t8+3.0*t1*t2*t30*t36*t55*t8-6.0*t1*t2*t29*t35*
754         t36*t39*t8-2.0*t1*t13*t14*t2*t64+2.0*t1*t14*t2*t21*t30*t62-
755         4.0*t1*t14*t2*t21*t39*t59+4.0*t1*t14*t2*t30*t35*t36);
756     res[8] = 0.5*(1.5*t1*t2*t21*t30*t8*t99+9.0*t1*t2*t21*
757         t67*t8*t96+4.5*t1*t2*t30*t35*t64*t8-9.0*t1*t2*t21*t35*t39*
758         t62*t8+4.5*t1*t2*t30*t36*t62*t8-9.0*t1*t2*t36*t39*t59*t8-1.5*
759         t1*t101*t13*t2*t8);
760     res[9] = 0.5*(-1.5*t1*t13*t2*t8*(-3447.645176040237*t113*
761         t19*t3*pow(grada,10.0)/pow(rhoa,17.33333333333333)+64.64197530864197*
762         t110*t15*t16*t17*t20+283.6343209876543*t3*t4*t5*t108*t19-1325.447187928669*
763         t48*t49*t106*t19+10136.0768175583*t15*t78*t79*t104*t19+t112+
764         t111+t109+t107+t105)+6.0*t1*t2*t29*t30*t8*t81-8.0*t1*t13*t14*
765         t2*t81-12.0*t1*t2*t21*t29*t39*t77*t8+6.0*t1*t2*t30*t31*t77*
766         t8+36.0*t1*t2*t31*t66*t67*t8+54.0*t1*t2*t21*t38*t47*t67*t8+
767         9.0*t1*t2*t30*t47*t51*t8-18.0*t1*t2*t38*t39*t51*t8-9.0*t1*
768         t2*t21*t39*pow(t47,2.0)*t8-36.0*t1*t2*t29*t31*t39*t47*t8+1.5*
769         t1*(t112+t111+t109+0.086008230452675*t3*t4*t5*t108+t107+t105)*
770         t2*t21*t30*t8-36.0*t1*t102*t2*t21*pow(t29,4.0)*t8+8.0*t1*t14*
771         t2*t21*t30*t77+48.0*t1*t14*t2*t21*t66*t67+1.777777777777778*
772         t1*t13*t2*t31*t65-1.777777777777778*t1*t2*t21*t29*t30*t65-
773         4.0*t1*t13*t2*t37*t51+24.0*t1*t14*t2*t29*t30*t51-48.0*t1*t14*
774         t2*t21*t29*t39*t47+4.0*t1*t2*t21*t30*t37*t47+24.0*t1*t14*t2*
775         t30*t31*t47-8.0*t1*t2*t21*t37*t38*t39-48.0*t1*t14*t2*t31*t38*
776         t39+8.0*t1*t2*t29*t30*t31*t37-0.740740740740741*t1*t13*t18*
777         t2*t21);
778     res[10] = 0.5*(-1.5*t1*t13*t2*t8*(2585.733882030178*t113*
779         t19*t3*pow(grada,9.0)/pow(rhoa,16.33333333333333)-22.81481481481481*
780         t15*t16*t20*t73*grada-6903.909465020575*t15*t78*t87*t80*t19-
781         141.3051851851852*t3*t4*t32*t71*t19+795.1131687242797*t48*
782         t56*t69*t19+t118+t117+t116+t115+t114)+4.5*t1*t2*t29*t30*t8*
783         t89-6.0*t1*t13*t14*t2*t89-9.0*t1*t2*t21*t29*t39*t8*t86+4.5*
784         t1*t2*t30*t31*t8*t86+6.0*t1*t14*t2*t21*t30*t86+1.5*t1*t2*t30*
785         t35*t8*t81-3.0*t1*t2*t21*t35*t39*t77*t8+1.5*t1*t2*t30*t36*
786         t77*t8+9.0*t1*t2*t36*t66*t67*t8+27.0*t1*t2*t21*t38*t55*t67*
787         t8+27.0*t1*t2*t21*t29*t35*t47*t67*t8+27.0*t1*t2*t31*t35*t38*
788         t67*t8-36.0*t1*t102*t2*t21*t35*t66*t8+4.5*t1*t2*t30*t47*t58*
789         t8-9.0*t1*t2*t38*t39*t58*t8+4.5*t1*t2*t30*t51*t55*t8-9.0*t1*
790         t2*t21*t39*t47*t55*t8-18.0*t1*t2*t29*t31*t39*t55*t8-9.0*t1*
791         t2*t29*t35*t39*t51*t8-9.0*t1*t2*t29*t36*t39*t47*t8-9.0*t1*
792         t2*t31*t35*t39*t47*t8+1.5*t1*(t118+t117+t116-0.041283950617284*
793         t3*t4*t32*t71+t115+t114)*t2*t21*t30*t8+36.0*t1*t14*t2*t21*
794         t35*t38*t67+0.444444444444444*t1*t13*t2*t36*t65-0.444444444444444*
795         t1*t2*t21*t30*t35*t65-2.0*t1*t13*t2*t37*t58+12.0*t1*t14*t2*
796         t29*t30*t58-24.0*t1*t14*t2*t21*t29*t39*t55+2.0*t1*t2*t21*t30*
797         t37*t55+12.0*t1*t14*t2*t30*t31*t55+6.0*t1*t14*t2*t30*t35*t51-
798         12.0*t1*t14*t2*t21*t35*t39*t47+6.0*t1*t14*t2*t30*t36*t47-12.0*
799         t1*t14*t2*t36*t38*t39-4.0*t1*t2*t21*t29*t35*t37*t39-24.0*t1*
800         t14*t2*t29*t31*t35*t39+2.0*t1*t2*t29*t30*t36*t37+2.0*t1*t2*
801         t30*t31*t35*t37);
802     res[11] = 0.5*(-1.5*t1*t13*t2*t8*(-1939.300411522634*
803         t113*t19*t3*t79/pow(rhoa,15.33333333333333)+4.88888888888889*
804         t15*t16*t20*t43+4596.141975308642*t15*t78*t49*t88*t19-447.9783950617283*
805         t48*t5*t50*t19+60.59925925925925*t3*t4*t17*t40*t19+t122+t121+
806         t120+t119)+3.0*t1*t2*t29*t30*t8*t95-4.0*t1*t13*t14*t2*t95-
807         6.0*t1*t2*t21*t29*t39*t8*t93+3.0*t1*t2*t30*t31*t8*t93+4.0*
808         t1*t14*t2*t21*t30*t93+3.0*t1*t2*t30*t35*t8*t89-6.0*t1*t2*t21*
809         t35*t39*t8*t86+3.0*t1*t2*t30*t36*t8*t86+9.0*t1*t2*t21*t38*
810         t62*t67*t8+9.0*t1*t2*t21*t47*t59*t67*t8+18.0*t1*t2*t29*t31*
811         t59*t67*t8+36.0*t1*t2*t21*t29*t35*t55*t67*t8+18.0*t1*t2*t35*
812         t36*t38*t67*t8+1.5*t1*t2*t30*t47*t64*t8-3.0*t1*t2*t38*t39*
813         t64*t8+1.5*t1*t2*t30*t51*t62*t8-3.0*t1*t2*t21*t39*t47*t62*
814         t8-6.0*t1*t2*t29*t31*t39*t62*t8-3.0*t1*t2*t39*t51*t59*t8-36.0*
815         t1*t102*t2*t21*t38*t59*t8+6.0*t1*t2*t30*t55*t58*t8-12.0*t1*
816         t2*t29*t35*t39*t58*t8-6.0*t1*t2*t21*t39*pow(t55,2.0)*t8-12.0*
817         t1*t2*t29*t36*t39*t55*t8-12.0*t1*t2*t31*t35*t39*t55*t8-6.0*
818         t1*t2*t35*t36*t39*t47*t8+1.5*t1*(t122+t121+0.016888888888889*
819         t3*t4*t17*t40+t120+t119)*t2*t21*t30*t8+24.0*t1*t14*t2*t21*
820         t29*t59*t67-0.666666666666667*t1*t13*t2*t37*t64+4.0*t1*t14*
821         t2*t29*t30*t64-8.0*t1*t14*t2*t21*t29*t39*t62+0.666666666666667*
822         t1*t2*t21*t30*t37*t62+4.0*t1*t14*t2*t30*t31*t62-1.333333333333333*
823         t1*t2*t21*t37*t39*t59-8.0*t1*t14*t2*t31*t39*t59+8.0*t1*t14*
824         t2*t30*t35*t58-16.0*t1*t14*t2*t21*t35*t39*t55+8.0*t1*t14*t2*
825         t30*t36*t55-16.0*t1*t14*t2*t29*t35*t36*t39+1.333333333333333*
826         t1*t2*t30*t35*t36*t37);
827     res[12] = 0.5*(-1.5*t1*t13*t2*t8*(1454.475308641976*t113*
828         t19*t3*t87/pow(rhoa,14.33333333333333)-2967.12962962963*t15*
829         t78*t56*t94*t19+230.3888888888889*t48*t32*t57*t19-20.10666666666666*
830         t3*t4*grada*t22*t19+t125+t124+t123)-3.0*t1*t2*t21*t29*t39*
831         t8*t99+1.5*t1*t2*t30*t31*t8*t99+2.0*t1*t14*t2*t21*t30*t99+
832         9.0*t1*t2*t31*t67*t8*t96-36.0*t1*t102*t2*t21*t29*t8*t96+12.0*
833         t1*t14*t2*t21*t67*t96+4.5*t1*t2*t30*t35*t8*t95-9.0*t1*t2*t21*
834         t35*t39*t8*t93+4.5*t1*t2*t30*t36*t8*t93+27.0*t1*t2*t21*t29*
835         t35*t62*t67*t8+27.0*t1*t2*t21*t55*t59*t67*t8+27.0*t1*t2*t29*
836         t36*t59*t67*t8+4.5*t1*t2*t30*t55*t64*t8-9.0*t1*t2*t29*t35*
837         t39*t64*t8+4.5*t1*t2*t30*t58*t62*t8-9.0*t1*t2*t21*t39*t55*
838         t62*t8-9.0*t1*t2*t29*t36*t39*t62*t8-9.0*t1*t2*t31*t35*t39*
839         t62*t8-9.0*t1*t2*t39*t58*t59*t8-18.0*t1*t2*t35*t36*t39*t55*
840         t8+1.5*t1*t101*t2*t29*t30*t8+1.5*t1*(t125-0.005333333333333*
841         t3*t4*grada*t22+t124+t123)*t2*t21*t30*t8+6.0*t1*t14*t2*t30*
842         t35*t64-12.0*t1*t14*t2*t21*t35*t39*t62+6.0*t1*t14*t2*t30*t36*
843         t62-12.0*t1*t14*t2*t36*t39*t59-2.0*t1*t101*t13*t14*t2);
844     res[13] = 0.5*(-1.5*t1*t13*t2*t8*(-1090.856481481482*
845         t113*t19*t3*t49/pow(rhoa,13.33333333333333)-102.1041666666666*
846         t48*t17*t63*t19+3.769999999999999*t3*t4*t6*t19+1832.638888888889*
847         t15*t78*t5*t100*t19+t128+t127+t126)-12.0*t1*t2*t21*t35*t39*
848         t8*t99+6.0*t1*t2*t30*t36*t8*t99+36.0*t1*t2*t36*t67*t8*t96+
849         54.0*t1*t2*t21*t59*t62*t67*t8+9.0*t1*t2*t30*t62*t64*t8-18.0*
850         t1*t2*t39*t59*t64*t8-9.0*t1*t2*t21*t39*pow(t62,2.0)*t8-36.0*
851         t1*t2*t35*t36*t39*t62*t8-36.0*t1*t102*t2*t21*pow(t35,4.0)*
852         t8+6.0*t1*t101*t2*t30*t35*t8+1.5*t1*(t128+0.001*t3*t4*t6+t127+
853         t126)*t2*t21*t30*t8);
854 
855 }
856 
857 static void
pw91x_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)858 pw91x_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
859 {
860     real res[14];
861 
862     pw91x_fourth_helper(dp->rhoa, dp->grada, res);
863 
864     ds->df1000 += factor*res[0];
865     ds->df0010 += factor*res[1];
866 
867     ds->df2000 += factor*res[2];
868     ds->df1010 += factor*res[3];
869     ds->df0020 += factor*res[4];
870 
871     ds->df3000 += factor*res[5];
872     ds->df2010 += factor*res[6];
873     ds->df1020 += factor*res[7];
874     ds->df0030 += factor*res[8];
875 
876     ds->df4000 += factor*res[9];
877     ds->df3010 += factor*res[10];
878     ds->df2020 += factor*res[11];
879     ds->df1030 += factor*res[12];
880     ds->df0040 += factor*res[13];
881 
882 
883     if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
884        fabs(dp->grada-dp->gradb)>1e-13)
885         pw91x_fourth_helper(dp->rhob, dp->gradb, res);
886 
887     ds->df0100 += factor*res[0];
888     ds->df0001 += factor*res[1];
889 
890     ds->df0200 += factor*res[2];
891     ds->df0101 += factor*res[3];
892     ds->df0002 += factor*res[4];
893 
894     ds->df0300 += factor*res[5];
895     ds->df0201 += factor*res[6];
896     ds->df0102 += factor*res[7];
897     ds->df0003 += factor*res[8];
898 
899     ds->df0400 += factor*res[9];
900     ds->df0301 += factor*res[10];
901     ds->df0202 += factor*res[11];
902     ds->df0103 += factor*res[12];
903     ds->df0004 += factor*res[13];
904 
905 }
906