1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // SCANFunctional.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18
19 #include "SCANFunctional.h"
20 #include <cmath>
21 #include <cassert>
22 #include <iostream>
23 #include <vector>
24 using namespace std;
25
SCANFunctional(const vector<vector<double>> & rhoe,double x_coeff,double c_coeff)26 SCANFunctional::SCANFunctional(const vector<vector<double> > &rhoe,
27 double x_coeff, double c_coeff)
28 {
29 x_coeff_ = x_coeff;
30 c_coeff_ = c_coeff;
31 _nspin = rhoe.size();
32 if ( _nspin > 1 ) assert(rhoe[0].size() == rhoe[1].size());
33 _np = rhoe[0].size();
34
35 if ( _nspin == 1 )
36 {
37 exc_.resize(_np);
38 vxc1_.resize(_np);
39 vxc2_.resize(_np);
40 vxc3_.resize(_np);
41 grad_rho_[0].resize(_np);
42 grad_rho_[1].resize(_np);
43 grad_rho_[2].resize(_np);
44 tau_.resize(_np);
45 rho = &rhoe[0][0];
46 grad_rho[0] = &grad_rho_[0][0];
47 grad_rho[1] = &grad_rho_[1][0];
48 grad_rho[2] = &grad_rho_[2][0];
49 tau = &tau_[0];
50 exc = &exc_[0];
51 vxc1 = &vxc1_[0];
52 vxc2 = &vxc2_[0];
53 vxc3 = &vxc3_[0];
54 }
55 else
56 {
57 exc_up_.resize(_np);
58 exc_dn_.resize(_np);
59 vxc1_up_.resize(_np);
60 vxc1_dn_.resize(_np);
61 vxc2_upup_.resize(_np);
62 vxc2_updn_.resize(_np);
63 vxc2_dnup_.resize(_np);
64 vxc2_dndn_.resize(_np);
65 vxc3_up_.resize(_np);
66 vxc3_dn_.resize(_np);
67 grad_rho_up_[0].resize(_np);
68 grad_rho_up_[1].resize(_np);
69 grad_rho_up_[2].resize(_np);
70 grad_rho_dn_[0].resize(_np);
71 grad_rho_dn_[1].resize(_np);
72 grad_rho_dn_[2].resize(_np);
73 tau_up_.resize(_np);
74 tau_dn_.resize(_np);
75
76 rho_up = &rhoe[0][0];
77 rho_dn = &rhoe[1][0];
78 grad_rho_up[0] = &grad_rho_up_[0][0];
79 grad_rho_up[1] = &grad_rho_up_[1][0];
80 grad_rho_up[2] = &grad_rho_up_[2][0];
81 grad_rho_dn[0] = &grad_rho_dn_[0][0];
82 grad_rho_dn[1] = &grad_rho_dn_[1][0];
83 grad_rho_dn[2] = &grad_rho_dn_[2][0];
84 tau_up = &tau_up_[0];
85 tau_dn = &tau_dn_[0];
86 exc_up = &exc_up_[0];
87 exc_dn = &exc_dn_[0];
88 vxc1_up = &vxc1_up_[0];
89 vxc1_dn = &vxc1_dn_[0];
90 vxc2_upup = &vxc2_upup_[0];
91 vxc2_updn = &vxc2_updn_[0];
92 vxc2_dnup = &vxc2_dnup_[0];
93 vxc2_dndn = &vxc2_dndn_[0];
94 vxc3_up = &vxc3_up_[0];
95 vxc3_dn = &vxc3_dn_[0];
96 }
97 }
98
setxc(void)99 void SCANFunctional::setxc(void)
100 {
101 if ( _np == 0 ) return;
102 if ( _nspin == 1 )
103 {
104 assert( rho != 0 );
105 assert( grad_rho[0] != 0 && grad_rho[1] != 0 && grad_rho[2] != 0 );
106 assert( tau != 0 );
107 assert( exc != 0 );
108 assert( vxc1 != 0 );
109 assert( vxc2 != 0 );
110 assert( vxc3 != 0 );
111 #pragma omp parallel for
112 for ( int i = 0; i < _np; i++ )
113 {
114 double grad = sqrt(grad_rho[0][i]*grad_rho[0][i] +
115 grad_rho[1][i]*grad_rho[1][i] +
116 grad_rho[2][i]*grad_rho[2][i]);
117 excSCAN(rho[i], grad, tau[i], &exc[i], &vxc1[i], &vxc2[i], &vxc3[i]);
118 }
119 }
120 else
121 {
122 assert( rho_up != 0 );
123 assert( rho_dn != 0 );
124 assert( grad_rho_up[0] != 0 && grad_rho_up[1] != 0 && grad_rho_up[2] != 0 );
125 assert( grad_rho_dn[0] != 0 && grad_rho_dn[1] != 0 && grad_rho_dn[2] != 0 );
126 assert( tau_up != 0 );
127 assert( tau_dn != 0 );
128 assert( exc_up != 0 );
129 assert( exc_dn != 0 );
130 assert( vxc1_up != 0 );
131 assert( vxc1_dn != 0 );
132 assert( vxc2_upup != 0 );
133 assert( vxc2_updn != 0 );
134 assert( vxc2_dnup != 0 );
135 assert( vxc2_dndn != 0 );
136 assert( vxc3_up != 0 );
137 assert( vxc3_dn != 0 );
138 #pragma omp parallel for
139 for ( int i = 0; i < _np; i++ )
140 {
141 double grx_up = grad_rho_up[0][i];
142 double gry_up = grad_rho_up[1][i];
143 double grz_up = grad_rho_up[2][i];
144 double grx_dn = grad_rho_dn[0][i];
145 double gry_dn = grad_rho_dn[1][i];
146 double grz_dn = grad_rho_dn[2][i];
147 double grx = grx_up + grx_dn;
148 double gry = gry_up + gry_dn;
149 double grz = grz_up + grz_dn;
150 double grad_up = sqrt(grx_up*grx_up + gry_up*gry_up + grz_up*grz_up);
151 double grad_dn = sqrt(grx_dn*grx_dn + gry_dn*gry_dn + grz_dn*grz_dn);
152 double grad = sqrt(grx*grx + gry*gry + grz*grz);
153 excSCAN_sp(rho_up[i],rho_dn[i],grad_up,grad_dn,grad,tau_up[i],tau_dn[i],
154 &exc_up[i],&exc_dn[i],&vxc1_up[i],&vxc1_dn[i],&vxc2_upup[i],
155 &vxc2_dndn[i],&vxc2_updn[i],&vxc2_dnup[i],&vxc3_up[i],
156 &vxc3_dn[i]);
157 }
158 }
159 }
160
161 ////////////////////////////////////////////////////////////////////////////////
162 //
163 // excSCAN: SCAN exchange-correlation
164 //
165 // input:
166 // rho: density
167 // grad: abs(grad(rho))
168 // tau: kinetic energy term
169 // output:
170 // exc: exchange-correlation energy per electron
171 // vxc1, vxc2 : quantities such that the total exchange potential is:
172 //
173 // vxc = vxc1 + div ( vxc2 * grad(n) )
174 //
175 // vxc3: additional potential energy term defined as change in energy with
176 // respect to kinetic energy tau
177 //
178 // References:
179 // [a] S. Jianwei, A. Ruzsinsky, and J.P.Perdew,
180 // "Strongly Constrained and Appropriately Normed Semilocal Density
181 // Functional", Phys.Rev.Lett. 115, 036402-1, (2015).
182 //
183 ////////////////////////////////////////////////////////////////////////////////
excSCAN(double rho,double grad,double tau,double * exc,double * vxc1,double * vxc2,double * vxc3)184 void SCANFunctional::excSCAN(double rho, double grad, double tau, double *exc,
185 double *vxc1, double *vxc2, double *vxc3)
186 {
187 const double pi = M_PI;
188 const double hx0 = 1.174;
189 const double a1 = 4.9479;
190 const double cx1 = 0.667;
191 const double cx2 = 0.8;
192 const double dx = 1.24;
193 const double k1 = 0.065;
194 const double mu_AK = 10.0/81.0;
195 const double b2 = sqrt(5913.0/405000.0);
196 const double b1 = 511.0/(13500.0 * 2.0 * b2);
197 const double b3 = 0.5;
198 const double b4 = (mu_AK * mu_AK / k1) - 1606.0/18225.0 - b1 * b1;
199 const double bc0 = 0.0285764;
200 const double bc1 = 0.0889;
201 const double bc2 = 0.125541;
202 const double cc1 = 0.64;
203 const double cc2 = 1.5;
204 const double dc = 0.7;
205 const double alpha0 = 0.2137;
206 const double beta00 = 0.0310907;
207 const double beta10 = 7.5957;
208 const double beta20 = 3.5876;
209 const double beta30 = 1.6382;
210 const double beta40 = 0.49294;
211 const double gamma = 0.03109069086965489; //gamma = (1-ln2)/pi^2
212 const double chi_inf = 0.12802585262625815;
213 const double z1 = 0.066724550603149220;
214 const double z2 = 0.1;
215 const double z3 = 0.1778;
216
217 double tmp0, tmp1, tmp2;
218
219 double rs, rtrs, kF, s, s2, tau_W, tau_unif, XCalpha, oneMalpha;
220 double exunif, x, gx, hx1, fx, FXSCAN;
221 double dxds, dxdalpha, dgxds, dhx1dx, dfxdalpha;
222 double dFXdalpha, dFXds;
223 double ex,vx1,vx2,vx3,ec,vc1,vc2,vc3;
224
225 double fc, ec0, ec1, H0, H1, w0, w1;
226 double beta1, ginf, A1, t1, g1, g5;
227 double ecLDA, decLDAdrs, ecLSDA, decLSDAdrs;
228
229 double dec0drs, dw0drs, dbeta1drs, dw1drs, dA1drs, dt1drs, dH1drs, dec1drs;
230 double dt1ds, dg1ds, dec1ds, dginfds, dec0ds;
231
232 double drsdn, dsdn, dalphadn, decdn;
233 double dsdgrad, dalphadgrad, decdgrad, dFXdsdsdgrad;
234 double dfcdalpha, dalphadtau, decdtau;
235
236 *exc = 0.0;
237 *vxc1 = 0.0;
238 *vxc2 = 0.0;
239 *vxc3 = 0.0;
240
241 if ( rho < 1.e-18 )
242 {
243 return;
244 }
245
246 rs = pow(4.0 * pi * rho / 3.0, -1.0/3.0);
247 rtrs = sqrt(rs);
248 kF = pow(3.0 * pi * pi * rho, 1.0/3.0);
249 s = grad / ( 2.0 * kF * rho );
250 s2 = s * s;
251 tau_W = grad * grad / (8.0 * rho);
252 tau_unif = 0.3 * pow(3.0 * pi * pi, 2.0 / 3.0) * pow(rho, 5.0 / 3.0);
253 //!! abs value in next line
254 XCalpha = fabs(tau - tau_W) / tau_unif;
255 oneMalpha = 1.0 - XCalpha;
256
257 // exchange
258
259 exunif = -3.0 / 4.0 * pow(3.0 * rho / pi, 1.0/3.0);
260
261 // SCAN exchange enhancement factor
262 tmp0 = (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha));
263 x = mu_AK * s2 * (1.0 + (b4 * s2 / mu_AK) * exp(-b4 * s2 / mu_AK)) +
264 tmp0 * tmp0;
265 gx = 1.0 - exp(-a1 / sqrt(s));
266 hx1 = 1.0 + k1 - k1/(1.0 + x / k1);
267 if ( XCalpha < 1.0 )
268 {
269 fx = exp(-cx1 * XCalpha / oneMalpha);
270 }
271 else if ( XCalpha > 1.0 )
272 {
273 fx = -dx * exp(cx2 / oneMalpha);
274 }
275 else
276 {
277 fx = 0.0;
278 }
279
280 FXSCAN = gx * (hx1 + fx * (hx0 - hx1));
281
282 //exchange energy
283 ex = exunif * FXSCAN;
284
285 // energy done, now the potential
286 dxds = 2.0 * s * mu_AK * (1.0 + 2.0 * (b4 * s2/ mu_AK) *
287 exp(-b4 * s2 / mu_AK) - (b4 * s2/ mu_AK) * (b4 * s2/ mu_AK) *
288 exp(-b4 * s2 / mu_AK)) + 4.0 * b1 * s *
289 (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha));
290 dxdalpha = (2.0 * b3 * oneMalpha * oneMalpha - 1.0) *
291 (2.0 * b2 * exp(-b3 * oneMalpha * oneMalpha)) *
292 (b1 * s2 + b2 * (oneMalpha) * exp(-b3 * oneMalpha * oneMalpha));
293 if ( s == 0 )
294 {
295 dgxds = 0;
296 }
297 else
298 {
299 dgxds = -a1 / (2.0 * pow(s , 1.5)) * exp(-a1 / sqrt(s));
300 }
301 dhx1dx = (k1 / (k1 + x)) * (k1 / (k1 + x));
302
303 if ( XCalpha < 1.0 )
304 {
305 dfxdalpha = -cx1 / oneMalpha / oneMalpha * exp(-cx1 * XCalpha / oneMalpha);
306 }
307 else if ( XCalpha > 1.0 )
308 {
309 dfxdalpha = -cx2 * dx / oneMalpha / oneMalpha * exp(cx2 / oneMalpha);
310 }
311 else
312 {
313 dfxdalpha = 0.0;
314 }
315
316 dsdn = -4.0 * s / (3.0 * rho);
317 dalphadn = 1.0 / rho * (tau_W / tau_unif - 5.0 / 3.0 * XCalpha);
318
319 dsdgrad = 1.0 / ( 2.0 * kF * rho ) / grad;
320 dalphadgrad = -1.0 / (4.0 * rho * tau_unif);
321 dalphadtau = 1.0 / tau_unif;
322
323 dFXds = dgxds * (hx1 + fx * (hx0 - hx1)) + gx * (1.0 - fx) * dhx1dx * dxds;
324 dFXdalpha = gx * (dhx1dx * dxdalpha + dfxdalpha * (hx0 - hx1) -
325 fx * dhx1dx * dxdalpha);
326
327 if ( s == 0 )
328 {
329 dFXdsdsdgrad = gx * (1.0 - fx) * dhx1dx *
330 (2.0 * mu_AK * (1.0 + 2.0 * (b4 * s2/ mu_AK) *
331 exp(-b4 * s2 / mu_AK) - (b4 * s2/ mu_AK) * (b4 * s2/ mu_AK) *
332 exp(-b4 * s2 / mu_AK)) + 4.0 * b1 *
333 (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha))) *
334 (1.0 / ( 2.0 * kF * rho ) / ( 2.0 * kF * rho ));
335 }
336 else
337 {
338 dFXdsdsdgrad = dFXds * dsdgrad;
339 }
340
341 //Note: 1/sigma from vx2 distributed to dsdgrad and dalphadgrad to
342 //prevent problems with 1/(grad=0)
343 vx1 = exunif * ((4.0/3.0) * FXSCAN + rho * (dFXds * dsdn +
344 dFXdalpha * dalphadn));
345 vx2 = -rho * exunif * (dFXdsdsdgrad + dFXdalpha * dalphadgrad);
346 vx3 = rho * exunif * dFXdalpha * dalphadtau;
347
348 // correlation
349 // SCAN Correlation energy
350 ecLDA = -bc0/(1.0 + bc1 * rtrs + bc2 * rs);
351 ginf = pow(1.0 + 4.0 * chi_inf * s2, -0.25);
352 w0 = exp(-ecLDA/bc0) - 1.0;
353 H0 = bc0 * log(1.0 + w0 * (1.0 - ginf));
354 ec0 = ecLDA + H0;
355
356 beta1 = z1 * ( 1.0 + z2 * rs) / (1.0 + z3 * rs);
357 gPW92(alpha0, beta00, beta10, beta20, beta30, beta40, rtrs,
358 &ecLSDA, &decLSDAdrs);
359 w1 = exp(-ecLSDA/gamma) - 1.0;
360 A1 = beta1 / (gamma * w1);
361 t1 = pow(3.0 * pi * pi / 16.0, 1.0/3.0) * s / rtrs;
362 g1 = pow(1.0 + 4.0 * A1 * t1 * t1, -0.25);
363 g5 = g1 * g1 * g1 * g1 * g1;
364 H1 = gamma * log(1.0 + w1 * (1.0 - g1));
365 ec1 = ecLSDA + H1;
366 if ( XCalpha < 1.0 )
367 {
368 fc = exp(-cc1 * XCalpha / oneMalpha);
369 }
370 else if ( XCalpha > 1.0 )
371 {
372 fc = -dc * exp(cc2 / oneMalpha);
373 }
374 else
375 {
376 fc = 0.0;
377 }
378
379 ec = ec1 + fc * (ec0 - ec1);
380
381 // rs derivatives
382 tmp1 = 1.0 + bc1 * rtrs + bc2 * rs;
383 decLDAdrs = bc0 * (bc1 / rtrs + 2.0 * bc2) / (2.0 * tmp1 * tmp1);
384 dw0drs = -1.0 * (1.0 + w0) / bc0 * decLDAdrs;
385 dec0drs = decLDAdrs + (bc0 * (1.0 - ginf))/(1.0 + w0 * (1.0 - ginf)) * dw0drs;
386
387 tmp2 = (1.0 + z3 * rs);
388 dbeta1drs = z1 * (z2 - z3) / tmp2 / tmp2;
389 dw1drs = - 1.0 * (1.0 + w1) * decLSDAdrs / gamma;
390 dA1drs = dbeta1drs / (gamma * w1) - beta1 * dw1drs / (gamma * w1 * w1);
391 dt1drs = -1.0 * pow(3.0 * pi * pi / 16.0, 1.0/3.0) * s / (2.0 * rtrs * rtrs *
392 rtrs);
393 dH1drs = (1.0 - g1) * gamma / (1.0 + w1 * (1.0 - g1)) * dw1drs + w1 * gamma /
394 (1.0 + w1 * (1.0 - g1)) * (t1 * t1 * g5 *
395 dA1drs + 2.0 * A1 * t1 * g5 * dt1drs);
396 dec1drs = decLSDAdrs + dH1drs;
397
398 // s derivatives
399
400 dginfds = -2.0 * chi_inf * s * ginf * ginf * ginf * ginf * ginf;
401 dec0ds = -bc0 * w0 * dginfds / (1.0 + w0 * (1.0 - ginf));
402
403 dt1ds = pow(3.0 * pi * pi / 16.0, 1.0/3.0) / rtrs;
404 dg1ds = -2.0 * A1 * t1 * g5 * dt1ds;
405 dec1ds = -gamma * w1 / (1.0 + w1 * (1.0 - g1)) * dg1ds;
406
407 // alpha derivatives
408 if ( XCalpha < 1.0 )
409 {
410 dfcdalpha = -cc1 / oneMalpha / oneMalpha * exp(-cc1 * XCalpha/ oneMalpha);
411 }
412 else if ( XCalpha > 1.0 )
413 {
414 dfcdalpha = -cc2 * dc / oneMalpha / oneMalpha * exp(cc2 / oneMalpha);
415 }
416 else
417 {
418 dfcdalpha = 0.0;
419 }
420
421 // V1
422 drsdn = -rs / (3.0 * rho);
423 dsdn = -4.0 * s / (3.0 * rho);
424 dalphadn = (tau_W / tau_unif - 5.0 * XCalpha / 3.0) / rho;
425 decdn = (dec1drs * drsdn + dec1ds * dsdn) +
426 dfcdalpha * dalphadn * (ec0 - ec1) + fc *
427 ((dec0drs * drsdn + dec0ds * dsdn) -
428 (dec1drs * drsdn + dec1ds * dsdn));
429 vc1 = ec + rho * decdn;
430
431 // V2
432 dsdgrad = 1.0 / ( 2.0 * kF * rho );
433 dalphadgrad = -2.0 * tau_W / (grad * tau_unif);
434 decdgrad = dec1ds * dsdgrad + dfcdalpha * dalphadgrad *
435 (ec0 - ec1) + fc * ( dec0ds * dsdgrad - dec1ds * dsdgrad );
436
437 // V2 grad=0
438 if ( s == 0 )
439 {
440 double dg1ds2 = -A1 * g5 * dt1ds /
441 (pow(16.0 * rho * rho * rho * rho,1.0/3.0) * rtrs);
442 double dec1ds2 = -gamma * w1 / (1.0 + w1 * (1.0 - g1)) * dg1ds2;
443 double dginfds2 = -chi_inf * ginf * ginf * ginf * ginf * ginf /
444 pow(3.0 * pi * pi * rho * rho * rho * rho,1.0/3.0);
445 double dec0ds2 = -bc0 * w0 * dginfds2 / (1.0 + w0 * (1.0 - ginf));
446 double dsdgrad2 = pow(24.0 * pi * pi * rho * rho * rho * rho,-1.0/3.0);
447 double dalphadgrad2 = -1.0 / (4.0 * rho * tau_unif);
448 double decdgrad2 = dec1ds2 * dsdgrad2 + dfcdalpha * dalphadgrad2 *
449 (ec0 - ec1) + fc * (dec0ds2 * dsdgrad2 - dec1ds2 * dsdgrad2);
450 vc2 = -rho * decdgrad2;
451 }
452 else
453 {
454 vc2 = -rho * decdgrad / grad;
455 }
456
457 // V3
458 dalphadtau = 1.0 / tau_unif;
459 decdtau = dfcdalpha * dalphadtau * (ec0 - ec1);
460 vc3 = rho * decdtau;
461
462 // XC
463 *exc = x_coeff_ * ex + c_coeff_ * ec;
464 *vxc1 = x_coeff_ * vx1 + c_coeff_ * vc1;
465 *vxc2 = x_coeff_ * vx2 + c_coeff_ * vc2;
466 *vxc3 = x_coeff_ * vx3 + c_coeff_ * vc3;
467 }
468
469 //////////////////////////////////////////////////////////////////////////////
excSCAN_sp(double rho_up,double rho_dn,double grad_up,double grad_dn,double grad,double tau_up,double tau_dn,double * exc_up,double * exc_dn,double * vxc1_up,double * vxc1_dn,double * vxc2_upup,double * vxc2_dndn,double * vxc2_updn,double * vxc2_dnup,double * vxc3_up,double * vxc3_dn)470 void SCANFunctional::excSCAN_sp(double rho_up, double rho_dn, double grad_up,
471 double grad_dn, double grad, double tau_up, double tau_dn, double *exc_up,
472 double *exc_dn, double *vxc1_up, double *vxc1_dn, double *vxc2_upup,
473 double *vxc2_dndn, double *vxc2_updn, double *vxc2_dnup, double *vxc3_up,
474 double *vxc3_dn)
475 {
476 const double pi = M_PI;
477 const double hx0 = 1.174;
478 const double a1 = 4.9479;
479 const double cx1 = 0.667;
480 const double cx2 = 0.8;
481 const double dx = 1.24;
482 const double k1 = 0.065;
483 const double mu_AK = 10.0/81.0;
484 const double b2 = sqrt(5913.0/405000.0);
485 const double b1 = 511.0/(13500.0 * 2.0 * b2);
486 const double b3 = 0.5;
487 const double b4 = (mu_AK * mu_AK / k1) - 1606.0/18225.0 - b1 * b1;
488 const double bc0 = 0.0285764;
489 const double bc1 = 0.0889;
490 const double bc2 = 0.125541;
491 const double cc1 = 0.64;
492 const double cc2 = 1.5;
493 const double dc = 0.7;
494 const double GC1 = 2.3631;
495 const double alpha0 = 0.2137;
496 const double beta00 = 0.0310907;
497 const double beta10 = 7.5957;
498 const double beta20 = 3.5876;
499 const double beta30 = 1.6382;
500 const double beta40 = 0.49294;
501 const double alpha1 = 0.20548;
502 const double beta01 = 0.01554535;
503 const double beta11 = 14.1189;
504 const double beta21 = 6.1977;
505 const double beta31 = 3.3662;
506 const double beta41 = 0.62517;
507 const double alpha2 = 0.11125;
508 const double beta02 = 0.0168869;
509 const double beta12 = 10.357;
510 const double beta22 = 3.6231;
511 const double beta32 = 0.88026;
512 const double beta42 = 0.49671;
513 const double gamma = 0.03109069086965489; //gamma = (1-ln2)/pi^2
514 const double gamma2 = 0.5198420997897463295344212145565;
515 const double chi_inf = 0.12802585262625815;
516 const double z1 = 0.066724550603149220;
517 const double z2 = 0.1;
518 const double z3 = 0.1778;
519
520 double tmp0, tmp1, tmp2;
521
522 double rs, rtrs, kF, s, s2, tau_W, tau_unif, ds, XCalpha, oneMalpha;
523 double exunif, x, gx, hx1, fx, FXSCAN;
524 double dxds, dxdalpha, dgxds, dhx1dx, dfxdalpha;
525 double dFXdalpha, dFXds, dFXdsdsdgrad;
526 double ex_up,ex_dn,vx1_up,vx1_dn,vx2_up,vx2_dn,vx3_up,vx3_dn,
527 ec,vc1_up,vc1_dn,vc2,vc3;
528 double rhotot, tautot = tau_up + tau_dn;
529
530 double fc, ec0, ec1, H0, H1, w0, w1, GC;
531 double Dx, phi, phi3, zeta, beta1, F, ginf, A1, t1, g1, g5;
532 double ecLDA, decLDAdrs, ecLSDA, ecLSDA0, ecLSDA1, ecLSDA2, decLSDAdrs,
533 decLSDAdrs0, decLSDAdrs1, decLSDAdrs2;
534
535 double dw0drs, dbeta1drs, dw1drs, dA1drs, dt1drs, dH1drs;
536 double dt1ds, dg1ds, dec1ds, dginfds, dH0ds, dec0ds;
537
538 double drsdn, dsdn, dalphadn;
539 double dsdgrad, dalphadgrad, decdgrad;
540 double dfcdalpha, dalphadtau, decdtau;
541 double dalphadzeta, dDxdzeta, dFdzeta, dw1dzeta, decLSDAdzeta,
542 dtdzeta, dAdzeta, dgdzeta, dH1dzeta, dphidzeta;
543 double dzetadnup, dec1dnup, dH0dnup, dGCdnup, dec0dnup, dfcdnup, decdnup;
544 double dzetadndn, dec1dndn, dH0dndn, dGCdndn, dec0dndn, dfcdndn, decdndn;
545
546 *exc_up = 0.0;
547 *exc_dn = 0.0;
548 *vxc1_up = 0.0;
549 *vxc1_dn = 0.0;
550 *vxc2_upup = 0.0;
551 *vxc2_updn = 0.0;
552 *vxc2_dnup = 0.0;
553 *vxc2_dndn = 0.0;
554 *vxc3_up = 0.0;
555 *vxc3_dn = 0.0;
556
557 if ( rho_up < 1.e-18 && rho_dn < 1.e-18 )
558 {
559 return;
560 }
561
562 ex_up = 0.0;
563 ex_dn = 0.0;
564 vx1_up = 0.0;
565 vx1_dn = 0.0;
566 vx2_up = 0.0;
567 vx2_dn = 0.0;
568 vx3_up = 0.0;
569 vx3_dn = 0.0;
570
571 //exchange up
572
573 if ( rho_up > 1.e-18 )
574 {
575 double tworhoup = 2.0 * rho_up;
576 double twogradup = 2.0 * grad_up;
577 rs = pow(4.0 * pi * tworhoup / 3.0, -1.0/3.0);
578 rtrs = sqrt(rs);
579 kF = pow(3.0 * pi * pi * tworhoup, 1.0/3.0);
580 s = twogradup / ( 2.0 * kF * tworhoup );
581 s2 = s * s;
582 tau_W = twogradup * twogradup / (8.0 * tworhoup);
583 tau_unif = 0.3 * pow(3.0 * pi * pi, 2.0 / 3.0) * pow(tworhoup, 5.0 / 3.0);
584 //!! abs value in next line
585 XCalpha = fabs(2.0 * tau_up - tau_W) / tau_unif;
586 oneMalpha = 1.0 - XCalpha;
587
588 exunif = -3.0 / 4.0 * pow(3.0 * tworhoup / pi, 1.0/3.0);
589
590 // SCAN exchange enhancement factor
591 tmp0 = (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha));
592 x = mu_AK * s2 * (1.0 + (b4 * s2 / mu_AK) * exp(-b4 * s2 / mu_AK)) +
593 tmp0 * tmp0;
594 gx = 1.0 - exp(-a1 / sqrt(s));
595 hx1 = 1.0 + k1 - k1/(1.0 + x / k1);
596 if ( XCalpha < 1.0 )
597 {
598 fx = exp(-cx1 * XCalpha / oneMalpha);
599 }
600 else if ( XCalpha > 1.0 )
601 {
602 fx = -dx * exp(cx2 / oneMalpha);
603 }
604 else
605 {
606 fx = 0.0;
607 }
608
609 FXSCAN = gx * (hx1 + fx * (hx0 - hx1));
610
611 //exchange energy
612 ex_up = exunif * FXSCAN;
613
614 // energy done, now the potential
615 dxds = 2.0 * s * mu_AK * (1.0 + 2.0 * (b4 * s2/ mu_AK) *
616 exp(-b4 * s2 / mu_AK) - (b4 * s2/ mu_AK) * (b4 * s2/ mu_AK) *
617 exp(-b4 * s2 / mu_AK)) + 4.0 * b1 * s *
618 (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha));
619 dxdalpha = (2.0 * b3 * oneMalpha * oneMalpha - 1.0) *
620 (2.0 * b2 * exp(-b3 * oneMalpha * oneMalpha)) *
621 (b1 * s2 + b2 * (oneMalpha) * exp(-b3 * oneMalpha * oneMalpha));
622 if ( s == 0 )
623 {
624 dgxds = 0;
625 }
626 else
627 {
628 dgxds = -a1 / (2.0 * pow(s , 1.5)) * exp(-a1 / sqrt(s));
629 }
630 dhx1dx = (k1 / (k1 + x)) * (k1 / (k1 + x));
631
632 if ( XCalpha < 1.0 )
633 {
634 dfxdalpha = -cx1 / oneMalpha / oneMalpha * exp(-cx1 * XCalpha / oneMalpha);
635 }
636 else if ( XCalpha > 1.0 )
637 {
638 dfxdalpha = -cx2 * dx / oneMalpha / oneMalpha * exp(cx2 / oneMalpha);
639 }
640 else
641 {
642 dfxdalpha = 0.0;
643 }
644
645 dsdn = -4.0 * s / (3.0 * tworhoup);
646 dalphadn = 1.0 / tworhoup * (tau_W / tau_unif - 5.0 / 3.0 * XCalpha);
647
648 dsdgrad = 1.0 / ( 2.0 * kF * tworhoup) / twogradup;
649 dalphadgrad = -1.0 / (4.0 * tworhoup * tau_unif);
650 dalphadtau = 1.0 / tau_unif;
651
652 dFXds = dgxds * (hx1 + fx * (hx0 - hx1)) + gx * (1.0 - fx) * dhx1dx * dxds;
653 dFXdalpha = gx * (dhx1dx * dxdalpha + dfxdalpha * (hx0 - hx1) -
654 fx * dhx1dx * dxdalpha);
655 if ( s == 0 )
656 {
657 dFXdsdsdgrad = gx * (1.0 - fx) * dhx1dx *
658 (2.0 * mu_AK * (1.0 + 2.0 * (b4 * s2/ mu_AK) *
659 exp(-b4 * s2 / mu_AK) - (b4 * s2/ mu_AK) * (b4 * s2/ mu_AK) *
660 exp(-b4 * s2 / mu_AK)) + 4.0 * b1 *
661 (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha))) *
662 (1.0 / ( 2.0 * kF * tworhoup ) / ( 2.0 * kF * tworhoup ));
663 }
664 else
665 {
666 dFXdsdsdgrad = dFXds * dsdgrad;
667 }
668
669 //Note: 1/sigma_up from vx2 distributed to dsdgrad and dalphadgrad to
670 //prevent problems with 1/(grad=0)
671 vx1_up = exunif * ((4.0/3.0) * FXSCAN + tworhoup * (dFXds * dsdn +
672 dFXdalpha * dalphadn));
673 vx2_up = -tworhoup * exunif * (dFXdsdsdgrad + dFXdalpha * dalphadgrad);
674 vx3_up = tworhoup * exunif * dFXdalpha * dalphadtau;
675 }
676
677 // exchange dn
678
679 if ( rho_dn > 1.e-18 )
680 {
681 double tworhodn = 2.0 * rho_dn;
682 double twograddn = 2.0 * grad_dn;
683 rs = pow(4.0 * pi * tworhodn / 3.0, -1.0/3.0);
684 rtrs = sqrt(rs);
685 kF = pow(3.0 * pi * pi * tworhodn, 1.0/3.0);
686 s = twograddn / ( 2.0 * kF * tworhodn );
687 s2 = s * s;
688 tau_W = twograddn * twograddn / (8.0 * tworhodn);
689 tau_unif = 0.3 * pow(3.0 * pi * pi, 2.0 / 3.0) * pow(tworhodn, 5.0 / 3.0);
690 //!! abs value in next line
691 XCalpha = fabs(2.0 * tau_dn - tau_W) / tau_unif;
692 oneMalpha = 1.0 - XCalpha;
693
694 exunif = -3.0 / 4.0 * pow(3.0 * tworhodn / pi, 1.0/3.0);
695
696 // SCAN exchange enhancement factor
697 tmp0 = (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha));
698 x = mu_AK * s2 * (1.0 + (b4 * s2 / mu_AK) * exp(-b4 * s2 / mu_AK)) +
699 tmp0 * tmp0;
700 gx = 1.0 - exp(-a1 / sqrt(s));
701 hx1 = 1.0 + k1 - k1/(1.0 + x / k1);
702 if ( XCalpha < 1.0 )
703 {
704 fx = exp(-cx1 * XCalpha / oneMalpha);
705 }
706 else if ( XCalpha > 1.0 )
707 {
708 fx = -dx * exp(cx2 / oneMalpha);
709 }
710 else
711 {
712 fx = 0.0;
713 }
714
715 FXSCAN = gx * (hx1 + fx * (hx0 - hx1));
716
717 //exchange energy
718 ex_dn = exunif * FXSCAN;
719
720 // energy done, now the potential
721 dxds = 2.0 * s * mu_AK * (1.0 + 2.0 * (b4 * s2/ mu_AK) *
722 exp(-b4 * s2 / mu_AK) - (b4 * s2/ mu_AK) * (b4 * s2/ mu_AK) *
723 exp(-b4 * s2 / mu_AK)) + 4.0 * b1 * s *
724 (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha));
725 dxdalpha = (2.0 * b3 * oneMalpha * oneMalpha - 1.0) *
726 (2.0 * b2 * exp(-b3 * oneMalpha * oneMalpha)) *
727 (b1 * s2 + b2 * (oneMalpha) * exp(-b3 * oneMalpha * oneMalpha));
728 if ( s == 0 )
729 {
730 dgxds = 0;
731 }
732 else
733 {
734 dgxds = -a1 / (2.0 * pow(s , 1.5)) * exp(-a1 / sqrt(s));
735 }
736 dhx1dx = (k1 / (k1 + x)) * (k1 / (k1 + x));
737
738 if ( XCalpha < 1.0 )
739 {
740 dfxdalpha = -cx1 / oneMalpha / oneMalpha * exp(-cx1 * XCalpha / oneMalpha);
741 }
742 else if ( XCalpha > 1.0 )
743 {
744 dfxdalpha = -cx2 * dx / oneMalpha / oneMalpha * exp(cx2 / oneMalpha);
745 }
746 else
747 {
748 dfxdalpha = 0.0;
749 }
750
751 dsdn = -4.0 * s / (3.0 * tworhodn);
752 dalphadn = 1.0 / tworhodn * (tau_W / tau_unif - 5.0 / 3.0 * XCalpha);
753
754 dsdgrad = 1.0 / ( 2.0 * kF * tworhodn ) / twograddn;
755 dalphadgrad = -1.0 / (4.0 * tworhodn * tau_unif);
756 dalphadtau = 1.0 / tau_unif;
757
758 dFXds = dgxds * (hx1 + fx * (hx0 - hx1)) + gx * (1.0 - fx) * dhx1dx * dxds;
759 dFXdalpha = gx * (dhx1dx * dxdalpha + dfxdalpha * (hx0 - hx1) -
760 fx * dhx1dx * dxdalpha);
761 if ( s == 0 )
762 {
763 dFXdsdsdgrad = gx * (1.0 - fx) * dhx1dx *
764 (2.0 * mu_AK * (1.0 + 2.0 * (b4 * s2/ mu_AK) *
765 exp(-b4 * s2 / mu_AK) - (b4 * s2/ mu_AK) * (b4 * s2/ mu_AK) *
766 exp(-b4 * s2 / mu_AK)) + 4.0 * b1 *
767 (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha))) *
768 (1.0 / ( 2.0 * kF * tworhodn ) / ( 2.0 * kF * tworhodn ));
769 }
770 else
771 {
772 dFXdsdsdgrad = dFXds * dsdgrad;
773 }
774
775 vx1_dn = exunif * ((4.0/3.0) * FXSCAN + tworhodn * (dFXds * dsdn +
776 dFXdalpha * dalphadn));
777 vx2_dn = -tworhodn * exunif * (dFXdsdsdgrad + dFXdalpha * dalphadgrad);
778 vx3_dn = tworhodn * exunif * dFXdalpha * dalphadtau;
779 }
780
781 // set negative densities to positive for correlation part
782 if ( rho_up < 0.0 ) rho_up = -rho_up;
783 if ( rho_dn < 0.0 ) rho_dn = -rho_dn;
784 // set small densities to 0 for correlation part
785 if (rho_up < 1.e-18 ) rho_up = 0.0;
786 if (rho_dn < 1.e-18 ) rho_dn = 0.0;
787
788 ec = 0.0;
789 vc1_up = 0.0;
790 vc1_dn = 0.0;
791 vc2 = 0.0;
792 vc3 = 0.0;
793
794 rhotot = rho_up + rho_dn;
795
796 // correlation
797 if ( rhotot > 1.e-18 )
798 {
799 zeta = (rho_up - rho_dn) / rhotot;
800 if ( zeta == 1.0 )
801 zeta = 1.0 - 1e-9;
802 if ( zeta == -1.0 )
803 zeta = -1.0 + 1e-9;
804 rs = pow(4.0 * pi * rhotot / 3.0, -1.0/3.0);
805 rtrs = sqrt(rs);
806 kF = pow(3.0 * pi * pi * rhotot, 1.0/3.0);
807 s = grad / ( 2.0 * kF * rhotot );
808 s2 = s * s;
809 ds = ((pow(1.0 + zeta, 5.0 / 3.0)) + (pow(1.0 - zeta, 5.0 / 3.0))) / 2.0;
810 tau_W = grad * grad / (8.0 * rhotot);
811 tau_unif = 0.3 * pow(3.0 * pi * pi, 2.0 / 3.0) * pow(rhotot, 5.0 / 3.0) *ds;
812 phi = ((pow(1.0 + zeta, 2.0 / 3.0)) + (pow(1.0 - zeta, 2.0 / 3.0))) / 2.0;
813 phi3 = phi * phi * phi;
814 //!! abs value in next line
815 XCalpha = fabs(tautot - tau_W) / tau_unif;
816 oneMalpha = 1.0 - XCalpha;
817
818 ecLDA = -bc0/(1.0 + bc1 * rtrs + bc2 * rs);
819 ginf = pow(1.0 + 4.0 * chi_inf * s2, -0.25);
820 w0 = exp(-ecLDA/bc0) - 1.0;
821 H0 = bc0 * log(1.0 + w0 * (1.0 - ginf));
822 Dx = ((pow(1.0 + zeta, 4.0 / 3.0)) + (pow(1.0 - zeta, 4.0 / 3.0))) / 2.0;
823 GC = (1.0 - GC1 * (Dx - 1.0)) * (1.0 - pow(zeta,12.0));
824 ec0 = (ecLDA + H0) * GC;
825 beta1 = z1 * ( 1.0 + z2 * rs) / (1.0 + z3 * rs);
826 F = (pow(1.0 + zeta,4.0/3.0) + pow(1.0 - zeta,4.0/3.0) - 2.0) / gamma2;
827 gPW92(alpha0, beta00, beta10, beta20, beta30, beta40, rtrs,
828 &ecLSDA0, &decLSDAdrs0);
829 gPW92(alpha1, beta01, beta11, beta21, beta31, beta41, rtrs,
830 &ecLSDA1, &decLSDAdrs1);
831 gPW92(alpha2, beta02, beta12, beta22, beta32, beta42, rtrs,
832 &ecLSDA2, &decLSDAdrs2);
833 ecLSDA = ecLSDA0 * (1.0 - F * pow(zeta,4.0)) + ecLSDA1 * F * pow(zeta,4.0)
834 - ecLSDA2 * F * (1.0 - pow(zeta,4.0))/(8.0/(9.0 * gamma2));
835 w1 = exp(-ecLSDA / (gamma * phi3)) - 1.0;
836 A1 = beta1 / (gamma * w1);
837 t1 = pow(3.0 * pi * pi / 16.0, 1.0/3.0) * s / (phi * rtrs);
838 g1 = pow(1.0 + 4.0 * A1 * t1 * t1, -0.25);
839 g5 = g1 * g1 * g1 * g1 * g1;
840 H1 = gamma * phi3 * log(1.0 + w1 * (1.0 - g1));
841 ec1 = ecLSDA + H1;
842 if ( XCalpha < 1.0 )
843 {
844 fc = exp(-cc1 * XCalpha / oneMalpha);
845 }
846 else if ( XCalpha > 1.0 )
847 {
848 fc = -dc * exp(cc2 / oneMalpha);
849 }
850 else
851 {
852 fc = 0.0;
853 }
854
855 ec = ec1 + fc * (ec0 - ec1);
856
857 // rs derivatives
858 tmp1 = 1.0 + bc1 * rtrs + bc2 * rs;
859 decLDAdrs = bc0 * (bc1 / rtrs + 2.0 * bc2) / (2.0 * tmp1 * tmp1);
860 dw0drs = -1.0 * (1.0 + w0) * decLDAdrs / bc0;
861
862 tmp2 = (1.0 + z3 * rs);
863 dbeta1drs = z1 * (z2 - z3) / tmp2 / tmp2;
864 decLSDAdrs = decLSDAdrs0 * (1.0 - F * pow(zeta,4.0)) + decLSDAdrs1 * F *
865 pow(zeta,4.0) - decLSDAdrs2 * F * (1.0 - pow(zeta,4.0)) /
866 (8.0 / (9.0 * gamma2));
867 dw1drs = - 1.0 * (1.0 + w1) * decLSDAdrs / (gamma * phi3);
868 dA1drs = dbeta1drs / (gamma * w1) - beta1 * dw1drs / (gamma * w1 * w1);
869 dt1drs = -1.0 * pow(3.0 * pi * pi / 16.0, 1.0/3.0) * s /
870 (2.0 * phi * rtrs * rtrs * rtrs);
871 dH1drs = (1.0 - g1) * gamma * phi3 * dw1drs / (1.0 + w1 * (1.0 - g1)) +
872 w1 * gamma * phi3 / (1.0 + w1 * (1.0 - g1)) *
873 (t1 * t1 * g5 * dA1drs + 2.0 * A1 * t1 * g5 * dt1drs);
874 //dec1drs = decLSDAdrs + dH1drs;
875
876 // s derivatives
877
878 dginfds = -2.0 * chi_inf * s * ginf * ginf * ginf * ginf * ginf;
879 dt1ds = pow(3.0 * pi * pi / 16.0, 1.0/3.0) / phi / rtrs;
880 dg1ds = -2.0 * A1 * t1 * g5 * dt1ds;
881 dec1ds = -gamma * phi3 * w1 / (1.0 + w1 * (1.0 - g1)) * dg1ds;
882 dH0ds = -bc0 * w0 * dginfds / (1.0 + w0 * (1.0 - ginf));
883 dec0ds = dH0ds * GC;
884
885 // alpha derivatives
886 if ( XCalpha < 1.0 )
887 {
888 dfcdalpha = -cc1 / oneMalpha / oneMalpha * exp(-cc1 * XCalpha/ oneMalpha);
889 }
890 else if ( XCalpha > 1.0 )
891 {
892 dfcdalpha = -cc2 * dc / oneMalpha / oneMalpha * exp(cc2 / oneMalpha);
893 }
894 else
895 {
896 dfcdalpha = 0.0;
897 }
898
899 // zeta derivatives
900 dphidzeta = (pow(1.0 + zeta, -1.0/3.0) - pow(1.0 - zeta, -1.0/3.0)) / 3.0;
901 dDxdzeta = 2.0 / 3.0 * (pow(1.0 + zeta, 1.0/3.0) -
902 pow(1.0 - zeta, 1.0/3.0));
903 dalphadzeta = 5.0 * XCalpha *
904 (pow(1.0 - zeta, 2.0/3.0) - pow(1.0 + zeta, 2.0/3.0)) /
905 (3.0 * (pow(1.0 - zeta, 5.0/3.0) + pow(1.0 + zeta, 5.0/3.0)));
906 dFdzeta = 4.0 * (pow(1.0 + zeta, 1.0/3.0) - pow(1.0 - zeta, 1.0/3.0)) /
907 (3.0 * gamma2);
908 decLSDAdzeta = ecLSDA0 * (-1.0 * dFdzeta * pow(zeta,4.0) -
909 4.0 * F * pow(zeta,3.0)) + ecLSDA1 * (dFdzeta*pow(zeta,4.) +
910 4.0 * F * pow(zeta,3.0)) - ecLSDA2 * (dFdzeta *
911 (1.0 - pow(zeta,4.0)) - 4.0 * F * pow(zeta,3.0)) /
912 (8.0 / (9.0 * gamma2));
913 dw1dzeta = (-decLSDAdzeta / (gamma * phi3) + ((3.0 * ecLSDA * dphidzeta) /
914 (gamma * phi * phi3))) * exp(-ecLSDA / (gamma * phi3));
915 dtdzeta = pow((3.0 * pi * pi) / 16.0, 1.0 / 3.0) * (-s * dphidzeta) /
916 (phi * phi * rtrs);
917 dAdzeta = (-beta1 * dw1dzeta)/(gamma * w1 * w1);
918 dgdzeta = (-dAdzeta * t1 * t1) / pow(1.0 + 4.0 * A1 * t1 * t1, 5.0 / 4.0) +
919 (-2.0 * A1 * t1 * dtdzeta) /
920 pow(1.0 + 4.0 * A1 * t1 * t1, 5.0 / 4.0);
921 dH1dzeta = 3.0 * dphidzeta * H1 / phi + (gamma * phi3 * (dw1dzeta *
922 (1.0 - g1))) / (1.0 + w1 * (1.0 - g1)) + (gamma * phi3 *
923 (w1 * (-dgdzeta))) / (1.0 + w1 * (1.0 - g1));
924
925 // n derivatives
926 drsdn = -rs / (3.0 * rhotot);
927 dsdn = -4.0 * s / (3.0 * rhotot);
928 dalphadn = (tau_W / tau_unif - 5.0 * XCalpha / 3.0) / rhotot;
929
930 // V1C_up
931 dzetadnup = 2.0 * rho_dn / rhotot / rhotot;
932 dec1dnup = decLSDAdrs * drsdn + decLSDAdzeta * dzetadnup +
933 dH1dzeta * dzetadnup + dH1drs * drsdn + dec1ds * dsdn;
934 dH0dnup = bc0 * ((1.0 - ginf) * dw0drs * drsdn - w0 * dginfds * dsdn) /
935 (1.0 + w0 * (1.0 - ginf));
936 dGCdnup = -dzetadnup * (GC1 * dDxdzeta * (1.0 - pow(zeta,12.0)) +
937 12.0 * (1.0 - GC1 * (Dx - 1.0)) * pow(zeta,11.0));
938 dec0dnup = (decLDAdrs * drsdn + dH0dnup) * GC + (ecLDA + H0) * dGCdnup;
939 dfcdnup = dfcdalpha * (dalphadn + dalphadzeta * dzetadnup);
940 decdnup = dec1dnup + dfcdnup * (ec0 - ec1) + fc * (dec0dnup - dec1dnup);
941 vc1_up = ec + rhotot * decdnup;
942
943 // V1C_dn
944 dzetadndn = -2.0 * rho_up / rhotot / rhotot;
945 dec1dndn = decLSDAdrs * drsdn + decLSDAdzeta * dzetadndn +
946 dH1dzeta * dzetadndn + dH1drs * drsdn + dec1ds * dsdn;
947 dH0dndn = bc0 * ((1.0 - ginf) * dw0drs * drsdn - w0 * dginfds * dsdn) /
948 (1.0 + w0 * (1.0 - ginf));
949 dGCdndn = -dzetadndn * (GC1 * dDxdzeta * (1.0 - pow(zeta,12.0)) +
950 12.0 * (1.0 - GC1 * (Dx - 1.0)) * pow(zeta,11.0));
951 dec0dndn = (decLDAdrs * drsdn + dH0dndn) * GC + (ecLDA + H0) * dGCdndn;
952 dfcdndn = dfcdalpha * (dalphadn + dalphadzeta * dzetadndn);
953 decdndn = dec1dndn + dfcdndn * (ec0 - ec1) + fc * (dec0dndn - dec1dndn);
954 vc1_dn = ec + rhotot * decdndn;
955
956 // VC2
957 dsdgrad = s / grad;
958 dalphadgrad = -2.0 * tau_W / (grad * tau_unif);
959 decdgrad = dec1ds * dsdgrad + dfcdalpha * dalphadgrad *
960 (ec0 - ec1) + fc * ( dec0ds * dsdgrad - dec1ds * dsdgrad );
961
962 // V2 grad=0
963 if ( s == 0 )
964 {
965 double dg1ds2 = -A1 * g5 * dt1ds /
966 (pow(16.0 * pow(rhotot,4.0),1.0/3.0) * phi * rtrs);
967 double dec1ds2 = -gamma * w1 / (1.0 + w1 * (1.0 - g1)) * dg1ds2;
968 double dginfds2 = -chi_inf * ginf * ginf * ginf * ginf * ginf /
969 pow(3.0 * pi * pi * pow(rhotot,4.0),1.0/3.0);
970 double dec0ds2 = -bc0 * w0 * dginfds2 * GC / (1.0 + w0 * (1.0 - ginf));
971 double dsdgrad2 = pow(24.0 * pi * pi * pow(rhotot,5.0),-1.0/3.0);
972 double dalphadgrad2 = -1.0 / (4.0 * rhotot * tau_unif);
973 double decdgrad2 = dec1ds2 * dsdgrad2 + dfcdalpha * dalphadgrad2 *
974 (ec0 - ec1) + fc * (dec0ds2 * dsdgrad2 - dec1ds2 * dsdgrad2);
975 vc2 = -rhotot * decdgrad2;
976 }
977 else
978 {
979 vc2 = -rhotot * decdgrad / grad;
980 }
981
982 // VC3
983 dalphadtau = 1.0 / tau_unif;
984 decdtau = dfcdalpha * dalphadtau * (ec0 - ec1);
985 vc3 = rhotot * decdtau;
986 }
987 *exc_up = x_coeff_ * ex_up + c_coeff_ * ec;
988 *exc_dn = x_coeff_ * ex_dn + c_coeff_ * ec;
989 *vxc1_up = x_coeff_ * vx1_up + c_coeff_ * vc1_up;
990 *vxc1_dn = x_coeff_ * vx1_dn + c_coeff_ * vc1_dn;
991 *vxc2_upup = x_coeff_ * 2.0 * vx2_up + c_coeff_ * vc2;
992 *vxc2_dndn = x_coeff_ * 2.0 * vx2_dn + c_coeff_ * vc2;
993 *vxc2_updn = c_coeff_ * vc2;
994 *vxc2_dnup = c_coeff_ * vc2;
995 *vxc3_up = x_coeff_ * vx3_up + c_coeff_ * vc3;
996 *vxc3_dn = x_coeff_ * vx3_dn + c_coeff_ * vc3;
997 }
998
999 ////////////////////////////////////////////////////////////////////////////////
1000 //
1001 // gPW92.c: Interpolate LSD correlation energy
1002 // as given by (10) of Perdew & Wang, Phys Rev B45 13244 (1992)
1003 // Translated into C Dec 9, 1996
1004 //
1005 ////////////////////////////////////////////////////////////////////////////////
1006
gPW92(double alpha,double beta0,double beta1,double beta2,double beta3,double beta4,double rtrs,double * gg,double * dgdrs)1007 void SCANFunctional::gPW92(double alpha, double beta0, double beta1,
1008 double beta2, double beta3, double beta4, double rtrs, double *gg,
1009 double *dgdrs)
1010 {
1011 double q0,q1,q2,q3;
1012 q0 = -2.0 * beta0 * ( 1.0 + alpha * rtrs * rtrs );
1013 q1 = 2.0 * beta0 * rtrs * ( beta1 + rtrs * ( beta2 + rtrs *
1014 ( beta3 + rtrs * beta4 ) ) );
1015 q2 = log ( 1.0 + 1.0 / q1 );
1016 q3 = beta0 * ( beta1 / rtrs + 2.0 * beta2 + rtrs *
1017 ( 3.0 * beta3 + 4.0 * beta4 * rtrs ));
1018 *gg = q0 * q2;
1019 *dgdrs = -2.0 * alpha * beta0 * q2 - q0 * q3 / ( q1 * ( 1.0 + q1 ));
1020 }
1021