1 /*
2
3
4 !
5 ! Dalton, a molecular electronic structure program
6 ! Copyright (C) by the authors of Dalton.
7 !
8 ! This program is free software; you can redistribute it and/or
9 ! modify it under the terms of the GNU Lesser General Public
10 ! License version 2.1 as published by the Free Software Foundation.
11 !
12 ! This program is distributed in the hope that it will be useful,
13 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
14 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 ! Lesser General Public License for more details.
16 !
17 ! If a copy of the GNU LGPL v2.1 was not distributed with this
18 ! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
19 !
20
21 !
22
23 */
24 /*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
25 /* fun-ft97bx.c:
26
27 Automatically generated code implementing FT97BX 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 ft97bxFunctional;" to 'functionals.h'
34 2. add "&ft97bxFunctional," to 'functionals.c'
35 3. add "fun-ft97bx.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 bet0: 0.002913644;
44 bet1: 0.000947417;
45 bet2: 2501.149;
46 beta: bet0+bet1*grada^2/(bet2+grada^2);
47 betb: bet0+bet1*gradb^2/(bet2+gradb^2);
48
49 xa: grada*rhoa^(-4/3);
50 xb: gradb*rhob^(-4/3);
51
52 sa: asinh(xa^2);
53 sb: asinh(xb^2);
54
55 da: 1+9*beta^2*xa^2*sa^2;
56 db: 1+9*betb^2*xb^2*sb^2;
57
58 Exa: -rhoa^(4/3)*beta*xa^2/(da^(1/2));
59 Exb: -rhob^(4/3)*betb*xb^2/(db^(1/2));
60
61 K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb);
62
63
64 ------ cut here -------
65 */
66
67
68 /* strictly conform to XOPEN ANSI C standard */
69 #if !defined(SYS_DEC)
70 /* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
71 * is not specified. */
72 #define _XOPEN_SOURCE 500
73 #define _XOPEN_SOURCE_EXTENDED 1
74 #endif
75 #include <math.h>
76 #include <stddef.h>
77 #include "general.h"
78
79 #define __CVERSION__
80
81 #include "functionals.h"
82
83 /* INTERFACE PART */
ft97bx_isgga(void)84 static integer ft97bx_isgga(void) { return 1; } /* FIXME: detect! */
85 static integer ft97bx_read(const char *conf_line);
86 static real ft97bx_energy(const FunDensProp* dp);
87 static void ft97bx_first(FunFirstFuncDrv *ds, real factor,
88 const FunDensProp* dp);
89 static void ft97bx_second(FunSecondFuncDrv *ds, real factor,
90 const FunDensProp* dp);
91 static void ft97bx_third(FunThirdFuncDrv *ds, real factor,
92 const FunDensProp* dp);
93 static void ft97bx_fourth(FunFourthFuncDrv *ds, real factor,
94 const FunDensProp* dp);
95
96 Functional FT97bxFunctional = {
97 "FT97bx", /* name */
98 ft97bx_isgga, /* gga-corrected */
99 1,
100 ft97bx_read,
101 NULL,
102 ft97bx_energy,
103 ft97bx_first,
104 ft97bx_second,
105 ft97bx_third,
106 ft97bx_fourth
107 };
108
109 /* IMPLEMENTATION PART */
110 static integer
ft97bx_read(const char * conf_line)111 ft97bx_read(const char *conf_line)
112 {
113 fun_set_hf_weight(0);
114 return 1;
115 }
116
117 static real
ft97bx_energy(const FunDensProp * dp)118 ft97bx_energy(const FunDensProp *dp)
119 {
120 real res;
121 real rhoa = dp->rhoa, rhob = dp->rhob;
122 real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
123
124 real t1, t2, t3, t4, t5, t6;
125
126 t1 = pow(grada,2.0);
127 t2 = 9.47417E-4*t1/(t1+2501.149)+0.002913644;
128 t3 = 1/pow(rhoa,2.666666666666667);
129 t4 = pow(gradb,2.0);
130 t5 = 9.47417E-4*t4/(t4+2501.149)+0.002913644;
131 t6 = 1/pow(rhob,2.666666666666667);
132
133 /* code */
134 res = -1.0*t4*t5/(sqrt(9.0*t4*pow(t5,2.0)*t6*pow(asinh(t4*
135 t6),2.0)+1.0)*pow(rhob,1.333333333333333))-1.0*t1*t2/(sqrt(9.0*
136 t1*pow(t2,2.0)*t3*pow(asinh(t1*t3),2.0)+1.0)*pow(rhoa,1.333333333333333));
137
138 return res;
139 }
140
141 static void
ft97bx_first_helper(real rhoa,real grada,real * res)142 ft97bx_first_helper(real rhoa, real grada, real *res)
143 { real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
144 real t11, t12, t13, t14, t15, t16, t17;
145
146 t1 = pow(grada,2.0);
147 t2 = t1+2501.149;
148 t3 = 1/t2;
149 t4 = 9.47417E-4*t1*t3+0.002913644;
150 t5 = pow(t4,2.0);
151 t6 = 1/pow(rhoa,2.666666666666667);
152 t7 = asinh(t1*t6);
153 t8 = pow(t7,2.0);
154 t9 = sqrt(9.0*t1*t5*t6*t8+1.0);
155 t10 = 1/t9;
156 t11 = pow(grada,4.0);
157 t12 = 1/pow(rhoa,5.333333333333333);
158 t13 = 1/sqrt(t11*t12+1.0);
159 t14 = 1/pow(t9,3.0);
160 t15 = 1/pow(rhoa,1.333333333333333);
161 t16 = pow(grada,3.0);
162 t17 = 0.001894834*grada*t3-0.001894834*t16/pow(t2,2.0);
163
164 /*
165 code */
166 res[0] = 0.5*t1*t14*t15*t4*(-48.0*t11*t13*t5*t7/pow(rhoa,
167 6.333333333333333)-24.0*t1*t5*t8/pow(rhoa,3.666666666666667))+
168 1.333333333333333*t1*t10*t4/pow(rhoa,2.333333333333333);
169 res[1] = 0.5*t1*t14*t15*t4*(18.0*t5*t6*t8*grada+18.0*
170 t1*t17*t4*t6*t8+36.0*t12*t13*t16*t5*t7)-2.0*t10*t15*t4*grada-
171 1.0*t1*t10*t15*t17;
172 }
173
174 static void
ft97bx_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)175 ft97bx_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
176 {
177 real res[2];
178
179 ft97bx_first_helper(dp->rhoa, dp->grada, res);
180 /* Final assignment */
181 ds->df1000 += factor*res[0];
182 ds->df0010 += factor*res[1];
183
184
185 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
186 fabs(dp->grada-dp->gradb)>1e-13)
187 ft97bx_first_helper(dp->rhob, dp->gradb, res);
188 ds->df0100 += factor*res[0];
189 ds->df0001 += factor*res[1];
190
191 }
192
193 static void
ft97bx_second_helper(real rhoa,real grada,real * res)194 ft97bx_second_helper(real rhoa, real grada, real *res)
195 {
196 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
197 real t11, t12, t13, t14, t15, t16, t17, t18;
198 real t19, t20, t21, t22, t23, t24, t25, t26;
199 real t27, t28, t29, t30;
200
201 t1 = pow(grada,2.0);
202 t2 = t1+2501.149;
203 t3 = 1/t2;
204 t4 = 9.47417E-4*t1*t3+0.002913644;
205 t5 = pow(t4,2.0);
206 t6 = 1/pow(rhoa,2.666666666666667);
207 t7 = asinh(t1*t6);
208 t8 = pow(t7,2.0);
209 t9 = sqrt(9.0*t1*t5*t6*t8+1.0);
210 t10 = 1/t9;
211 t11 = 1/pow(rhoa,2.333333333333333);
212 t12 = pow(grada,4.0);
213 t13 = 1/pow(rhoa,5.333333333333333);
214 t14 = t12*t13+1.0;
215 t15 = sqrt(t14);
216 t16 = 1/t15;
217 t17 = 1/pow(rhoa,6.333333333333333);
218 t18 = 1/pow(rhoa,3.666666666666667);
219 t19 = -24.0*t1*t18*t5*t8-48.0*t12*t16*t17*t5*t7;
220 t20 = 1/pow(t9,3.0);
221 t21 = 1/pow(rhoa,1.333333333333333);
222 t22 = pow(grada,3.0);
223 t23 = 1/pow(t2,2.0);
224 t24 = 0.001894834*grada*t3-0.001894834*t22*t23;
225 t25 = 18.0*t5*t6*t8*grada+18.0*t1*t24*t4*t6*t8+36.0*t13*
226 t16*t22*t5*t7;
227 t26 = 1/pow(t9,5.0);
228 t27 = 1/pow(t15,3.0);
229 t28 = pow(grada,6.0);
230 t29 = 1/t14;
231 t30 = 0.001894834*t3-0.00947417*t1*t23+0.007579336*t12/
232 pow(t2,3.0);
233
234 /* code */
235 res[0] = 0.5*t1*t19*t20*t21*t4+1.333333333333333*t1*t10*
236 t11*t4;
237 res[1] = -2.0*t10*t21*t4*grada+0.5*t1*t20*t21*t25*t4-
238 1.0*t1*t10*t21*t24;
239 res[2] = 0.5*t1*t20*t21*t4*(-128.0*t27*t5*t7*pow(grada,
240 8.0)/pow(rhoa,12.66666666666667)+128.0*t28*t29*t5/pow(rhoa,
241 10.0)+432.0*t12*t16*t5*t7/pow(rhoa,7.333333333333333)+88.0*
242 t1*t5*t8/pow(rhoa,4.666666666666667))-3.111111111111111*t1*
243 t10*t4/pow(rhoa,3.333333333333333)-0.75*t1*pow(t19,2.0)*t21*
244 t26*t4-1.333333333333333*t1*t11*t19*t20*t4;
245 res[3] = 0.5*t1*t20*t21*t4*(96.0*t27*t5*t7*pow(grada,
246 7.0)/pow(rhoa,11.66666666666667)-96.0*t29*t5*pow(grada,5.0)/
247 pow(rhoa,9.0)-48.0*t18*t5*t8*grada-48.0*t1*t18*t24*t4*t8-288.0*
248 t16*t17*t22*t5*t7-96.0*t12*t16*t17*t24*t4*t7)+2.666666666666667*
249 t10*t11*t4*grada-0.75*t1*t19*t21*t25*t26*t4-0.666666666666667*
250 t1*t11*t20*t25*t4+0.5*t1*t19*t20*t21*t24+1.333333333333333*
251 t1*t10*t11*t24+grada*t4*t19*t20*t21;
252 res[4] = 0.5*t1*t20*t21*t4*(-72.0*t27*t28*t5*t7/pow(rhoa,
253 10.66666666666667)+72.0*t12*t29*t5/pow(rhoa,8.0)+72.0*t24*
254 t4*t6*t8*grada+18.0*t5*t6*t8+18.0*t1*t30*t4*t6*t8+18.0*t1*
255 pow(t24,2.0)*t6*t8+180.0*t1*t13*t16*t5*t7+144.0*t13*t16*t22*
256 t24*t4*t7)+2.0*t20*t21*t25*t4*grada-4.0*t10*t21*t24*grada-
257 0.75*t1*t21*pow(t25,2.0)*t26*t4-2.0*t10*t21*t4-1.0*t1*t10*
258 t21*t30+t1*t24*t25*t20*t21;
259
260 }
261
262 static void
ft97bx_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)263 ft97bx_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
264 {
265 real res[5];
266
267 ft97bx_second_helper(dp->rhoa, dp->grada, res);
268
269 ds->df1000 += factor*res[0];
270 ds->df0010 += factor*res[1];
271
272 ds->df2000 += factor*res[2];
273 ds->df1010 += factor*res[3];
274 ds->df0020 += factor*res[4];
275
276
277 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
278 fabs(dp->grada-dp->gradb)>1e-13)
279 ft97bx_second_helper(dp->rhob, dp->gradb, res);
280 ds->df0100 += factor*res[0];
281 ds->df0001 += factor*res[1];
282
283 ds->df0200 += factor*res[2];
284 ds->df0101 += factor*res[3];
285 ds->df0002 += factor*res[4];
286
287 }
288
289 static void
ft97bx_third_helper(real rhoa,real grada,real * res)290 ft97bx_third_helper(real rhoa, real grada, real *res)
291 {
292 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
293 real t11, t12, t13, t14, t15, t16, t17, t18;
294 real t19, t20, t21, t22, t23, t24, t25, t26;
295 real t27, t28, t29, t30, t31, t32, t33, t34;
296 real t35, t36, t37, t38, t39, t40, t41, t42;
297 real t43, t44, t45, t46, t47, t48, t49, t50;
298 real t51, t52, t53, t54, t55;
299
300 t1 = pow(grada,2.0);
301 t2 = t1+2501.149;
302 t3 = 1/t2;
303 t4 = 9.47417E-4*t1*t3+0.002913644;
304 t5 = pow(t4,2.0);
305 t6 = 1/pow(rhoa,2.666666666666667);
306 t7 = asinh(t1*t6);
307 t8 = pow(t7,2.0);
308 t9 = sqrt(9.0*t1*t5*t6*t8+1.0);
309 t10 = 1/t9;
310 t11 = 1/pow(rhoa,2.333333333333333);
311 t12 = pow(grada,4.0);
312 t13 = 1/pow(rhoa,5.333333333333333);
313 t14 = t12*t13+1.0;
314 t15 = sqrt(t14);
315 t16 = 1/t15;
316 t17 = 1/pow(rhoa,6.333333333333333);
317 t18 = 1/pow(rhoa,3.666666666666667);
318 t19 = -24.0*t1*t18*t5*t8-48.0*t12*t16*t17*t5*t7;
319 t20 = 1/pow(t9,3.0);
320 t21 = 1/pow(rhoa,1.333333333333333);
321 t22 = pow(grada,3.0);
322 t23 = 1/pow(t2,2.0);
323 t24 = 0.001894834*grada*t3-0.001894834*t22*t23;
324 t25 = 18.0*t5*t6*t8*grada+18.0*t1*t24*t4*t6*t8+36.0*t13*
325 t16*t22*t5*t7;
326 t26 = 1/pow(rhoa,3.333333333333333);
327 t27 = pow(t19,2.0);
328 t28 = 1/pow(t9,5.0);
329 t29 = pow(grada,8.0);
330 t30 = 1/pow(t15,3.0);
331 t31 = 1/pow(rhoa,12.66666666666667);
332 t32 = pow(grada,6.0);
333 t33 = 1/t14;
334 t34 = 1/pow(rhoa,10.0);
335 t35 = 1/pow(rhoa,7.333333333333333);
336 t36 = 1/pow(rhoa,4.666666666666667);
337 t37 = 88.0*t1*t36*t5*t8+432.0*t12*t16*t35*t5*t7-128.0*
338 t29*t30*t31*t5*t7+128.0*t32*t33*t34*t5;
339 t38 = pow(grada,7.0);
340 t39 = 1/pow(rhoa,11.66666666666667);
341 t40 = pow(grada,5.0);
342 t41 = 1/pow(rhoa,9.0);
343 t42 = -48.0*t18*t5*t8*grada-48.0*t1*t18*t24*t4*t8+96.0*
344 t30*t38*t39*t5*t7-288.0*t16*t17*t22*t5*t7-96.0*t12*t16*t17*
345 t24*t4*t7-96.0*t33*t40*t41*t5;
346 t43 = pow(t25,2.0);
347 t44 = 1/pow(rhoa,10.66666666666667);
348 t45 = 1/pow(rhoa,8.0);
349 t46 = pow(t24,2.0);
350 t47 = 1/pow(t2,3.0);
351 t48 = 0.001894834*t3-0.00947417*t1*t23+0.007579336*t12*
352 t47;
353 t49 = 72.0*t24*t4*t6*t8*grada+18.0*t5*t6*t8+18.0*t1*t4*
354 t48*t6*t8+18.0*t1*t46*t6*t8-72.0*t30*t32*t44*t5*t7+180.0*t1*
355 t13*t16*t5*t7+144.0*t13*t16*t22*t24*t4*t7+72.0*t12*t33*t45*
356 t5;
357 t50 = 1/pow(t9,7.0);
358 t51 = 1/pow(t15,5.0);
359 t52 = pow(grada,10.0);
360 t53 = 1/pow(t14,2.0);
361 t54 = pow(grada,9.0);
362 t55 = 0.068214024*t22*t47-0.045476016*t40/pow(t2,4.0)-
363 0.022738008*grada*t23;
364
365 /* code */
366 res[0] = 0.5*t1*t19*t20*t21*t4+1.333333333333333*t1*t10*
367 t11*t4;
368 res[1] = -2.0*t10*t21*t4*grada+0.5*t1*t20*t21*t25*t4-
369 1.0*t1*t10*t21*t24;
370 res[2] = 0.5*t1*t20*t21*t37*t4-0.75*t1*t21*t27*t28*t4-
371 3.111111111111111*t1*t10*t26*t4-1.333333333333333*t1*t11*t19*
372 t20*t4;
373 res[3] = 2.666666666666667*t10*t11*t4*grada+0.5*t1*t20*
374 t21*t4*t42-0.75*t1*t19*t21*t25*t28*t4-0.666666666666667*t1*
375 t11*t20*t25*t4+0.5*t1*t19*t20*t21*t24+1.333333333333333*t1*
376 t10*t11*t24+grada*t4*t19*t20*t21;
377 res[4] = 2.0*t20*t21*t25*t4*grada-4.0*t10*t21*t24*grada+
378 0.5*t1*t20*t21*t4*t49-1.0*t1*t10*t21*t48-0.75*t1*t21*t28*t4*
379 t43-2.0*t10*t21*t4+t1*t24*t25*t20*t21;
380 res[5] = 0.5*t1*t20*t21*t4*(-1024.0*t5*t51*t7*pow(grada,
381 12.0)/pow(rhoa,19.0)+1024.0*t5*t52*t53/pow(rhoa,16.33333333333333)+
382 2773.333333333333*t29*t30*t5*t7/pow(rhoa,13.66666666666667)-
383 2432.0*t32*t33*t5/pow(rhoa,11.0)-3637.333333333333*t12*t16*
384 t5*t7/pow(rhoa,8.333333333333334)-410.6666666666667*t1*t5*
385 t8/pow(rhoa,5.666666666666667))+10.37037037037037*t1*t10*t4/
386 pow(rhoa,4.333333333333333)+1.875*t1*pow(t19,3.0)*t21*t4*t50-
387 2.25*t1*t19*t21*t28*t37*t4-2.0*t1*t11*t20*t37*t4+3.0*t1*t11*
388 t27*t28*t4+4.666666666666667*t1*t19*t20*t26*t4;
389 res[6] = 0.5*t1*t20*t21*t4*(768.0*t5*t51*t7*pow(grada,
390 11.0)/pow(rhoa,18.0)-768.0*t5*t53*t54/pow(rhoa,15.33333333333333)+
391 176.0*t36*t5*t8*grada+176.0*t1*t24*t36*t4*t8-1888.0*t30*t31*
392 t38*t5*t7+2080.0*t16*t22*t35*t5*t7+864.0*t12*t16*t24*t35*t4*
393 t7-256.0*t24*t29*t30*t31*t4*t7+1632.0*t33*t34*t40*t5+256.0*
394 t24*t32*t33*t34*t4)-1.5*t21*t27*t28*t4*grada-6.222222222222222*
395 t10*t26*t4*grada-2.666666666666667*t11*t19*t20*t4*grada+1.875*
396 t1*t21*t25*t27*t4*t50-1.5*t1*t19*t21*t28*t4*t42-1.333333333333333*
397 t1*t11*t20*t4*t42-0.75*t1*t21*t25*t28*t37*t4+2.0*t1*t11*t19*
398 t25*t28*t4+1.555555555555556*t1*t20*t25*t26*t4+0.5*t1*t20*
399 t21*t24*t37-0.75*t1*t21*t24*t27*t28-3.111111111111111*t1*t10*
400 t24*t26-1.333333333333333*t1*t11*t19*t20*t24+grada*t4*t37*
401 t20*t21;
402 res[7] = 0.5*t1*t20*t21*t4*(-576.0*t5*t51*t52*t7/pow(rhoa,
403 17.0)+576.0*t29*t5*t53/pow(rhoa,14.33333333333333)-192.0*t18*
404 t24*t4*t8*grada-48.0*t18*t5*t8-48.0*t1*t18*t4*t48*t8-48.0*
405 t1*t18*t46*t8+1248.0*t30*t32*t39*t5*t7-1056.0*t1*t16*t17*t5*
406 t7-96.0*t12*t16*t17*t4*t48*t7-96.0*t12*t16*t17*t46*t7+384.0*
407 t24*t30*t38*t39*t4*t7-1152.0*t16*t17*t22*t24*t4*t7-1056.0*
408 t12*t33*t41*t5-384.0*t24*t33*t4*t40*t41)+2.0*t20*t21*t4*t42*
409 grada-3.0*t19*t21*t25*t28*t4*grada-2.666666666666667*t11*t20*
410 t25*t4*grada+2.0*t19*t20*t21*t24*grada+5.333333333333333*t10*
411 t11*t24*grada+1.875*t1*t19*t21*t4*t43*t50-0.75*t1*t19*t21*
412 t28*t4*t49-0.666666666666667*t1*t11*t20*t4*t49+0.5*t1*t19*
413 t20*t21*t48+1.333333333333333*t1*t10*t11*t48-1.5*t1*t21*t25*
414 t28*t4*t42+2.666666666666667*t10*t11*t4-1.5*t1*t19*t21*t24*
415 t25*t28-1.333333333333333*t1*t11*t20*t24*t25+t1*t24*t42*t20*
416 t21+t4*t19*t20*t21+t1*t4*t43*t28*t11;
417 res[8] = 0.5*t1*t20*t21*t4*(432.0*t5*t51*t54*t7/pow(rhoa,
418 16.0)-432.0*t38*t5*t53/pow(rhoa,13.33333333333333)+108.0*t4*
419 t48*t6*t8*grada+108.0*t46*t6*t8*grada+432.0*t13*t16*t5*t7*
420 grada+18.0*t1*t4*t55*t6*t8+54.0*t1*t24*t48*t6*t8+108.0*t24*
421 t4*t6*t8-792.0*t30*t40*t44*t5*t7+216.0*t13*t16*t22*t4*t48*
422 t7+216.0*t13*t16*t22*t46*t7-432.0*t24*t30*t32*t4*t44*t7+1080.0*
423 t1*t13*t16*t24*t4*t7+648.0*t22*t33*t45*t5+432.0*t12*t24*t33*
424 t4*t45)+3.0*t20*t21*t4*t49*grada-6.0*t10*t21*t48*grada-4.5*
425 t21*t28*t4*t43*grada+6.0*t20*t21*t24*t25*grada-1.0*t1*t10*
426 t21*t55+1.875*t1*t21*pow(t25,3.0)*t4*t50-2.25*t1*t21*t25*t28*
427 t4*t49+1.5*t1*t20*t21*t24*t49+1.5*t1*t20*t21*t25*t48-2.25*
428 t1*t21*t24*t28*t43+3.0*t20*t21*t25*t4-6.0*t10*t21*t24;
429
430 }
431
432 static void
ft97bx_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)433 ft97bx_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
434 {
435 real res[9];
436
437 ft97bx_third_helper(dp->rhoa, dp->grada, res);
438
439 ds->df1000 += factor*res[0];
440 ds->df0010 += factor*res[1];
441
442 ds->df2000 += factor*res[2];
443 ds->df1010 += factor*res[3];
444 ds->df0020 += factor*res[4];
445
446 ds->df3000 += factor*res[5];
447 ds->df2010 += factor*res[6];
448 ds->df1020 += factor*res[7];
449 ds->df0030 += factor*res[8];
450
451
452 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
453 fabs(dp->grada-dp->gradb)>1e-13)
454 ft97bx_third_helper(dp->rhob, dp->gradb, res);
455
456 ds->df0100 += factor*res[0];
457 ds->df0001 += factor*res[1];
458
459 ds->df0200 += factor*res[2];
460 ds->df0101 += factor*res[3];
461 ds->df0002 += factor*res[4];
462
463 ds->df0300 += factor*res[5];
464 ds->df0201 += factor*res[6];
465 ds->df0102 += factor*res[7];
466 ds->df0003 += factor*res[8];
467
468 }
469
470 static void
ft97bx_fourth_helper(real rhoa,real grada,real * res)471 ft97bx_fourth_helper(real rhoa, real grada, real *res)
472 {
473 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
474 real t11, t12, t13, t14, t15, t16, t17, t18;
475 real t19, t20, t21, t22, t23, t24, t25, t26;
476 real t27, t28, t29, t30, t31, t32, t33, t34;
477 real t35, t36, t37, t38, t39, t40, t41, t42;
478 real t43, t44, t45, t46, t47, t48, t49, t50;
479 real t51, t52, t53, t54, t55, t56, t57, t58;
480 real t59, t60, t61, t62, t63, t64, t65, t66;
481 real t67, t68, t69, t70, t71, t72, t73, t74;
482 real t75, t76, t77, t78, t79, t80, t81, t82;
483 real t83;
484
485 t1 = pow(grada,2.0);
486 t2 = t1+2501.149;
487 t3 = 1/t2;
488 t4 = 9.47417E-4*t1*t3+0.002913644;
489 t5 = pow(t4,2.0);
490 t6 = 1/pow(rhoa,2.666666666666667);
491 t7 = asinh(t1*t6);
492 t8 = pow(t7,2.0);
493 t9 = sqrt(9.0*t1*t5*t6*t8+1.0);
494 t10 = 1/t9;
495 t11 = 1/pow(rhoa,2.333333333333333);
496 t12 = pow(grada,4.0);
497 t13 = 1/pow(rhoa,5.333333333333333);
498 t14 = t12*t13+1.0;
499 t15 = sqrt(t14);
500 t16 = 1/t15;
501 t17 = 1/pow(rhoa,6.333333333333333);
502 t18 = 1/pow(rhoa,3.666666666666667);
503 t19 = -24.0*t1*t18*t5*t8-48.0*t12*t16*t17*t5*t7;
504 t20 = 1/pow(t9,3.0);
505 t21 = 1/pow(rhoa,1.333333333333333);
506 t22 = pow(grada,3.0);
507 t23 = 1/pow(t2,2.0);
508 t24 = 0.001894834*grada*t3-0.001894834*t22*t23;
509 t25 = 18.0*t5*t6*t8*grada+18.0*t1*t24*t4*t6*t8+36.0*t13*
510 t16*t22*t5*t7;
511 t26 = 1/pow(rhoa,3.333333333333333);
512 t27 = pow(t19,2.0);
513 t28 = 1/pow(t9,5.0);
514 t29 = pow(grada,8.0);
515 t30 = 1/pow(t15,3.0);
516 t31 = 1/pow(rhoa,12.66666666666667);
517 t32 = pow(grada,6.0);
518 t33 = 1/t14;
519 t34 = 1/pow(rhoa,10.0);
520 t35 = 1/pow(rhoa,7.333333333333333);
521 t36 = 1/pow(rhoa,4.666666666666667);
522 t37 = 88.0*t1*t36*t5*t8+432.0*t12*t16*t35*t5*t7-128.0*
523 t29*t30*t31*t5*t7+128.0*t32*t33*t34*t5;
524 t38 = pow(grada,7.0);
525 t39 = 1/pow(rhoa,11.66666666666667);
526 t40 = pow(grada,5.0);
527 t41 = 1/pow(rhoa,9.0);
528 t42 = -48.0*t18*t5*t8*grada-48.0*t1*t18*t24*t4*t8+96.0*
529 t30*t38*t39*t5*t7-288.0*t16*t17*t22*t5*t7-96.0*t12*t16*t17*
530 t24*t4*t7-96.0*t33*t40*t41*t5;
531 t43 = pow(t25,2.0);
532 t44 = 1/pow(rhoa,10.66666666666667);
533 t45 = 1/pow(rhoa,8.0);
534 t46 = pow(t24,2.0);
535 t47 = 1/pow(t2,3.0);
536 t48 = 0.001894834*t3-0.00947417*t1*t23+0.007579336*t12*
537 t47;
538 t49 = 72.0*t24*t4*t6*t8*grada+18.0*t5*t6*t8+18.0*t1*t4*
539 t48*t6*t8+18.0*t1*t46*t6*t8-72.0*t30*t32*t44*t5*t7+180.0*t1*
540 t13*t16*t5*t7+144.0*t13*t16*t22*t24*t4*t7+72.0*t12*t33*t45*
541 t5;
542 t50 = 1/pow(rhoa,4.333333333333333);
543 t51 = pow(t19,3.0);
544 t52 = 1/pow(t9,7.0);
545 t53 = pow(grada,12.0);
546 t54 = 1/pow(t15,5.0);
547 t55 = 1/pow(rhoa,19.0);
548 t56 = pow(grada,10.0);
549 t57 = 1/pow(t14,2.0);
550 t58 = 1/pow(rhoa,16.33333333333333);
551 t59 = 1/pow(rhoa,13.66666666666667);
552 t60 = 1/pow(rhoa,11.0);
553 t61 = 1/pow(rhoa,8.333333333333334);
554 t62 = 1/pow(rhoa,5.666666666666667);
555 t63 = -410.6666666666667*t1*t5*t62*t8-3637.333333333333*
556 t12*t16*t5*t61*t7+2773.333333333333*t29*t30*t5*t59*t7-1024.0*
557 t5*t53*t54*t55*t7-2432.0*t32*t33*t5*t60+1024.0*t5*t56*t57*
558 t58;
559 t64 = pow(grada,11.0);
560 t65 = 1/pow(rhoa,18.0);
561 t66 = pow(grada,9.0);
562 t67 = 1/pow(rhoa,15.33333333333333);
563 t68 = 176.0*t36*t5*t8*grada+176.0*t1*t24*t36*t4*t8+768.0*
564 t5*t54*t64*t65*t7-1888.0*t30*t31*t38*t5*t7+2080.0*t16*t22*
565 t35*t5*t7+864.0*t12*t16*t24*t35*t4*t7-256.0*t24*t29*t30*t31*
566 t4*t7-768.0*t5*t57*t66*t67+1632.0*t33*t34*t40*t5+256.0*t24*
567 t32*t33*t34*t4;
568 t69 = 1/pow(rhoa,17.0);
569 t70 = 1/pow(rhoa,14.33333333333333);
570 t71 = -192.0*t18*t24*t4*t8*grada-48.0*t18*t5*t8-48.0*
571 t1*t18*t4*t48*t8-48.0*t1*t18*t46*t8+576.0*t29*t5*t57*t70-576.0*
572 t5*t54*t56*t69*t7+1248.0*t30*t32*t39*t5*t7-1056.0*t1*t16*t17*
573 t5*t7-96.0*t12*t16*t17*t4*t48*t7-96.0*t12*t16*t17*t46*t7+384.0*
574 t24*t30*t38*t39*t4*t7-1152.0*t16*t17*t22*t24*t4*t7-1056.0*
575 t12*t33*t41*t5-384.0*t24*t33*t4*t40*t41;
576 t72 = pow(t25,3.0);
577 t73 = 1/pow(rhoa,16.0);
578 t74 = 1/pow(rhoa,13.33333333333333);
579 t75 = 1/pow(t2,4.0);
580 t76 = -0.022738008*grada*t23+0.068214024*t22*t47-0.045476016*
581 t40*t75;
582 t77 = 108.0*t4*t48*t6*t8*grada+108.0*t46*t6*t8*grada+
583 432.0*t13*t16*t5*t7*grada+18.0*t1*t4*t6*t76*t8+54.0*t1*t24*
584 t48*t6*t8+108.0*t24*t4*t6*t8-432.0*t38*t5*t57*t74+432.0*t5*
585 t54*t66*t7*t73-792.0*t30*t40*t44*t5*t7+216.0*t13*t16*t22*t4*
586 t48*t7+216.0*t13*t16*t22*t46*t7-432.0*t24*t30*t32*t4*t44*t7+
587 1080.0*t1*t13*t16*t24*t4*t7+648.0*t22*t33*t45*t5+432.0*t12*
588 t24*t33*t4*t45;
589 t78 = 1/pow(t9,9.0);
590 t79 = 1/pow(t15,7.0);
591 t80 = pow(grada,14.0);
592 t81 = 1/pow(t14,3.0);
593 t82 = pow(grada,13.0);
594 t83 = -0.636664224*t12*t75+0.295594104*t1*t47+0.363808128*
595 t32/pow(t2,5.0)-0.022738008*t23;
596
597 /* code */
598 res[0] = 0.5*t1*t19*t20*t21*t4+1.333333333333333*t1*t10*
599 t11*t4;
600 res[1] = -2.0*t10*t21*t4*grada+0.5*t1*t20*t21*t25*t4-
601 1.0*t1*t10*t21*t24;
602 res[2] = 0.5*t1*t20*t21*t37*t4-0.75*t1*t21*t27*t28*t4-
603 3.111111111111111*t1*t10*t26*t4-1.333333333333333*t1*t11*t19*
604 t20*t4;
605 res[3] = 2.666666666666667*t10*t11*t4*grada+0.5*t1*t20*
606 t21*t4*t42-0.75*t1*t19*t21*t25*t28*t4-0.666666666666667*t1*
607 t11*t20*t25*t4+0.5*t1*t19*t20*t21*t24+1.333333333333333*t1*
608 t10*t11*t24+grada*t4*t19*t20*t21;
609 res[4] = 2.0*t20*t21*t25*t4*grada-4.0*t10*t21*t24*grada+
610 0.5*t1*t20*t21*t4*t49-1.0*t1*t10*t21*t48-0.75*t1*t21*t28*t4*
611 t43-2.0*t10*t21*t4+t1*t24*t25*t20*t21;
612 res[5] = 0.5*t1*t20*t21*t4*t63+1.875*t1*t21*t4*t51*t52+
613 10.37037037037037*t1*t10*t4*t50-2.25*t1*t19*t21*t28*t37*t4-
614 2.0*t1*t11*t20*t37*t4+3.0*t1*t11*t27*t28*t4+4.666666666666667*
615 t1*t19*t20*t26*t4;
616 res[6] = -1.5*t21*t27*t28*t4*grada-6.222222222222222*
617 t10*t26*t4*grada-2.666666666666667*t11*t19*t20*t4*grada+0.5*
618 t1*t20*t21*t4*t68+1.875*t1*t21*t25*t27*t4*t52-1.5*t1*t19*t21*
619 t28*t4*t42-1.333333333333333*t1*t11*t20*t4*t42-0.75*t1*t21*
620 t25*t28*t37*t4+2.0*t1*t11*t19*t25*t28*t4+1.555555555555556*
621 t1*t20*t25*t26*t4+0.5*t1*t20*t21*t24*t37-0.75*t1*t21*t24*t27*
622 t28-3.111111111111111*t1*t10*t24*t26-1.333333333333333*t1*
623 t11*t19*t20*t24+grada*t4*t37*t20*t21;
624 res[7] = 2.0*t20*t21*t4*t42*grada-3.0*t19*t21*t25*t28*
625 t4*grada-2.666666666666667*t11*t20*t25*t4*grada+2.0*t19*t20*
626 t21*t24*grada+5.333333333333333*t10*t11*t24*grada+0.5*t1*t20*
627 t21*t4*t71+1.875*t1*t19*t21*t4*t43*t52-0.75*t1*t19*t21*t28*
628 t4*t49-0.666666666666667*t1*t11*t20*t4*t49+0.5*t1*t19*t20*
629 t21*t48+1.333333333333333*t1*t10*t11*t48-1.5*t1*t21*t25*t28*
630 t4*t42+2.666666666666667*t10*t11*t4-1.5*t1*t19*t21*t24*t25*
631 t28-1.333333333333333*t1*t11*t20*t24*t25+t1*t24*t42*t20*t21+
632 t4*t19*t20*t21+t1*t4*t43*t28*t11;
633 res[8] = 3.0*t20*t21*t4*t49*grada-6.0*t10*t21*t48*grada-
634 4.5*t21*t28*t4*t43*grada+6.0*t20*t21*t24*t25*grada+0.5*t1*
635 t20*t21*t4*t77-1.0*t1*t10*t21*t76+1.875*t1*t21*t4*t52*t72-
636 2.25*t1*t21*t25*t28*t4*t49+1.5*t1*t20*t21*t24*t49+1.5*t1*t20*
637 t21*t25*t48-2.25*t1*t21*t24*t28*t43+3.0*t20*t21*t25*t4-6.0*
638 t10*t21*t24;
639 res[9] = 0.5*t1*t20*t21*t4*(-13653.33333333333*t5*t7*
640 t79*pow(grada,16.0)/pow(rhoa,25.33333333333333)+13653.33333333333*
641 t5*t80*t81/pow(rhoa,22.66666666666667)+41642.66666666666*t5*
642 t53*t54*t7/pow(rhoa,20.0)-37091.55555555555*t5*t56*t57/pow(rhoa,
643 17.33333333333333)-47601.77777777778*t29*t30*t5*t7/pow(rhoa,
644 14.66666666666667)+36451.55555555555*t32*t33*t5/pow(rhoa,12.0)+
645 32501.33333333333*t12*t16*t5*t7/pow(rhoa,9.333333333333334)+
646 2327.111111111111*t1*t5*t8/pow(rhoa,6.666666666666667))-6.5625*
647 t1*pow(t19,4.0)*t21*t4*t78-3.0*t1*t19*t21*t28*t4*t63-2.666666666666667*
648 t1*t11*t20*t4*t63-10.0*t1*t11*t4*t51*t52+11.25*t1*t21*t27*
649 t37*t4*t52-20.74074074074074*t1*t19*t20*t4*t50-2.25*t1*t21*
650 t28*pow(t37,2.0)*t4+12.0*t1*t11*t19*t28*t37*t4+9.333333333333334*
651 t1*t20*t26*t37*t4-14.0*t1*t26*t27*t28*t4-44.93827160493827*
652 t1*t10*t13*t4;
653 res[10] = 0.5*t1*t20*t21*t4*(10240.0*t5*t7*t79*pow(grada,
654 15.0)/pow(rhoa,24.33333333333333)-10240.0*t5*t81*t82/pow(rhoa,
655 21.66666666666667)-821.3333333333334*t5*t62*t8*grada-821.3333333333334*
656 t1*t24*t4*t62*t8-28928.0*t5*t54*t55*t64*t7-16192.0*t16*t22*
657 t5*t61*t7-7274.666666666667*t12*t16*t24*t4*t61*t7+29461.33333333333*
658 t30*t38*t5*t59*t7+5546.666666666667*t24*t29*t30*t4*t59*t7-
659 2048.0*t24*t4*t53*t54*t55*t7+25514.66666666667*t5*t57*t58*
660 t66-21866.66666666667*t33*t40*t5*t60-4864.0*t24*t32*t33*t4*
661 t60+2048.0*t24*t4*t56*t57*t58)+3.75*t21*t4*t51*t52*grada+20.74074074074074*
662 t10*t4*t50*grada-4.5*t19*t21*t28*t37*t4*grada-4.0*t11*t20*
663 t37*t4*grada+6.0*t11*t27*t28*t4*grada+9.333333333333334*t19*
664 t20*t26*t4*grada-6.5625*t1*t21*t25*t4*t51*t78-2.25*t1*t19*
665 t21*t28*t4*t68-2.0*t1*t11*t20*t4*t68-0.75*t1*t21*t25*t28*t4*
666 t63+0.5*t1*t20*t21*t24*t63+1.875*t1*t21*t24*t51*t52+5.625*
667 t1*t21*t27*t4*t42*t52+5.625*t1*t19*t21*t25*t37*t4*t52-7.5*
668 t1*t11*t25*t27*t4*t52-5.185185185185185*t1*t20*t25*t4*t50+
669 10.37037037037037*t1*t10*t24*t50-2.25*t1*t21*t28*t37*t4*t42+
670 6.0*t1*t11*t19*t28*t4*t42+4.666666666666667*t1*t20*t26*t4*
671 t42+3.0*t1*t11*t25*t28*t37*t4-7.0*t1*t19*t25*t26*t28*t4-2.25*
672 t1*t19*t21*t24*t28*t37-2.0*t1*t11*t20*t24*t37+3.0*t1*t11*t24*
673 t27*t28+4.666666666666667*t1*t19*t20*t24*t26+grada*t4*t63*
674 t20*t21;
675 res[11] = 0.5*t1*t20*t21*t4*(-7680.0*t5*t7*t79*t80/pow(rhoa,
676 23.33333333333333)+7680.0*t5*t53*t81/pow(rhoa,20.66666666666667)+
677 704.0*t24*t36*t4*t8*grada+176.0*t36*t5*t8+176.0*t1*t36*t4*
678 t48*t8+176.0*t1*t36*t46*t8+3072.0*t24*t4*t54*t64*t65*t7+19776.0*
679 t5*t54*t56*t65*t7+6944.0*t1*t16*t35*t5*t7-17376.0*t30*t31*
680 t32*t5*t7+864.0*t12*t16*t35*t4*t48*t7-256.0*t29*t30*t31*t4*
681 t48*t7+864.0*t12*t16*t35*t46*t7-256.0*t29*t30*t31*t46*t7-7552.0*
682 t24*t30*t31*t38*t4*t7+8320.0*t16*t22*t24*t35*t4*t7-3072.0*
683 t24*t4*t57*t66*t67-17216.0*t29*t5*t57*t67+12320.0*t12*t33*
684 t34*t5+256.0*t32*t33*t34*t4*t48+256.0*t32*t33*t34*t46+6528.0*
685 t24*t33*t34*t4*t40)+2.0*t20*t21*t4*t68*grada+7.5*t21*t25*t27*
686 t4*t52*grada-6.0*t19*t21*t28*t4*t42*grada-5.333333333333333*
687 t11*t20*t4*t42*grada-3.0*t21*t25*t28*t37*t4*grada+8.0*t11*
688 t19*t25*t28*t4*grada+6.222222222222222*t20*t25*t26*t4*grada+
689 2.0*t20*t21*t24*t37*grada-3.0*t21*t24*t27*t28*grada-12.44444444444444*
690 t10*t24*t26*grada-5.333333333333333*t11*t19*t20*t24*grada-
691 6.5625*t1*t21*t27*t4*t43*t78-1.5*t1*t19*t21*t28*t4*t71-1.333333333333333*
692 t1*t11*t20*t4*t71-1.5*t1*t21*t25*t28*t4*t68+1.875*t1*t21*t27*
693 t4*t49*t52+1.875*t1*t21*t37*t4*t43*t52-5.0*t1*t11*t19*t4*t43*
694 t52+7.5*t1*t19*t21*t25*t4*t42*t52+3.75*t1*t21*t24*t25*t27*
695 t52-0.75*t1*t21*t28*t37*t4*t49+2.0*t1*t11*t19*t28*t4*t49+1.555555555555556*
696 t1*t20*t26*t4*t49+0.5*t1*t20*t21*t37*t48-0.75*t1*t21*t27*t28*
697 t48-3.111111111111111*t1*t10*t26*t48-1.333333333333333*t1*
698 t11*t19*t20*t48-2.333333333333333*t1*t26*t28*t4*t43-1.5*t1*
699 t21*t28*t4*pow(t42,2.0)+4.0*t1*t11*t25*t28*t4*t42-3.0*t1*t19*
700 t21*t24*t28*t42-2.666666666666667*t1*t11*t20*t24*t42-1.5*t21*
701 t27*t28*t4-6.222222222222222*t10*t26*t4-2.666666666666667*
702 t11*t19*t20*t4-1.5*t1*t21*t24*t25*t28*t37+4.0*t1*t11*t19*t24*
703 t25*t28+3.111111111111111*t1*t20*t24*t25*t26+t1*t24*t68*t20*
704 t21+t4*t37*t20*t21;
705 res[12] = 0.5*t1*t20*t21*t4*(5760.0*t5*t7*t79*t82/pow(rhoa,
706 22.33333333333333)-5760.0*t5*t64*t81/pow(rhoa,19.66666666666667)-
707 288.0*t18*t4*t48*t8*grada-288.0*t18*t46*t8*grada-2304.0*t16*
708 t17*t5*t7*grada-48.0*t1*t18*t4*t76*t8-144.0*t1*t18*t24*t48*
709 t8-288.0*t18*t24*t4*t8-96.0*t12*t16*t17*t4*t7*t76+11328.0*
710 t38*t5*t57*t70+3456.0*t24*t29*t4*t57*t70-13248.0*t5*t54*t66*
711 t69*t7-3456.0*t24*t4*t54*t56*t69*t7+9600.0*t30*t39*t40*t5*
712 t7+576.0*t30*t38*t39*t4*t48*t7-1728.0*t16*t17*t22*t4*t48*t7-
713 288.0*t12*t16*t17*t24*t48*t7+576.0*t30*t38*t39*t46*t7-1728.0*
714 t16*t17*t22*t46*t7+7488.0*t24*t30*t32*t39*t4*t7-6336.0*t1*
715 t16*t17*t24*t4*t7-6336.0*t22*t33*t41*t5-576.0*t33*t4*t40*t41*
716 t48-576.0*t33*t40*t41*t46-6336.0*t12*t24*t33*t4*t41)+3.0*t20*
717 t21*t4*t71*grada+11.25*t19*t21*t4*t43*t52*grada-4.5*t19*t21*
718 t28*t4*t49*grada-4.0*t11*t20*t4*t49*grada+3.0*t19*t20*t21*
719 t48*grada+8.0*t10*t11*t48*grada+6.0*t11*t28*t4*t43*grada-9.0*
720 t21*t25*t28*t4*t42*grada+6.0*t20*t21*t24*t42*grada-9.0*t19*
721 t21*t24*t25*t28*grada-8.0*t11*t20*t24*t25*grada-6.5625*t1*
722 t19*t21*t4*t72*t78-0.75*t1*t19*t21*t28*t4*t77-0.666666666666667*
723 t1*t11*t20*t4*t77+0.5*t1*t19*t20*t21*t76+1.333333333333333*
724 t1*t10*t11*t76-2.5*t1*t11*t4*t52*t72-2.25*t1*t21*t25*t28*t4*
725 t71+1.5*t1*t20*t21*t24*t71+5.625*t1*t19*t21*t25*t4*t49*t52+
726 5.625*t1*t21*t4*t42*t43*t52+5.625*t1*t19*t21*t24*t43*t52-2.25*
727 t1*t21*t28*t4*t42*t49+3.0*t1*t11*t25*t28*t4*t49-2.25*t1*t19*
728 t21*t24*t28*t49-2.0*t1*t11*t20*t24*t49+1.5*t1*t20*t21*t42*
729 t48-2.25*t1*t19*t21*t25*t28*t48-2.0*t1*t11*t20*t25*t48+3.0*
730 t1*t11*t24*t28*t43+3.0*t20*t21*t4*t42-4.5*t1*t21*t24*t25*t28*
731 t42-4.5*t19*t21*t25*t28*t4-4.0*t11*t20*t25*t4+3.0*t19*t20*
732 t21*t24+8.0*t10*t11*t24;
733 res[13] = 0.5*t1*t20*t21*t4*(-4320.0*t5*t53*t7*t79/pow(rhoa,
734 21.33333333333333)+4320.0*t5*t56*t81/pow(rhoa,18.66666666666667)+
735 144.0*t4*t6*t76*t8*grada+432.0*t24*t48*t6*t8*grada+3456.0*
736 t13*t16*t24*t4*t7*grada+18.0*t1*t4*t6*t8*t83+72.0*t1*t24*t6*
737 t76*t8+54.0*t1*pow(t48,2.0)*t6*t8+216.0*t4*t48*t6*t8+216.0*
738 t46*t6*t8+288.0*t13*t16*t22*t4*t7*t76-7200.0*t32*t5*t57*t74-
739 3456.0*t24*t38*t4*t57*t74+3456.0*t24*t4*t54*t66*t7*t73+8640.0*
740 t29*t5*t54*t7*t73-4824.0*t12*t30*t44*t5*t7+432.0*t13*t16*t5*
741 t7-864.0*t30*t32*t4*t44*t48*t7+2160.0*t1*t13*t16*t4*t48*t7+
742 864.0*t13*t16*t22*t24*t48*t7-864.0*t30*t32*t44*t46*t7+2160.0*
743 t1*t13*t16*t46*t7-6336.0*t24*t30*t4*t40*t44*t7+2808.0*t1*t33*
744 t45*t5+864.0*t12*t33*t4*t45*t48+864.0*t12*t33*t45*t46+5184.0*
745 t22*t24*t33*t4*t45)+4.0*t20*t21*t4*t77*grada-8.0*t10*t21*t76*
746 grada+15.0*t21*t4*t52*t72*grada-18.0*t21*t25*t28*t4*t49*grada+
747 12.0*t20*t21*t24*t49*grada+12.0*t20*t21*t25*t48*grada-18.0*
748 t21*t24*t28*t43*grada-1.0*t1*t10*t21*t83-6.5625*t1*t21*pow(t25,
749 4.0)*t4*t78-3.0*t1*t21*t25*t28*t4*t77+2.0*t1*t20*t21*t24*t77+
750 2.0*t1*t20*t21*t25*t76+7.5*t1*t21*t24*t52*t72+11.25*t1*t21*
751 t4*t43*t49*t52-2.25*t1*t21*t28*t4*pow(t49,2.0)+3.0*t1*t20*
752 t21*t48*t49+6.0*t20*t21*t4*t49-9.0*t1*t21*t24*t25*t28*t49-
753 4.5*t1*t21*t28*t43*t48-12.0*t10*t21*t48-9.0*t21*t28*t4*t43+
754 12.0*t20*t21*t24*t25;
755
756 }
757
758 static void
ft97bx_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)759 ft97bx_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
760 {
761 real res[14];
762
763 ft97bx_fourth_helper(dp->rhoa, dp->grada, res);
764
765 ds->df1000 += factor*res[0];
766 ds->df0010 += factor*res[1];
767
768 ds->df2000 += factor*res[2];
769 ds->df1010 += factor*res[3];
770 ds->df0020 += factor*res[4];
771
772 ds->df3000 += factor*res[5];
773 ds->df2010 += factor*res[6];
774 ds->df1020 += factor*res[7];
775 ds->df0030 += factor*res[8];
776
777 ds->df4000 += factor*res[9];
778 ds->df3010 += factor*res[10];
779 ds->df2020 += factor*res[11];
780 ds->df1030 += factor*res[12];
781 ds->df0040 += factor*res[13];
782
783
784 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
785 fabs(dp->grada-dp->gradb)>1e-13)
786 ft97bx_fourth_helper(dp->rhob, dp->gradb, res);
787
788 ds->df0100 += factor*res[0];
789 ds->df0001 += factor*res[1];
790
791 ds->df0200 += factor*res[2];
792 ds->df0101 += factor*res[3];
793 ds->df0002 += factor*res[4];
794
795 ds->df0300 += factor*res[5];
796 ds->df0201 += factor*res[6];
797 ds->df0102 += factor*res[7];
798 ds->df0003 += factor*res[8];
799
800 ds->df0400 += factor*res[9];
801 ds->df0301 += factor*res[10];
802 ds->df0202 += factor*res[11];
803 ds->df0103 += factor*res[12];
804 ds->df0004 += factor*res[13];
805
806 }
807