1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
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
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
30 /*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
31 /** @file fun-p86c.c P86c implementation.
32 
33    Automatically generated code implementing p86c functional and its
34    derivatives. Generated by func-codegen.pl being a part of a
35    "Automatic code generation framework for analytical functional
36     derivative evaluation", Pawel Salek, 2004
37 
38     This functional is connected by making following changes:
39     1. add "extern Functional p86cFunctional;" to 'functionals.h'
40     2. add "&p86cFunctional," to 'functionals.c'
41     3. add "fun-p86c}.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.
42 
43     This functional has been generated from following input:
44     ------ cut here -------
45 c0:    1.667e-3;
46 cn1:   2.568e-3;
47 cn2:   2.3266e-2;
48 cn3:   7.389e-6;
49 cd1:   1.0;
50 cd2:   8.723;
51 cd3:   0.472;
52 cd4:   1.0e4*cn3;
53 cinf:  c0+cn1/cd1;
54 phi1:  (9.0*%PI)^(1/6);
55 phi2:  0.11;
56 
57 rho:  rhoa + rhob;
58 grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
59 zeta: (rhoa-rhob)/(rhoa+rhob);
60 dzet: 2^(1/3)*sqrt(((1+zeta)/2)^(5/3) + ((1-zeta)/2)^(5/3));
61 
62 rs:    (3/(4*%PI*rho))^(1/3);
63 crho:  c0 + (cn1+cn2*rs+cn3*rs^2)/(1+cd2*rs+cd3*(rs^2)+cd4*(rs^3));
64 phi:   phi1*phi2*(cinf/crho)*grad*(rho^(-7.0/6.0));
65 
66 K(rhoa,rhob,grada,gradb,gradab):=exp(-phi)*crho*(grad^2)*(rho^(-4.0/3.0))/dzet;
67 
68     ------ cut here -------
69 */
70 
71 #include <math.h>
72 #include <stddef.h>
73 
74 #define __CVERSION__
75 
76 #include "functionals.h"
77 
78 /* INTERFACE PART */
p86c_isgga(void)79 static int p86c_isgga(void) { return 1; } /* FIXME: detect! */
80 static int p86c_read(const char *conf_line);
81 static real p86c_energy(const FunDensProp *dp);
82 static void p86c_first(FunFirstFuncDrv *ds,   real factor,
83                          const FunDensProp *dp);
84 static void p86c_second(FunSecondFuncDrv *ds, real factor,
85                           const FunDensProp *dp);
86 static void p86c_third(FunThirdFuncDrv *ds,   real factor,
87                          const FunDensProp *dp);
88 
89 Functional P86cFunctional = {
90   "P86c",       /* name */
91   p86c_isgga,   /* gga-corrected */
92   p86c_read,
93   NULL,
94   p86c_energy,
95   p86c_first,
96   p86c_second,
97   p86c_third
98 };
99 
100 /* IMPLEMENTATION PART */
101 static int
p86c_read(const char * conf_line)102 p86c_read(const char *conf_line)
103 {
104     fun_set_hf_weight(0);
105     return 1;
106 }
107 
108 static real
p86c_energy(const FunDensProp * dp)109 p86c_energy(const FunDensProp *dp)
110 {
111     real res;
112     real rhoa = dp->rhoa, rhob = dp->rhob;
113     real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
114 
115     real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
116     real t11, t12, t13, t14;
117 
118     t1 = POW(gradb,2.0)+2.0*gradab+POW(grada,2.0);
119     t2 = rhob+rhoa;
120     t3 = 1/POW(2.0,0.66666666666667);
121     t4 = rhoa-1.0*rhob;
122     t5 = 1/t2;
123     t6 = POW(3.0,0.66666666666667);
124     t7 = 1/POW(4.0,0.66666666666667);
125     t8 = 1/POW(3.141592653589793,0.66666666666667);
126     t9 = 1/POW(t2,0.66666666666667);
127     t10 = POW(3.0,0.33333333333333);
128     t11 = 1/POW(4.0,0.33333333333333);
129     t12 = 1/POW(3.141592653589793,0.33333333333333);
130     t13 = 1/POW(t2,0.33333333333333);
131     t14 = (0.023266*t10*t11*t12*t13+7.389E-6*t6*t7*t8*t9+
132         0.002568)/(0.472*t6*t7*t8*t9+0.01763993811759*t5+8.723000000000001*
133         t10*t11*t12*t13+1.0)+0.001667;
134 
135    /* code */
136     res = t1*t14/(POW(2.0,0.33333333333333)*POW(2.718281828459045,
137         6.718719623277062E-4*POW(3.141592653589793,0.16666666666667)*
138         SQRT(t1)/(t14*POW(t2,1.166666666666667)))*POW(t2,1.333333333333333)*
139         SQRT(0.5*t3*POW(t4*t5+1.0,1.666666666666667)+0.5*t3*POW(1.0-
140         1.0*t4*t5,1.666666666666667)));
141 
142     return res;
143 }
144 
145 static void
p86c_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)146 p86c_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
147 {
148     real dfdra, dfdrb, dfdga, dfdgb, dfdgab;
149     real rhoa = dp->rhoa, rhob = dp->rhob;
150     real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
151 
152     real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
153     real t11, t12, t13, t14, t15, t16, t17, t18;
154     real t19, t20, t21, t22, t23, t24, t25, t26;
155     real t27, t28, t29, t30, t31, t32, t33, t34;
156     real t35, t36, t37, t38, t39, t40, t41, t42;
157     real t43;
158 
159     t1 = 1/POW(2.0,0.33333333333333);
160     t2 = POW(gradb,2.0)+2.0*gradab+POW(grada,2.0);
161     t3 = rhob+rhoa;
162     t4 = 1/POW(t3,1.333333333333333);
163     t5 = 1/POW(2.0,0.66666666666667);
164     t6 = rhoa-1.0*rhob;
165     t7 = 1/t3;
166     t8 = 1.0-1.0*t6*t7;
167     t9 = t6*t7+1.0;
168     t10 = SQRT(0.5*t5*POW(t9,1.666666666666667)+0.5*t5*POW(t8,
169         1.666666666666667));
170     t11 = 1/t10;
171     t12 = 0.31830988618379;
172     t13 = 1/POW(t3,2.0);
173     t14 = POW(3.0,0.66666666666667);
174     t15 = 1/POW(4.0,0.66666666666667);
175     t16 = 1/POW(3.141592653589793,0.66666666666667);
176     t17 = 1/POW(t3,1.666666666666667);
177     t18 = POW(3.0,0.33333333333333);
178     t19 = 1/POW(4.0,0.33333333333333);
179     t20 = 1/POW(3.141592653589793,0.33333333333333);
180     t21 = 1/POW(t3,1.333333333333333);
181     t22 = 1/POW(t3,0.66666666666667);
182     t23 = 1/POW(t3,0.33333333333333);
183     t24 = 0.023266*t18*t19*t20*t23+7.389E-6*t14*t15*t16*t22+
184         0.002568;
185     t25 = 0.0554175*t12*t7+8.723000000000001*t18*t19*t20*
186         t23+0.472*t14*t15*t16*t22+1.0;
187     t26 = 1/t25;
188     t27 = (-0.00775533333333*t18*t19*t20*t21-4.926E-6*t14*
189         t15*t16*t17)*t26-1.0*(-2.907666666666667*t18*t19*t20*t21-0.31466666666667*
190         t14*t15*t16*t17-0.0554175*t12*t13)*t24/POW(t25,2.0);
191     t28 = POW(3.141592653589793,0.16666666666667);
192     t29 = SQRT(t2);
193     t30 = 1/POW(t3,1.166666666666667);
194     t31 = t24*t26+0.001667;
195     t32 = 1/t31;
196     t33 = 1/POW(2.718281828459045,6.718719623277062E-4*t28*
197         t29*t30*t32);
198     t34 = t1*t2*t4*t11*t27*t33;
199     t35 = t6*t13;
200     t36 = -1.0*t7;
201     t37 = POW(t8,0.66666666666667);
202     t38 = -1.0*t13*t6;
203     t39 = POW(t9,0.66666666666667);
204     t40 = 1/POW(t10,3.0);
205     t41 = -1.333333333333333*t1*t11*t2*t31*t33/POW(t3,2.333333333333333);
206     t42 = t1*
207         t11*t2*t31*(7.838506227156574E-4*t28*t29*t32/POW(t3,2.166666666666667)+
208         6.718719623277062E-4*t27*t28*t29*t30/POW(t31,2.0))*t33*t4;
209     t43 = 1/
210         POW(t3,2.5);
211 
212    /* code */
213     dfdra = -0.5*t1*t2*t31*t33*t4*t40*(0.83333333333333*t39*
214         t5*(t7+t38)+0.83333333333333*(t36+t35)*t37*t5)+t42+t41+t34;
215     dfdrb = -
216         0.5*t1*t2*t31*t33*t4*t40*(0.83333333333333*t37*t5*(t7+t35)+
217         0.83333333333333*(t36+t38)*t39*t5)+t42+t41+t34;
218     dfdga = 2.0*t1*t11*t31*t33*t4*grada-6.718719623277062E-4*
219         t1*t28*grada*t29*t43*t11*t33;
220     dfdgb = 2.0*t1*t11*t31*t33*t4*gradb-6.718719623277062E-4*
221         t1*t28*gradb*t29*t43*t11*t33;
222     dfdgab = 2.0*t1*t11*t31*t33*t4-6.718719623277062E-4*t1*
223         t28*t29*t43*t11*t33;
224 
225 
226    /* Final assignment */
227     ds->df1000 += factor*dfdra;
228     ds->df0100 += factor*dfdrb;
229     ds->df0010 += factor*dfdga;
230     ds->df0001 += factor*dfdgb;
231     ds->df00001+= factor*dfdgab;
232 }
233 
234 static void
p86c_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)235 p86c_second(FunSecondFuncDrv *ds, real factor, const FunDensProp *dp)
236 {
237     real dfdra, dfdrb, dfdga, dfdgb, dfdgab;
238     real d2fdrara, d2fdrarb, d2fdraga, d2fdragb, d2fdrbrb, d2fdraab,
239          d2fdgaga, d2fdgbgb, d2fdrbga, d2fdrbgb,
240          d2fdrbgab, d2fdgagb, d2fdgagab, d2fdgbgab, d2fdgabgab;
241     real rhoa = dp->rhoa, rhob = dp->rhob;
242     real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
243 
244     real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
245     real t11, t12, t13, t14, t15, t16, t17, t18;
246     real t19, t20, t21, t22, t23, t24, t25, t26;
247     real t27, t28, t29, t30, t31, t32, t33, t34;
248     real t35, t36, t37, t38, t39, t40, t41, t42;
249     real t43, t44, t45, t46, t47, t48, t49, t50;
250     real t51, t52, t53, t54, t55, t56, t57, t58;
251     real t59, t60, t61, t62, t63, t64, t65, t66;
252     real t67, t68, t69, t70, t71, t72, t73, t74;
253     real t75, t76, t77, t78, t79, t80, t81, t82;
254     real t83, t84, t85, t86, t87, t88, t89, t90;
255     real t91, t92, t93, t94, t95, t96, t97, t98;
256     real t99, t100, t101, t102, t103, t104, t105;
257 
258     t1 = 1/POW(2.0,0.33333333333333);
259     t2 = POW(grada,2.0);
260     t3 = POW(gradb,2.0);
261     t4 = 2.0*gradab+t3+t2;
262     t5 = rhob+rhoa;
263     t6 = 1/POW(t5,1.333333333333333);
264     t7 = 1/POW(2.0,0.66666666666667);
265     t8 = rhoa-1.0*rhob;
266     t9 = 1/t5;
267     t10 = 1.0-1.0*t8*t9;
268     t11 = t8*t9+1.0;
269     t12 = SQRT(0.5*POW(t11,1.666666666666667)*t7+0.5*POW(t10,
270         1.666666666666667)*t7);
271     t13 = 1/t12;
272     t14 = 0.31830988618379;
273     t15 = 1/POW(t5,2.0);
274     t16 = POW(3.0,0.66666666666667);
275     t17 = 1/POW(4.0,0.66666666666667);
276     t18 = 1/POW(3.141592653589793,0.66666666666667);
277     t19 = 1/POW(t5,1.666666666666667);
278     t20 = POW(3.0,0.33333333333333);
279     t21 = 1/POW(4.0,0.33333333333333);
280     t22 = POW(3.141592653589793,0.33333333333333);
281     t23 = 1/t22;
282     t24 = 1/POW(t5,1.333333333333333);
283     t25 = -2.907666666666667*t20*t21*t23*t24-0.31466666666667*
284         t16*t17*t18*t19-0.0554175*t14*t15;
285     t26 = 1/POW(t5,0.66666666666667);
286     t27 = 1/POW(t5,0.33333333333333);
287     t28 = 0.023266*t20*t21*t23*t27+7.389E-6*t16*t17*t18*t26+
288         0.002568;
289     t29 = 0.0554175*t14*t9+8.723000000000001*t20*t21*t23*
290         t27+0.472*t16*t17*t18*t26+1.0;
291     t30 = 1/POW(t29,2.0);
292     t31 = -0.00775533333333*t20*t21*t23*t24-4.926E-6*t16*
293         t17*t18*t19;
294     t32 = 1/t29;
295     t33 = t31*t32-1.0*t25*t28*t30;
296     t34 = POW(3.141592653589793,0.16666666666667);
297     t35 = SQRT(t4);
298     t36 = 1/POW(t5,1.166666666666667);
299     t37 = t28*t32+0.001667;
300     t38 = 1/t37;
301     t39 = 1/POW(2.718281828459045,6.718719623277062E-4*t34*
302         t35*t36*t38);
303     t40 = t1*t4*t6*t13*t33*t39;
304     t41 = t8*t15;
305     t42 = -1.0*t9;
306     t43 = t42+t41;
307     t44 = POW(t10,0.66666666666667);
308     t45 = -1.0*t15*t8;
309     t46 = t9+t45;
310     t47 = POW(t11,0.66666666666667);
311     t48 = 0.83333333333333*t46*t47*t7+0.83333333333333*t43*
312         t44*t7;
313     t49 = 1/POW(t12,3.0);
314     t50 = 1/POW(t5,2.333333333333333);
315     t51 = -1.333333333333333*t1*t4*t50*t13*t37*t39;
316     t52 = 1/POW(t37,2.0);
317     t53 = 1/POW(t5,2.166666666666667);
318     t54 = 7.838506227156574E-4*t34*t35*t53*t38+6.718719623277062E-4*
319         t34*t35*t36*t33*t52;
320     t55 = t1*t4*t6*t13*t37*t54*t39;
321     t56 = t9+t41;
322     t57 = t42+t45;
323     t58 = 0.83333333333333*t47*t57*t7+0.83333333333333*t44*
324         t56*t7;
325     t59 = 1/POW(t5,2.5);
326     t60 = -6.718719623277062E-4*t1*t34*t35*t59*t13*t39;
327     t61 = 2.0*
328         t1*t13*t37*t39*t6;
329     t62 = 1/POW(t5,3.0);
330     t63 = 1/POW(t5,2.666666666666667);
331     t64 = 1/POW(t5,2.333333333333334);
332     t65 = -1.0*t28*t30*(3.876888888888889*t20*t21*t23*t64+
333         0.52444444444444*t16*t17*t18*t63+0.110835*t14*t62)+(0.01034044444444*
334         t20*t21*t23*t64+8.210000000000001E-6*t16*t17*t18*t63)*t32-
335         2.0*t25*t30*t31+2.0*POW(t25,2.0)*t28/POW(t29,3.0);
336     t66 = t1*t4*t6*t13*t65*t39;
337     t67 = -2.666666666666667*t1*t4*t50*t13*t33*t39;
338     t68 = 1/POW(t12,5.0);
339     t69 = 1/POW(t10,0.33333333333333);
340     t70 = -2.0*t62*t8;
341     t71 = 2.0*t15;
342     t72 = 1/POW(t11,0.33333333333333);
343     t73 = 2.0*t62*t8;
344     t74 = -2.0*t15;
345     t75 = 3.111111111111111*t1*t13*t37*t39*t4/POW(t5,3.333333333333333);
346     t76 = t1*
347         t13*t37*t39*t4*(6.718719623277062E-4*t34*t35*t36*t65*t52-0.00156770124543*
348         t34*t35*t53*t33*t52-0.00169834301588*t34*t35*t38/POW(t5,3.166666666666667)-
349         0.00134374392466*POW(t33,2.0)*t34*t35*t36/POW(t37,3.0))*t6;
350     t77 = 2.0*
351         t1*t13*t33*t39*t4*t54*t6;
352     t78 = -2.666666666666667*t1*t4*t50*t13*t37*t54*t39;
353     t79 = t1*
354         t13*t37*t39*t4*POW(t54,2.0)*t6;
355     t80 = 1/POW(t5,3.5);
356     t81 = 8.958292831036083E-4*t1*t34*grada*t35*t80*t13*t39;
357     t82 = 2.0*
358         t1*t13*t33*t39*t6*grada;
359     t83 = -6.718719623277062E-4*t1*t34*grada*t35*t59*t13*
360         t33*t38*t39;
361     t84 = -2.666666666666667*t1*grada*t50*t13*t37*t39;
362     t85 = 1/t35;
363     t86 = t1*t4*t6*t13*t37*(7.838506227156574E-4*t34*grada*
364         t85*t53*t38+6.718719623277062E-4*t34*grada*t85*t36*t33*t52)*
365         t39;
366     t87 = -6.718719623277062E-4*t1*t34*grada*t35*t59*t13*
367         t54*t39;
368     t88 = 2.0*t1*t13*t37*t39*t54*t6*grada;
369     t89 = 8.958292831036083E-4*t1*t34*gradb*t35*t80*t13*t39;
370     t90 = 2.0*
371         t1*t13*t33*t39*t6*gradb;
372     t91 = -6.718719623277062E-4*t1*t34*gradb*t35*t59*t13*
373         t33*t38*t39;
374     t92 = -2.666666666666667*t1*gradb*t50*t13*t37*t39;
375     t93 = t1*t4*t6*t13*t37*(7.838506227156574E-4*t34*gradb*
376         t85*t53*t38+6.718719623277062E-4*t34*gradb*t85*t36*t33*t52)*
377         t39;
378     t94 = -6.718719623277062E-4*t1*t34*gradb*t35*t59*t13*
379         t54*t39;
380     t95 = 2.0*t1*t13*t37*t39*t54*t6*gradb;
381     t96 = 8.958292831036083E-4*t1*t34*t35*t80*t13*t39;
382     t97 = 2.0*t1*t13*t33*t39*t6;
383     t98 = -6.718719623277062E-4*t1*t34*t35*t59*t13*t33*t38*
384         t39;
385     t99 = -2.666666666666667*t1*t50*t13*t37*t39;
386     t100 = t1*t4*t6*t13*t37*(7.838506227156574E-4*t34*t85*
387         t53*t38+6.718719623277062E-4*t34*t85*t36*t33*t52)*t39;
388     t101 = -6.718719623277062E-4*t1*t34*t35*t59*t13*t54*t39;
389     t102 = 2.0*
390         t1*t13*t37*t39*t54*t6;
391     t103 = t99+t98+t97+t96-1.0*t1*t37*t39*t48*t49*t6+3.359359811638531E-4*
392         t1*t34*t35*t59*t48*t49*t39+t102+t101+t100;
393     t104 = t99+t98+t97+t96-1.0*t1*t37*t39*t49*t58*t6+3.359359811638531E-4*
394         t1*t34*t35*t59*t58*t49*t39+t102+t101+t100;
395     t105 = 1/POW(t5,3.666666666666667);
396 
397    /* code */
398     dfdra = -0.5*t1*t37*t39*t4*t48*t49*t6+t55+t51+t40;
399     dfdrb = -0.5*t1*t37*t39*t4*t49*t58*t6+t55+t51+t40;
400     dfdga = 2.0*t1*t13*t37*t39*t6*grada-6.718719623277062E-4*
401         t1*t34*grada*t35*t59*t13*t39;
402     dfdgb = 2.0*t1*t13*t37*t39*t6*gradb-6.718719623277062E-4*
403         t1*t34*gradb*t35*t59*t13*t39;
404     dfdgab = t61+t60;
405     d2fdrara = t79+t78+t77+t76+t75-0.5*t1*t37*t39*t4*t49*
406         t6*(0.83333333333333*t47*t7*(t74+t73)+0.55555555555556*POW(t46,
407         2.0)*t7*t72+0.83333333333333*t44*t7*(t71+t70)+0.55555555555556*
408         POW(t43,2.0)*t69*t7)+0.75*t1*t37*t39*t4*POW(t48,2.0)*t6*t68+
409         t67+t66-1.0*t1*t37*t39*t4*t48*t49*t54*t6-1.0*t1*t33*t39*t4*
410         t48*t49*t6+1.333333333333333*t1*t4*t50*t48*t49*t37*t39;
411     d2fdrarb = -0.5*t1*t37*t39*t4*t49*t6*(1.666666666666667*
412         t47*t62*t7*t8-1.666666666666667*t44*t62*t7*t8+0.55555555555556*
413         t46*t57*t7*t72+0.55555555555556*t43*t56*t69*t7)+t79+t78+t77+
414         t76+t75+0.75*t1*t37*t39*t4*t48*t58*t6*t68+t67+t66-0.5*t1*t37*
415         t39*t4*t49*t54*t58*t6-0.5*t1*t33*t39*t4*t49*t58*t6-0.5*t1*
416         t37*t39*t4*t48*t49*t54*t6-0.5*t1*t33*t39*t4*t48*t49*t6+0.66666666666667*
417         t1*t4*t50*t58*t49*t37*t39+0.66666666666667*t1*t4*t50*t48*t49*
418         t37*t39;
419     d2fdraga = -1.0*t1*t37*t39*t48*t49*t6*grada+t88+t87+t86+
420         t84+t83+t82+t81+3.359359811638531E-4*t1*t34*grada*t35*t59*
421         t48*t49*t39;
422     d2fdragb = -1.0*t1*t37*t39*t48*t49*t6*gradb+t95+t94+t93+
423         t92+t91+t90+t89+3.359359811638531E-4*t1*t34*gradb*t35*t59*
424         t48*t49*t39;
425     d2fdrbrb = t79+t78+t77+t76+t75-0.5*t1*t37*t39*t4*t49*
426         t6*(0.83333333333333*t44*t7*(t74+t70)+0.55555555555556*POW(t57,
427         2.0)*t7*t72+0.83333333333333*t47*t7*(t71+t73)+0.55555555555556*
428         POW(t56,2.0)*t69*t7)+0.75*t1*t37*t39*t4*POW(t58,2.0)*t6*t68+
429         t67+t66-1.0*t1*t37*t39*t4*t49*t54*t58*t6-1.0*t1*t33*t39*t4*
430         t49*t58*t6+1.333333333333333*t1*t4*t50*t58*t49*t37*t39;
431     d2fdraab = t103;
432     d2fdgaga = t61+4.514119337620827E-7*t1*t22*t2*t105*t13*
433         t38*t39+t60-0.00201561588698*t1*t34*t2*t85*t59*t13*t39;
434     d2fdgbgb = t61+4.514119337620827E-7*t1*t22*t3*t105*t13*
435         t38*t39+t60-0.00201561588698*t1*t34*t3*t85*t59*t13*t39;
436     d2fdrbga = -1.0*t1*t37*t39*t49*t58*t6*grada+t88+t87+t86+
437         t84+t83+t82+t81+3.359359811638531E-4*t1*t34*grada*t35*t59*
438         t58*t49*t39;
439     d2fdrbgb = -1.0*t1*t37*t39*t49*t58*t6*gradb+t95+t94+t93+
440         t92+t91+t90+t89+3.359359811638531E-4*t1*t34*gradb*t35*t59*
441         t58*t49*t39;
442     d2fdgagb = 4.514119337620827E-7*t1*t22*grada*gradb*t105*
443         t13*t38*t39-0.00201561588698*t1*t34*grada*gradb*t85*t59*t13*
444         t39;
445     d2fdrbgab = t104;
446     d2fdgagab = 4.514119337620827E-7*t1*t22*grada*t105*t13*
447         t38*t39-0.00201561588698*t1*t34*grada*t85*t59*t13*t39;
448     d2fdgbgab = 4.514119337620827E-7*t1*t22*gradb*t105*t13*
449         t38*t39-0.00201561588698*t1*t34*gradb*t85*t59*t13*t39;
450     d2fdgabgab = 4.514119337620827E-7*t1*t22*t105*t13*t38*
451         t39-0.00201561588698*t1*t34*t85*t59*t13*t39;
452 
453 
454     ds->df1000 += factor*dfdra;
455     ds->df0100 += factor*dfdrb;
456     ds->df0010 += factor*dfdga;
457     ds->df0001 += factor*dfdgb;
458     ds->df00001+= factor*dfdgab;
459 
460     ds->df2000 += factor*d2fdrara;
461     ds->df1100 += factor*d2fdrarb;
462     ds->df1010 += factor*d2fdraga;
463     ds->df1001 += factor*d2fdragb;
464     ds->df10001+= factor*d2fdraab;
465     ds->df0200 += factor*d2fdrbrb;
466     ds->df0110 += factor*d2fdrbga;
467     ds->df0101 += factor*d2fdrbgb;
468     ds->df01001+= factor*d2fdrbgab;
469     ds->df0020 += factor*d2fdgaga;
470     ds->df0011 += factor*d2fdgagb;
471     ds->df00101+= factor*d2fdgagab;
472     ds->df0002 += factor*d2fdgbgb;
473     ds->df00011+= factor*d2fdgbgab;
474     ds->df00002+= factor*d2fdgabgab;
475 }
476 
477 static void
p86c_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)478 p86c_third(FunThirdFuncDrv *ds, real factor, const FunDensProp *dp)
479 {
480     real dfdra, dfdrb, dfdga, dfdgb, dfdgab;
481     real d2fdrara, d2fdrarb, d2fdraga, d2fdragb, d2fdrbrb, d2fdraab,
482          d2fdrbab, d2fdgaga, d2fdgbgb, d2fdrbga, d2fdrbgb,
483          d2fdgagb;
484     real d3fdrararb, d3fdraraga, d3fdraragb, d3fdrbrbab, d3fdraraab,
485          d3fdrarbrb, d3fdrarbga, d3fdrarbgb, d3fdrarbab, d3fdragaga,
486          d3fdragbgb, d3fdrarara, d3fdrbrbrb, d3fdrbrbga, d3fdrbrbgb,
487          d3fdrbgaga, d3fdrbgbgb, d3fdgagaga,
488          d3fdragagb, d3fdrbgagb, d3fdgagagb, d3fdgagbgb, d3fdgbgbgb;
489     real rhoa = dp->rhoa, rhob = dp->rhob;
490     real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
491 
492     real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
493     real t11, t12, t13, t14, t15, t16, t17, t18;
494     real t19, t20, t21, t22, t23, t24, t25, t26;
495     real t27, t28, t29, t30, t31, t32, t33, t34;
496     real t35, t36, t37, t38, t39, t40, t41, t42;
497     real t43, t44, t45, t46, t47, t48, t49, t50;
498     real t51, t52, t53, t54, t55, t56, t57, t58;
499     real t59, t60, t61, t62, t63, t64, t65, t66;
500     real t67, t68, t69, t70, t71, t72, t73, t74;
501     real t75, t76, t77, t78, t79, t80, t81, t82;
502     real t83, t84, t85, t86, t87, t88, t89, t90;
503     real t91, t92, t93, t94, t95, t96, t97, t98;
504     real t99, t100, t101, t102, t103, t104, t105;
505     real t106, t107, t108, t109, t110, t111, t112;
506     real t113, t114, t115, t116, t117, t118, t119;
507     real t120, t121, t122, t123, t124, t125, t126;
508     real t127, t128, t129, t130, t131, t132, t133;
509     real t134, t135, t136, t137, t138, t139, t140;
510     real t141, t142, t143, t144, t145, t146, t147;
511     real t148, t149, t150, t151, t152, t153, t154;
512     real t155, t156, t157, t158, t159, t160, t161;
513     real t162, t163, t164, t165, t166, t167, t168;
514     real t169, t170, t171, t172, t173, t174, t175;
515     real t176, t177, t178, t179, t180, t181, t182;
516     real t183, t184, t185, t186, t187, t188, t189;
517     real t190, t191, t192, t193, t194, t195, t196;
518     real t197, t198, t199, t200, t201, t202, t203;
519     real t204, t205, t206, t207, t208, t209, t210;
520     real t211, t212, t213, t214, t215, t216, t217;
521     real t218, t219, t220, t221, t222, t223, t224;
522     real t225, t226, t227, t228, t229, t230, t231;
523     real t232, t233, t234, t235, t236, t237, t238;
524     real t239, t240, t241, t242, t243, t244, t245;
525     real t246, t247, t248, t249, t250, t251, t252;
526     real t253, t254, t255, t256, t257, t258, t259;
527     real t260, t261, t262, t263, t264;
528 
529     t1 = 1/POW(2.0,0.33333333333333);
530     t2 = POW(grada,2.0);
531     t3 = POW(gradb,2.0);
532     t4 = 2.0*gradab+t3+t2;
533     t5 = rhob+rhoa;
534     t6 = 1/POW(t5,1.333333333333333);
535     t7 = 1/POW(2.0,0.66666666666667);
536     t8 = rhoa-1.0*rhob;
537     t9 = 1/t5;
538     t10 = 1.0-1.0*t8*t9;
539     t11 = t8*t9+1.0;
540     t12 = SQRT(0.5*POW(t11,1.666666666666667)*t7+0.5*POW(t10,
541         1.666666666666667)*t7);
542     t13 = 1/t12;
543     t14 = 0.31830988618379;
544     t15 = 1/POW(t5,2.0);
545     t16 = POW(3.0,0.66666666666667);
546     t17 = 1/POW(4.0,0.66666666666667);
547     t18 = 1/POW(3.141592653589793,0.66666666666667);
548     t19 = 1/POW(t5,1.666666666666667);
549     t20 = POW(3.0,0.33333333333333);
550     t21 = 1/POW(4.0,0.33333333333333);
551     t22 = POW(3.141592653589793,0.33333333333333);
552     t23 = 1/t22;
553     t24 = 1/POW(t5,1.333333333333333);
554     t25 = -2.907666666666667*t20*t21*t23*t24-0.31466666666667*
555         t16*t17*t18*t19-0.0554175*t14*t15;
556     t26 = 1/POW(t5,0.66666666666667);
557     t27 = 1/POW(t5,0.33333333333333);
558     t28 = 0.023266*t20*t21*t23*t27+7.389E-6*t16*t17*t18*t26+
559         0.002568;
560     t29 = 0.0554175*t14*t9+8.723000000000001*t20*t21*t23*
561         t27+0.472*t16*t17*t18*t26+1.0;
562     t30 = 1/POW(t29,2.0);
563     t31 = -0.00775533333333*t20*t21*t23*t24-4.926E-6*t16*
564         t17*t18*t19;
565     t32 = 1/t29;
566     t33 = t31*t32-1.0*t25*t28*t30;
567     t34 = POW(3.141592653589793,0.16666666666667);
568     t35 = SQRT(t4);
569     t36 = 1/POW(t5,1.166666666666667);
570     t37 = t28*t32+0.001667;
571     t38 = 1/t37;
572     t39 = 1/POW(2.718281828459045,6.718719623277062E-4*t34*
573         t35*t36*t38);
574     t40 = t1*t4*t6*t13*t33*t39;
575     t41 = t8*t15;
576     t42 = -1.0*t9;
577     t43 = t42+t41;
578     t44 = POW(t10,0.66666666666667);
579     t45 = -1.0*t15*t8;
580     t46 = t9+t45;
581     t47 = POW(t11,0.66666666666667);
582     t48 = 0.83333333333333*t46*t47*t7+0.83333333333333*t43*
583         t44*t7;
584     t49 = 1/POW(t12,3.0);
585     t50 = 1/POW(t5,2.333333333333333);
586     t51 = -1.333333333333333*t1*t4*t50*t13*t37*t39;
587     t52 = 1/POW(t37,2.0);
588     t53 = 1/POW(t5,2.166666666666667);
589     t54 = 7.838506227156574E-4*t34*t35*t53*t38+6.718719623277062E-4*
590         t34*t35*t36*t33*t52;
591     t55 = t1*t4*t6*t13*t37*t54*t39;
592     t56 = t9+t41;
593     t57 = t42+t45;
594     t58 = 0.83333333333333*t47*t57*t7+0.83333333333333*t44*
595         t56*t7;
596     t59 = 1/POW(t5,2.5);
597     t60 = -6.718719623277062E-4*t1*t34*t35*t59*t13*t39;
598     t61 = 2.0*
599         t1*t13*t37*t39*t6;
600     t62 = POW(t25,2.0);
601     t63 = 1/POW(t29,3.0);
602     t64 = 1/POW(t5,3.0);
603     t65 = 1/POW(t5,2.666666666666667);
604     t66 = 1/POW(t5,2.333333333333334);
605     t67 = 3.876888888888889*t20*t21*t23*t66+0.52444444444444*
606         t16*t17*t18*t65+0.110835*t14*t64;
607     t68 = 0.01034044444444*t20*t21*t23*t66+8.210000000000001E-6*
608         t16*t17*t18*t65;
609     t69 = -1.0*t28*t30*t67+2.0*t28*t62*t63+t68*t32-2.0*t25*
610         t30*t31;
611     t70 = t1*t4*t6*t13*t69*t39;
612     t71 = -2.666666666666667*t1*t4*t50*t13*t33*t39;
613     t72 = POW(t48,2.0);
614     t73 = 1/POW(t12,5.0);
615     t74 = POW(t43,2.0);
616     t75 = 1/POW(t10,0.33333333333333);
617     t76 = -2.0*t64*t8;
618     t77 = 2.0*t15;
619     t78 = t77+t76;
620     t79 = POW(t46,2.0);
621     t80 = 1/POW(t11,0.33333333333333);
622     t81 = 2.0*t64*t8;
623     t82 = -2.0*t15;
624     t83 = t82+t81;
625     t84 = 0.83333333333333*t47*t7*t83+0.55555555555556*t7*
626         t79*t80+0.83333333333333*t44*t7*t78+0.55555555555556*t7*t74*
627         t75;
628     t85 = 1/POW(t5,3.333333333333333);
629     t86 = 3.111111111111111*t1*t4*t85*t13*t37*t39;
630     t87 = POW(t33,2.0);
631     t88 = 1/POW(t37,3.0);
632     t89 = 1/POW(t5,3.166666666666667);
633     t90 = -0.00169834301588*t34*t35*t89*t38-0.00156770124543*
634         t34*t35*t53*t33*t52+6.718719623277062E-4*t34*t35*t36*t69*t52-
635         0.00134374392466*t34*t35*t36*t87*t88;
636     t91 = t1*t4*t6*t13*t37*t90*t39;
637     t92 = 2.0*t1*t13*t33*t39*t4*t54*t6;
638     t93 = -2.666666666666667*t1*t4*t50*t13*t37*t54*t39;
639     t94 = POW(t54,
640         2.0);
641     t95 = t1*t4*t6*t13*t37*t94*t39;
642     t96 = 0.55555555555556*t46*t57*t7*t80+1.666666666666667*
643         t47*t64*t7*t8-1.666666666666667*t44*t64*t7*t8+0.55555555555556*
644         t43*t56*t7*t75;
645     t97 = 1/POW(t5,3.5);
646     t98 = 8.958292831036083E-4*t1*t34*grada*t35*t97*t13*t39;
647     t99 = 2.0*
648         t1*t13*t33*t39*t6*grada;
649     t100 = -6.718719623277062E-4*t1*t34*grada*t35*t59*t13*
650         t33*t38*t39;
651     t101 = -2.666666666666667*t1*grada*t50*t13*t37*t39;
652     t102 = 1/
653         t35;
654     t103 = 7.838506227156574E-4*t34*grada*t102*t53*t38+6.718719623277062E-4*
655         t34*grada*t102*t36*t33*t52;
656     t104 = t1*t4*t6*t13*t37*t103*t39;
657     t105 = -6.718719623277062E-4*t1*t34*grada*t35*t59*t13*
658         t54*t39;
659     t106 = 2.0*t1*t13*t37*t39*t54*t6*grada;
660     t107 = 8.958292831036083E-4*t1*t34*gradb*t35*t97*t13*
661         t39;
662     t108 = 2.0*t1*t13*t33*t39*t6*gradb;
663     t109 = -6.718719623277062E-4*t1*t34*gradb*t35*t59*t13*
664         t33*t38*t39;
665     t110 = -2.666666666666667*t1*gradb*t50*t13*t37*t39;
666     t111 = 7.838506227156574E-4*t34*gradb*t102*t53*t38+6.718719623277062E-4*
667         t34*gradb*t102*t36*t33*t52;
668     t112 = t1*t4*t6*t13*t37*t111*t39;
669     t113 = -6.718719623277062E-4*t1*t34*gradb*t35*t59*t13*
670         t54*t39;
671     t114 = 2.0*t1*t13*t37*t39*t54*t6*gradb;
672     t115 = POW(t58,2.0);
673     t116 = POW(t56,2.0);
674     t117 = t82+t76;
675     t118 = POW(t57,2.0);
676     t119 = t77+t81;
677     t120 = 0.55555555555556*t118*t7*t80+0.55555555555556*
678         t116*t7*t75+0.83333333333333*t119*t47*t7+0.83333333333333*
679         t117*t44*t7;
680     t121 = 3.359359811638531E-4*t1*t34*t35*t59*t48*t49*t39;
681     t122 = 8.958292831036083E-4*t1*t34*t35*t97*t13*t39;
682     t123 = 2.0*
683         t1*t13*t33*t39*t6;
684     t124 = -6.718719623277062E-4*t1*t34*t35*t59*t13*t33*t38*
685         t39;
686     t125 = -1.0*t1*t37*t39*t48*t49*t6;
687     t126 = -2.666666666666667*t1*t50*t13*t37*t39;
688     t127 = 6.718719623277062E-4*t34*t102*t36*t33*t52;
689     t128 = 7.838506227156574E-4*t34*t102*t53*t38;
690     t129 = t128+t127;
691     t130 = t1*t4*t6*t13*t37*t129*t39;
692     t131 = -6.718719623277062E-4*t1*t34*t35*t59*t13*t54*t39;
693     t132 = 2.0*
694         t1*t13*t37*t39*t54*t6;
695     t133 = t132+t131+t130+t126+t125+t124+t123+t122+t121;
696     t134 = 3.359359811638531E-4*t1*t34*t35*t59*t58*t49*t39;
697     t135 = -
698         1.0*t1*t37*t39*t49*t58*t6;
699     t136 = t132+t131+t130+t126+t135+t124+t123+t122+t134;
700     t137 = 1/
701         POW(t5,3.666666666666667);
702     t138 = -0.00201561588698*t1*t34*grada*t102*t59*t13*t39;
703     t139 = 4.514119337620827E-7*t1*t22*grada*t137*t13*t38*
704         t39;
705     t140 = -0.00201561588698*t1*t34*gradb*t102*t59*t13*t39;
706     t141 = 4.514119337620827E-7*t1*t22*gradb*t137*t13*t38*
707         t39;
708     t142 = 1/POW(t5,4.0);
709     t143 = 1/POW(t5,3.666666666666667);
710     t144 = 1/POW(t5,3.333333333333334);
711     t145 = -3.0*t25*t30*t68+6.0*t25*t28*t63*t67-3.0*t30*t31*
712         t67+6.0*t31*t62*t63+(-0.0241277037037*t20*t21*t23*t144-2.189333333333334E-5*
713         t16*t17*t18*t143)*t32-1.0*(-9.046074074074074*t20*t21*t23*
714         t144-1.398518518518519*t16*t17*t18*t143-0.332505*t14*t142)*
715         t28*t30-6.0*POW(t25,3.0)*t28/POW(t29,4.0);
716     t146 = t1*t4*t6*t13*t145*t39;
717     t147 = -4.0*t1*t4*t50*t13*t69*t39;
718     t148 = 1.5*t1*t33*t39*t4*t48*t58*t6*t73;
719     t149 = -1.0*t1*t33*t39*t4*t49*t6*t96;
720     t150 = 9.333333333333332*t1*t4*t85*t13*t33*t39;
721     t151 = 1/POW(t12,7.0);
722     t152 = -2.0*t1*t4*t50*t58*t48*t73*t37*t39;
723     t153 = 1.333333333333333*t1*t4*t50*t96*t49*t37*t39;
724     t154 = 1/
725         POW(t10,1.333333333333333);
726     t155 = 6.0*t142*t8;
727     t156 = 1/POW(t11,1.333333333333333);
728     t157 = -6.0*t142*t8;
729     t158 = -10.37037037037037*t1*t13*t37*t39*t4/POW(t5,4.333333333333333);
730     t159 = t1*
731         t13*t37*t39*t4*t6*(0.00470310373629*t34*t35*t53*t87*t88-0.00403123177397*
732         t34*t35*t36*t69*t33*t88-0.00235155186815*t34*t35*t53*t69*t52+
733         0.00509502904765*t34*t35*t89*t33*t52+6.718719623277062E-4*
734         t34*t35*t36*t145*t52+0.00537808621697*t34*t35*t38/POW(t5,4.166666666666667)+
735         0.00403123177397*POW(t33,3.0)*t34*t35*t36/POW(t37,4.0));
736     t160 = 3.0*
737         t1*t13*t33*t39*t4*t6*t90;
738     t161 = -4.0*t1*t4*t50*t13*t37*t90*t39;
739     t162 = 3.0*t1*t13*t39*t4*t54*t6*t69;
740     t163 = -8.0*t1*t4*t50*t13*t33*t54*t39;
741     t164 = 1.5*t1*t37*t39*t4*t48*t54*t58*t6*t73;
742     t165 = -1.0*t1*t37*t39*t4*t49*t54*t6*t96;
743     t166 = 9.333333333333332*t1*t4*t85*t13*t37*t54*t39;
744     t167 = 3.0*
745         t1*t13*t37*t39*t4*t54*t6*t90;
746     t168 = 3.0*t1*t13*t33*t39*t4*t6*t94;
747     t169 = -4.0*t1*t4*t50*t13*t37*t94*t39;
748     t170 = t1*t13*t37*t39*t4*POW(t54,3.0)*t6;
749     t171 = 1/POW(t5,4.5);
750     t172 = -0.00209026832724*t1*t34*grada*t35*t171*t13*t39;
751     t173 = 2.0*
752         t1*t13*t39*t6*t69*grada;
753     t174 = -5.333333333333333*t1*grada*t50*t13*t33*t39;
754     t175 = -
755         6.718719623277062E-4*t1*t34*grada*t35*t59*t13*t69*t38*t39;
756     t176 = 0.00179165856621*
757         t1*t34*grada*t35*t97*t13*t33*t38*t39;
758     t177 = 6.222222222222221*t1*grada*t85*t13*t37*t39;
759     t178 = t1*t4*t6*t13*t37*(-0.00169834301588*t34*grada*
760         t102*t89*t38-0.00156770124543*t34*grada*t102*t53*t33*t52+6.718719623277062E-4*
761         t34*grada*t102*t36*t69*t52-0.00134374392466*t34*grada*t102*
762         t36*t87*t88)*t39;
763     t179 = -6.718719623277062E-4*t1*t34*grada*t35*t59*t13*
764         t90*t39;
765     t180 = 2.0*t1*t13*t37*t39*t6*t90*grada;
766     t181 = 2.0*t1*t103*t13*t33*t39*t4*t6;
767     t182 = -2.666666666666667*t1*t4*t50*t13*t37*t103*t39;
768     t183 = 0.00179165856621*
769         t1*t34*grada*t35*t97*t13*t54*t39;
770     t184 = 4.0*t1*t13*t33*t39*t54*t6*grada;
771     t185 = -0.00134374392466*t1*t34*grada*t35*t59*t13*t33*
772         t38*t54*t39;
773     t186 = -5.333333333333333*t1*grada*t50*t13*t37*t54*t39;
774     t187 = 2.0*
775         t1*t103*t13*t37*t39*t4*t54*t6;
776     t188 = -6.718719623277062E-4*t1*t34*grada*t35*t59*t13*
777         t94*t39;
778     t189 = 2.0*t1*t13*t37*t39*t6*t94*grada;
779     t190 = -0.00209026832724*t1*t34*gradb*t35*t171*t13*t39;
780     t191 = 2.0*
781         t1*t13*t39*t6*t69*gradb;
782     t192 = -5.333333333333333*t1*gradb*t50*t13*t33*t39;
783     t193 = -
784         6.718719623277062E-4*t1*t34*gradb*t35*t59*t13*t69*t38*t39;
785     t194 = 0.00179165856621*
786         t1*t34*gradb*t35*t97*t13*t33*t38*t39;
787     t195 = 6.222222222222221*t1*gradb*t85*t13*t37*t39;
788     t196 = t1*t4*t6*t13*t37*(-0.00169834301588*t34*gradb*
789         t102*t89*t38-0.00156770124543*t34*gradb*t102*t53*t33*t52+6.718719623277062E-4*
790         t34*gradb*t102*t36*t69*t52-0.00134374392466*t34*gradb*t102*
791         t36*t87*t88)*t39;
792     t197 = -6.718719623277062E-4*t1*t34*gradb*t35*t59*t13*
793         t90*t39;
794     t198 = 2.0*t1*t13*t37*t39*t6*t90*gradb;
795     t199 = 2.0*t1*t111*t13*t33*t39*t4*t6;
796     t200 = -2.666666666666667*t1*t4*t50*t13*t37*t111*t39;
797     t201 = 0.00179165856621*
798         t1*t34*gradb*t35*t97*t13*t54*t39;
799     t202 = 4.0*t1*t13*t33*t39*t54*t6*gradb;
800     t203 = -0.00134374392466*t1*t34*gradb*t35*t59*t13*t33*
801         t38*t54*t39;
802     t204 = -5.333333333333333*t1*gradb*t50*t13*t37*t54*t39;
803     t205 = 2.0*
804         t1*t111*t13*t37*t39*t4*t54*t6;
805     t206 = -6.718719623277062E-4*t1*t34*gradb*t35*t59*t13*
806         t94*t39;
807     t207 = 2.0*t1*t13*t37*t39*t6*t94*gradb;
808     t208 = -0.00209026832724*t1*t34*t35*t171*t13*t39;
809     t209 = 2.0*t1*t13*t39*t6*t69;
810     t210 = -5.333333333333333*t1*t50*t13*t33*t39;
811     t211 = -6.718719623277062E-4*t1*t34*t35*t59*t13*t69*t38*
812         t39;
813     t212 = 0.00179165856621*t1*t34*t35*t97*t13*t33*t38*t39;
814     t213 = 6.222222222222221*
815         t1*t85*t13*t37*t39;
816     t214 = t1*t4*t6*t13*t37*(-0.00169834301588*t34*t102*t89*
817         t38-0.00156770124543*t34*t102*t53*t33*t52+6.718719623277062E-4*
818         t34*t102*t36*t69*t52-0.00134374392466*t34*t102*t36*t87*t88)*
819         t39;
820     t215 = -6.718719623277062E-4*t1*t34*t35*t59*t13*t90*t39;
821     t216 = 2.0*
822         t1*t13*t37*t39*t6*t90;
823     t217 = 2.0*t1*t129*t13*t33*t39*t4*t6;
824     t218 = -2.666666666666667*t1*t4*t50*t13*t37*t129*t39;
825     t219 = 0.00179165856621*
826         t1*t34*t35*t97*t13*t54*t39;
827     t220 = 4.0*t1*t13*t33*t39*t54*t6;
828     t221 = -0.00134374392466*t1*t34*t35*t59*t13*t33*t38*t54*
829         t39;
830     t222 = -5.333333333333333*t1*t50*t13*t37*t54*t39;
831     t223 = 2.0*t1*t129*t13*t37*t39*t4*t54*t6;
832     t224 = -6.718719623277062E-4*t1*t34*t35*t59*t13*t94*t39;
833     t225 = 2.0*
834         t1*t13*t37*t39*t6*t94;
835     t226 = 0.00268748784931*t1*t34*t2*t102*t97*t13*t39;
836     t227 = 4.514119337620827E-7*t1*t22*t2*t137*t13*t33*t52*
837         t39;
838     t228 = 1/POW(t5,4.666666666666667);
839     t229 = -6.018825783494436E-7*t1*t22*t2*t228*t13*t38*t39;
840     t230 = -
841         0.00201561588698*t1*t34*t2*t102*t59*t13*t33*t38*t39;
842     t231 = 1/POW(t35,3.0);
843     t232 = t1*t4*t6*t13*t37*(t128-7.838506227156574E-4*t34*
844         t2*t231*t53*t38+t127-6.718719623277062E-4*t34*t2*t231*t36*
845         t33*t52)*t39;
846     t233 = -0.00134374392466*t1*t34*grada*t35*t59*t13*t103*
847         t39;
848     t234 = 4.0*t1*t103*t13*t37*t39*t6*grada;
849     t235 = -0.00201561588698*t1*t34*t2*t102*t59*t13*t54*t39;
850     t236 = 4.514119337620827E-7*t1*t22*t2*t137*t13*t38*t54*
851         t39;
852     t237 = 0.00268748784931*t1*t34*grada*gradb*t102*t97*t13*
853         t39;
854     t238 = 4.514119337620827E-7*t1*t22*grada*gradb*t137*t13*
855         t33*t52*t39;
856     t239 = -6.018825783494436E-7*t1*t22*grada*gradb*t228*
857         t13*t38*t39;
858     t240 = -0.00201561588698*t1*t34*grada*gradb*t102*t59*
859         t13*t33*t38*t39;
860     t241 = t1*t4*t6*t13*t37*(-7.838506227156574E-4*t34*grada*
861         gradb*t231*t53*t38-6.718719623277062E-4*t34*grada*gradb*t231*
862         t36*t33*t52)*t39;
863     t242 = -6.718719623277062E-4*t1*t34*gradb*t35*t59*t13*
864         t103*t39;
865     t243 = 2.0*t1*t103*t13*t37*t39*t6*gradb;
866     t244 = -6.718719623277062E-4*t1*t34*grada*t35*t59*t13*
867         t111*t39;
868     t245 = 2.0*t1*t111*t13*t37*t39*t6*grada;
869     t246 = -0.00201561588698*t1*t34*grada*gradb*t102*t59*
870         t13*t54*t39;
871     t247 = 4.514119337620827E-7*t1*t22*grada*gradb*t137*t13*
872         t38*t54*t39;
873     t248 = 0.00268748784931*t1*t34*t3*t102*t97*t13*t39;
874     t249 = 4.514119337620827E-7*t1*t22*t3*t137*t13*t33*t52*
875         t39;
876     t250 = -6.018825783494436E-7*t1*t22*t3*t228*t13*t38*t39;
877     t251 = -
878         0.00201561588698*t1*t34*t3*t102*t59*t13*t33*t38*t39;
879     t252 = t1*t4*t6*t13*t37*(t128-7.838506227156574E-4*t34*
880         t3*t231*t53*t38+t127-6.718719623277062E-4*t34*t3*t231*t36*
881         t33*t52)*t39;
882     t253 = -0.00134374392466*t1*t34*gradb*t35*t59*t13*t111*
883         t39;
884     t254 = 4.0*t1*t111*t13*t37*t39*t6*gradb;
885     t255 = -0.00201561588698*t1*t34*t3*t102*t59*t13*t54*t39;
886     t256 = 4.514119337620827E-7*t1*t22*t3*t137*t13*t38*t54*
887         t39;
888     t257 = -6.0*t64;
889     t258 = 6.0*t64;
890     t259 = t132+t236+t131+t235+t234+t233+t232+t126+t135+t124+
891         t230+t229-2.257059668810414E-7*t1*t22*t2*t137*t58*t49*t38*
892         t39+t227+t123+t122+t226+t134+0.00100780794349*t1*t34*t2*t102*
893         t59*t58*t49*t39;
894     t260 = POW(grada,3.0);
895     t261 = 1.772453850905516;
896     t262 = 1/POW(t5,4.833333333333334);
897     t263 = 1/t4;
898     t264 = POW(gradb,3.0);
899 
900    /* code */
901     dfdra = -0.5*t1*t37*t39*t4*t48*t49*t6+t55+t51+t40;
902     dfdrb = -0.5*t1*t37*t39*t4*t49*t58*t6+t55+t51+t40;
903     dfdga = 2.0*t1*t13*t37*t39*t6*grada-6.718719623277062E-4*
904         t1*t34*grada*t35*t59*t13*t39;
905     dfdgb = 2.0*t1*t13*t37*t39*t6*gradb-6.718719623277062E-4*
906         t1*t34*gradb*t35*t59*t13*t39;
907     dfdgab = t61+t60;
908     d2fdrara = t95+t93+t92+t91+t86-0.5*t1*t37*t39*t4*t49*
909         t6*t84+0.75*t1*t37*t39*t4*t6*t72*t73+t71+t70-1.0*t1*t37*t39*
910         t4*t48*t49*t54*t6-1.0*t1*t33*t39*t4*t48*t49*t6+1.333333333333333*
911         t1*t4*t50*t48*t49*t37*t39;
912     d2fdrarb = -0.5*t1*t37*t39*t4*t49*t6*t96+t95+t93+t92+
913         t91+t86+0.75*t1*t37*t39*t4*t48*t58*t6*t73+t71+t70-0.5*t1*t37*
914         t39*t4*t49*t54*t58*t6-0.5*t1*t33*t39*t4*t49*t58*t6-0.5*t1*
915         t37*t39*t4*t48*t49*t54*t6-0.5*t1*t33*t39*t4*t48*t49*t6+0.66666666666667*
916         t1*t4*t50*t58*t49*t37*t39+0.66666666666667*t1*t4*t50*t48*t49*
917         t37*t39;
918     d2fdraga = -1.0*t1*t37*t39*t48*t49*t6*grada+t99+t98+3.359359811638531E-4*
919         t1*t34*grada*t35*t59*t48*t49*t39+t106+t105+t104+t101+t100;
920     d2fdragb = -
921         1.0*t1*t37*t39*t48*t49*t6*gradb+3.359359811638531E-4*t1*t34*
922         gradb*t35*t59*t48*t49*t39+t114+t113+t112+t110+t109+t108+t107;
923     d2fdrbrb = t95+
924         t93+t92+t91+t86+0.75*t1*t115*t37*t39*t4*t6*t73+t71+t70-1.0*
925         t1*t37*t39*t4*t49*t54*t58*t6-1.0*t1*t33*t39*t4*t49*t58*t6-
926         0.5*t1*t120*t37*t39*t4*t49*t6+1.333333333333333*t1*t4*t50*
927         t58*t49*t37*t39;
928     d2fdraab = t133;
929     d2fdrbab = t136;
930     d2fdgaga = t61+4.514119337620827E-7*t1*t22*t2*t137*t13*
931         t38*t39+t60-0.00201561588698*t1*t34*t2*t102*t59*t13*t39;
932     d2fdgbgb = t61+
933         4.514119337620827E-7*t1*t22*t3*t137*t13*t38*t39+t60-0.00201561588698*
934         t1*t34*t3*t102*t59*t13*t39;
935     d2fdrbga = -1.0*t1*t37*t39*t49*t58*t6*grada+t99+t98+3.359359811638531E-4*
936         t1*t34*grada*t35*t59*t58*t49*t39+t106+t105+t104+t101+t100;
937     d2fdrbgb = -
938         1.0*t1*t37*t39*t49*t58*t6*gradb+3.359359811638531E-4*t1*t34*
939         gradb*t35*t59*t58*t49*t39+t114+t113+t112+t110+t109+t108+t107;
940     d2fdgagb = 4.514119337620827E-7*t1*t22*grada*gradb*t137*
941         t13*t38*t39-0.00201561588698*t1*t34*grada*gradb*t102*t59*t13*
942         t39;
943     d3fdrararb = 1.5*t1*t37*t39*t4*t48*t6*t73*t96-0.5*t1*
944         t37*t39*t4*t49*t58*t6*t94-1.0*t1*t37*t39*t4*t48*t49*t6*t94-
945         0.5*t1*t37*t39*t4*t49*t58*t6*t90-1.0*t1*t37*t39*t4*t48*t49*
946         t6*t90+0.75*t1*t37*t39*t4*t58*t6*t73*t84-0.5*t1*t37*t39*t4*
947         t49*t54*t6*t84-0.5*t1*t33*t39*t4*t49*t6*t84-0.5*t1*t37*t39*
948         t4*t49*t6*(0.55555555555556*t57*t7*t80*t83+2.222222222222222*
949         t46*t64*t7*t8*t80-2.222222222222222*t43*t64*t7*t75*t8-0.18518518518519*
950         t156*t57*t7*t79+0.55555555555556*t56*t7*t75*t78-0.18518518518519*
951         t154*t56*t7*t74+0.83333333333333*t47*(2.0*t64+t157)*t7+0.83333333333333*
952         t44*(t155-2.0*t64)*t7)+0.75*t1*t37*t39*t4*t54*t6*t72*t73+0.75*
953         t1*t33*t39*t4*t6*t72*t73-1.875*t1*t151*t37*t39*t4*t58*t6*t72-
954         0.5*t1*t39*t4*t49*t58*t6*t69-1.0*t1*t39*t4*t48*t49*t6*t69-
955         1.0*t1*t33*t39*t4*t49*t54*t58*t6-2.0*t1*t33*t39*t4*t48*t49*
956         t54*t6+1.333333333333333*t1*t4*t50*t58*t49*t37*t54*t39+2.666666666666667*
957         t1*t4*t50*t48*t49*t37*t54*t39-1.0*t1*t4*t50*t72*t73*t37*t39+
958         0.66666666666667*t1*t4*t50*t84*t49*t37*t39-1.555555555555555*
959         t1*t4*t85*t58*t49*t37*t39-3.111111111111111*t1*t4*t85*t48*
960         t49*t37*t39+1.333333333333333*t1*t4*t50*t58*t49*t33*t39+2.666666666666667*
961         t1*t4*t50*t48*t49*t33*t39+t170+t169+t168+t167+t166+t165+t164+
962         t163+t162+t161+t160+t159+t158+t153+t152+t150+t149+t148+t147+
963         t146;
964     d3fdraraga = -1.0*t1*t37*t39*t49*t6*t84*grada+1.5*t1*
965         t37*t39*t6*t72*t73*grada-2.0*t1*t37*t39*t48*t49*t54*t6*grada-
966         2.0*t1*t33*t39*t48*t49*t6*grada-1.0*t1*t103*t37*t39*t4*t48*
967         t49*t6-5.039039717457796E-4*t1*t34*grada*t35*t59*t72*t73*t39+
968         6.718719623277062E-4*t1*t34*grada*t35*t59*t48*t49*t54*t39+
969         3.359359811638531E-4*t1*t34*grada*t35*t59*t84*t49*t39-8.958292831036083E-4*
970         t1*t34*grada*t35*t97*t48*t49*t39+6.718719623277062E-4*t1*t34*
971         grada*t35*t59*t48*t49*t33*t38*t39+2.666666666666667*t1*grada*
972         t50*t48*t49*t37*t39+t189+t188+t187+t186+t185+t184+t183+t182+
973         t181+t180+t179+t178+t177+t176+t175+t174+t173+t172;
974     d3fdraragb = -1.0*t1*t37*t39*t49*t6*t84*gradb+1.5*t1*
975         t37*t39*t6*t72*t73*gradb-2.0*t1*t37*t39*t48*t49*t54*t6*gradb-
976         2.0*t1*t33*t39*t48*t49*t6*gradb-1.0*t1*t111*t37*t39*t4*t48*
977         t49*t6-5.039039717457796E-4*t1*t34*gradb*t35*t59*t72*t73*t39+
978         6.718719623277062E-4*t1*t34*gradb*t35*t59*t48*t49*t54*t39+
979         3.359359811638531E-4*t1*t34*gradb*t35*t59*t84*t49*t39-8.958292831036083E-4*
980         t1*t34*gradb*t35*t97*t48*t49*t39+6.718719623277062E-4*t1*t34*
981         gradb*t35*t59*t48*t49*t33*t38*t39+2.666666666666667*t1*gradb*
982         t50*t48*t49*t37*t39+t207+t206+t205+t204+t203+t202+t201+t200+
983         t199+t198+t197+t196+t195+t194+t193+t192+t191+t190;
984     d3fdrbrbab = 1.5*t1*t115*t37*t39*t6*t73-2.0*t1*t37*t39*
985         t49*t54*t58*t6-1.0*t1*t129*t37*t39*t4*t49*t58*t6-2.0*t1*t33*
986         t39*t49*t58*t6-1.0*t1*t120*t37*t39*t49*t6-5.039039717457796E-4*
987         t1*t34*t35*t59*t115*t73*t39+6.718719623277062E-4*t1*t34*t35*
988         t59*t58*t49*t54*t39-8.958292831036083E-4*t1*t34*t35*t97*t58*
989         t49*t39+3.359359811638531E-4*t1*t34*t35*t59*t120*t49*t39+6.718719623277062E-4*
990         t1*t34*t35*t59*t58*t49*t33*t38*t39+2.666666666666667*t1*t50*
991         t58*t49*t37*t39+t225+t224+t223+t222+t221+t220+t219+t218+t217+
992         t216+t215+t214+t213+t212+t211+t210+t209+t208;
993     d3fdraraab = -1.0*t1*t37*t39*t49*t6*t84+1.5*t1*t37*t39*
994         t6*t72*t73-2.0*t1*t37*t39*t48*t49*t54*t6-1.0*t1*t129*t37*t39*
995         t4*t48*t49*t6-2.0*t1*t33*t39*t48*t49*t6-5.039039717457796E-4*
996         t1*t34*t35*t59*t72*t73*t39+6.718719623277062E-4*t1*t34*t35*
997         t59*t48*t49*t54*t39+3.359359811638531E-4*t1*t34*t35*t59*t84*
998         t49*t39-8.958292831036083E-4*t1*t34*t35*t97*t48*t49*t39+6.718719623277062E-4*
999         t1*t34*t35*t59*t48*t49*t33*t38*t39+2.666666666666667*t1*t50*
1000         t48*t49*t37*t39+t225+t224+t223+t222+t221+t220+t219+t218+t217+
1001         t216+t215+t214+t213+t212+t211+t210+t209+t208;
1002     d3fdrarbrb = 1.5*t1*t37*t39*t4*t58*t6*t73*t96-1.0*t1*
1003         t37*t39*t4*t49*t58*t6*t94-0.5*t1*t37*t39*t4*t48*t49*t6*t94-
1004         1.0*t1*t37*t39*t4*t49*t58*t6*t90-0.5*t1*t37*t39*t4*t48*t49*
1005         t6*t90-0.5*t1*t37*t39*t4*t49*t6*(2.222222222222222*t57*t64*
1006         t7*t8*t80+0.55555555555556*t119*t46*t7*t80-2.222222222222222*
1007         t56*t64*t7*t75*t8-5.0*t142*t47*t7*t8+5.0*t142*t44*t7*t8+0.55555555555556*
1008         t117*t43*t7*t75-1.666666666666667*t47*t64*t7+1.666666666666667*
1009         t44*t64*t7-0.18518518518519*t118*t156*t46*t7-0.18518518518519*
1010         t116*t154*t43*t7)+0.75*t1*t115*t37*t39*t4*t54*t6*t73+0.75*
1011         t1*t120*t37*t39*t4*t48*t6*t73+0.75*t1*t115*t33*t39*t4*t6*t73-
1012         1.0*t1*t39*t4*t49*t58*t6*t69-0.5*t1*t39*t4*t48*t49*t6*t69-
1013         2.0*t1*t33*t39*t4*t49*t54*t58*t6-1.0*t1*t33*t39*t4*t48*t49*
1014         t54*t6-0.5*t1*t120*t37*t39*t4*t49*t54*t6-0.5*t1*t120*t33*t39*
1015         t4*t49*t6-1.875*t1*t115*t151*t37*t39*t4*t48*t6+2.666666666666667*
1016         t1*t4*t50*t58*t49*t37*t54*t39+1.333333333333333*t1*t4*t50*
1017         t48*t49*t37*t54*t39-1.0*t1*t4*t50*t115*t73*t37*t39-3.111111111111111*
1018         t1*t4*t85*t58*t49*t37*t39-1.555555555555555*t1*t4*t85*t48*
1019         t49*t37*t39+0.66666666666667*t1*t4*t50*t120*t49*t37*t39+2.666666666666667*
1020         t1*t4*t50*t58*t49*t33*t39+1.333333333333333*t1*t4*t50*t48*
1021         t49*t33*t39+t170+t169+t168+t167+t166+t165+t164+t163+t162+t161+
1022         t160+t159+t158+t153+t152+t150+t149+t148+t147+t146;
1023     d3fdrarbga = -1.0*t1*t37*t39*t49*t6*t96*grada+1.5*t1*
1024         t37*t39*t48*t58*t6*t73*grada-1.0*t1*t37*t39*t49*t54*t58*t6*
1025         grada-1.0*t1*t33*t39*t49*t58*t6*grada-1.0*t1*t37*t39*t48*t49*
1026         t54*t6*grada-1.0*t1*t33*t39*t48*t49*t6*grada-0.5*t1*t103*t37*
1027         t39*t4*t49*t58*t6-0.5*t1*t103*t37*t39*t4*t48*t49*t6-5.039039717457796E-4*
1028         t1*t34*grada*t35*t59*t58*t48*t73*t39+3.359359811638531E-4*
1029         t1*t34*grada*t35*t59*t58*t49*t54*t39+3.359359811638531E-4*
1030         t1*t34*grada*t35*t59*t48*t49*t54*t39+3.359359811638531E-4*
1031         t1*t34*grada*t35*t59*t96*t49*t39-4.479146415518041E-4*t1*t34*
1032         grada*t35*t97*t58*t49*t39-4.479146415518041E-4*t1*t34*grada*
1033         t35*t97*t48*t49*t39+3.359359811638531E-4*t1*t34*grada*t35*
1034         t59*t58*t49*t33*t38*t39+3.359359811638531E-4*t1*t34*grada*
1035         t35*t59*t48*t49*t33*t38*t39+1.333333333333333*t1*grada*t50*
1036         t58*t49*t37*t39+1.333333333333333*t1*grada*t50*t48*t49*t37*
1037         t39+t189+t188+t187+t186+t185+t184+t183+t182+t181+t180+t179+
1038         t178+t177+t176+t175+t174+t173+t172;
1039     d3fdrarbgb = -1.0*t1*t37*t39*t49*t6*t96*gradb+1.5*t1*
1040         t37*t39*t48*t58*t6*t73*gradb-1.0*t1*t37*t39*t49*t54*t58*t6*
1041         gradb-1.0*t1*t33*t39*t49*t58*t6*gradb-1.0*t1*t37*t39*t48*t49*
1042         t54*t6*gradb-1.0*t1*t33*t39*t48*t49*t6*gradb-0.5*t1*t111*t37*
1043         t39*t4*t49*t58*t6-0.5*t1*t111*t37*t39*t4*t48*t49*t6-5.039039717457796E-4*
1044         t1*t34*gradb*t35*t59*t58*t48*t73*t39+3.359359811638531E-4*
1045         t1*t34*gradb*t35*t59*t58*t49*t54*t39+3.359359811638531E-4*
1046         t1*t34*gradb*t35*t59*t48*t49*t54*t39+3.359359811638531E-4*
1047         t1*t34*gradb*t35*t59*t96*t49*t39-4.479146415518041E-4*t1*t34*
1048         gradb*t35*t97*t58*t49*t39-4.479146415518041E-4*t1*t34*gradb*
1049         t35*t97*t48*t49*t39+3.359359811638531E-4*t1*t34*gradb*t35*
1050         t59*t58*t49*t33*t38*t39+3.359359811638531E-4*t1*t34*gradb*
1051         t35*t59*t48*t49*t33*t38*t39+1.333333333333333*t1*gradb*t50*
1052         t58*t49*t37*t39+1.333333333333333*t1*gradb*t50*t48*t49*t37*
1053         t39+t207+t206+t205+t204+t203+t202+t201+t200+t199+t198+t197+
1054         t196+t195+t194+t193+t192+t191+t190;
1055     d3fdrarbab = -1.0*t1*t37*t39*t49*t6*t96+1.5*t1*t37*t39*
1056         t48*t58*t6*t73-1.0*t1*t37*t39*t49*t54*t58*t6-0.5*t1*t129*t37*
1057         t39*t4*t49*t58*t6-1.0*t1*t33*t39*t49*t58*t6-1.0*t1*t37*t39*
1058         t48*t49*t54*t6-0.5*t1*t129*t37*t39*t4*t48*t49*t6-1.0*t1*t33*
1059         t39*t48*t49*t6-5.039039717457796E-4*t1*t34*t35*t59*t58*t48*
1060         t73*t39+3.359359811638531E-4*t1*t34*t35*t59*t58*t49*t54*t39+
1061         3.359359811638531E-4*t1*t34*t35*t59*t48*t49*t54*t39+3.359359811638531E-4*
1062         t1*t34*t35*t59*t96*t49*t39-4.479146415518041E-4*t1*t34*t35*
1063         t97*t58*t49*t39-4.479146415518041E-4*t1*t34*t35*t97*t48*t49*
1064         t39+3.359359811638531E-4*t1*t34*t35*t59*t58*t49*t33*t38*t39+
1065         3.359359811638531E-4*t1*t34*t35*t59*t48*t49*t33*t38*t39+1.333333333333333*
1066         t1*t50*t58*t49*t37*t39+1.333333333333333*t1*t50*t48*t49*t37*
1067         t39+t225+t224+t223+t222+t221+t220+t219+t218+t217+t216+t215+
1068         t214+t213+t212+t211+t210+t209+t208;
1069     d3fdragaga = t132+t236+t131+t235+t234+t233+t232+t126+
1070         t125+t124+t230+t229-2.257059668810414E-7*t1*t22*t2*t137*t48*
1071         t49*t38*t39+t227+t123+t122+t226+t121+0.00100780794349*t1*t34*
1072         t2*t102*t59*t48*t49*t39;
1073     d3fdragagb = t247+t246+t245+t244+t243+t242+t241+t240+
1074         t239-2.257059668810414E-7*t1*t22*grada*gradb*t137*t48*t49*
1075         t38*t39+t238+t237+0.00100780794349*t1*t34*grada*gradb*t102*
1076         t59*t48*t49*t39;
1077     d3fdrbgagb = t247+t246+t245+t244+t243+t242+t241+t240+
1078         t239-2.257059668810414E-7*t1*t22*grada*gradb*t137*t58*t49*
1079         t38*t39+t238+t237+0.00100780794349*t1*t34*grada*gradb*t102*
1080         t59*t58*t49*t39;
1081     d3fdragbgb = t132+t256+t131+t255+t254+t253+t252+t126+
1082         t125+t124+t251+t250-2.257059668810414E-7*t1*t22*t3*t137*t48*
1083         t49*t38*t39+t249+t123+t122+t248+t121+0.00100780794349*t1*t34*
1084         t3*t102*t59*t48*t49*t39;
1085     d3fdrarara = -1.5*t1*t37*t39*t4*t48*t49*t6*t94-1.5*t1*
1086         t37*t39*t4*t48*t49*t6*t90+2.25*t1*t37*t39*t4*t48*t6*t73*t84-
1087         1.5*t1*t37*t39*t4*t49*t54*t6*t84-1.5*t1*t33*t39*t4*t49*t6*
1088         t84-0.5*t1*t37*t39*t4*t49*t6*(1.666666666666667*t46*t7*t80*
1089         t83+1.666666666666667*t43*t7*t75*t78+0.83333333333333*(t258+
1090         t157)*t47*t7-0.18518518518519*t156*POW(t46,3.0)*t7+0.83333333333333*
1091         (t257+t155)*t44*t7-0.18518518518519*t154*POW(t43,3.0)*t7)+
1092         2.25*t1*t37*t39*t4*t54*t6*t72*t73+2.25*t1*t33*t39*t4*t6*t72*
1093         t73-1.5*t1*t39*t4*t48*t49*t6*t69-3.0*t1*t33*t39*t4*t48*t49*
1094         t54*t6-1.875*t1*t151*t37*t39*t4*POW(t48,3.0)*t6+4.0*t1*t4*
1095         t50*t48*t49*t37*t54*t39-3.0*t1*t4*t50*t72*t73*t37*t39+2.0*
1096         t1*t4*t50*t84*t49*t37*t39-4.666666666666666*t1*t4*t85*t48*
1097         t49*t37*t39+4.0*t1*t4*t50*t48*t49*t33*t39+t170+t169+t168+t167+
1098         t166+t163+t162+t161+t160+t159+t158+t150+t147+t146;
1099     d3fdrbrbrb = -1.5*t1*t37*t39*t4*t49*t58*t6*t94-1.5*t1*
1100         t37*t39*t4*t49*t58*t6*t90-0.5*t1*t37*t39*t4*t49*t6*(1.666666666666667*
1101         t119*t57*t7*t80+1.666666666666667*t117*t56*t7*t75-0.18518518518519*
1102         t156*POW(t57,3.0)*t7-0.18518518518519*t154*POW(t56,3.0)*t7+
1103         0.83333333333333*(t257+t157)*t47*t7+0.83333333333333*(t258+
1104         t155)*t44*t7)+2.25*t1*t120*t37*t39*t4*t58*t6*t73+2.25*t1*t115*
1105         t37*t39*t4*t54*t6*t73+2.25*t1*t115*t33*t39*t4*t6*t73-1.5*t1*
1106         t39*t4*t49*t58*t6*t69-1.875*t1*t151*t37*t39*t4*POW(t58,3.0)*
1107         t6-3.0*t1*t33*t39*t4*t49*t54*t58*t6-1.5*t1*t120*t37*t39*t4*
1108         t49*t54*t6-1.5*t1*t120*t33*t39*t4*t49*t6+4.0*t1*t4*t50*t58*
1109         t49*t37*t54*t39-3.0*t1*t4*t50*t115*t73*t37*t39-4.666666666666666*
1110         t1*t4*t85*t58*t49*t37*t39+2.0*t1*t4*t50*t120*t49*t37*t39+4.0*
1111         t1*t4*t50*t58*t49*t33*t39+t170+t169+t168+t167+t166+t163+t162+
1112         t161+t160+t159+t158+t150+t147+t146;
1113     d3fdrbrbga = 1.5*t1*t115*t37*t39*t6*t73*grada-2.0*t1*
1114         t37*t39*t49*t54*t58*t6*grada-2.0*t1*t33*t39*t49*t58*t6*grada-
1115         1.0*t1*t120*t37*t39*t49*t6*grada-1.0*t1*t103*t37*t39*t4*t49*
1116         t58*t6-5.039039717457796E-4*t1*t34*grada*t35*t59*t115*t73*
1117         t39+6.718719623277062E-4*t1*t34*grada*t35*t59*t58*t49*t54*
1118         t39-8.958292831036083E-4*t1*t34*grada*t35*t97*t58*t49*t39+
1119         3.359359811638531E-4*t1*t34*grada*t35*t59*t120*t49*t39+6.718719623277062E-4*
1120         t1*t34*grada*t35*t59*t58*t49*t33*t38*t39+2.666666666666667*
1121         t1*grada*t50*t58*t49*t37*t39+t189+t188+t187+t186+t185+t184+
1122         t183+t182+t181+t180+t179+t178+t177+t176+t175+t174+t173+t172;
1123     d3fdrbrbgb = 1.5*
1124         t1*t115*t37*t39*t6*t73*gradb-2.0*t1*t37*t39*t49*t54*t58*t6*
1125         gradb-2.0*t1*t33*t39*t49*t58*t6*gradb-1.0*t1*t120*t37*t39*
1126         t49*t6*gradb-1.0*t1*t111*t37*t39*t4*t49*t58*t6-5.039039717457796E-4*
1127         t1*t34*gradb*t35*t59*t115*t73*t39+6.718719623277062E-4*t1*
1128         t34*gradb*t35*t59*t58*t49*t54*t39-8.958292831036083E-4*t1*
1129         t34*gradb*t35*t97*t58*t49*t39+3.359359811638531E-4*t1*t34*
1130         gradb*t35*t59*t120*t49*t39+6.718719623277062E-4*t1*t34*gradb*
1131         t35*t59*t58*t49*t33*t38*t39+2.666666666666667*t1*gradb*t50*
1132         t58*t49*t37*t39+t207+t206+t205+t204+t203+t202+t201+t200+t199+
1133         t198+t197+t196+t195+t194+t193+t192+t191+t190;
1134     d3fdrbgaga = t259;
1135     d3fdrbgbgb = t132+t256+t131+t255+t254+t253+t252+t126+
1136         t135+t124+t251+t250-2.257059668810414E-7*t1*t22*t3*t137*t58*
1137         t49*t38*t39+t249+t123+t122+t248+t134+0.00100780794349*t1*t34*
1138         t3*t102*t59*t58*t49*t39;
1139     d3fdgagaga = 1.354235801286248E-6*t1*t22*t260*t263*t137*
1140         t13*t38*t39+1.354235801286248E-6*t1*t22*grada*t137*t13*t38*
1141         t39-3.03291021754875E-10*t1*t261*t260*t102*t262*t13*t52*t39-
1142         0.00604684766095*t1*t34*grada*t102*t59*t13*t39+0.00201561588698*
1143         t1*t34*t260*t231*t59*t13*t39;
1144     d3fdgagagb = 1.354235801286248E-6*t1*t22*t2*gradb*t263*
1145         t137*t13*t38*t39+t141-3.03291021754875E-10*t1*t261*t2*gradb*
1146         t102*t262*t13*t52*t39+t140+0.00201561588698*t1*t34*t2*gradb*
1147         t231*t59*t13*t39;
1148     d3fdgagbgb = 1.354235801286248E-6*t1*t22*grada*t3*t263*
1149         t137*t13*t38*t39+t139-3.03291021754875E-10*t1*t261*grada*t3*
1150         t102*t262*t13*t52*t39+t138+0.00201561588698*t1*t34*grada*t3*
1151         t231*t59*t13*t39;
1152     d3fdgbgbgb = 1.354235801286248E-6*t1*t22*t264*t263*t137*
1153         t13*t38*t39+1.354235801286248E-6*t1*t22*gradb*t137*t13*t38*
1154         t39-3.03291021754875E-10*t1*t261*t264*t102*t262*t13*t52*t39-
1155         0.00604684766095*t1*t34*gradb*t102*t59*t13*t39+0.00201561588698*
1156         t1*t34*t264*t231*t59*t13*t39;
1157 
1158 
1159     ds->df1000 += factor*dfdra;
1160     ds->df0100 += factor*dfdrb;
1161     ds->df0010 += factor*dfdga;
1162     ds->df0001 += factor*dfdgb;
1163     ds->df00001+= factor*dfdgab;
1164 
1165     ds->df2000 += factor*d2fdrara;
1166     ds->df1100 += factor*d2fdrarb;
1167     ds->df1010 += factor*d2fdraga;
1168     ds->df1001 += factor*d2fdragb;
1169     ds->df0200 += factor*d2fdrbrb;
1170     ds->df10001+= factor*d2fdraab;
1171     ds->df01001+= factor*d2fdrbab;
1172     ds->df0020 += factor*d2fdgaga;
1173     ds->df0011 += factor*d2fdgagb;
1174     ds->df0002 += factor*d2fdgbgb;
1175     ds->df0110 += factor*d2fdrbga;
1176     ds->df0101 += factor*d2fdrbgb;
1177 
1178     ds->df2010 += factor*d3fdraraga;
1179     ds->df2001 += factor*d3fdraragb;
1180     ds->df1101 += factor*d3fdrarbgb;
1181     ds->df11001 += factor*d3fdrarbab;
1182     ds->df1020 += factor*d3fdragaga;
1183     ds->df1011 += factor*d3fdragagb;
1184     ds->df0111 += factor*d3fdrbgagb;
1185     ds->df1002 += factor*d3fdragbgb;
1186     ds->df3000 += factor*d3fdrarara;
1187     ds->df2100 += factor*d3fdrararb;
1188     ds->df20001 += factor*d3fdraraab;
1189     ds->df02001 += factor*d3fdrbrbab;
1190     ds->df1200 += factor*d3fdrarbrb;
1191     ds->df1110 += factor*d3fdrarbga;
1192     ds->df0300 += factor*d3fdrbrbrb;
1193     ds->df0210 += factor*d3fdrbrbga;
1194     ds->df0201 += factor*d3fdrbrbgb;
1195     ds->df0120 += factor*d3fdrbgaga;
1196     ds->df0102 += factor*d3fdrbgbgb;
1197     ds->df0030 += factor*d3fdgagaga;
1198     ds->df0021 += factor*d3fdgagagb;
1199     ds->df0012 += factor*d3fdgagbgb;
1200     ds->df0003 += factor*d3fdgbgbgb;
1201 }
1202