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