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