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 /* Partially automatically generated PZ81 functional.
25    Reference: J.P. Perdew and A. Zunger, Phys. Rev. B, 23, 5048 (1981).
26    Implemented and tested by: Pawel Salek.
27 */
28 
29 #if !defined(SYS_DEC)
30 /* XOPEN compliance is missing on old Tru64 4.0E Alphas */
31 #define _XOPEN_SOURCE          500
32 #define _XOPEN_SOURCE_EXTENDED 1
33 #endif
34 #include <math.h>
35 #include <stddef.h>
36 #include "general.h"
37 
38 #define __CVERSION__
39 
40 #include "functionals.h"
41 #define LOG log
42 #define ABS fabs
43 #define ASINH asinh
44 #define SQRT sqrt
45 
46 /* INTERFACE PART */
47 static integer pz81_read(const char* conf_line);
48 static real pz81_energy(const FunDensProp* dp);
49 static void pz81_first (FunFirstFuncDrv *ds,  real factor, const FunDensProp* dp);
50 static void pz81_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp);
51 static void pz81_third (FunThirdFuncDrv *ds,  real factor, const FunDensProp* dp);
52 
53 Functional PZ81Functional = {
54   "PZ81",
55   fun_false,
56   1,
57   pz81_read,
58   NULL,
59   pz81_energy,
60   pz81_first,
61   pz81_second,
62   pz81_third
63 };
64 
65 /* IMPLEMENTATION PART */
66 
67 /* the constants, as defined in the article */
68 static const real Au = 0.0311,  Bu = -0.048;
69 static const real Ap = 0.01555, Bp = -0.0269;
70 
71 /* These follow table XII, Ceperley-Alder parameters used */
72 /* correct values */
73 static const real gu = -0.1423, b1u = 1.0529, b2u = 0.3334,
74     Cu = 0.0020, Du = -0.0116;
75 static const real gp = -0.0843, b1p = 1.3981, b2p = 0.2611,
76     Cp = 0.0007, Dp = -0.0048;
77 
78 
79 static integer
pz81_read(const char * conf_line)80 pz81_read(const char* conf_line)
81 {
82     fun_set_hf_weight(0);
83     return 1;
84 }
85 
86 
87 /* ******************************************************************* */
88 /*                 Low density (rs>=1) part                            */
89 /* ******************************************************************* */
90 
91 
92 static real
pz81a_energy(const FunDensProp * dp)93 pz81a_energy(const FunDensProp* dp)
94 {
95     real t[7],zk;
96     real rhoa = dp->rhoa;
97     real rhob = dp->rhob;
98 
99     t[1] = rhob+rhoa;
100     t[2] = 1/pow(t[1],0.33333333333333);
101     t[3] = 1/pow(t[1],0.16666666666667);
102     t[4] = 1/(0.78762331789974*b1u*t[3]+0.6203504908994*b2u*t[2]+1.0);
103     t[5] = rhoa-1.0*rhob;
104     t[6] = 1/t[1];
105     zk = t[1]*((pow(t[5]*t[6]+1.0,1.333333333333333)+pow(1.0-1.0*t[5]*t[6],1.333333333333333)-2.0)*(gp/(0.78762331789974*b1p*t[3]+0.6203504908994*b2p*t[2]+1.0)-1.0*t[4]*gu)/(2.0*pow(2.0,0.33333333333333)-2.0)+gu*t[4]);
106     return zk;
107 }
108 
109 static void
pz81a_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)110 pz81a_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp)
111 {
112     real t[28];
113     real dfdra, dfdrb;
114     real rhoa = dp->rhoa;
115     real rhob = dp->rhob;
116 
117     t[1] = rhob+rhoa;
118     t[2] = 1/pow(t[1],0.33333333333333);
119     t[3] = 1/pow(t[1],0.16666666666667);
120     t[4] = 0.78762331789974*b1u*t[3]+0.6203504908994*b2u*t[2]+1.0;
121     t[5] = 1/t[4];
122     t[6] = gu*t[5];
123     t[7] = 1/(2.0*pow(2.0,0.33333333333333)-2.0);
124     t[8] = rhoa-1.0*rhob;
125     t[9] = 1/t[1];
126     t[10] = 1.0-1.0*t[8]*t[9];
127     t[11] = t[8]*t[9]+1.0;
128     t[12] = pow(t[11],1.333333333333333)+pow(t[10],1.333333333333333)-2.0;
129     t[13] = 0.78762331789974*b1p*t[3]+0.6203504908994*b2p*t[2]+1.0;
130     t[14] = gp/t[13]-1.0*t[5]*gu;
131     t[15] = t[7]*t[12]*t[14];
132     t[16] = 1/pow(t[1],1.333333333333333);
133     t[17] = 1/pow(t[1],1.166666666666667);
134     t[18] = -0.13127055298329*b1u*t[17]-0.20678349696647*b2u*t[16];
135     t[19] = 1/pow(t[4],2.0);
136     t[20] = -1.0*t[18]*t[19]*gu;
137     t[21] = t[12]*t[7]*(gu*t[18]*t[19]-1.0*(-0.13127055298329*b1p*t[17]-0.20678349696647*b2p*t[16])*gp/pow(t[13],2.0));
138     t[22] = 1/pow(t[1],2.0);
139     t[23] = t[8]*t[22];
140     t[24] = -1.0*t[9];
141     t[25] = pow(t[10],0.33333333333333);
142     t[26] = -1.0*t[22]*t[8];
143     t[27] = pow(t[11],0.33333333333333);
144     dfdra = t[1]*(t[14]*t[7]*(1.333333333333333*t[27]*(t[9]+t[26])+1.333333333333333*(t[24]+t[23])*t[25])+t[21]+t[20])+t[6]+t[15];
145     dfdrb = t[1]*(t[14]*t[7]*(1.333333333333333*t[25]*(t[9]+t[23])+1.333333333333333*(t[24]+t[26])*t[27])+t[21]+t[20])+t[6]+t[15];
146     ds->df1000 += factor*dfdra;
147     ds->df0100 += factor*dfdrb;
148 }
149 
150 static void
pz81a_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)151 pz81a_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
152 {
153     real t[56];
154     real dfdra, dfdrb;
155     real d2fdrara, d2fdrarb, d2fdrbrb;
156     real rhoa = dp->rhoa;
157     real rhob = dp->rhob;
158 
159     t[1] = rhob+rhoa;
160     t[2] = 1/pow(t[1],0.33333333333333);
161     t[3] = 1/pow(t[1],0.16666666666667);
162     t[4] = 0.78762331789974*b1u*t[3]+0.6203504908994*b2u*t[2]+1.0;
163     t[5] = 1/t[4];
164     t[6] = gu*t[5];
165     t[7] = 1/(2.0*pow(2.0,0.33333333333333)-2.0);
166     t[8] = rhoa-1.0*rhob;
167     t[9] = 1/t[1];
168     t[10] = 1.0-1.0*t[8]*t[9];
169     t[11] = t[8]*t[9]+1.0;
170     t[12] = pow(t[11],1.333333333333333)+pow(t[10],1.333333333333333)-2.0;
171     t[13] = 0.78762331789974*b1p*t[3]+0.6203504908994*b2p*t[2]+1.0;
172     t[14] = gp/t[13]-1.0*t[5]*gu;
173     t[15] = t[7]*t[12]*t[14];
174     t[16] = 1/pow(t[1],1.333333333333333);
175     t[17] = 1/pow(t[1],1.166666666666667);
176     t[18] = -0.13127055298329*b1u*t[17]-0.20678349696647*b2u*t[16];
177     t[19] = 1/pow(t[4],2.0);
178     t[20] = -1.0*t[18]*t[19]*gu;
179     t[21] = -0.13127055298329*b1p*t[17]-0.20678349696647*b2p*t[16];
180     t[22] = 1/pow(t[13],2.0);
181     t[23] = gu*t[18]*t[19]-1.0*t[21]*t[22]*gp;
182     t[24] = t[7]*t[12]*t[23];
183     t[25] = 1/pow(t[1],2.0);
184     t[26] = t[8]*t[25];
185     t[27] = -1.0*t[9];
186     t[28] = t[27]+t[26];
187     t[29] = pow(t[10],0.33333333333333);
188     t[30] = -1.0*t[25]*t[8];
189     t[31] = t[9]+t[30];
190     t[32] = pow(t[11],0.33333333333333);
191     t[33] = 1.333333333333333*t[31]*t[32]+1.333333333333333*t[28]*t[29];
192     t[34] = t[7]*t[33]*t[14];
193     t[35] = t[9]+t[26];
194     t[36] = t[27]+t[30];
195     t[37] = 1.333333333333333*t[32]*t[36]+1.333333333333333*t[29]*t[35];
196     t[38] = t[7]*t[37]*t[14];
197     t[39] = -2.0*t[18]*t[19]*gu;
198     t[40] = 2.0*t[12]*t[23]*t[7];
199     t[41] = pow(t[18],2.0);
200     t[42] = 1/pow(t[4],3.0);
201     t[43] = 2.0*t[41]*t[42]*gu;
202     t[44] = 1/pow(t[1],2.333333333333334);
203     t[45] = 1/pow(t[1],2.166666666666667);
204     t[46] = 0.15314897848051*b1u*t[45]+0.27571132928862*b2u*t[44];
205     t[47] = -1.0*t[19]*t[46]*gu;
206     t[48] = t[12]*t[7]*(-2.0*t[41]*t[42]*gu-1.0*t[22]*(0.15314897848051*b1p*t[45]+0.27571132928862*b2p*t[44])*gp+2.0*pow(t[21],2.0)*gp/pow(t[13],3.0)+gu*t[46]*t[19]);
207     t[49] = 1/pow(t[10],0.66666666666667);
208     t[50] = 1/pow(t[1],3.0);
209     t[51] = -2.0*t[50]*t[8];
210     t[52] = 2.0*t[25];
211     t[53] = 1/pow(t[11],0.66666666666667);
212     t[54] = 2.0*t[50]*t[8];
213     t[55] = -2.0*t[25];
214     dfdra = t[1]*(t[34]+t[24]+t[20])+t[15]+t[6];
215     dfdrb = t[1]*(t[38]+t[24]+t[20])+t[15]+t[6];
216     d2fdrara = t[1]*(t[14]*(1.333333333333333*t[32]*(t[55]+t[54])+0.44444444444444*pow(t[31],2.0)*t[53]+1.333333333333333*t[29]*(t[52]+t[51])+0.44444444444444*pow(t[28],2.0)*t[49])*t[7]+2.0*t[23]*t[33]*t[7]+t[48]+t[47]+t[43])+2.0*t[14]*t[33]*t[7]+t[40]+t[39];
217     d2fdrarb = t[1]*(t[14]*t[7]*(2.666666666666667*t[32]*t[50]*t[8]-2.666666666666667*t[29]*t[50]*t[8]+0.44444444444444*t[31]*t[36]*t[53]+0.44444444444444*t[28]*t[35]*t[49])+t[48]+t[47]+t[43]+t[7]*t[37]*t[23]+t[7]*t[33]*t[23])+t[40]+t[39]+t[38]+t[34];
218     d2fdrbrb = t[1]*(t[14]*(1.333333333333333*t[29]*(t[55]+t[51])+0.44444444444444*pow(t[36],2.0)*t[53]+1.333333333333333*t[32]*(t[52]+t[54])+0.44444444444444*pow(t[35],2.0)*t[49])*t[7]+2.0*t[23]*t[37]*t[7]+t[48]+t[47]+t[43])+2.0*t[14]*t[37]*t[7]+t[40]+t[39];
219     ds->df1000 += factor*dfdra;
220     ds->df0100 += factor*dfdrb;
221     ds->df2000 += factor*d2fdrara;
222     ds->df1100 += factor*d2fdrarb;
223     ds->df0200 += factor*d2fdrbrb;
224 }
225 
226 static void
pz81a_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)227 pz81a_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
228 {
229     real t[95];
230     real dfdra, dfdrb;
231     real d2fdrara, d2fdrarb, d2fdrbrb;
232     real d3fdrarara, d3fdrararb, d3fdrarbrb, d3fdrbrbrb;
233     real rhoa = dp->rhoa;
234     real rhob = dp->rhob;
235 
236     t[1] = rhob+rhoa;
237     t[2] = 1/pow(t[1],0.33333333333333);
238     t[3] = 1/pow(t[1],0.16666666666667);
239     t[4] = 0.78762331789974*b1u*t[3]+0.6203504908994*b2u*t[2]+1.0;
240     t[5] = 1/t[4];
241     t[6] = gu*t[5];
242     t[7] = 1/(2.0*pow(2.0,0.33333333333333)-2.0);
243     t[8] = rhoa-1.0*rhob;
244     t[9] = 1/t[1];
245     t[10] = 1.0-1.0*t[8]*t[9];
246     t[11] = t[8]*t[9]+1.0;
247     t[12] = pow(t[11],1.333333333333333)+pow(t[10],1.333333333333333)-2.0;
248     t[13] = 0.78762331789974*b1p*t[3]+0.6203504908994*b2p*t[2]+1.0;
249     t[14] = gp/t[13]-1.0*t[5]*gu;
250     t[15] = t[7]*t[12]*t[14];
251     t[16] = 1/pow(t[1],1.333333333333333);
252     t[17] = 1/pow(t[1],1.166666666666667);
253     t[18] = -0.13127055298329*b1u*t[17]-0.20678349696647*b2u*t[16];
254     t[19] = 1/pow(t[4],2.0);
255     t[20] = -1.0*t[18]*t[19]*gu;
256     t[21] = -0.13127055298329*b1p*t[17]-0.20678349696647*b2p*t[16];
257     t[22] = 1/pow(t[13],2.0);
258     t[23] = gu*t[18]*t[19]-1.0*t[21]*t[22]*gp;
259     t[24] = t[7]*t[12]*t[23];
260     t[25] = 1/pow(t[1],2.0);
261     t[26] = t[8]*t[25];
262     t[27] = -1.0*t[9];
263     t[28] = t[27]+t[26];
264     t[29] = pow(t[10],0.33333333333333);
265     t[30] = -1.0*t[25]*t[8];
266     t[31] = t[9]+t[30];
267     t[32] = pow(t[11],0.33333333333333);
268     t[33] = 1.333333333333333*t[31]*t[32]+1.333333333333333*t[28]*t[29];
269     t[34] = t[7]*t[33]*t[14];
270     t[35] = t[9]+t[26];
271     t[36] = t[27]+t[30];
272     t[37] = 1.333333333333333*t[32]*t[36]+1.333333333333333*t[29]*t[35];
273     t[38] = t[7]*t[37]*t[14];
274     t[39] = -2.0*t[18]*t[19]*gu;
275     t[40] = 2.0*t[12]*t[23]*t[7];
276     t[41] = pow(t[18],2.0);
277     t[42] = 1/pow(t[4],3.0);
278     t[43] = 2.0*t[41]*t[42]*gu;
279     t[44] = 1/pow(t[1],2.333333333333334);
280     t[45] = 1/pow(t[1],2.166666666666667);
281     t[46] = 0.15314897848051*b1u*t[45]+0.27571132928862*b2u*t[44];
282     t[47] = -1.0*t[19]*t[46]*gu;
283     t[48] = 1/pow(t[13],3.0);
284     t[49] = 0.15314897848051*b1p*t[45]+0.27571132928862*b2p*t[44];
285     t[50] = -2.0*t[41]*t[42]*gu-1.0*t[22]*t[49]*gp+2.0*pow(t[21],2.0)*t[48]*gp+gu*t[46]*t[19];
286     t[51] = t[7]*t[12]*t[50];
287     t[52] = 2.0*t[23]*t[33]*t[7];
288     t[53] = pow(t[28],2.0);
289     t[54] = 1/pow(t[10],0.66666666666667);
290     t[55] = 1/pow(t[1],3.0);
291     t[56] = -2.0*t[55]*t[8];
292     t[57] = 2.0*t[25];
293     t[58] = t[57]+t[56];
294     t[59] = pow(t[31],2.0);
295     t[60] = 1/pow(t[11],0.66666666666667);
296     t[61] = 2.0*t[55]*t[8];
297     t[62] = -2.0*t[25];
298     t[63] = t[62]+t[61];
299     t[64] = 1.333333333333333*t[32]*t[63]+0.44444444444444*t[59]*t[60]+1.333333333333333*t[29]*t[58]+0.44444444444444*t[53]*t[54];
300     t[65] = t[7]*t[64]*t[14];
301     t[66] = 2.666666666666667*t[32]*t[55]*t[8]-2.666666666666667*t[29]*t[55]*t[8]+0.44444444444444*t[31]*t[36]*t[60]+0.44444444444444*t[28]*t[35]*t[54];
302     t[67] = 2.0*t[23]*t[37]*t[7];
303     t[68] = pow(t[35],2.0);
304     t[69] = t[62]+t[56];
305     t[70] = pow(t[36],2.0);
306     t[71] = t[57]+t[61];
307     t[72] = 1.333333333333333*t[32]*t[71]+0.44444444444444*t[60]*t[70]+1.333333333333333*t[29]*t[69]+0.44444444444444*t[54]*t[68];
308     t[73] = t[7]*t[72]*t[14];
309     t[74] = 6.0*t[41]*t[42]*gu;
310     t[75] = -3.0*t[19]*t[46]*gu;
311     t[76] = 3.0*t[12]*t[50]*t[7];
312     t[77] = 2.0*t[14]*t[66]*t[7];
313     t[78] = pow(t[18],3.0);
314     t[79] = 1/pow(t[4],4.0);
315     t[80] = -6.0*t[78]*t[79]*gu;
316     t[81] = 6.0*t[18]*t[42]*t[46]*gu;
317     t[82] = 1/pow(t[1],3.333333333333334);
318     t[83] = 1/pow(t[1],3.166666666666667);
319     t[84] = -0.33182278670776*b1u*t[83]-0.64332643500679*b2u*t[82];
320     t[85] = -1.0*t[19]*t[84]*gu;
321     t[86] = t[12]*t[7]*(6.0*t[78]*t[79]*gu-6.0*t[18]*t[42]*t[46]*gu-1.0*t[22]*(-0.33182278670776*b1p*t[83]-0.64332643500679*b2p*t[82])*gp+6.0*t[21]*t[48]*t[49]*gp-6.0*pow(t[21],3.0)*gp/pow(t[13],4.0)+gu*t[84]*t[19]);
322     t[87] = 2.0*t[23]*t[66]*t[7];
323     t[88] = 1/pow(t[10],1.666666666666667);
324     t[89] = 1/pow(t[1],4.0);
325     t[90] = 6.0*t[8]*t[89];
326     t[91] = 1/pow(t[11],1.666666666666667);
327     t[92] = -6.0*t[8]*t[89];
328     t[93] = -6.0*t[55];
329     t[94] = 6.0*t[55];
330     dfdra = t[1]*(t[34]+t[24]+t[20])+t[15]+t[6];
331     dfdrb = t[1]*(t[38]+t[24]+t[20])+t[15]+t[6];
332     d2fdrara = 2.0*t[14]*t[33]*t[7]+t[1]*(t[65]+t[52]+t[51]+t[47]+t[43])+t[40]+t[39];
333     d2fdrarb = t[1]*(t[7]*t[66]*t[14]+t[7]*t[33]*t[23]+t[7]*t[37]*t[23]+t[51]+t[47]+t[43])+t[34]+t[38]+t[40]+t[39];
334     d2fdrbrb = t[1]*(t[73]+t[67]+t[51]+t[47]+t[43])+2.0*t[14]*t[37]*t[7]+t[40]+t[39];
335     d3fdrararb = t[1]*(t[14]*t[7]*(1.333333333333333*t[32]*(t[92]+2.0*t[55])-0.2962962962963*t[36]*t[59]*t[91]+1.333333333333333*t[29]*(t[90]-2.0*t[55])-0.2962962962963*t[35]*t[53]*t[88]+1.777777777777778*t[31]*t[55]*t[60]*t[8]-1.777777777777778*t[28]*t[54]*t[55]*t[8]+0.44444444444444*t[36]*t[60]*t[63]+0.44444444444444*t[35]*t[54]*t[58])+t[87]+t[86]+t[85]+t[81]+t[80]+2.0*t[33]*t[50]*t[7]+t[7]*t[37]*t[50]+t[7]*t[64]*t[23])+t[77]+t[76]+t[75]+t[74]+4.0*t[23]*t[33]*t[7]+t[67]+t[65];
336     d3fdrarbrb = t[1]*(t[14]*t[7]*(-0.2962962962963*t[31]*t[70]*t[91]-8.0*t[32]*t[8]*t[89]+8.0*t[29]*t[8]*t[89]-0.2962962962963*t[28]*t[68]*t[88]+1.777777777777778*t[36]*t[55]*t[60]*t[8]-1.777777777777778*t[35]*t[54]*t[55]*t[8]+0.44444444444444*t[31]*t[60]*t[71]+0.44444444444444*t[28]*t[54]*t[69]-2.666666666666667*t[32]*t[55]+2.666666666666667*t[29]*t[55])+t[87]+t[86]+t[85]+t[81]+t[80]+2.0*t[37]*t[50]*t[7]+t[7]*t[33]*t[50]+t[7]*t[72]*t[23])+t[77]+t[76]+t[75]+t[74]+t[73]+4.0*t[23]*t[37]*t[7]+t[52];
337     d3fdrarara = t[1]*(t[14]*t[7]*(1.333333333333333*t[32]*(t[94]+t[92])+1.333333333333333*t[29]*(t[93]+t[90])-0.2962962962963*pow(t[31],3.0)*t[91]-0.2962962962963*pow(t[28],3.0)*t[88]+1.333333333333333*t[31]*t[60]*t[63]+1.333333333333333*t[28]*t[54]*t[58])+t[86]+t[85]+t[81]+t[80]+3.0*t[23]*t[64]*t[7]+3.0*t[33]*t[50]*t[7])+t[76]+t[75]+t[74]+3.0*t[14]*t[64]*t[7]+6.0*t[23]*t[33]*t[7];
338     d3fdrbrbrb = t[1]*(t[14]*t[7]*(1.333333333333333*t[29]*(t[94]+t[90])+1.333333333333333*t[32]*(t[93]+t[92])-0.2962962962963*pow(t[36],3.0)*t[91]-0.2962962962963*pow(t[35],3.0)*t[88]+1.333333333333333*t[36]*t[60]*t[71]+1.333333333333333*t[35]*t[54]*t[69])+t[86]+t[85]+t[81]+t[80]+3.0*t[23]*t[7]*t[72]+3.0*t[37]*t[50]*t[7])+t[76]+t[75]+t[74]+3.0*t[14]*t[7]*t[72]+6.0*t[23]*t[37]*t[7];
339     ds->df1000 += factor*dfdra;
340     ds->df0100 += factor*dfdrb;
341     ds->df2000 += factor*d2fdrara;
342     ds->df1100 += factor*d2fdrarb;
343     ds->df0200 += factor*d2fdrbrb;
344     ds->df3000 += factor*d3fdrarara;
345     ds->df2100 += factor*d3fdrararb;
346     ds->df1200 += factor*d3fdrarbrb;
347     ds->df0300 += factor*d3fdrbrbrb;
348 }
349 
350 /* ******************************************************************* */
351 /*                 High density (rs<1) part                            */
352 /* ******************************************************************* */
353 
354 
355 
356 static real
pz81b_energy(const FunDensProp * dp)357 pz81b_energy(const FunDensProp* dp)
358 {
359     real t[6],zk;
360     real rhoa = dp->rhoa;
361     real rhob = dp->rhob;
362 
363     t[1] = rhob+rhoa;
364     t[2] = 1/pow(t[1],0.33333333333333);
365     t[3] = LOG(0.6203504908994*t[2]);
366     t[4] = rhoa-1.0*rhob;
367     t[5] = 1/t[1];
368     zk = t[1]*(Bu+(pow(t[4]*t[5]+1.0,1.333333333333333)+pow(1.0-1.0*t[4]*t[5],1.333333333333333)-2.0)*(-1.0*Bu+Bp-1.0*t[3]*Au+Ap*t[3]-0.6203504908994*Cu*t[2]*t[3]+0.6203504908994*Cp*t[2]*t[3]-0.6203504908994*Du*t[2]+0.6203504908994*Dp*t[2])/(2.0*pow(2.0,0.33333333333333)-2.0)+Au*t[3]+0.6203504908994*Cu*t[2]*t[3]+0.6203504908994*Du*t[2]);
369     return zk;
370 }
371 
372 static void
pz81b_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)373 pz81b_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp)
374 {
375     real t[27];
376     real dfdra, dfdrb;
377     real rhoa = dp->rhoa;
378     real rhob = dp->rhob;
379 
380     t[1] = rhob+rhoa;
381     t[2] = 1/pow(t[1],0.33333333333333);
382     t[3] = 0.6203504908994*Du*t[2];
383     t[4] = LOG(0.6203504908994*t[2]);
384     t[5] = Au*t[4];
385     t[6] = 0.6203504908994*Cu*t[2]*t[4];
386     t[7] = 1/(2.0*pow(2.0,0.33333333333333)-2.0);
387     t[8] = rhoa-1.0*rhob;
388     t[9] = 1/t[1];
389     t[10] = 1.0-1.0*t[8]*t[9];
390     t[11] = t[8]*t[9]+1.0;
391     t[12] = pow(t[11],1.333333333333333)+pow(t[10],1.333333333333333)-2.0;
392     t[13] = -1.0*Bu+Bp-1.0*t[4]*Au+Ap*t[4]-0.6203504908994*Cu*t[2]*t[4]+0.6203504908994*Cp*t[2]*t[4]-0.6203504908994*Du*t[2]+0.6203504908994*Dp*t[2];
393     t[14] = t[7]*t[12]*t[13];
394     t[15] = 1/pow(t[1],1.333333333333333);
395     t[16] = -0.20678349696647*Cu*t[15];
396     t[17] = -0.20678349696647*Du*t[15];
397     t[18] = -0.33333333333333*Au*t[9];
398     t[19] = -0.20678349696647*Cu*t[15]*t[4];
399     t[20] = t[7]*t[12]*(0.20678349696647*Cu*t[15]*t[4]-0.20678349696647*Cp*t[15]*t[4]+0.33333333333333*Au*t[9]-0.33333333333333*Ap*t[9]+0.20678349696647*Du*t[15]-0.20678349696647*Dp*t[15]+0.20678349696647*Cu*t[15]-0.20678349696647*Cp*t[15]);
400     t[21] = 1/pow(t[1],2.0);
401     t[22] = t[8]*t[21];
402     t[23] = -1.0*t[9];
403     t[24] = pow(t[10],0.33333333333333);
404     t[25] = -1.0*t[21]*t[8];
405     t[26] = pow(t[11],0.33333333333333);
406     dfdra = Bu+t[1]*(t[13]*t[7]*(1.333333333333333*t[26]*(t[9]+t[25])+1.333333333333333*(t[23]+t[22])*t[24])+t[20]+t[19]+t[18]+t[17]+t[16])+t[6]+t[5]+t[3]+t[14];
407     dfdrb = Bu+t[1]*(t[13]*t[7]*(1.333333333333333*t[24]*(t[9]+t[22])+1.333333333333333*(t[23]+t[25])*t[26])+t[20]+t[19]+t[18]+t[17]+t[16])+t[6]+t[5]+t[3]+t[14];
408     ds->df1000 += factor*dfdra;
409     ds->df0100 += factor*dfdrb;
410 }
411 
412 static void
pz81b_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)413 pz81b_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
414 {
415     real t[54];
416     real dfdra, dfdrb;
417     real d2fdrara, d2fdrarb, d2fdrbrb;
418     real rhoa = dp->rhoa;
419     real rhob = dp->rhob;
420 
421     t[1] = rhob+rhoa;
422     t[2] = 1/pow(t[1],0.33333333333333);
423     t[3] = 0.6203504908994*Du*t[2];
424     t[4] = LOG(0.6203504908994*t[2]);
425     t[5] = Au*t[4];
426     t[6] = 0.6203504908994*Cu*t[2]*t[4];
427     t[7] = 1/(2.0*pow(2.0,0.33333333333333)-2.0);
428     t[8] = rhoa-1.0*rhob;
429     t[9] = 1/t[1];
430     t[10] = 1.0-1.0*t[8]*t[9];
431     t[11] = t[8]*t[9]+1.0;
432     t[12] = pow(t[11],1.333333333333333)+pow(t[10],1.333333333333333)-2.0;
433     t[13] = -1.0*Bu+Bp-1.0*t[4]*Au+Ap*t[4]-0.6203504908994*Cu*t[2]*t[4]+0.6203504908994*Cp*t[2]*t[4]-0.6203504908994*Du*t[2]+0.6203504908994*Dp*t[2];
434     t[14] = t[7]*t[12]*t[13];
435     t[15] = 1/pow(t[1],1.333333333333333);
436     t[16] = -0.20678349696647*Cu*t[15];
437     t[17] = -0.20678349696647*Du*t[15];
438     t[18] = -0.33333333333333*Au*t[9];
439     t[19] = -0.20678349696647*Cu*t[15]*t[4];
440     t[20] = 0.20678349696647*Cu*t[15]*t[4]-0.20678349696647*Cp*t[15]*t[4]+0.33333333333333*Au*t[9]-0.33333333333333*Ap*t[9]+0.20678349696647*Du*t[15]-0.20678349696647*Dp*t[15]+0.20678349696647*Cu*t[15]-0.20678349696647*Cp*t[15];
441     t[21] = t[7]*t[12]*t[20];
442     t[22] = 1/pow(t[1],2.0);
443     t[23] = t[8]*t[22];
444     t[24] = -1.0*t[9];
445     t[25] = t[24]+t[23];
446     t[26] = pow(t[10],0.33333333333333);
447     t[27] = -1.0*t[22]*t[8];
448     t[28] = t[9]+t[27];
449     t[29] = pow(t[11],0.33333333333333);
450     t[30] = 1.333333333333333*t[28]*t[29]+1.333333333333333*t[25]*t[26];
451     t[31] = t[7]*t[30]*t[13];
452     t[32] = t[9]+t[23];
453     t[33] = t[24]+t[27];
454     t[34] = 1.333333333333333*t[29]*t[33]+1.333333333333333*t[26]*t[32];
455     t[35] = t[7]*t[34]*t[13];
456     t[36] = -0.41356699393293*Cu*t[15];
457     t[37] = -0.41356699393293*Du*t[15];
458     t[38] = -0.66666666666667*Au*t[9];
459     t[39] = -0.41356699393293*Cu*t[15]*t[4];
460     t[40] = 2.0*t[12]*t[20]*t[7];
461     t[41] = 1/pow(t[1],2.333333333333334);
462     t[42] = 0.34463916161078*Cu*t[41];
463     t[43] = 0.27571132928862*Du*t[41];
464     t[44] = 0.33333333333333*Au*t[22];
465     t[45] = 0.27571132928862*Cu*t[41]*t[4];
466     t[46] = t[7]*t[12]*(-0.27571132928862*Cu*t[41]*t[4]+0.27571132928862*Cp*t[41]*t[4]-0.33333333333333*Au*t[22]+0.33333333333333*Ap*t[22]-0.27571132928862*Du*t[41]+0.27571132928862*Dp*t[41]-0.34463916161078*Cu*t[41]+0.34463916161078*Cp*t[41]);
467     t[47] = 1/pow(t[10],0.66666666666667);
468     t[48] = 1/pow(t[1],3.0);
469     t[49] = -2.0*t[48]*t[8];
470     t[50] = 2.0*t[22];
471     t[51] = 1/pow(t[11],0.66666666666667);
472     t[52] = 2.0*t[48]*t[8];
473     t[53] = -2.0*t[22];
474     dfdra = t[1]*(t[31]+t[21]+t[19]+t[18]+t[17]+t[16])+t[14]+t[6]+t[5]+t[3]+Bu;
475     dfdrb = t[1]*(t[35]+t[21]+t[19]+t[18]+t[17]+t[16])+t[14]+t[6]+t[5]+t[3]+Bu;
476     d2fdrara =t[1]*(t[13]*(1.333333333333333*t[29]*(t[53]+t[52])+0.44444444444444*pow(t[28],2.0)*t[51]+1.333333333333333*t[26]*(t[50]+t[49])+0.44444444444444*pow(t[25],2.0)*t[47])*t[7]+2.0*t[20]*t[30]*t[7]+t[46]+t[45]+t[44]+t[43]+t[42])+2.0*t[13]*t[30]*t[7]+t[40]+t[39]+t[38]+t[37]+t[36];
477     d2fdrarb = t[1]*(t[13]*t[7]*(2.666666666666667*t[29]*t[48]*t[8]-2.666666666666667*t[26]*t[48]*t[8]+0.44444444444444*t[28]*t[33]*t[51]+0.44444444444444*t[25]*t[32]*t[47])+t[46]+t[45]+t[44]+t[43]+t[42]+t[7]*t[34]*t[20]+t[7]*t[30]*t[20])+t[40]+t[39]+t[38]+t[37]+t[36]+t[35]+t[31];
478     d2fdrbrb =t[1]*(t[13]*(1.333333333333333*t[26]*(t[53]+t[49])+0.44444444444444*pow(t[33],2.0)*t[51]+1.333333333333333*t[29]*(t[50]+t[52])+0.44444444444444*pow(t[32],2.0)*t[47])*t[7]+2.0*t[20]*t[34]*t[7]+t[46]+t[45]+t[44]+t[43]+t[42])+2.0*t[13]*t[34]*t[7]+t[40]+t[39]+t[38]+t[37]+t[36];
479     ds->df1000 += factor*dfdra;
480     ds->df0100 += factor*dfdrb;
481     ds->df2000 += factor*d2fdrara;
482     ds->df1100 += factor*d2fdrarb;
483     ds->df0200 += factor*d2fdrbrb;
484 }
485 
486 static void
pz81b_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)487 pz81b_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
488 {
489     real t[90];
490     real dfdra, dfdrb;
491     real d2fdrara, d2fdrarb, d2fdrbrb;
492     real d3fdrarara, d3fdrararb, d3fdrarbrb, d3fdrbrbrb;
493     real rhoa = dp->rhoa;
494     real rhob = dp->rhob;
495 
496     t[1] = rhob+rhoa;
497     t[2] = 1/pow(t[1],0.33333333333333);
498     t[3] = 0.6203504908994*Du*t[2];
499     t[4] = LOG(0.6203504908994*t[2]);
500     t[5] = Au*t[4];
501     t[6] = 0.6203504908994*Cu*t[2]*t[4];
502     t[7] = 1/(2.0*pow(2.0,0.33333333333333)-2.0);
503     t[8] = rhoa-1.0*rhob;
504     t[9] = 1/t[1];
505     t[10] = 1.0-1.0*t[8]*t[9];
506     t[11] = t[8]*t[9]+1.0;
507     t[12] = pow(t[11],1.333333333333333)+pow(t[10],1.333333333333333)-2.0;
508     t[13] = -1.0*Bu+Bp-1.0*t[4]*Au+Ap*t[4]-0.6203504908994*Cu*t[2]*t[4]+0.6203504908994*Cp*t[2]*t[4]-0.6203504908994*Du*t[2]+0.6203504908994*Dp*t[2];
509     t[14] = t[7]*t[12]*t[13];
510     t[15] = 1/pow(t[1],1.333333333333333);
511     t[16] = -0.20678349696647*Cu*t[15];
512     t[17] = -0.20678349696647*Du*t[15];
513     t[18] = -0.33333333333333*Au*t[9];
514     t[19] = -0.20678349696647*Cu*t[15]*t[4];
515     t[20] = 0.20678349696647*Cu*t[15]*t[4]-0.20678349696647*Cp*t[15]*t[4]+0.33333333333333*Au*t[9]-0.33333333333333*Ap*t[9]+0.20678349696647*Du*t[15]-0.20678349696647*Dp*t[15]+0.20678349696647*Cu*t[15]-0.20678349696647*Cp*t[15];
516     t[21] = t[7]*t[12]*t[20];
517     t[22] = 1/pow(t[1],2.0);
518     t[23] = t[8]*t[22];
519     t[24] = -1.0*t[9];
520     t[25] = t[24]+t[23];
521     t[26] = pow(t[10],0.33333333333333);
522     t[27] = -1.0*t[22]*t[8];
523     t[28] = t[9]+t[27];
524     t[29] = pow(t[11],0.33333333333333);
525     t[30] = 1.333333333333333*t[28]*t[29]+1.333333333333333*t[25]*t[26];
526     t[31] = t[7]*t[30]*t[13];
527     t[32] = t[9]+t[23];
528     t[33] = t[24]+t[27];
529     t[34] = 1.333333333333333*t[29]*t[33]+1.333333333333333*t[26]*t[32];
530     t[35] = t[7]*t[34]*t[13];
531     t[36] = -0.41356699393293*Cu*t[15];
532     t[37] = -0.41356699393293*Du*t[15];
533     t[38] = -0.66666666666667*Au*t[9];
534     t[39] = -0.41356699393293*Cu*t[15]*t[4];
535     t[40] = 2.0*t[12]*t[20]*t[7];
536     t[41] = 1/pow(t[1],2.333333333333334);
537     t[42] = 0.34463916161078*Cu*t[41];
538     t[43] = 0.27571132928862*Du*t[41];
539     t[44] = 0.33333333333333*Au*t[22];
540     t[45] = 0.27571132928862*Cu*t[41]*t[4];
541     t[46] = -0.27571132928862*Cu*t[41]*t[4]+0.27571132928862*Cp*t[41]*t[4]-0.33333333333333*Au*t[22]+0.33333333333333*Ap*t[22]-0.27571132928862*Du*t[41]+0.27571132928862*Dp*t[41]-0.34463916161078*Cu*t[41]+0.34463916161078*Cp*t[41];
542     t[47] = t[7]*t[12]*t[46];
543     t[48] = 2.0*t[20]*t[30]*t[7];
544     t[49] = pow(t[25],2.0);
545     t[50] = 1/pow(t[10],0.66666666666667);
546     t[51] = 1/pow(t[1],3.0);
547     t[52] = -2.0*t[51]*t[8];
548     t[53] = 2.0*t[22];
549     t[54] = t[53]+t[52];
550     t[55] = pow(t[28],2.0);
551     t[56] = 1/pow(t[11],0.66666666666667);
552     t[57] = 2.0*t[51]*t[8];
553     t[58] = -2.0*t[22];
554     t[59] = t[58]+t[57];
555     t[60] = 1.333333333333333*t[29]*t[59]+0.44444444444444*t[55]*t[56]+1.333333333333333*t[26]*t[54]+0.44444444444444*t[49]*t[50];
556     t[61] = t[7]*t[60]*t[13];
557     t[62] = 2.666666666666667*t[29]*t[51]*t[8]-2.666666666666667*t[26]*t[51]*t[8]+0.44444444444444*t[28]*t[33]*t[56]+0.44444444444444*t[25]*t[32]*t[50];
558     t[63] = 2.0*t[20]*t[34]*t[7];
559     t[64] = pow(t[32],2.0);
560     t[65] = t[58]+t[52];
561     t[66] = pow(t[33],2.0);
562     t[67] = t[53]+t[57];
563     t[68] = 1.333333333333333*t[29]*t[67]+0.44444444444444*t[56]*t[66]+1.333333333333333*t[26]*t[65]+0.44444444444444*t[50]*t[64];
564     t[69] = t[7]*t[68]*t[13];
565     t[70] = 1.033917484832333*Cu*t[41];
566     t[71] = 0.82713398786587*Du*t[41];
567     t[72] = Au*t[22];
568     t[73] = 0.82713398786587*Cu*t[41]*t[4];
569     t[74] = 3.0*t[12]*t[46]*t[7];
570     t[75] = 2.0*t[13]*t[62]*t[7];
571     t[76] = 1/pow(t[1],3.333333333333334);
572     t[77] = -0.89606182018802*Cu*t[76];
573     t[78] = -0.64332643500679*Du*t[76];
574     t[79] = -0.66666666666667*Au*t[51];
575     t[80] = -0.64332643500679*Cu*t[76]*t[4];
576     t[81] = t[7]*t[12]*(0.64332643500679*Cu*t[76]*t[4]-0.64332643500679*Cp*t[76]*t[4]+0.66666666666667*Au*t[51]-0.66666666666667*Ap*t[51]+0.64332643500679*Du*t[76]-0.64332643500679*Dp*t[76]+0.89606182018802*Cu*t[76]-0.89606182018802*Cp*t[76]);
577     t[82] = 2.0*t[20]*t[62]*t[7];
578     t[83] = 1/pow(t[10],1.666666666666667);
579     t[84] = 1/pow(t[1],4.0);
580     t[85] = 6.0*t[8]*t[84];
581     t[86] = 1/pow(t[11],1.666666666666667);
582     t[87] = -6.0*t[8]*t[84];
583     t[88] = -6.0*t[51];
584     t[89] = 6.0*t[51];
585     dfdra = t[1]*(t[31]+t[21]+t[19]+t[18]+t[17]+t[16])+t[14]+t[6]+t[5]+t[3]+Bu;
586     dfdrb = t[1]*(t[35]+t[21]+t[19]+t[18]+t[17]+t[16])+t[14]+t[6]+t[5]+t[3]+Bu;
587     d2fdrara = 2.0*t[13]*t[30]*t[7]+t[1]*(t[61]+t[48]+t[47]+t[45]+t[44]+t[43]+t[42])+t[40]+t[39]+t[38]+t[37]+t[36];
588     d2fdrarb = t[1]*(t[7]*t[62]*t[13]+t[7]*t[30]*t[20]+t[7]*t[34]*t[20]+t[47]+t[45]+t[44]+t[43]+t[42])+t[31]+t[35]+t[40]+t[39]+t[38]+t[37]+t[36];
589     d2fdrbrb = 2.0*t[13]*t[34]*t[7]+t[1]*(t[69]+t[63]+t[47]+t[45]+t[44]+t[43]+t[42])+t[40]+t[39]+t[38]+t[37]+t[36];
590     d3fdrararb = t[1]*(t[13]*t[7]*(1.333333333333333*t[29]*(t[87]+2.0*t[51])-0.2962962962963*t[33]*t[55]*t[86]+1.333333333333333*t[26]*(t[85]-2.0*t[51])-0.2962962962963*t[32]*t[49]*t[83]+1.777777777777778*t[28]*t[51]*t[56]*t[8]-1.777777777777778*t[25]*t[50]*t[51]*t[8]+0.44444444444444*t[33]*t[56]*t[59]+0.44444444444444*t[32]*t[50]*t[54])+t[82]+t[81]+t[80]+t[79]+t[78]+t[77]+2.0*t[30]*t[46]*t[7]+t[7]*t[34]*t[46]+t[7]*t[60]*t[20])+t[75]+t[74]+t[73]+t[72]+t[71]+t[70]+4.0*t[20]*t[30]*t[7]+t[63]+t[61];
591     d3fdrarbrb = t[1]*(t[13]*t[7]*(-0.2962962962963*t[28]*t[66]*t[86]-8.0*t[29]*t[8]*t[84]+8.0*t[26]*t[8]*t[84]-0.2962962962963*t[25]*t[64]*t[83]+1.777777777777778*t[33]*t[51]*t[56]*t[8]-1.777777777777778*t[32]*t[50]*t[51]*t[8]+0.44444444444444*t[28]*t[56]*t[67]+0.44444444444444*t[25]*t[50]*t[65]-2.666666666666667*t[29]*t[51]+2.666666666666667*t[26]*t[51])+t[82]+t[81]+t[80]+t[79]+t[78]+t[77]+2.0*t[34]*t[46]*t[7]+t[7]*t[30]*t[46]+t[7]*t[68]*t[20])+t[75]+t[74]+t[73]+t[72]+t[71]+t[70]+4.0*t[20]*t[34]*t[7]+t[69]+t[48];
592     d3fdrarara = t[1]*(t[13]*t[7]*(1.333333333333333*t[29]*(t[89]+t[87])+1.333333333333333*t[26]*(t[88]+t[85])-0.2962962962963*pow(t[28],3.0)*t[86]-0.2962962962963*pow(t[25],3.0)*t[83]+1.333333333333333*t[28]*t[56]*t[59]+1.333333333333333*t[25]*t[50]*t[54])+t[81]+t[80]+t[79]+t[78]+t[77]+3.0*t[20]*t[60]*t[7]+3.0*t[30]*t[46]*t[7])+t[74]+t[73]+t[72]+t[71]+t[70]+3.0*t[13]*t[60]*t[7]+6.0*t[20]*t[30]*t[7];
593     d3fdrbrbrb = t[1]*(t[13]*t[7]*(1.333333333333333*t[26]*(t[89]+t[85])+1.333333333333333*t[29]*(t[88]+t[87])-0.2962962962963*pow(t[33],3.0)*t[86]-0.2962962962963*pow(t[32],3.0)*t[83]+1.333333333333333*t[33]*t[56]*t[67]+1.333333333333333*t[32]*t[50]*t[65])+t[81]+t[80]+t[79]+t[78]+t[77]+3.0*t[20]*t[68]*t[7]+3.0*t[34]*t[46]*t[7])+t[74]+t[73]+t[72]+t[71]+t[70]+3.0*t[13]*t[68]*t[7]+6.0*t[20]*t[34]*t[7];
594     ds->df1000 += factor*dfdra;
595     ds->df0100 += factor*dfdrb;
596     ds->df2000 += factor*d2fdrara;
597     ds->df1100 += factor*d2fdrarb;
598     ds->df0200 += factor*d2fdrbrb;
599     ds->df3000 += factor*d3fdrarara;
600     ds->df2100 += factor*d3fdrararb;
601     ds->df1200 += factor*d3fdrarbrb;
602     ds->df0300 += factor*d3fdrbrbrb;
603 }
604 
605 /* ******************************************************************* */
606 /*                        Dispatch part                                */
607 /* ******************************************************************* */
608 static real
pz81_energy(const FunDensProp * dp)609 pz81_energy(const FunDensProp* dp)
610 {
611     real rs3 = 3/(4*M_PI*(dp->rhoa+dp->rhob));
612     if(rs3>=1)
613         return pz81a_energy(dp);
614     else
615         return pz81b_energy(dp);
616 }
617 
618 static void
pz81_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)619 pz81_first(FunFirstFuncDrv *ds,  real factor, const FunDensProp* dp)
620 {
621     real rs3 = 3/(4*M_PI*(dp->rhoa+dp->rhob));
622     if(rs3>=1)
623         pz81a_first(ds, factor, dp);
624     else
625         pz81b_first(ds, factor, dp);
626 }
627 
628 static void
pz81_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)629 pz81_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
630 {
631     real rs3 = 3/(4*M_PI*(dp->rhoa+dp->rhob));
632     if(rs3>=1)
633         pz81a_second(ds, factor, dp);
634     else
635         pz81b_second(ds, factor, dp);
636 }
637 
638 static void
pz81_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)639 pz81_third (FunThirdFuncDrv *ds,  real factor, const FunDensProp* dp)
640 {
641     real rs3 = 3/(4*M_PI*(dp->rhoa+dp->rhob));
642     if(rs3>=1)
643         pz81a_third(ds, factor, dp);
644     else
645         pz81b_third(ds, factor, dp);
646 }
647