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