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 /* -*-mode:c; c-style:bsd; c-basic-offset:4;indent-tabs-mode:nil; -*- */
25 /* fun-vwn.c:
26 implementation of VWN functional and its derivatives
27 (c), Pawel Salek, pawsa@theochem.kth.se, sep 2001, nov 2002
28 */
29
30 /* strictly conform to XOPEN ANSI C standard */
31 #if !defined(SYS_DEC)
32 /* XOPEN compliance is missing on old Tru64 4.0E Alphas */
33 #define _XOPEN_SOURCE 500
34 #define _XOPEN_SOURCE_EXTENDED 1
35 #endif
36 #include <math.h>
37 #include <stddef.h>
38 #include <stdlib.h>
39 #include "general.h"
40
41 #define __CVERSION__
42
43 #include "functionals.h"
44
45 /* INTERFACE PART */
vwn_isgga(void)46 static integer vwn_isgga(void) { return 0; }
47 static integer vwn_read(const char* conf_line);
48 static real vwn3_energy(const FunDensProp* dp);
49 static void vwn3_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp);
50 static void vwn3_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp);
51 static void vwn3_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp);
52 static void vwn3_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp);
53 static real vwn_energy(const FunDensProp* dp);
54 static void vwn_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp);
55 static void vwn_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp);
56 static void vwn_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp);
57 static void vwn_fourth(FunFourthFuncDrv *ds, real factor,
58 const FunDensProp* dp);
59 static real vwni_energy(const FunDensProp* dp);
60 static void vwni_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp);
61 static void vwni_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp);
62 static void vwni_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp);
63 static void vwni_fourth(FunFourthFuncDrv *ds, real factor,
64 const FunDensProp* dp);
65 static real vwn3i_energy(const FunDensProp* dp);
66 static void vwn3i_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp);
67 static void vwn3i_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp);
68 static void vwn3i_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp);
69 static void vwn3i_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp);
70
71 /* VWN3 is a Gaussian version of the VWN functional based on suboptimal
72 * set of parameters */
73 Functional VWN3Functional = {
74 "VWN3", /* name */
75 vwn_isgga, /* gga-corrected */
76 1,
77 vwn_read, /* no extra input expected, just set the common block */
78 NULL,
79 vwn3_energy,
80 vwn3_first,
81 vwn3_second,
82 vwn3_third,
83 vwn3_fourth
84 };
85
86 Functional VWN5Functional = {
87 "VWN5", /* name */
88 vwn_isgga, /* gga-corrected */
89 3,
90 vwn_read, /* no extra input expected, just set the common block */
91 NULL,
92 vwn_energy,
93 vwn_first,
94 vwn_second,
95 vwn_third,
96 vwn_fourth
97 };
98
99 /* VWN is used for backward compatibility only */
100 Functional VWNFunctional = {
101 "VWN", /* name */
102 vwn_isgga, /* gga-corrected */
103 1,
104 vwn_read, /* no extra input expected, just set the common block */
105 NULL,
106 vwn_energy,
107 vwn_first,
108 vwn_second,
109 vwn_third,
110 vwn_fourth
111 };
112
113 /* VWNIFunctional is a variant of VWN5 functional with another spin
114 polarization dependence:
115
116 F(r,zeta) = (E_p + f(zeta)*(E_f - E_p))*rho
117
118 */
119 Functional VWNIFunctional = {
120 "VWNI", /* name */
121 vwn_isgga, /* gga-corrected */
122 1,
123 vwn_read, /* no extra input expected, just set the common block */
124 NULL,
125 vwni_energy,
126 vwni_first,
127 vwni_second,
128 vwni_third,
129 vwni_fourth
130 };
131 Functional VWN3IFunctional = {
132 "VWN3I", /* name */
133 vwn_isgga, /* gga-corrected */
134 1,
135 vwn_read, /* no extra input expected, just set the common block */
136 NULL,
137 vwn3i_energy,
138 vwn3i_first,
139 vwn3i_second,
140 vwn3i_third,
141 vwn3i_fourth
142 };
143
144
145 /* IMPLEMENTATION PART */
146 #define VWN_ZERO 1e-40
147
148 static integer
vwn_read(const char * conf_line)149 vwn_read(const char* conf_line)
150 {
151 fun_set_hf_weight(0);
152 return 1;
153 }
154
155 /* vwn_params contains two sets of parameters for paramagnetic and
156 ferromagnetic cases. See Table 5 in VWN paper.
157 */
158
159
160 static const struct vwn_params {
161 real X0, A, B, C;
162 } vwn_paramagnetic = { -0.1049800, 0.0621814, 3.72744, 12.9352 },
163 vwn_ferromagnetic = { -0.3250000, 0.0310907, 7.06042, 18.0578 },
164 vwn_interp = { -0.0047584,-0.0337737, 1.13107, 13.0045 },
165 vwn3_paramagnetic = { -0.4092860, 0.0621814, 13.0720, 42.7198 },
166 vwn3_ferromagnetic= { -0.7432940, 0.0310907, 20.1231, 101.578 };
167
168
169 static const real SPINPOLF = 1.92366105093154; /* 1/(2(2^(1/3)-1)) */
170 static const real THREEFTHRD2 = 0.584822305543806;/* hm? 4.5/(4*SPINPOLF) */
171 static const real FOURTHREE = 1.333333333333333;
172
173 /* vwn_en_pot:
174 returns "order" numbers in enpot array
175 enpot[0]: energy of given type p E = rho F - THIS IS AN EXCEPTION!
176 DO NOT BLAME IT ON ME.
177 enpot[1]: E'=enpot[0] + rho*F'
178 enpot[2]: E''
179 enpot[3]: E'''
180 */
181 static void
vwn_en_pot(real * enpot,real rho,integer order,const struct vwn_params * p)182 vwn_en_pot(real* enpot, real rho, integer order, const struct vwn_params* p)
183 {
184 const real
185 AI = p->A,
186 BI = p->B,
187 CI = p->C,
188 X0I = p->X0;
189 const real
190 Q = sqrt(4*CI - BI*BI),
191 XF0I = X0I*X0I + BI*X0I + CI,
192 YF0I = Q/(BI + 2*X0I),
193 DCRS = pow(3.0/(4*M_PI),1.0/6.0),
194 B = X0I/XF0I,
195 C = XF0I*YF0I,
196 ACON = B*BI - 1.0,
197 BCON = 2*ACON + 2.0,
198 CCON = 2*BI*(1.0/Q - X0I/C);
199
200 real rho13 = pow(rho,1.0/3.0);
201 real x = DCRS/sqrt(rho13);
202 real xrho = -DCRS*pow(rho,-7.0/6.0)/6.0;
203 real xxrho = +DCRS*pow(rho,-13.0/6.0)*7.0/36.0;
204 real xf = x*x + BI*x+CI;
205 real xfx = 2*x + BI;
206 real yf = Q/xfx;
207 real e1, ex1, exx1, exxx1;
208 real xf_p2, xf_p3, xf_p4, xfx_p2, xfx_p4;
209 real xrho_p2, xrho_p3, xrho_p4, xxrho_p2, xxxrho, xxxxrho;
210 real Q_p2, exxxx1;
211 real x_p4;
212
213 e1 = 2*log(x)
214 + ACON*log(xf)
215 - BCON*log(x - X0I)
216 + CCON*atan(yf);
217 enpot[0] = 0.5*AI*e1;
218 if(order<1) return;
219
220 ex1 = 2.0/x
221 + ACON*xfx/xf
222 - BCON/(x - X0I)
223 - CCON*(2*yf/xfx)/(1.0 + yf*yf);
224 enpot[1] = 0.5*AI*(e1 + rho*ex1*xrho);
225 if(order<2) return;
226
227 exx1= -2.0/(x*x)
228 + ACON*(2.0/xf - xfx*xfx/(xf*xf))
229 + BCON/((x - X0I)*(x - X0I))
230 + CCON*8*Q*xfx/((Q*Q + xfx*xfx)*(Q*Q + xfx*xfx));
231 enpot[2] = 0.5*AI*xrho*(2*ex1 + rho*exx1*xrho - ex1*7.0/6.0);
232 if(order<3) return;
233
234 exxx1= 4.0/(x*x*x)
235 + ACON*(2.0*xfx/(xf*xf))*(xfx*xfx/xf-3.0)
236 - BCON*2.0/pow(x - X0I,3.0)
237 + CCON*16.0*Q*(Q*Q-3.0*xfx*xfx)/pow(Q*Q + xfx*xfx,3.0);
238 enpot[3] = 0.5*AI*xxrho*(2*ex1 + rho*exx1*xrho - 7.0/6.0*ex1)
239 + 0.5*AI*xrho*(2*exx1*xrho
240 +exx1*xrho + rho*(exxx1*xrho*xrho+exx1*xxrho)
241 -7.0/6.0*exx1*xrho);
242
243 if(order<4) return;
244 x_p4 = pow(x,4.0);
245 xf_p2 = xf*xf;
246 xf_p3 = xf_p2*xf;
247 xf_p4 = xf_p2*xf_p2;
248 xfx_p2 = xfx*xfx;
249 xfx_p4 = xfx_p2*xfx_p2;
250 xrho_p2 = xrho*xrho;
251 xrho_p3 = xrho_p2*xrho;
252 xrho_p4 = xrho_p2*xrho_p2;
253 xxrho_p2 = xxrho*xxrho;
254 xxxrho = (-91.0/216.0)*DCRS/pow(rho,19.0/6.0);
255 xxxxrho = (1729.0/1296.0)*DCRS/pow(rho,25.0/6.0);
256 Q_p2 = Q*Q;
257
258 exxxx1 = 6*((-2*ACON)/xf_p2 + (4*xfx_p2*ACON)/xf_p3 -
259 (xfx_p4*ACON)/xf_p4 + (64*xfx*CCON*(xfx - Q)*Q*(xfx + Q))/
260 pow((xfx_p2 + Q_p2),4.0) - 2/x_p4 + BCON/pow((x - X0I),4.0));
261
262 enpot[4] = 0.5*AI*(4*exxx1*xrho_p3 + exxxx1*rho*xrho_p4 +
263 6*exxx1*rho*xrho_p2*xxrho +
264 3*exx1*rho*xxrho_p2 +
265 4*exx1*xrho*(3*xxrho + rho*xxxrho) +
266 ex1*(4*xxxrho + rho*xxxxrho));
267 /* for fourth derivatives we need also
268 first derivative of enpot[0] with respect to rho.
269 It is stored in enpot[5] (what a mess!) */
270 enpot[5] = 0.5*AI*ex1*xrho;
271
272 }
273
274 static real
par_energy(const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)275 par_energy(const FunDensProp* dp, const struct vwn_params* para,
276 const struct vwn_params* ferro)
277 {
278 real ep_p[2], ep_f[2], ep_i[2], zeta, zeta4, f_zeta, delta;
279 real rhoa = dp->rhoa, rhob = dp->rhob, rho;
280
281 if(rhoa<VWN_ZERO) rhoa = VWN_ZERO;
282 if(rhob<VWN_ZERO) rhob = VWN_ZERO;
283 rho = rhoa + rhob;
284 vwn_en_pot(ep_p, rho, 0, para);
285
286 if(fabs(dp->rhoa-dp->rhob)<VWN_ZERO) return ep_p[0]*rho;
287 vwn_en_pot(ep_f, rho, 0, ferro);
288 vwn_en_pot(ep_i, rho, 0, &vwn_interp);
289
290 zeta = (dp->rhoa-dp->rhob)/rho;
291 zeta4 = pow(zeta,4.0);
292 f_zeta = SPINPOLF*(pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
293 delta = f_zeta*((ep_f[0]-ep_p[0])*zeta4 + ep_i[0]*(1-zeta4)*THREEFTHRD2);
294
295 return (ep_p[0]+ delta)*rho;
296 }
297
298 static void
par_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)299 par_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp,
300 const struct vwn_params* para, const struct vwn_params* ferro)
301 {
302 real zeta, zeta3,zeta4, f_zeta, f_zet1, e_f,acp, dacp, vcfp, g_f;
303 real delta, ep_p[2], ep_f[2], ep_i[2];
304 real rhoa = dp->rhoa, rhob = dp->rhob, rho;
305
306 if(rhoa<VWN_ZERO) rhoa = VWN_ZERO;
307 if(rhob<VWN_ZERO) rhob = VWN_ZERO;
308 rho = rhoa + rhob;
309 vwn_en_pot(ep_p, rho, 1, para);
310
311 ds->df1000 += ep_p[1]*factor;
312 ds->df0100 += ep_p[1]*factor;
313
314 if(fabs(dp->rhoa-dp->rhob)<VWN_ZERO) return;
315
316 /* contribution from spin-polarized case; first order */
317 zeta = (dp->rhoa-dp->rhob)/rho;
318 zeta3 = pow(zeta,3.0);
319 zeta4 = zeta3*zeta;
320 f_zeta = SPINPOLF*(pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
321 f_zet1 = SPINPOLF*4.0/3.0*(pow(1+zeta,1.0/3.0)-pow(1-zeta,1.0/3.0));
322 vwn_en_pot(ep_f, rho, 1, ferro);
323 e_f = ep_f[0] - ep_p[0];
324 g_f = ep_f[1] - ep_p[1];
325 vwn_en_pot(ep_i, rho, 1, &vwn_interp);
326 acp = ep_i[0]*THREEFTHRD2;
327 dacp = ep_i[1]*THREEFTHRD2;
328
329 vcfp = f_zeta*(g_f*zeta4 + dacp*(1-zeta4));
330 delta= (f_zet1*(e_f*zeta4 + acp*(1-zeta4)) +
331 4*f_zeta*(e_f - acp)*zeta3);
332
333 /* the final section: begin */
334 ds->df1000 += (vcfp + delta*(1-zeta))*factor;
335 ds->df0100 += (vcfp - delta*(1+zeta))*factor;
336 /* the final section: end */
337 }
338
339 /* vwn_second:
340 CAUTION: may raise zeros to a negative power!
341 */
342 static void
par_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)343 par_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp,
344 const struct vwn_params* para, const struct vwn_params* ferro)
345 {
346 real zeta, zeta2, zeta3,zeta4, f_zeta, f_zet1, f_zet2, vcfp;
347 real delta, ep_p[3], ep_f[3], ep_i[3];
348 real vcfp1, ef0, ef1, ef2, ei0, ei1, ei2, bterm, cterm, dterm;
349 real spA, spB;
350 real rhoa = dp->rhoa, rhob = dp->rhob, rho = dp->rhoa + dp->rhob;
351 real rho2 = rho*rho;
352 real dAA, dAB, dBB;
353
354 vwn_en_pot(ep_p, rho, 2, para);
355
356 ds->df1000 += ep_p[1]*factor;
357 ds->df0100 += ep_p[1]*factor;
358 ds->df2000 += ep_p[2]*factor;
359 ds->df0200 += ep_p[2]*factor;
360 ds->df1100 += ep_p[2]*factor;
361
362 /* if(0&&dp->rhoa==dp->rhob) return; */
363 /* contribution from spin-polarized case; second order */
364 zeta = (dp->rhoa-dp->rhob)/rho;
365 zeta2 = zeta*zeta;
366 zeta3 = zeta2*zeta;
367 zeta4 = zeta3*zeta;
368 f_zeta = SPINPOLF*(pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
369 f_zet1 = SPINPOLF*4.0/3.0*(pow(1+zeta,1.0/3.0)-pow(1-zeta,1.0/3.0));
370 /* CAUTION: may raise 0 to negative power! */
371 f_zet2 = SPINPOLF*4.0/9.0*(pow(1+zeta,-2.0/3.0)+pow(1-zeta,-2.0/3.0));
372 vwn_en_pot(ep_f, rho, 2, ferro);
373 ef0 = ep_f[0] - ep_p[0];
374 ef1 = ep_f[1] - ep_p[1];
375 ef2 = ep_f[2] - ep_p[2];
376 vwn_en_pot(ep_i, rho, 2, &vwn_interp);
377 ei0 = ep_i[0]*THREEFTHRD2;
378 ei1 = ep_i[1]*THREEFTHRD2;
379 ei2 = ep_i[2]*THREEFTHRD2;
380
381 bterm = ef1*zeta4 + ei1*(1-zeta4);
382 vcfp = f_zeta*bterm;
383 delta = (f_zet1*(ef0*zeta4 + ei0*(1-zeta4))
384 + 4*f_zeta*(ef0 - ei0)*zeta3);
385
386 spA = 2*rhob/rho2; /* = 2(1-zeta)/rho */
387 spB =-2*rhoa/rho2; /* = -2(1+zeta)/rho */
388 /* contribution from spin-polarized case; second order */
389 /* spin independent part of vcfp */
390 vcfp1 = f_zeta*(ef2*zeta4 + ei2*(1-zeta4));
391 /* spin dependent part of vcfp */
392 cterm = 4*f_zeta*(ef1-ei1)*zeta3 + bterm*f_zet1 - delta;
393
394 /* spin dependent part of delta */
395 dterm = (f_zet2*(ef0*zeta4+ei0*(1-zeta4))
396 +8*f_zet1*(ef0-ei0)*zeta3
397 +12*f_zeta*(ef0-ei0)*zeta2)*rho;
398
399 dAA = dterm*spA*spA;
400 dAB = dterm*spA*spB;
401 dBB = dterm*spB*spB;
402
403 /* the final section: begin */
404 ds->df1000 += (vcfp + delta*(1-zeta))*factor;
405 ds->df0100 += (vcfp - delta*(1+zeta))*factor;
406
407 ds->df2000 += (vcfp1+ cterm*(spA+spA) + dAA)*factor;
408 ds->df1100 += (vcfp1+ cterm*(spA+spB) + dAB)*factor;
409 ds->df0200 += (vcfp1+ cterm*(spB+spB) + dBB)*factor;
410 /* the final section: end */
411 }
412
413 /* third not tested for open-shell! */
414 static void
par_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)415 par_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp,
416 const struct vwn_params* para, const struct vwn_params* ferro)
417 {
418 real zeta, zeta2, zeta3,zeta4, f_zeta, f_zet1, f_zet2, f_zet3, vcfp;
419 real delta, ep_p[4], ep_f[4], ep_i[4];
420 real vcfp1, vcfp2, ef0, ef1, ef2, ef3, ei0, ei1, ei2, ei3;
421 real bterm, cterm, dterm, eterm, ctrm1, ctrm2, dtrm1, dtrm2;
422 real ef2bi, ei2bi;
423 real spA, spB, spAA, spAB, spBB;
424 real rhoa = dp->rhoa, rhob = dp->rhob, rho, rho2, rho3;
425
426 if(rhoa<VWN_ZERO) rhoa = VWN_ZERO;
427 if(rhob<VWN_ZERO) rhob = VWN_ZERO;
428 rho = rhoa + rhob;
429 rho2 = rho*rho; rho3 = rho2*rho;
430
431 vwn_en_pot(ep_p, rho, 3, para);
432
433 ds->df1000 += ep_p[1]*factor;
434 ds->df0100 += ep_p[1]*factor;
435 ds->df2000 += ep_p[2]*factor;
436 ds->df0200 += ep_p[2]*factor;
437 ds->df1100 += ep_p[2]*factor;
438
439 ds->df3000 += ep_p[3]*factor;
440 ds->df2100 += ep_p[3]*factor;
441 ds->df1200 += ep_p[3]*factor;
442 ds->df0300 += ep_p[3]*factor;
443
444 /* if(0&&dp->rhoa==dp->rhob) return; */
445 /* contribution from spin-polarized case; second order */
446 zeta = (dp->rhoa-dp->rhob)/rho;
447 zeta2 = zeta*zeta;
448 zeta3 = zeta2*zeta;
449 zeta4 = zeta3*zeta;
450 f_zeta = SPINPOLF* (pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
451 f_zet1 = SPINPOLF*4.0/3.0 *(pow(1+zeta, 1.0/3.0)-pow(1-zeta, 1.0/3.0));
452 f_zet2 = SPINPOLF*4.0/9.0 *(pow(1+zeta,-2.0/3.0)+pow(1-zeta,-2.0/3.0));
453 f_zet3 =-SPINPOLF*8.0/27.0*(pow(1+zeta,-5.0/3.0)-pow(1-zeta,-5.0/3.0));
454
455 vwn_en_pot(ep_f, rho, 3, ferro);
456 ef0 = ep_f[0] - ep_p[0];
457 ef1 = ep_f[1] - ep_p[1];
458 ef2 = ep_f[2] - ep_p[2];
459 ef3 = ep_f[3] - ep_p[3];
460 vwn_en_pot(ep_i, rho, 3,&vwn_interp);
461 ei0 = ep_i[0]*THREEFTHRD2;
462 ei1 = ep_i[1]*THREEFTHRD2;
463 ei2 = ep_i[2]*THREEFTHRD2;
464 ei3 = ep_i[3]*THREEFTHRD2;
465
466 bterm = ef1*zeta4 + ei1*(1-zeta4);
467 vcfp = f_zeta*bterm;
468 delta = (f_zet1*(ef0*zeta4 + ei0*(1-zeta4))
469 + 4*f_zeta*(ef0 - ei0)*zeta3);
470
471 spA = 2*rhob/rho2; /* = 2(1-zeta)/rho */
472 spB =-2*rhoa/rho2; /* = -2(1+zeta)/rho */
473 spAA = -4*rhob/rho3;
474 spAB = 2*(rho-2*rhob)/rho3;
475 spBB = 4*rhoa/rho3;
476 /* contribution from spin-polarized case; second order */
477 /* spin independent part of vcfp */
478 vcfp1 = f_zeta*(ef2*zeta4 + ei2*(1-zeta4));
479 /* spin dependent part of vcfp */
480 cterm = 4*f_zeta*(ef1-ei1)*zeta3 + bterm*f_zet1 - delta;
481
482 /* spin dependent part of delta */
483 dterm = (f_zet2*(ef0*zeta4+ei0*(1-zeta4))
484 +8*f_zet1*(ef0-ei0)*zeta3
485 +12*f_zeta*(ef0-ei0)*zeta2)*rho;
486
487 /* third order terms */
488 vcfp2 = f_zeta*(ef3*zeta4+ei3*(1-zeta4));
489 eterm = f_zet1*(ef2*zeta4 + ei2*(1-zeta4)) + f_zeta*(ef2-ei2)*4*zeta3;
490 ef2bi = ef2-(ef1-ef0)/rho;
491 ei2bi = ei2-(ei1-ei0)/rho;
492 ctrm1 = 4*f_zeta*(ef2bi-ei2bi)*zeta3 +f_zet1*(ef2bi*zeta4+ei2bi*(1-zeta4));
493
494 ctrm2 = (ef1-ei1-ef0+ei0)*(8*f_zet1*zeta3+12*f_zeta*zeta2)
495 +f_zet2*(bterm-(ef0*zeta4 + ei0*(1-zeta4)));
496
497 dtrm1 = f_zet2*((ef1-ef0)*zeta4+(ei1-ei0)*(1-zeta4))
498 +(8*f_zet1*zeta3+12*f_zeta*zeta2)*(ef1-ei1-ef0+ei0)
499 +dterm/rho;
500 dtrm2 = ((12*f_zet2*zeta3 + 36*f_zet1*zeta2 + 24*f_zeta*zeta)*(ef0-ei0)+
501 f_zet3*(ef0*zeta4+ei0*(1-zeta4)))*rho;
502
503 /* the final section: begin */
504 ds->df1000 += (vcfp + delta*(1-zeta))*factor;
505 ds->df0100 += (vcfp - delta*(1+zeta))*factor;
506
507 ds->df2000 += (vcfp1+ cterm*(spA+spA) + dterm*spA*spA)*factor;
508 ds->df1100 += (vcfp1+ cterm*(spA+spB) + dterm*spA*spB)*factor;
509 ds->df0200 += (vcfp1+ cterm*(spB+spB) + dterm*spB*spB)*factor;
510 ds->df3000 += (vcfp2+ eterm*spA +
511 ctrm1*(spA+spA)+ ctrm2*spA*(spA+spA) + cterm*(spAA+spAA) +
512 dtrm1*(spA*spA)+ dtrm2*spA*spA*spA + dterm*(2*spAA*spA)
513 )*factor;
514 ds->df2100 += (vcfp2+ eterm*spB +
515 ctrm1*(spA+spA)+ ctrm2*spB*(spA+spA) + cterm*(spAB+spAB) +
516 dtrm1*(spA*spA)+ dtrm2*spA*spA*spB + dterm*(2*spAB*spA)
517 )*factor;
518 ds->df1200 += (vcfp2+ eterm*spA +
519 ctrm1*(spB+spB)+ ctrm2*spA*(spB+spB) + cterm*(spAB+spAB) +
520 dtrm1*(spB*spB)+ dtrm2*spB*spB*spA + dterm*(2*spAB*spB)
521 )*factor;
522 ds->df0300 += (vcfp2+ eterm*spB +
523 ctrm1*(spB+spB)+ ctrm2*spB*(spB+spB) + cterm*(spBB+spBB) +
524 dtrm1*(spB*spB)+ dtrm2*spB*spB*spB + dterm*(2*spBB*spB)
525 )*factor;
526 /* the final section: end */
527 }
528
529 static void
par_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)530 par_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp,
531 const struct vwn_params* para, const struct vwn_params* ferro)
532 {
533 real zeta, zeta2, zeta3,zeta4, f_zeta, f_zet1, f_zet2, f_zet3, vcfp;
534 real delta, ep_p[6], ep_f[6], ep_i[6], d_ef0, d_ei0;
535 real vcfp1, vcfp2, vcfp3, ef0, ef1, ef2, ef3, ei0, ei1, ei2, ei3;
536 real bterm, cterm, dterm, eterm, ctrm1, ctrm2, dtrm1, dtrm2;
537 real ef2bi, ei2bi;
538 real spA, spB, spAA, spAB, spBB;
539 real rhoa = dp->rhoa, rhob = dp->rhob, rho = dp->rhoa + dp->rhob;
540 real rhoa_p2 = rhoa*rhoa, rhob_p2 = rhob*rhob;
541 real rho2 = rho*rho;
542 real rho3 = rho2*rho;
543 real rho4 = rho2*rho2;
544 real zeta_p2, zeta_p3, zeta_p4, zeta_p5;
545 real mzeta, pzeta, mzeta_p2, pzeta_p2;
546 real f_zet4;
547 real ef4, ei4, ef_ei, ef1_ei1, ef2_ei2, ef3_ei3, ef4_ei4, ef, ei;
548 real spAAA, spBBB, spAAB, spABB;
549 real spA_p2, spA_p3, spB_p2, spB_p3;
550 real d_ef2bi, d_ei2bi, d_bterm_indep, d_bterm_dep, d_delta_indep,
551 d_delta_dep, d_vcfp2_indep, d_vcfp2_dep, d_eterm_indep,
552 d_eterm_dep, d_ctrm1_indep, d_ctrm1_dep, d_ctrm2_indep,
553 d_ctrm2_dep, d_cterm_indep, d_cterm_dep, d_dterm_indep,
554 d_dterm_dep, d_dtrm1_indep, d_dtrm1_dep, d_dtrm2_indep,
555 d_dtrm2_dep;
556
557 vwn_en_pot(ep_p, rho, 4, para);
558
559 ds->df1000 += ep_p[1]*factor;
560 ds->df0100 += ep_p[1]*factor;
561
562 ds->df2000 += ep_p[2]*factor;
563 ds->df0200 += ep_p[2]*factor;
564 ds->df1100 += ep_p[2]*factor;
565
566 ds->df3000 += ep_p[3]*factor;
567 ds->df2100 += ep_p[3]*factor;
568 ds->df1200 += ep_p[3]*factor;
569 ds->df0300 += ep_p[3]*factor;
570
571 ds->df4000 += ep_p[4]*factor;
572 ds->df3100 += ep_p[4]*factor;
573 ds->df2200 += ep_p[4]*factor;
574 ds->df1300 += ep_p[4]*factor;
575 ds->df0400 += ep_p[4]*factor;
576
577 //if(dp->rhoa==dp->rhob) return;
578 /* contribution from spin-polarized case; second order */
579 zeta = (dp->rhoa-dp->rhob)/rho;
580 zeta_p2 = zeta2 = zeta*zeta;
581 zeta_p3 = zeta3 = zeta2*zeta;
582 zeta_p4 = zeta4 = zeta3*zeta;
583 zeta_p5 = zeta4*zeta;
584 mzeta = -1.0+zeta;
585 pzeta = 1.0+zeta;
586 mzeta_p2 = mzeta*mzeta;
587 pzeta_p2 = pzeta*pzeta;
588 f_zeta = SPINPOLF* (pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
589 f_zet1 = SPINPOLF*4.0/3.0 *(pow(1+zeta, 1.0/3.0)-pow(1-zeta, 1.0/3.0));
590 f_zet2 = SPINPOLF*4.0/9.0 *(pow(1+zeta,-2.0/3.0)+pow(1-zeta,-2.0/3.0));
591 f_zet3 =-SPINPOLF*8.0/27.0*(pow(1+zeta,-5.0/3.0)-pow(1-zeta,-5.0/3.0));
592 f_zet4 = SPINPOLF*(40.0/81.0)*(pow(1-zeta,-8.0/3.0)+pow(1+zeta,-8.0/3.0));
593
594 vwn_en_pot(ep_f, rho, 4, ferro);
595 ef0 = ep_f[0] - ep_p[0];
596 ef1 = ep_f[1] - ep_p[1];
597 ef2 = ep_f[2] - ep_p[2];
598 ef3 = ep_f[3] - ep_p[3];
599 ef4 = ep_f[4] - ep_p[4];
600 d_ef0 = ep_f[5] - ep_p[5];
601 vwn_en_pot(ep_i, rho, 4,&vwn_interp);
602 ei0 = ep_i[0]*THREEFTHRD2;
603 ei1 = ep_i[1]*THREEFTHRD2;
604 ei2 = ep_i[2]*THREEFTHRD2;
605 ei3 = ep_i[3]*THREEFTHRD2;
606 ei4 = ep_i[4]*THREEFTHRD2;
607 d_ei0 = ep_i[5]*THREEFTHRD2;
608
609 ef = ef0;
610 ei = ei0;
611
612 ef_ei = ef - ei;
613
614 ef1_ei1 = ef1 - ei1;
615 ef2_ei2 = ef2 - ei2;
616 ef3_ei3 = ef3 - ei3;
617 ef4_ei4 = ef4 - ei4;
618
619
620 bterm = ef1*zeta4 + ei1*(1-zeta4);
621 vcfp = f_zeta*bterm;
622 delta = (f_zet1*(ef0*zeta4 + ei0*(1-zeta4))
623 + 4*f_zeta*(ef0 - ei0)*zeta3);
624
625 spA = 2*rhob/rho2; /* = 2(1-zeta)/rho */
626 spB =-2*rhoa/rho2; /* = -2(1+zeta)/rho */
627 spAA = -4*rhob/rho3;
628 spAB = 2*(rho-2*rhob)/rho3;
629 spBB = 4*rhoa/rho3;
630 spAAA = 12*rhob/rho4;
631 spAAB = -4*(rhoa - 2* rhob)/rho4;
632 spABB = (-8*rhoa + 4*rhob)/rho4;
633 spBBB = -12*rhoa/rho4;
634
635 /* Powers of sp */
636 spA_p2 = spA*spA;
637 spA_p3 = spA_p2*spA;
638
639 spB_p2 = spB*spB;
640 spB_p3 = spB_p2*spB;
641
642
643 /* contribution from spin-polarized case; second order */
644 /* spin independent part of vcfp */
645 vcfp1 = f_zeta*(ef2*zeta4 + ei2*(1-zeta4));
646 /* spin dependent part of vcfp */
647 cterm = 4*f_zeta*(ef1-ei1)*zeta3 + bterm*f_zet1 - delta;
648
649 /* spin dependent part of delta */
650 dterm = (f_zet2*(ef0*zeta4+ei0*(1-zeta4))
651 +8*f_zet1*(ef0-ei0)*zeta3
652 +12*f_zeta*(ef0-ei0)*zeta2)*rho;
653
654 /* third order terms */
655 vcfp2 = f_zeta*(ef3*zeta4+ei3*(1-zeta4));
656
657 eterm = f_zet1*(ef2*zeta4 + ei2*(1-zeta4)) + f_zeta*(ef2-ei2)*4*zeta3;
658
659 ef2bi = ef2-(ef1-ef0)/rho;
660 ei2bi = ei2-(ei1-ei0)/rho;
661 ctrm1 = 4*f_zeta*(ef2bi-ei2bi)*zeta3 +f_zet1*(ef2bi*zeta4+ei2bi*(1-zeta4));
662
663 ctrm2 = (ef1-ei1-ef0+ei0)*(8*f_zet1*zeta3+12*f_zeta*zeta2)
664 +f_zet2*(bterm-(ef0*zeta4 + ei0*(1-zeta4)));
665
666 dtrm1 = f_zet2*((ef1-ef0)*zeta4+(ei1-ei0)*(1-zeta4))
667 +(8*f_zet1*zeta3+12*f_zeta*zeta2)*(ef1-ei1-ef0+ei0)
668 +dterm/rho;
669 dtrm2 = ((12*f_zet2*zeta3 + 36*f_zet1*zeta2 + 24*f_zeta*zeta)*(ef0-ei0)+
670 f_zet3*(ef0*zeta4+ei0*(1-zeta4)))*rho;
671
672 /* the final section: begin */
673 ds->df1000 += (vcfp + delta*(1-zeta))*factor;
674 ds->df0100 += (vcfp - delta*(1+zeta))*factor;
675
676 ds->df2000 += (vcfp1+ cterm*(spA+spA) + dterm*spA*spA)*factor;
677 ds->df1100 += (vcfp1+ cterm*(spA+spB) + dterm*spA*spB)*factor;
678 ds->df0200 += (vcfp1+ cterm*(spB+spB) + dterm*spB*spB)*factor;
679
680 ds->df3000 += (vcfp2+ eterm*spA +
681 ctrm1*(spA+spA)+ ctrm2*spA*(spA+spA) + cterm*(spAA+spAA) +
682 dtrm1*(spA*spA)+ dtrm2*spA*spA*spA + dterm*(2*spAA*spA)
683 )*factor;
684 ds->df2100 += (vcfp2+ eterm*spB +
685 ctrm1*(spA+spA)+ ctrm2*spB*(spA+spA) + cterm*(spAB+spAB) +
686 dtrm1*(spA*spA)+ dtrm2*spA*spA*spB + dterm*(2*spAB*spA)
687 )*factor;
688 ds->df1200 += (vcfp2+ eterm*spA +
689 ctrm1*(spB+spB)+ ctrm2*spA*(spB+spB) + cterm*(spAB+spAB) +
690 dtrm1*(spB*spB)+ dtrm2*spB*spB*spA + dterm*(2*spAB*spB)
691 )*factor;
692 ds->df0300 += (vcfp2+ eterm*spB +
693 ctrm1*(spB+spB)+ ctrm2*spB*(spB+spB) + cterm*(spBB+spBB) +
694 dtrm1*(spB*spB)+ dtrm2*spB*spB*spB + dterm*(2*spBB*spB)
695 )*factor;
696
697 d_ef2bi = ef3 + (-ef0 + ef1)/rho2 - (-d_ef0 + ef2)/rho;
698 d_ei2bi = ei3 + (-ei0 + ei1)/rho2 - (-d_ei0 + ei2)/rho;
699
700
701 d_bterm_indep = ei2*(1 - zeta4) + ef2*zeta4;
702 d_bterm_dep = (4*ef1*zeta3 - 4*ei1*zeta3);
703
704
705 d_delta_indep = 4*(d_ef0 - d_ei0)*f_zeta*zeta3 +
706 d_ei0*f_zet1*(1 - zeta4) + d_ef0*f_zet1*zeta4;
707
708 d_delta_dep = (12*(ef0 - ei0)*f_zeta*zeta2 + 4*ef0*f_zet1*zeta3 +
709 4*(ef0 - ei0)*f_zet1*zeta3 - 4*ei0*f_zet1*zeta3 +
710 f_zet2*(ei0*(1 - zeta4) + ef0*zeta4));
711
712 d_vcfp2_indep = ei4*f_zeta*(1 - zeta4) + ef4*f_zeta*zeta4;
713
714 d_vcfp2_dep = (4*ef3*f_zeta*zeta3 - 4*ei3*f_zeta*zeta3 +
715 f_zet1*(ei3*(1 - zeta4) + ef3*zeta4));
716
717
718 d_eterm_indep = 4*(ef3 - ei3)*f_zeta*zeta3 + ei3*f_zet1*(1 - zeta4) +
719 ef3*f_zet1*zeta4;
720
721 d_eterm_dep = (12*(ef2 - ei2)*f_zeta*zeta2 +
722 4*ef2*f_zet1*zeta3 + 4*(ef2 - ei2)*f_zet1*zeta3 -
723 4*ei2*f_zet1*zeta3 + f_zet2*(ei2*(1 - zeta4) +
724 ef2*zeta4));
725
726
727 d_ctrm1_indep = 4*(d_ef2bi - d_ei2bi)*f_zeta*zeta3 +
728 d_ei2bi*f_zet1*(1 - zeta4) + d_ef2bi*f_zet1*zeta4;
729
730 d_ctrm1_dep = (12*(ef2bi - ei2bi)*f_zeta*zeta2 +
731 4*ef2bi*f_zet1*zeta3 + 4*(ef2bi - ei2bi)*f_zet1*zeta3 -
732 4*ei2bi*f_zet1*zeta3 + f_zet2*(ei2bi*(1 - zeta4) +
733 ef2bi*zeta4));
734
735
736 d_ctrm2_indep = d_bterm_indep*f_zet2 + (-d_ef0 + d_ei0 + ef2 - ei2)*
737 (12*f_zeta*zeta2 + 8*f_zet1*zeta3) +
738 d_ei0*f_zet2*(-1 + zeta4) - d_ef0*f_zet2*zeta4;
739
740 d_ctrm2_dep = d_bterm_dep*f_zet2 +
741 (4*(-ef0 + ei0)*f_zet2*zeta3 -
742 4*(ef0 - ef1 - ei0 + ei1)*(6*f_zeta*zeta +
743 9*f_zet1*zeta2 + 2*f_zet2*zeta3) +
744 f_zet3*(bterm + ei0*(-1 + zeta4) - ef0*zeta4));
745
746
747 d_cterm_indep = -d_delta_indep + d_bterm_indep*f_zet1 +
748 4*(ef2 - ei2)*f_zeta*zeta3;
749
750 d_cterm_dep = -d_delta_dep + d_bterm_dep*f_zet1 +
751 (bterm*f_zet2 + 12*(ef1 - ei1)*f_zeta*
752 zeta2 + 4*(ef1 - ei1)*f_zet1*zeta3);
753
754
755 d_dterm_indep = 12*(ef0 - ei0)*f_zeta*zeta2 + 12*d_ef0*f_zeta*rho*
756 zeta2 + 8*(ef0 - ei0)*f_zet1*zeta3 +
757 8*d_ef0*f_zet1*rho*zeta3 + d_ef0*f_zet2*rho*zeta4 +
758 f_zet2*(ei0 + ef0*zeta4 - ei0*zeta4) +
759 d_ei0*rho*(f_zet2 - 12*f_zeta*zeta2 - 8*f_zet1*zeta3 -
760 f_zet2*zeta4);
761
762 d_dterm_dep = (24*ef0*f_zeta*rho*zeta +
763 36*ef0*f_zet1*rho*zeta2 + 12*ef0*f_zet2*rho*zeta3 -
764 ei0*rho*(12*(2*f_zeta*zeta + 3*f_zet1*zeta2 +
765 f_zet2*zeta3) + f_zet3*(-1 + zeta4)) +
766 ef0*f_zet3*rho*zeta4);
767
768
769 d_dtrm1_indep = -(d_ei0*f_zet2) + ei2*f_zet2 +
770 (-d_ef0 + d_ei0 + ef2 - ei2)*(12*f_zeta*zeta2 +
771 8*f_zet1*zeta3) - d_ef0*f_zet2*zeta4 +
772 (d_ei0 + ef2 - ei2)*f_zet2*zeta4 -
773 dterm/rho2 + d_dterm_indep/rho;
774
775 d_dtrm1_dep = (4*(-ef0 + ef1 + ei0 - ei1)*f_zet2*zeta3 -
776 4*(ef0 - ef1 - ei0 + ei1)*(6*f_zeta*zeta +
777 9*f_zet1*zeta2 + 2*f_zet2*zeta3)
778 +
779 f_zet3*(ei1 + ei0*(-1 + zeta4) - (ef0 - ef1 + ei1)*zeta4))
780 + d_dterm_dep/rho;
781
782
783 d_dtrm2_indep = d_ei0*f_zet3*rho + 12*(ef0 - ei0)*
784 (2*f_zeta*zeta + 3*f_zet1*zeta2 + f_zet2*zeta3) +
785 12*(d_ef0 - d_ei0)*rho*(2*f_zeta*zeta + 3*f_zet1*zeta2 +
786 f_zet2*zeta3) + d_ef0*f_zet3*rho*zeta4 -
787 d_ei0*f_zet3*rho*zeta4 + f_zet3*(ei0 + ef0*zeta4 -
788 ei0*zeta4);
789
790 d_dtrm2_dep = (4*(ef0 - ei0)*f_zet3*rho*zeta3 +
791 12*(ef0 - ei0)*rho*(2*f_zeta + 8*f_zet1*zeta +
792 6*f_zet2*zeta2 + f_zet3*zeta3)
793 + f_zet4*rho*(ei0 + ef0*zeta4 - ei0*zeta4));
794
795 ds->df4000 += ( d_vcfp2_indep + 2*d_ctrm1_indep*spA + d_eterm_indep*spA +
796 d_vcfp2_dep*spA + 2*d_ctrm1_dep*spA*spA +
797 2*d_ctrm2_indep*spA*spA + d_dtrm1_indep*spA*spA +
798 d_eterm_dep*spA*spA + 2*d_ctrm2_dep*spA*spA*spA +
799 d_dtrm1_dep*spA*spA*spA + d_dtrm2_indep*spA*spA*spA +
800 d_dtrm2_dep*spA*spA*spA*spA + 2*d_cterm_indep*spAA +
801 2*d_cterm_dep*spA*spAA + 2*d_dterm_indep*spA*spAA +
802 2*d_dterm_dep*spA*spA*spAA + 2*spAAA*cterm +
803 2*spAA*ctrm1 + 4*spA*spAA*ctrm2 +
804 2*spAA*spAA*dterm + 2*spA*spAAA*dterm +
805 2*spA*spAA*dtrm1 + 3*spA*spA*spAA*dtrm2 +
806 spAA*eterm
807 )*factor;
808
809 ds->df3100 += ( d_vcfp2_indep + 2*d_ctrm1_indep*spA + d_eterm_indep*spA +
810 2*d_ctrm2_indep*spA*spA + d_dtrm1_indep*spA*spA +
811 d_dtrm2_indep*spA*spA*spA + 2*d_cterm_indep*spAA +
812 2*d_dterm_indep*spA*spAA + d_vcfp2_dep*spB +
813 2*d_ctrm1_dep*spA*spB + d_eterm_dep*spA*spB +
814 2*d_ctrm2_dep*spA*spA*spB + d_dtrm1_dep*spA*spA*spB +
815 d_dtrm2_dep*spA*spA*spA*spB + 2*d_cterm_dep*spAA*spB +
816 2*d_dterm_dep*spA*spAA*spB + 2*spAAB*cterm +
817 2*spAB*ctrm1 + 4*spA*spAB*ctrm2 +
818 2*spA*spAAB*dterm + 2*spAA*spAB*dterm +
819 2*spA*spAB*dtrm1 + 3*spA*spA*spAB*dtrm2 +
820 spAB*eterm
821 )*factor;
822
823 ds->df2200 += (d_vcfp2_indep + 2*d_ctrm1_indep*spA +
824 d_dtrm1_indep*spA*spA + 2*d_cterm_indep*spAB +
825 2*d_dterm_indep*spA*spAB + d_eterm_indep*spB +
826 d_vcfp2_dep*spB + 2*d_ctrm1_dep*spA*spB +
827 2*d_ctrm2_indep*spA*spB + d_dtrm1_dep*spA*spA*spB +
828 d_dtrm2_indep*spA*spA*spB + 2*d_cterm_dep*spAB*spB +
829 2*d_dterm_dep*spA*spAB*spB + d_eterm_dep*spB*spB +
830 2*d_ctrm2_dep*spA*spB*spB + d_dtrm2_dep*spA*spA*spB*spB +
831 2*spABB*cterm + 2*spAB*ctrm1 +
832 2*spAB*spB*ctrm2 + 2*spA*spBB*ctrm2 +
833 2*spAB*spAB*dterm + 2*spA*spABB*dterm +
834 2*spA*spAB*dtrm1 + 2*spA*spAB*spB*dtrm2 +
835 spA*spA*spBB*dtrm2 + spBB*eterm
836 )*factor;
837
838 ds->df1300 += ( d_vcfp2_indep + 2*d_ctrm1_indep*spB + d_eterm_indep*spB +
839 2*d_ctrm2_indep*spB*spB + d_dtrm1_indep*spB*spB +
840 d_dtrm2_indep*spB*spB*spB + 2*d_cterm_indep*spBB +
841 2*d_dterm_indep*spB*spBB + d_vcfp2_dep*spA +
842 2*d_ctrm1_dep*spB*spA + d_eterm_dep*spB*spA +
843 2*d_ctrm2_dep*spB*spB*spA + d_dtrm1_dep*spB*spB*spA +
844 d_dtrm2_dep*spB*spB*spB*spA + 2*d_cterm_dep*spBB*spA +
845 2*d_dterm_dep*spB*spBB*spA + 2*spABB*cterm +
846 2*spAB*ctrm1 + 4*spB*spAB*ctrm2 +
847 2*spB*spABB*dterm + 2*spBB*spAB*dterm +
848 2*spB*spAB*dtrm1 + 3*spB*spB*spAB*dtrm2 +
849 spAB*eterm
850 )*factor;
851
852 ds->df0400 += ( d_vcfp2_indep + 2*d_ctrm1_indep*spB + d_eterm_indep*spB +
853 d_vcfp2_dep*spB + 2*d_ctrm1_dep*spB*spB +
854 2*d_ctrm2_indep*spB*spB + d_dtrm1_indep*spB*spB +
855 d_eterm_dep*spB*spB + 2*d_ctrm2_dep*spB*spB*spB +
856 d_dtrm1_dep*spB*spB*spB + d_dtrm2_indep*spB*spB*spB +
857 d_dtrm2_dep*spB*spB*spB*spB + 2*d_cterm_indep*spBB +
858 2*d_cterm_dep*spB*spBB + 2*d_dterm_indep*spB*spBB +
859 2*d_dterm_dep*spB*spB*spBB + 2*spBBB*cterm +
860 2*spBB*ctrm1 + 4*spB*spBB*ctrm2 +
861 2*spBB*spBB*dterm + 2*spB*spBBB*dterm +
862 2*spB*spBB*dtrm1 + 3*spB*spB*spBB*dtrm2 +
863 spBB*eterm
864 )*factor;
865
866
867 /* the final section: end */
868 }
869
870 /* The dispatch part of the functional implementation */
871 static real
vwn3_energy(const FunDensProp * dp)872 vwn3_energy(const FunDensProp* dp)
873 {
874 return par_energy(dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
875 }
876
877 static void
vwn3_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)878 vwn3_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp)
879 {
880 par_first(ds, factor, dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
881 }
882
883 static void
vwn3_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)884 vwn3_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
885 {
886 par_second(ds, factor, dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
887 }
888
889 static void
vwn3_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)890 vwn3_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
891 {
892 par_third(ds, factor, dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
893 }
894
895 static void
vwn3_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)896 vwn3_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
897 {
898 par_fourth(ds, factor, dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
899 }
900
901
902 static real
vwn_energy(const FunDensProp * dp)903 vwn_energy(const FunDensProp* dp)
904 {
905 return par_energy(dp, &vwn_paramagnetic, &vwn_ferromagnetic);
906 }
907
908 static void
vwn_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)909 vwn_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp)
910 {
911 par_first(ds, factor, dp, &vwn_paramagnetic, &vwn_ferromagnetic);
912 }
913
914 static void
vwn_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)915 vwn_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
916 {
917 par_second(ds, factor, dp, &vwn_paramagnetic, &vwn_ferromagnetic);
918 }
919
920 static void
vwn_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)921 vwn_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
922 {
923 par_third(ds, factor, dp, &vwn_paramagnetic, &vwn_ferromagnetic);
924 }
925
926 static void
vwn_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)927 vwn_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
928 {
929 par_fourth(ds, factor, dp, &vwn_paramagnetic, &vwn_ferromagnetic);
930 }
931
932
933
934 /* Other spin interpolation scheme */
935 static real
spni_energy(const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)936 spni_energy(const FunDensProp* dp, const struct vwn_params* para,
937 const struct vwn_params* ferro)
938 {
939 real ep_p[2], ep_f[2], ep_i[2], zeta, zeta4, f_zeta, delta;
940 real rhoa = dp->rhoa, rhob = dp->rhob, rho;
941
942 if(rhoa<VWN_ZERO) rhoa = 1e-40;
943 if(rhob<VWN_ZERO) rhob = 1e-40;
944 rho = rhoa + rhob;
945 vwn_en_pot(ep_p, rho, 0, para);
946
947 if( fabs(dp->rhoa - dp->rhob)<VWN_ZERO) return ep_p[0]*rho;
948 vwn_en_pot(ep_f, rho, 0, ferro);
949 vwn_en_pot(ep_i, rho, 0, &vwn_interp);
950
951 zeta = (dp->rhoa-dp->rhob)/rho;
952 zeta4 = pow(zeta,4.0);
953 f_zeta = SPINPOLF*(pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
954 delta = f_zeta*(ep_f[0]-ep_p[0]);
955
956 return (ep_p[0]+ delta)*rho;
957 }
958
959 static void
spni_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)960 spni_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp,
961 const struct vwn_params* para, const struct vwn_params* ferro)
962 {
963 real zeta, f_zeta, f_zet1, vcfp;
964 real delta, ep_p[2], ep_f[2];
965 real rhoa = dp->rhoa, rhob = dp->rhob, rho;
966
967 if(rhoa<VWN_ZERO) rhoa = 1e-40;
968 if(rhob<VWN_ZERO) rhob = 1e-40;
969 rho = rhoa + rhob;
970 vwn_en_pot(ep_p, rho, 1, para);
971
972 ds->df1000 += ep_p[1]*factor;
973 ds->df0100 += ep_p[1]*factor;
974
975 /* if(dp->rhoa==dp->rhob) return; */
976
977 /* contribution from spin-polarized case; first order */
978 zeta = (dp->rhoa-dp->rhob)/rho;
979 f_zeta = SPINPOLF*(pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
980 f_zet1 = SPINPOLF*4.0/3.0*(pow(1+zeta,1.0/3.0)-pow(1-zeta,1.0/3.0));
981 vwn_en_pot(ep_f, rho, 1, ferro);
982
983 vcfp = f_zeta*(ep_f[1] - ep_p[1]);
984 delta= f_zet1*(ep_f[0] - ep_p[0]);
985
986 /* the final section: begin */
987 ds->df1000 += (vcfp + delta*(1-zeta))*factor;
988 ds->df0100 += (vcfp - delta*(1+zeta))*factor;
989 /* the final section: end */
990 }
991
992 static void
spni_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)993 spni_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp,
994 const struct vwn_params* para, const struct vwn_params* ferro)
995 {
996 real zeta, f_zeta, f_zet1, f_zet2, vcfp;
997 real delta, ep_p[3], ep_f[3];
998 real rhoa = dp->rhoa, rhob = dp->rhob, rho = dp->rhoa + dp->rhob;
999 real rho2 = rho*rho;
1000 real rho3 = rho*rho2;
1001 real vcf2, fac2, vap2, del2, ef0, ef1, ef2;
1002 real zA, zB, zAAr, zABr, zBBr;
1003
1004 vwn_en_pot(ep_p, rho, 2, para);
1005
1006 ds->df1000 += ep_p[1]*factor;
1007 ds->df0100 += ep_p[1]*factor;
1008 ds->df2000 += ep_p[2]*factor;
1009 ds->df1100 += ep_p[2]*factor;
1010 ds->df0200 += ep_p[2]*factor;
1011
1012 /* if( fabs(rhoa - rhob)<VWN_ZERO) return; */
1013 /* contribution from spin-polarized case; first order */
1014 zeta = (rhoa - rhob)/rho;
1015 f_zeta = SPINPOLF*(pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
1016 f_zet1 = SPINPOLF*4.0/3.0*(pow(1+zeta, 1.0/3.0)-pow(1-zeta, 1.0/3.0));
1017 f_zet2 = SPINPOLF*4.0/9.0*(pow(1+zeta,-2.0/3.0)+pow(1-zeta,-2.0/3.0));
1018 vwn_en_pot(ep_f, rho, 2, ferro);
1019
1020 ef0 = ep_f[0] - ep_p[0];
1021 ef1 = ep_f[1] - ep_p[1];
1022 ef2 = ep_f[2] - ep_p[2];
1023 vcfp = f_zeta*ef1;
1024 delta= f_zet1*ef0;
1025
1026 vcf2 = f_zeta*ef2;
1027 vap2 = f_zet1*ef1;
1028 fac2 = f_zet2*ef0*rho;
1029 zA = 2*rhob/rho2;
1030 zB = -2*rhoa/rho2;
1031 zAAr = -4*rhob/rho2;
1032 zABr = 2*zeta/rho;
1033 zBBr = 4*rhoa/rho2;
1034
1035 /* the final section: begin */
1036 ds->df1000 += (vcfp + delta*rho*zA)*factor;
1037 ds->df0100 += (vcfp - delta*(1+zeta))*factor;
1038 /* first order */
1039 //ds->df1000 += (f_zeta*ef1 + f_zet1*zA*rho*ef0)*factor;
1040 //ds->df0100 += (f_zeta*ef1 + f_zet1*zB*rho*ef0)*factor;
1041
1042 ds->df2000 += (vcf2 + vap2*(zA+zA) + fac2*zA*zA + delta*zAAr)*factor;
1043 ds->df1100 += (vcf2 + vap2*(zA+zB) +fac2*zA*zB + delta*zABr)*factor;
1044 ds->df0200 += (vcf2 + vap2*(zB+zB) + fac2*zB*zB + delta*zBBr)*factor;
1045 /* second order */
1046 //ds->df2000 += (f_zeta*ef2 + f_zet1*((zA+zA)*ef1 + zAAr*ef0) + f_zet2*zA*zA*rho*ef0)*factor;
1047 //ds->df1100 += (f_zeta*ef2 + f_zet1*((zA+zB)*ef1 + zABr*ef0) + f_zet2*zA*zB*rho*ef0)*factor;
1048 //ds->df0200 += (f_zeta*ef2 + f_zet1*((zB+zB)*ef1 + zBBr*ef0) + f_zet2*zB*zB*rho*ef0)*factor;
1049 /* the final section: end */
1050 }
1051
1052 static void
spni_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)1053 spni_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp,
1054 const struct vwn_params* para, const struct vwn_params* ferro)
1055 {
1056 real zeta, f_zeta, f_zet1, f_zet2, f_zet3;
1057 real delta, ep_p[4], ep_f[4];
1058 real ef0, ef1, ef2, ef3;
1059 real zA, zB, zAA, zAB, zBB, zAAA, zAAB, zABB, zBBB;
1060 real rhoa = dp->rhoa, rhob = dp->rhob, rho, rho2, rho3, rho4;
1061
1062 if(rhoa<VWN_ZERO) rhoa = VWN_ZERO;
1063 if(rhob<VWN_ZERO) rhob = VWN_ZERO;
1064
1065 rho = rhoa + rhob;
1066 rho2 = rho*rho;
1067 rho3 = rho2*rho;
1068 rho4 = rho3*rho;
1069 vwn_en_pot(ep_p, rho, 3, para);
1070
1071 ds->df1000 += ep_p[1]*factor;
1072 ds->df0100 += ep_p[1]*factor;
1073
1074 ds->df2000 += ep_p[2]*factor;
1075 ds->df1100 += ep_p[2]*factor;
1076 ds->df0200 += ep_p[2]*factor;
1077
1078 ds->df3000 += ep_p[3]*factor;
1079 ds->df2100 += ep_p[3]*factor;
1080 ds->df1200 += ep_p[3]*factor;
1081 ds->df0300 += ep_p[3]*factor;
1082
1083 zeta = (dp->rhoa-dp->rhob)/rho;
1084 f_zeta = SPINPOLF* (pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
1085 f_zet1 = SPINPOLF*4.0/3.0 *(pow(1+zeta, 1.0/3.0)-pow(1-zeta, 1.0/3.0));
1086 f_zet2 = SPINPOLF*4.0/9.0 *(pow(1+zeta,-2.0/3.0)+pow(1-zeta,-2.0/3.0));
1087 f_zet3 =-SPINPOLF*8.0/27.0*(pow(1+zeta,-5.0/3.0)-pow(1-zeta,-5.0/3.0));
1088
1089 vwn_en_pot(ep_f, rho, 3, ferro);
1090 ef0 = ep_f[0] - ep_p[0];
1091 ef1 = ep_f[1] - ep_p[1];
1092 ef2 = ep_f[2] - ep_p[2];
1093 ef3 = ep_f[3] - ep_p[3];
1094
1095 zA = 2*rhob/rho2; /* = 2(1-zeta)/rho */
1096 zB =-2*rhoa/rho2; /* = -2(1+zeta)/rho */
1097 zAA = -4*rhob/rho3;
1098 zAB = 2*(rhoa-rhob)/rho3;
1099 zBB = 4*rhoa/rho3;
1100 zAAA = 12*rhob/rho4;
1101 zAAB = zAAA - 4/rho3;
1102 zBBB = -12*rhoa/rho4;
1103 zABB = zBBB + 4/rho3;
1104
1105 float f0e1, f0e2, f0e3, f1e0, f1e1, f1e2, f2e0, f2e1, f3e0;
1106 f0e1 = f_zeta*ef1;
1107 f0e2 = f_zeta*ef2;
1108 f0e3 = f_zeta*ef3;
1109 f1e0 = f_zet1*ef0*rho;
1110 f1e1 = f_zet1*ef1;
1111 f1e2 = f_zet1*ef2;
1112 f2e0 = f_zet2*ef0*rho;
1113 f2e1 = f_zet2*ef1;
1114 f3e0 = f_zet3*ef0*rho;
1115
1116
1117 /* first order */
1118 ds->df1000 += (zA*f1e0 + f0e1)*factor;
1119 ds->df0100 += (zB*f1e0 + f0e1)*factor;
1120
1121 /* second order */
1122 ds->df2000 += (zAA*f1e0 + zA*zA*f2e0 + (zA+zA)*f1e1 + f0e2)*factor;
1123 ds->df1100 += (zAB*f1e0 + zA*zB*f2e0 + (zA+zB)*f1e1 + f0e2)*factor;
1124 ds->df0200 += (zBB*f1e0 + zB*zB*f2e0 + (zB+zB)*f1e1 + f0e2)*factor;
1125
1126 /* third order terms */
1127
1128 ds->df3000 += (
1129 f0e3 +
1130 zAAA*f1e0 + (zAA + zAA + zAA)*f1e1 + (zA + zA + zA)*f1e2 +
1131 (zAA*zA + zA*zAA + zAA*zA)*f2e0 + (zA*zA + zA*zA +zA*zA)*f2e1 +
1132 zA*zA*zA*f3e0
1133 )*factor;
1134 ds->df2100 += (
1135 f0e3 +
1136 zAAB*f1e0 + (zAA + zAB + zAB)*f1e1 + (zA + zA + zB)*f1e2 +
1137 (zAA*zB +zAB*zA + zA*zAB)*f2e0 + (zA*zA + zA*zB + zA*zB)*f2e1 +
1138 zA*zA*zB*f3e0
1139 )*factor;
1140 ds->df1200 += (
1141 f0e3 +
1142 zABB*f1e0 + (zBB + zAB + zAB)*f1e1 + (zB + zB + zA)*f1e2 +
1143 (zBB*zA +zAB*zB + zB*zAB)*f2e0 + (zB*zB + zB*zA + zB*zA)*f2e1 +
1144 zB*zB*zA*f3e0
1145 )*factor;
1146 ds->df0300 += (
1147 f0e3 +
1148 zBBB*f1e0 + (zBB+zBB+zBB)*f1e1 + (zB + zB + zB)*f1e2 +
1149 (zBB*zB + zBB*zB + zB*zBB)*f2e0 + (zB*zB+ (zB+zB)*zB)*f2e1 +
1150 zB*zB*zB*f3e0
1151 )*factor;
1152
1153 /* the final section: end */
1154 }
1155
1156 static void
spni_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp,const struct vwn_params * para,const struct vwn_params * ferro)1157 spni_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp,
1158 const struct vwn_params* para, const struct vwn_params* ferro)
1159 {
1160 real zeta, f_zeta, f_zet1, f_zet2, f_zet3, f_zet4;
1161 real delta, ep_p[5], ep_f[5];
1162 real ef0, ef1, ef2, ef3, ef4;
1163 real zA, zB, zAA, zAB, zBB, zAAA, zAAB, zABB, zBBB,
1164 zAAAA, zAAAB, zAABB, zABBB, zBBBB;
1165 real rhoa = dp->rhoa, rhob = dp->rhob, rho, rho2, rho3, rho4, rho5;
1166
1167 if(rhoa<VWN_ZERO) rhoa = VWN_ZERO;
1168 if(rhob<VWN_ZERO) rhob = VWN_ZERO;
1169
1170 rho = rhoa + rhob;
1171 rho2 = rho*rho;
1172 rho3 = rho2*rho;
1173 rho4 = rho3*rho;
1174 rho5 = rho4*rho;
1175 vwn_en_pot(ep_p, rho, 4, para);
1176
1177 ds->df1000 += ep_p[1]*factor;
1178 ds->df0100 += ep_p[1]*factor;
1179
1180 ds->df2000 += ep_p[2]*factor;
1181 ds->df1100 += ep_p[2]*factor;
1182 ds->df0200 += ep_p[2]*factor;
1183
1184 ds->df3000 += ep_p[3]*factor;
1185 ds->df2100 += ep_p[3]*factor;
1186 ds->df1200 += ep_p[3]*factor;
1187 ds->df0300 += ep_p[3]*factor;
1188
1189 ds->df4000 += ep_p[4]*factor;
1190 ds->df3100 += ep_p[4]*factor;
1191 ds->df2200 += ep_p[4]*factor;
1192 ds->df1300 += ep_p[4]*factor;
1193 ds->df0400 += ep_p[4]*factor;
1194
1195 zeta = (dp->rhoa-dp->rhob)/rho;
1196 f_zeta = SPINPOLF* (pow(1+zeta,FOURTHREE)+pow(1-zeta,FOURTHREE)-2.0);
1197 f_zet1 = SPINPOLF*4.0/3.0 *(pow(1+zeta, 1.0/3.0)-pow(1-zeta, 1.0/3.0));
1198 f_zet2 = SPINPOLF*4.0/9.0 *(pow(1+zeta,-2.0/3.0)+pow(1-zeta,-2.0/3.0));
1199 f_zet3 =-SPINPOLF*8.0/27.0*(pow(1+zeta,-5.0/3.0)-pow(1-zeta,-5.0/3.0));
1200 f_zet4 = SPINPOLF*40.0/81.0*(pow(1+zeta,-8.0/3.0)-pow(1-zeta,-8.0/3.0));
1201
1202 vwn_en_pot(ep_f, rho, 4, ferro);
1203 ef0 = ep_f[0] - ep_p[0];
1204 ef1 = ep_f[1] - ep_p[1];
1205 ef2 = ep_f[2] - ep_p[2];
1206 ef3 = ep_f[3] - ep_p[3];
1207 ef4 = ep_f[4] - ep_p[4];
1208
1209 zA = 2*rhob/rho2; /* = 2(1-zeta)/rho */
1210 zB =-2*rhoa/rho2; /* = -2(1+zeta)/rho */
1211 zAA = -4*rhob/rho3;
1212 zAB = 2*(rhoa-rhob)/rho3;
1213 zBB = 4*rhoa/rho3;
1214 zAAA = 12*rhob/rho4;
1215 zAAB = zAAA - 4/rho3;
1216 zBBB = -12*rhoa/rho4;
1217 zABB = zBBB + 4/rho3;
1218 zAAAA = -48*rhob/rho5;
1219 zBBBB = 48*rhoa/rho5;
1220 zAAAB = zAAAA + 12/rho4;
1221 zABBB = zBBBB - 12/rho4;
1222 zAABB = zAAAB + 12/rho4;
1223
1224 float f0e1, f0e2, f0e3, f0e4, f1e0, f1e1, f1e2, f1e3, f2e0, f2e1, f2e2, f3e0, f3e1, f4e0;
1225 f0e1 = f_zeta*ef1;
1226 f0e2 = f_zeta*ef2;
1227 f0e3 = f_zeta*ef3;
1228 f0e4 = f_zeta*ef4;
1229 f1e0 = f_zet1*ef0*rho;
1230 f1e1 = f_zet1*ef1;
1231 f1e2 = f_zet1*ef2;
1232 f1e3 = f_zet1*ef3;
1233 f2e0 = f_zet2*ef0*rho;
1234 f2e1 = f_zet2*ef1;
1235 f2e2 = f_zet2*ef2;
1236 f3e0 = f_zet3*ef0*rho;
1237 f3e1 = f_zet3*ef1;
1238 f4e0 = f_zet4*ef0*rho;
1239
1240
1241 /* first order */
1242 ds->df1000 += (zA*f1e0 + f0e1)*factor;
1243 ds->df0100 += (zB*f1e0 + f0e1)*factor;
1244
1245 /* second order */
1246 ds->df2000 += (zAA*f1e0 + zA*zA*f2e0 + (zA+zA)*f1e1 + f0e2)*factor;
1247 ds->df1100 += (zAB*f1e0 + zA*zB*f2e0 + (zA+zB)*f1e1 + f0e2)*factor;
1248 ds->df0200 += (zBB*f1e0 + zB*zB*f2e0 + (zB+zB)*f1e1 + f0e2)*factor;
1249
1250 /* third order terms */
1251
1252 ds->df3000 += (
1253 f0e3 +
1254 zAAA*f1e0 + (zAA + zAA + zAA)*f1e1 + (zA + zA + zA)*f1e2 +
1255 (zAA*zA + zA*zAA + zAA*zA)*f2e0 + (zA*zA + zA*zA + zA*zA)*f2e1 +
1256 zA*zA*zA*f3e0
1257 )*factor;
1258 ds->df2100 += (
1259 f0e3 +
1260 zAAB*f1e0 + (zAA + zAB + zAB)*f1e1 + (zA + zA + zB)*f1e2 +
1261 (zAA*zB + zAB*zA + zA*zAB)*f2e0 + (zA*zA + zA*zB + zA*zB)*f2e1 +
1262 zA*zA*zB*f3e0
1263 )*factor;
1264 ds->df1200 += (
1265 f0e3 +
1266 zABB*f1e0 + (zBB + zAB + zAB)*f1e1 + (zB + zB + zA)*f1e2 +
1267 (zBB*zA + zAB*zB + zB*zAB)*f2e0 + (zB*zB + zB*zA + zB*zA)*f2e1 +
1268 zB*zB*zA*f3e0
1269 )*factor;
1270 ds->df0300 += (
1271 f0e3 +
1272 zBBB*f1e0 + (zBB + zBB + zBB)*f1e1 + (zB + zB + zB)*f1e2 +
1273 (zBB*zB + zBB*zB + zB*zBB)*f2e0 + (zB*zB+ (zB+zB)*zB)*f2e1 +
1274 zB*zB*zB*f3e0
1275 )*factor;
1276
1277 ds->df4000 += (
1278 f1e0*zAAAA +
1279 f1e1*(4*zAAA) + f1e2*(6*zAA) + f1e3*(4*zA) +
1280 f2e0*(4*zA*zAAA) +
1281 f2e1*(12*zA*zAA + 3*zAA*zAA) +
1282 f2e2*(6*zA*zA) +
1283 f3e0*(6*zAA*zA*zA) +
1284 f3e1*(4*zA*zA*zA) +
1285 f4e0*zA*zA*zA*zA
1286 )*factor;
1287 ds->df3100 += (
1288 f1e0*zAAAB +
1289 f1e1*(3*zAAB + zAAA) + f1e2*(3*zAA + 3*zAB) + f1e3*(3*zA + zB) +
1290 f2e0*(3*zA*zAAB + zB*zAAA) +
1291 f2e1*(3*zA*(zAA + 2*zAB) + zB*(3*zAA) + 3*zAA*zAB) +
1292 f2e2*(3*zA*zA + 3*zA*zB) +
1293 f3e0*(3*zAA*zA*zB + 3*zAB*zA*zA) +
1294 f3e1*(zA*zA*zA + zA*zA*zB) +
1295 f4e0*zA*zA*zA*zB
1296 )*factor;
1297 ds->df2200 += (
1298 f1e0*zAABB +
1299 f1e1*(2*zAAB + 2*zABB) + f1e2*(zAA + 4*zAB + zBB) + f1e3*(2*zA + 2*zB) +
1300 f2e0*(2*zA*zABB + 2*zB*zAAB) +
1301 f2e1*(2*zA*(2*zAB + zBB) + 2*zB*(zAA + 2*zAB) + (zAA*zBB + 2*zAB*zAB)) +
1302 f2e2*(zA*zA + 4*zA*zB + zB*zB) +
1303 f3e0*(zAA*zB*zB + 4*zAB*zA*zB + zBB*zA*zA) +
1304 f3e1*(2*zA*zA*zB + 2*zA*zB*zB) +
1305 f4e0*zA*zA*zB*zB
1306 )*factor;
1307 ds->df1300 += (
1308 f1e0*zABBB +
1309 f1e1*(3*zABB + zBBB) + f1e2*(3*zAB + 3*zBB) + f1e3*(zA + 3*zB) +
1310 f2e0*(zA*zBBB + 3*zB*zABB) +
1311 f2e1*(3*zA*zBB + 3*zB*(2*zAB + zBB) + 3*zAB*zBB) +
1312 f2e2*(3*zA*zB + 3*zB*zB) +
1313 f3e0*(3*zAB*zB*zB + 3*zBB*zA*zB) +
1314 f3e1*(3*zA*zB*zB + zB*zB*zB) +
1315 f4e0*zA*zB*zB*zB
1316 )*factor;
1317 ds->df0400 += (
1318 f1e0*zBBBB +
1319 f1e1*(4*zBBB) + f1e2*(6*zBB) + f1e3*(4*zB) +
1320 f2e0*(4*zB*zBBB) +
1321 f2e1*(12*zB*zBB + 3*zBB*zBB) +
1322 f2e2*(6*zB*zB) +
1323 f3e0*(6*zBB*zB*zB) +
1324 f3e1*(4*zB*zB*zB) +
1325 f4e0*zB*zB*zB*zB
1326 )*factor;
1327 /* the final section: end */
1328 }
1329
1330 static real
vwni_energy(const FunDensProp * dp)1331 vwni_energy(const FunDensProp* dp)
1332 {
1333 return spni_energy(dp, &vwn_paramagnetic, &vwn_ferromagnetic);
1334 }
1335
1336 static void
vwni_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)1337 vwni_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp)
1338 {
1339 spni_first(ds, factor, dp, &vwn_paramagnetic, &vwn_ferromagnetic);
1340 }
1341
1342 static void
vwni_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)1343 vwni_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
1344 {
1345 spni_second(ds, factor, dp, &vwn_paramagnetic, &vwn_ferromagnetic);
1346 }
1347
1348 static void
vwni_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)1349 vwni_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
1350 {
1351 spni_third(ds, factor, dp, &vwn_paramagnetic, &vwn_ferromagnetic);
1352 }
1353
1354 static void
vwni_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)1355 vwni_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
1356 {
1357 spni_fourth(ds, factor, dp, &vwn_paramagnetic, &vwn_ferromagnetic);
1358 }
1359
1360 static real
vwn3i_energy(const FunDensProp * dp)1361 vwn3i_energy(const FunDensProp* dp)
1362 {
1363 return spni_energy(dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
1364 }
1365
1366 static void
vwn3i_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)1367 vwn3i_first(FunFirstFuncDrv *ds, real factor, const FunDensProp* dp)
1368 {
1369 spni_first(ds, factor, dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
1370 }
1371
1372 static void
vwn3i_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)1373 vwn3i_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
1374 {
1375 spni_second(ds, factor, dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
1376 }
1377
1378 static void
vwn3i_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)1379 vwn3i_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
1380 {
1381 spni_third(ds, factor, dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
1382 }
1383
1384 static void
vwn3i_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)1385 vwn3i_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
1386 {
1387 spni_fourth(ds, factor, dp, &vwn3_paramagnetic, &vwn3_ferromagnetic);
1388 }
1389