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-pw86x.c:
24
25 Implemented by Dave Wilson (david.wilson@latrobe.edu.au), Jun 2005
26
27 Automatically generated code implementing PW86X 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 pw86xFunctional;" to 'functionals.h'
34 2. add "&pw86xFunctional," to 'functionals.c'
35 3. add "fun-pw86x.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 Fa: (1+1.296*sa^2+14.0*sa^4+0.2*sa^6)^(1/15);
51 Fb: (1+1.296*sb^2+14.0*sb^4+0.2*sb^6)^(1/15);
52
53 Exa: 2*rhoa*exa*Fa;
54 Exb: 2*rhob*exb*Fb;
55
56 K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb)/2;
57
58 ------ cut here -------
59 */
60
61
62 /* strictly conform to XOPEN ANSI C standard */
63 #if !defined(SYS_DEC)
64 /* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
65 * is not specified. */
66 #define _XOPEN_SOURCE 500
67 #define _XOPEN_SOURCE_EXTENDED 1
68 #endif
69 #include <math.h>
70 #include <stddef.h>
71 #include "general.h"
72
73 #define __CVERSION__
74
75 #include "functionals.h"
76
77 /* INTERFACE PART */
pw86x_isgga(void)78 static integer pw86x_isgga(void) { return 1; } /* FIXME: detect! */
79 static integer pw86x_read(const char *conf_line);
80 static real pw86x_energy(const FunDensProp* dp);
81 static void pw86x_first(FunFirstFuncDrv *ds, real factor,
82 const FunDensProp* dp);
83 static void pw86x_second(FunSecondFuncDrv *ds, real factor,
84 const FunDensProp* dp);
85 static void pw86x_third(FunThirdFuncDrv *ds, real factor,
86 const FunDensProp* dp);
87 static void pw86x_fourth(FunFourthFuncDrv *ds, real factor,
88 const FunDensProp* dp);
89
90 Functional PW86xFunctional = {
91 "PW86x", /* name */
92 pw86x_isgga, /* gga-corrected */
93 1,
94 pw86x_read,
95 NULL,
96 pw86x_energy,
97 pw86x_first,
98 pw86x_second,
99 pw86x_third,
100 pw86x_fourth
101 };
102
103 /* IMPLEMENTATION PART */
104 static integer
pw86x_read(const char * conf_line)105 pw86x_read(const char *conf_line)
106 {
107 fun_set_hf_weight(0);
108 return 1;
109 }
110
111 static real
pw86x_energy(const FunDensProp * dp)112 pw86x_energy(const FunDensProp *dp)
113 {
114 real res;
115 real rhoa = dp->rhoa, rhob = dp->rhob;
116 real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
117
118 real t1, t2, t3, t4, t5, t6, t7;
119
120 t1 = pow(6.0,0.333333333333333);
121 t2 = 1/pow(3.141592653589793,0.333333333333333);
122 t3 = 1/pow(3.141592653589793,4.0);
123 t4 = 1/t1;
124 t5 = 1/pow(3.141592653589793,2.666666666666667);
125 t6 = 1/pow(6.0,0.666666666666667);
126 t7 = 1/pow(3.141592653589793,1.333333333333333);
127
128 /* code */
129 res = 0.5*(-1.5*t1*t2*pow(rhob,1.333333333333333)*pow(8.680555555555556E-5*
130 t3*pow(gradb,6.0)/pow(rhob,8.0)+0.145833333333333*t4*t5*pow(gradb,
131 4.0)/pow(rhob,5.333333333333333)+0.324*t6*t7*pow(gradb,2.0)/
132 pow(rhob,2.666666666666667)+1.0,0.066666666666667)-1.5*t1*
133 t2*pow(rhoa,1.333333333333333)*pow(8.680555555555556E-5*t3*
134 pow(grada,6.0)/pow(rhoa,8.0)+0.145833333333333*t4*t5*pow(grada,
135 4.0)/pow(rhoa,5.333333333333333)+0.324*t6*t7*pow(grada,2.0)/
136 pow(rhoa,2.666666666666667)+1.0,0.066666666666667));
137
138 return res;
139 }
140
141 static void
pw86x_first_helper(real rhoa,real grada,real * res)142 pw86x_first_helper(real rhoa, real grada, real *res)
143 {
144 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
145 real t11, t12, t13, t14, t15, t16;
146
147 t1 = pow(6.0,0.333333333333333);
148 t2 = 1/pow(3.141592653589793,0.333333333333333);
149 t3 = 1/pow(3.141592653589793,4.0);
150 t4 = pow(grada,6.0);
151 t5 = 1/pow(rhoa,8.0);
152 t6 = 1/t1;
153 t7 = 1/pow(3.141592653589793,2.666666666666667);
154 t8 = pow(grada,4.0);
155 t9 = 1/pow(rhoa,5.333333333333333);
156 t10 = 1/pow(6.0,0.666666666666667);
157 t11 = 1/pow(3.141592653589793,1.333333333333333);
158 t12 = pow(grada,2.0);
159 t13 = 1/pow(rhoa,2.666666666666667);
160 t14 = 0.145833333333333*t6*t7*t8*t9+8.680555555555556E-5*
161 t3*t4*t5+0.324*t10*t11*t12*t13+1.0;
162 t15 = 1/pow(t14,0.933333333333333);
163 t16 = pow(rhoa,1.333333333333333);
164
165 /* code */
166 res[0] = 0.5*(-0.1*t1*t15*t16*t2*(-6.944444444444445E-4*
167 t3*t4/pow(rhoa,9.0)-0.777777777777778*t6*t7*t8/pow(rhoa,6.333333333333333)-
168 0.864*t10*t11*t12/pow(rhoa,3.666666666666667))-2.0*t1*pow(t14,
169 0.066666666666667)*t2*pow(rhoa,0.333333333333333));
170 res[1] = -0.05*t1*t15*t16*t2*(5.208333333333333E-4*t3*
171 t5*pow(grada,5.0)+0.583333333333333*t6*t7*t9*pow(grada,3.0)+
172 0.648*t10*t11*grada*t13);
173 }
174
175 static void
pw86x_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)176 pw86x_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
177 {
178 real res[2];
179
180 pw86x_first_helper(dp->rhoa, dp->grada, res);
181 /* Final assignment */
182 ds->df1000 += factor*res[0];
183 ds->df0010 += factor*res[1];
184
185
186 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
187 fabs(dp->grada-dp->gradb)>1e-13)
188 pw86x_first_helper(dp->rhob, dp->gradb, res);
189 ds->df0100 += factor*res[0];
190 ds->df0001 += factor*res[1];
191
192 }
193
194 static void
pw86x_second_helper(real rhoa,real grada,real * res)195 pw86x_second_helper(real rhoa, real grada, real *res)
196 {
197 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
198 real t11, t12, t13, t14, t15, t16, t17, t18;
199 real t19, t20, t21, t22, t23, t24, t25, t26;
200
201 t1 = pow(6.0,0.333333333333333);
202 t2 = 1/pow(3.141592653589793,0.333333333333333);
203 t3 = 1/pow(3.141592653589793,4.0);
204 t4 = pow(grada,6.0);
205 t5 = 1/pow(rhoa,8.0);
206 t6 = 1/t1;
207 t7 = 1/pow(3.141592653589793,2.666666666666667);
208 t8 = pow(grada,4.0);
209 t9 = 1/pow(rhoa,5.333333333333333);
210 t10 = 1/pow(6.0,0.666666666666667);
211 t11 = 1/pow(3.141592653589793,1.333333333333333);
212 t12 = pow(grada,2.0);
213 t13 = 1/pow(rhoa,2.666666666666667);
214 t14 = 0.145833333333333*t6*t7*t8*t9+8.680555555555556E-5*
215 t3*t4*t5+0.324*t10*t11*t12*t13+1.0;
216 t15 = pow(t14,0.066666666666667);
217 t16 = pow(rhoa,0.333333333333333);
218 t17 = 1/pow(rhoa,9.0);
219 t18 = 1/pow(rhoa,6.333333333333333);
220 t19 = 1/pow(rhoa,3.666666666666667);
221 t20 = -0.864*t10*t11*t12*t19-0.777777777777778*t6*t7*
222 t8*t18-6.944444444444445E-4*t3*t4*t17;
223 t21 = 1/pow(t14,0.933333333333333);
224 t22 = pow(rhoa,1.333333333333333);
225 t23 = pow(grada,5.0);
226 t24 = pow(grada,3.0);
227 t25 = 0.648*t10*t11*grada*t13+0.583333333333333*t6*t7*
228 t24*t9+5.208333333333333E-4*t3*t23*t5;
229 t26 = 1/pow(t14,1.933333333333333);
230
231 /* code */
232 res[0] = 0.5*(-0.1*t1*t2*t20*t21*t22-2.0*t1*t15*t16*t2);
233 res[1] = -0.05*t1*t2*t21*t22*t25;
234 res[2] = 0.5*(-0.1*t1*t2*t21*t22*(0.00625*t3*t4/pow(rhoa,
235 10.0)+4.925925925925925*t6*t7*t8/pow(rhoa,7.333333333333333)+
236 3.168*t10*t11*t12/pow(rhoa,4.666666666666667))-0.666666666666667*
237 t1*t15*t2/pow(rhoa,0.666666666666667)+0.093333333333333*t1*
238 t2*pow(t20,2.0)*t22*t26-0.266666666666667*t1*t16*t2*t20*t21);
239 res[3] = 0.5*(0.093333333333333*t1*t2*t20*t22*t25*t26-
240 0.133333333333333*t1*t16*t2*t21*t25-0.1*t1*(-1.728*t10*t11*
241 grada*t19-3.111111111111111*t6*t7*t24*t18-0.004166666666667*
242 t3*t23*t17)*t2*t21*t22);
243 res[4] = 0.046666666666667*t1*t2*t22*pow(t25,2.0)*t26-
244 0.05*t1*(0.648*t10*t11*t13+1.75*t6*t7*t12*t9+0.002604166666667*
245 t3*t8*t5)*t2*t21*t22;
246
247 }
248
249 static void
pw86x_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)250 pw86x_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
251 {
252 real res[5];
253
254 pw86x_second_helper(dp->rhoa, dp->grada, res);
255
256 ds->df1000 += factor*res[0];
257 ds->df0010 += factor*res[1];
258
259 ds->df2000 += factor*res[2];
260 ds->df1010 += factor*res[3];
261 ds->df0020 += factor*res[4];
262
263
264 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
265 fabs(dp->grada-dp->gradb)>1e-13)
266 pw86x_second_helper(dp->rhob, dp->gradb, res);
267 ds->df0100 += factor*res[0];
268 ds->df0001 += factor*res[1];
269
270 ds->df0200 += factor*res[2];
271 ds->df0101 += factor*res[3];
272 ds->df0002 += factor*res[4];
273
274 }
275
276 static void
pw86x_third_helper(real rhoa,real grada,real * res)277 pw86x_third_helper(real rhoa, real grada, real *res)
278 {
279 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
280 real t11, t12, t13, t14, t15, t16, t17, t18;
281 real t19, t20, t21, t22, t23, t24, t25, t26;
282 real t27, t28, t29, t30, t31, t32, t33, t34;
283 real t35, t36;
284
285 t1 = pow(6.0,0.333333333333333);
286 t2 = 1/pow(3.141592653589793,0.333333333333333);
287 t3 = 1/pow(3.141592653589793,4.0);
288 t4 = pow(grada,6.0);
289 t5 = 1/pow(rhoa,8.0);
290 t6 = 1/t1;
291 t7 = 1/pow(3.141592653589793,2.666666666666667);
292 t8 = pow(grada,4.0);
293 t9 = 1/pow(rhoa,5.333333333333333);
294 t10 = 1/pow(6.0,0.666666666666667);
295 t11 = 1/pow(3.141592653589793,1.333333333333333);
296 t12 = pow(grada,2.0);
297 t13 = 1/pow(rhoa,2.666666666666667);
298 t14 = 0.145833333333333*t6*t7*t8*t9+8.680555555555556E-5*
299 t3*t4*t5+0.324*t10*t11*t12*t13+1.0;
300 t15 = pow(t14,0.066666666666667);
301 t16 = pow(rhoa,0.333333333333333);
302 t17 = 1/pow(rhoa,9.0);
303 t18 = 1/pow(rhoa,6.333333333333333);
304 t19 = 1/pow(rhoa,3.666666666666667);
305 t20 = -0.864*t10*t11*t12*t19-0.777777777777778*t6*t7*
306 t8*t18-6.944444444444445E-4*t3*t4*t17;
307 t21 = 1/pow(t14,0.933333333333333);
308 t22 = pow(rhoa,1.333333333333333);
309 t23 = pow(grada,5.0);
310 t24 = pow(grada,3.0);
311 t25 = 0.648*t10*t11*grada*t13+0.583333333333333*t6*t7*
312 t24*t9+5.208333333333333E-4*t3*t23*t5;
313 t26 = 1/pow(rhoa,0.666666666666667);
314 t27 = pow(t20,2.0);
315 t28 = 1/pow(t14,1.933333333333333);
316 t29 = 1/pow(rhoa,10.0);
317 t30 = 1/pow(rhoa,7.333333333333333);
318 t31 = 1/pow(rhoa,4.666666666666667);
319 t32 = 3.168*t10*t11*t12*t31+4.925925925925925*t6*t7*t8*
320 t30+0.00625*t3*t4*t29;
321 t33 = -1.728*t10*t11*grada*t19-3.111111111111111*t6*t7*
322 t24*t18-0.004166666666667*t3*t23*t17;
323 t34 = pow(t25,2.0);
324 t35 = 0.648*t10*t11*t13+1.75*t6*t7*t12*t9+0.002604166666667*
325 t3*t8*t5;
326 t36 = 1/pow(t14,2.933333333333333);
327
328 /* code */
329 res[0] = 0.5*(-0.1*t1*t2*t20*t21*t22-2.0*t1*t15*t16*t2);
330 res[1] = -0.05*t1*t2*t21*t22*t25;
331 res[2] = 0.5*(-0.1*t1*t2*t21*t22*t32+0.093333333333333*
332 t1*t2*t22*t27*t28-0.666666666666667*t1*t15*t2*t26-0.266666666666667*
333 t1*t16*t2*t20*t21);
334 res[3] = 0.5*(-0.1*t1*t2*t21*t22*t33+0.093333333333333*
335 t1*t2*t20*t22*t25*t28-0.133333333333333*t1*t16*t2*t21*t25);
336 res[4] = 0.046666666666667*t1*t2*t22*t28*t34-0.05*t1*
337 t2*t21*t22*t35;
338 res[5] = 0.5*(-0.1*t1*t2*t21*t22*(-0.0625*t3*t4/pow(rhoa,
339 11.0)-36.12345679012345*t6*t7*t8/pow(rhoa,8.333333333333334)-
340 14.784*t10*t11*t12/pow(rhoa,5.666666666666667))+0.444444444444444*
341 t1*t15*t2/pow(rhoa,1.666666666666667)-0.180444444444444*t1*
342 t2*pow(t20,3.0)*t22*t36+0.28*t1*t2*t20*t22*t28*t32-0.4*t1*
343 t16*t2*t21*t32+0.373333333333333*t1*t16*t2*t27*t28-0.133333333333333*
344 t1*t2*t20*t21*t26);
345 res[6] = 0.5*(-0.180444444444444*t1*t2*t22*t25*t27*t36+
346 0.186666666666667*t1*t2*t20*t22*t28*t33-0.266666666666667*
347 t1*t16*t2*t21*t33+0.093333333333333*t1*t2*t22*t25*t28*t32-
348 0.1*t1*t2*t21*t22*(6.336*t10*t11*grada*t31+19.7037037037037*
349 t6*t7*t24*t30+0.0375*t3*t23*t29)+0.248888888888889*t1*t16*
350 t2*t20*t25*t28-0.044444444444444*t1*t2*t21*t25*t26);
351 res[7] = 0.5*(-0.180444444444444*t1*t2*t20*t22*t34*t36+
352 0.093333333333333*t1*t2*t20*t22*t28*t35-0.133333333333333*
353 t1*t16*t2*t21*t35+0.124444444444444*t1*t16*t2*t28*t34+0.186666666666667*
354 t1*t2*t22*t25*t28*t33-0.1*t1*(-1.728*t10*t11*t19-9.333333333333332*
355 t6*t7*t12*t18-0.020833333333333*t3*t8*t17)*t2*t21*t22);
356 res[8] = -0.05*t1*t2*t21*t22*(3.5*t6*t7*grada*t9+0.010416666666667*
357 t3*t24*t5)-0.090222222222222*t1*t2*t22*pow(t25,3.0)*t36+0.14*
358 t1*t2*t22*t25*t28*t35;
359
360 }
361
362 static void
pw86x_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)363 pw86x_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
364 {
365 real res[9];
366
367 pw86x_third_helper(dp->rhoa, dp->grada, res);
368
369 ds->df1000 += factor*res[0];
370 ds->df0010 += factor*res[1];
371
372 ds->df2000 += factor*res[2];
373 ds->df1010 += factor*res[3];
374 ds->df0020 += factor*res[4];
375
376 ds->df3000 += factor*res[5];
377 ds->df2010 += factor*res[6];
378 ds->df1020 += factor*res[7];
379 ds->df0030 += factor*res[8];
380
381
382 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
383 fabs(dp->grada-dp->gradb)>1e-13)
384 pw86x_third_helper(dp->rhob, dp->gradb, res);
385
386 ds->df0100 += factor*res[0];
387 ds->df0001 += factor*res[1];
388
389 ds->df0200 += factor*res[2];
390 ds->df0101 += factor*res[3];
391 ds->df0002 += factor*res[4];
392
393 ds->df0300 += factor*res[5];
394 ds->df0201 += factor*res[6];
395 ds->df0102 += factor*res[7];
396 ds->df0003 += factor*res[8];
397
398 }
399
400 static void
pw86x_fourth_helper(real rhoa,real grada,real * res)401 pw86x_fourth_helper(real rhoa, real grada, real *res)
402 {
403 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
404 real t11, t12, t13, t14, t15, t16, t17, t18;
405 real t19, t20, t21, t22, t23, t24, t25, t26;
406 real t27, t28, t29, t30, t31, t32, t33, t34;
407 real t35, t36, t37, t38, t39, t40, t41, t42;
408 real t43, t44, t45, t46, t47;
409
410 t1 = pow(6.0,0.333333333333333);
411 t2 = 1/pow(3.141592653589793,0.333333333333333);
412 t3 = 1/pow(3.141592653589793,4.0);
413 t4 = pow(grada,6.0);
414 t5 = 1/pow(rhoa,8.0);
415 t6 = 1/t1;
416 t7 = 1/pow(3.141592653589793,2.666666666666667);
417 t8 = pow(grada,4.0);
418 t9 = 1/pow(rhoa,5.333333333333333);
419 t10 = 1/pow(6.0,0.666666666666667);
420 t11 = 1/pow(3.141592653589793,1.333333333333333);
421 t12 = pow(grada,2.0);
422 t13 = 1/pow(rhoa,2.666666666666667);
423 t14 = 0.145833333333333*t6*t7*t8*t9+8.680555555555556E-5*
424 t3*t4*t5+0.324*t10*t11*t12*t13+1.0;
425 t15 = pow(t14,0.066666666666667);
426 t16 = pow(rhoa,0.333333333333333);
427 t17 = 1/pow(rhoa,9.0);
428 t18 = 1/pow(rhoa,6.333333333333333);
429 t19 = 1/pow(rhoa,3.666666666666667);
430 t20 = -0.864*t10*t11*t12*t19-0.777777777777778*t6*t7*
431 t8*t18-6.944444444444445E-4*t3*t4*t17;
432 t21 = 1/pow(t14,0.933333333333333);
433 t22 = pow(rhoa,1.333333333333333);
434 t23 = pow(grada,5.0);
435 t24 = pow(grada,3.0);
436 t25 = 0.648*t10*t11*grada*t13+0.583333333333333*t6*t7*
437 t24*t9+5.208333333333333E-4*t3*t23*t5;
438 t26 = 1/pow(rhoa,0.666666666666667);
439 t27 = pow(t20,2.0);
440 t28 = 1/pow(t14,1.933333333333333);
441 t29 = 1/pow(rhoa,10.0);
442 t30 = 1/pow(rhoa,7.333333333333333);
443 t31 = 1/pow(rhoa,4.666666666666667);
444 t32 = 3.168*t10*t11*t12*t31+4.925925925925925*t6*t7*t8*
445 t30+0.00625*t3*t4*t29;
446 t33 = -1.728*t10*t11*grada*t19-3.111111111111111*t6*t7*
447 t24*t18-0.004166666666667*t3*t23*t17;
448 t34 = pow(t25,2.0);
449 t35 = 0.648*t10*t11*t13+1.75*t6*t7*t12*t9+0.002604166666667*
450 t3*t8*t5;
451 t36 = 1/pow(rhoa,1.666666666666667);
452 t37 = pow(t20,3.0);
453 t38 = 1/pow(t14,2.933333333333333);
454 t39 = 1/pow(rhoa,11.0);
455 t40 = 1/pow(rhoa,8.333333333333334);
456 t41 = 1/pow(rhoa,5.666666666666667);
457 t42 = -14.784*t10*t11*t12*t41-36.12345679012345*t6*t7*
458 t8*t40-0.0625*t3*t4*t39;
459 t43 = 6.336*t10*t11*grada*t31+19.7037037037037*t6*t7*
460 t24*t30+0.0375*t3*t23*t29;
461 t44 = -1.728*t10*t11*t19-9.333333333333332*t6*t7*t12*
462 t18-0.020833333333333*t3*t8*t17;
463 t45 = pow(t25,3.0);
464 t46 = 3.5*t6*t7*grada*t9+0.010416666666667*t3*t24*t5;
465 t47 = 1/
466 pow(t14,3.933333333333333);
467
468 /* code */
469 res[0] = 0.5*(-0.1*t1*t2*t20*t21*t22-2.0*t1*t15*t16*t2);
470 res[1] = -0.05*t1*t2*t21*t22*t25;
471 res[2] = 0.5*(-0.1*t1*t2*t21*t22*t32+0.093333333333333*
472 t1*t2*t22*t27*t28-0.666666666666667*t1*t15*t2*t26-0.266666666666667*
473 t1*t16*t2*t20*t21);
474 res[3] = 0.5*(-0.1*t1*t2*t21*t22*t33+0.093333333333333*
475 t1*t2*t20*t22*t25*t28-0.133333333333333*t1*t16*t2*t21*t25);
476 res[4] = 0.046666666666667*t1*t2*t22*t28*t34-0.05*t1*
477 t2*t21*t22*t35;
478 res[5] = 0.5*(-0.1*t1*t2*t21*t22*t42-0.180444444444444*
479 t1*t2*t22*t37*t38+0.444444444444444*t1*t15*t2*t36+0.28*t1*
480 t2*t20*t22*t28*t32-0.4*t1*t16*t2*t21*t32+0.373333333333333*
481 t1*t16*t2*t27*t28-0.133333333333333*t1*t2*t20*t21*t26);
482 res[6] = 0.5*(-0.1*t1*t2*t21*t22*t43-0.180444444444444*
483 t1*t2*t22*t25*t27*t38+0.186666666666667*t1*t2*t20*t22*t28*
484 t33-0.266666666666667*t1*t16*t2*t21*t33+0.093333333333333*
485 t1*t2*t22*t25*t28*t32+0.248888888888889*t1*t16*t2*t20*t25*
486 t28-0.044444444444444*t1*t2*t21*t25*t26);
487 res[7] = 0.5*(-0.1*t1*t2*t21*t22*t44-0.180444444444444*
488 t1*t2*t20*t22*t34*t38+0.093333333333333*t1*t2*t20*t22*t28*
489 t35-0.133333333333333*t1*t16*t2*t21*t35+0.124444444444444*
490 t1*t16*t2*t28*t34+0.186666666666667*t1*t2*t22*t25*t28*t33);
491 res[8] = -0.05*t1*t2*t21*t22*t46-0.090222222222222*t1*
492 t2*t22*t38*t45+0.14*t1*t2*t22*t25*t28*t35;
493 res[9] = 0.5*(-0.1*t1*t2*t21*t22*(0.6875*t3*t4/pow(rhoa,
494 12.0)+301.0288065843621*t6*t7*t8/pow(rhoa,9.333333333333334)+
495 83.776*t10*t11*t12/pow(rhoa,6.666666666666667))+0.529303703703704*
496 t1*t2*pow(t20,4.0)*t22*t47+0.373333333333333*t1*t2*t20*t22*
497 t28*t42-0.533333333333333*t1*t16*t2*t21*t42-0.96237037037037*
498 t1*t16*t2*t37*t38-1.082666666666667*t1*t2*t22*t27*t32*t38+
499 0.118518518518519*t1*t2*t20*t21*t36+0.28*t1*t2*t22*t28*pow(t32,
500 2.0)+1.493333333333333*t1*t16*t2*t20*t28*t32-0.266666666666667*
501 t1*t2*t21*t26*t32+0.248888888888889*t1*t2*t26*t27*t28-0.740740740740741*
502 t1*t13*t15*t2);
503 res[10] = 0.5*(0.529303703703704*t1*t2*t22*t25*t37*t47+
504 0.28*t1*t2*t20*t22*t28*t43-0.4*t1*t16*t2*t21*t43+0.093333333333333*
505 t1*t2*t22*t25*t28*t42-0.1*t1*t2*t21*t22*(-29.568*t10*t11*grada*
506 t41-144.4938271604938*t6*t7*t24*t40-0.375*t3*t23*t39)-0.541333333333333*
507 t1*t2*t22*t27*t33*t38-0.541333333333333*t1*t2*t20*t22*t25*
508 t32*t38-0.721777777777778*t1*t16*t2*t25*t27*t38+0.02962962962963*
509 t1*t2*t21*t25*t36+0.28*t1*t2*t22*t28*t32*t33+0.746666666666667*
510 t1*t16*t2*t20*t28*t33-0.133333333333333*t1*t2*t21*t26*t33+
511 0.373333333333333*t1*t16*t2*t25*t28*t32+0.124444444444444*
512 t1*t2*t20*t25*t26*t28);
513 res[11] = 0.5*(0.529303703703704*t1*t2*t22*t27*t34*t47+
514 0.186666666666667*t1*t2*t20*t22*t28*t44-0.266666666666667*
515 t1*t16*t2*t21*t44+0.186666666666667*t1*t2*t22*t25*t28*t43-
516 0.180444444444444*t1*t2*t22*t27*t35*t38-0.180444444444444*
517 t1*t2*t22*t32*t34*t38-0.481185185185185*t1*t16*t2*t20*t34*
518 t38-0.721777777777778*t1*t2*t20*t22*t25*t33*t38+0.093333333333333*
519 t1*t2*t22*t28*t32*t35+0.248888888888889*t1*t16*t2*t20*t28*
520 t35-0.044444444444444*t1*t2*t21*t26*t35+0.041481481481481*
521 t1*t2*t26*t28*t34+0.186666666666667*t1*t2*t22*t28*pow(t33,
522 2.0)+0.497777777777778*t1*t16*t2*t25*t28*t33-0.1*t1*t2*t21*
523 t22*(6.336*t10*t11*t31+59.1111111111111*t6*t7*t12*t30+0.1875*
524 t3*t8*t29));
525 res[12] = 0.5*(0.529303703703704*t1*t2*t20*t22*t45*t47+
526 0.093333333333333*t1*t2*t20*t22*t28*t46-0.133333333333333*
527 t1*t16*t2*t21*t46-0.240592592592593*t1*t16*t2*t38*t45+0.28*
528 t1*t2*t22*t25*t28*t44-0.541333333333333*t1*t2*t20*t22*t25*
529 t35*t38-0.541333333333333*t1*t2*t22*t33*t34*t38+0.28*t1*t2*
530 t22*t28*t33*t35+0.373333333333333*t1*t16*t2*t25*t28*t35-0.1*
531 t1*(-18.66666666666666*t6*t7*grada*t18-0.083333333333333*t3*
532 t24*t17)*t2*t21*t22);
533 res[13] = -0.05*t1*t2*t21*t22*(3.5*t6*t7*t9+0.03125*t3*
534 t12*t5)+0.264651851851852*t1*t2*t22*pow(t25,4.0)*t47+0.186666666666667*
535 t1*t2*t22*t25*t28*t46-0.541333333333333*t1*t2*t22*t34*t35*
536 t38+0.14*t1*t2*t22*t28*pow(t35,2.0);
537
538 }
539
540 static void
pw86x_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)541 pw86x_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
542 {
543 real res[14];
544
545 pw86x_fourth_helper(dp->rhoa, dp->grada, res);
546
547 ds->df1000 += factor*res[0];
548 ds->df0010 += factor*res[1];
549
550 ds->df2000 += factor*res[2];
551 ds->df1010 += factor*res[3];
552 ds->df0020 += factor*res[4];
553
554 ds->df3000 += factor*res[5];
555 ds->df2010 += factor*res[6];
556 ds->df1020 += factor*res[7];
557 ds->df0030 += factor*res[8];
558
559 ds->df4000 += factor*res[9];
560 ds->df3010 += factor*res[10];
561 ds->df2020 += factor*res[11];
562 ds->df1030 += factor*res[12];
563 ds->df0040 += factor*res[13];
564
565
566 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
567 fabs(dp->grada-dp->gradb)>1e-13)
568 pw86x_fourth_helper(dp->rhob, dp->gradb, res);
569
570 ds->df0100 += factor*res[0];
571 ds->df0001 += factor*res[1];
572
573 ds->df0200 += factor*res[2];
574 ds->df0101 += factor*res[3];
575 ds->df0002 += factor*res[4];
576
577 ds->df0300 += factor*res[5];
578 ds->df0201 += factor*res[6];
579 ds->df0102 += factor*res[7];
580 ds->df0003 += factor*res[8];
581
582 ds->df0400 += factor*res[9];
583 ds->df0301 += factor*res[10];
584 ds->df0202 += factor*res[11];
585 ds->df0103 += factor*res[12];
586 ds->df0004 += factor*res[13];
587
588 }
589