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