1 /*
2 * XCFun, an arbitrary order exchange-correlation library
3 * Copyright (C) 2020 Ulf Ekström and contributors.
4 *
5 * This file is part of XCFun.
6 *
7 * This Source Code Form is subject to the terms of the Mozilla Public
8 * License, v. 2.0. If a copy of the MPL was not distributed with this
9 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 *
11 * For information on the complete list of contributors to the
12 * XCFun library, see: <https://xcfun.readthedocs.io/>
13 */
14
15 #include "config.hpp"
16 #include "constants.hpp"
17 #include "functional.hpp"
18 #include "pw92eps.hpp"
19
20 /*! Common code for SCAN and SCAN-like functionals
21 *
22 * Implemented by James Furness,
23 *
24 * SCAN - J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402
25 * (2015) rSCAN - A. P. Bartok and J. R. Yates, J. Chem. Phys. 150, 161101 (2019)
26 * r2SCAN - J Furness, A Kaplan, J Ning, J Perdew, & J Sun; J. Chem. Phys. Lett.;
27 * Accepted (DOI: 10.1021/acs.jpclett.0c02405)
28 *
29 * The r++SCAN (aka "rppSCAN") and r4SCAN functionals are relatives of r2SCAN. They
30 * are not recommended for general use.
31 *
32 * r++SCAN, r4SCAN - J Furness, A Kaplan, J Ning, J Perdew, & J Sun; in preparation.
33 */
34
35 namespace SCAN_eps {
36 template <class num>
37 static num get_SCAN_Fx(const num,
38 const num,
39 const num,
40 const int,
41 const int,
42 const int);
43 template <class num>
44 static num SCAN_X_Fx(const num, const num, const parameter, const int, const int);
45 template <class num>
46 static num r2SCAN_C(const densvars<num> &, const int, const int, const int);
47 template <class num>
48 static num scan_ec0(const num,
49 const num,
50 const num,
51 const parameter,
52 const parameter,
53 const parameter);
54 template <class num>
55 static num lda_0(const num, const parameter, const parameter, const parameter);
56 template <class num>
57 static num scan_ec1(const num,
58 const num,
59 const num,
60 const parameter[8],
61 const parameter,
62 const parameter,
63 const parameter,
64 const parameter,
65 const int);
66 template <class num>
67 static void get_lsda1(const num, const num, const num, num &, num &);
68 template <class num>
69 static void gcor2(const parameter[6], const num, const num, num &, num &);
70
fx_unif(const num & d)71 template <class num> static num fx_unif(const num & d) {
72 return (-0.75 * pow(3 / PI, 1.0 / 3.0)) * pow(d, 4.0 / 3.0);
73 }
74
75 template <class num>
get_SCAN_Fx(const num d_n,const num d_g,const num d_tau,const int IALPHA,const int IINTERP,const int IDELFX)76 static num get_SCAN_Fx(const num d_n,
77 const num d_g,
78 const num d_tau,
79 const int IALPHA,
80 const int IINTERP,
81 const int IDELFX) {
82 const parameter ETA = 1.0e-3;
83 const parameter TAU_R = 1.0e-4;
84 const parameter A_REG = 1.0e-3;
85
86 num tauw = d_g / (8.0 * d_n);
87
88 num tauUnif = 0.0;
89 if (IALPHA == 1) {
90 tauUnif = ((3.0 / 10.0) * pow(3 * PI2, 2.0 / 3.0) * pow(d_n, 5.0 / 3.0)) + TAU_R;
91 } else {
92 tauUnif = (0.3 * pow(3 * PI2, 2.0 / 3.0) * pow(d_n, 5.0 / 3.0));
93 }
94
95 num alpha = 0.0;
96 if (IALPHA == 0) {
97 // alpha (SCAN)
98 if (abs(d_tau - tauw) > 1.0e-14) {
99 alpha = (d_tau - tauw) / tauUnif;
100 }
101
102 } else if (IALPHA == 1) {
103 // alpha' (rSCAN)
104 alpha = (d_tau - tauw) / tauUnif;
105 alpha = pow(alpha, 3) / (pow2(alpha) + A_REG);
106
107 } else if (IALPHA == 2) {
108 // \bar{alpha} (r2SCAN, r4SCAN)
109 if (abs(d_tau - tauw) > 1.0e-14) {
110 alpha = (d_tau - tauw) / (tauUnif + ETA * tauw);
111 }
112
113 } else {
114 // Unknown IALPHA choice
115 printf("ERROR: Unknown IALPHA %d\n", IALPHA);
116 }
117
118 num p = 0.0;
119 if (abs(d_g) > 1.0e-16) {
120 p = d_g /
121 (4.0 * pow(3.0 * PI2, 2.0 / 3.0) * pow(d_n, 8.0 / 3.0));
122 } else {
123 p = 1e-16/
124 (4.0 * pow(3.0 * PI2, 2.0 / 3.0) * pow(d_n, 8.0 / 3.0));
125 }
126
127 num Fx = SCAN_X_Fx(p, alpha, ETA, IINTERP, IDELFX);
128
129 return Fx;
130 }
131
132 template <class num>
SCAN_X_Fx(const num p,const num alpha,const parameter ETA,const int IINTERP,const int IDELFX)133 static num SCAN_X_Fx(const num p,
134 const num alpha,
135 const parameter ETA,
136 const int IINTERP,
137 const int IDELFX) {
138 const parameter A1 = 4.9479;
139 const parameter K1 = 0.065;
140 const parameter K0 = 0.174;
141 const parameter MU = 10.0 / 81.0;
142
143 const parameter IE_PARAMS[8] = {1.0,
144 -0.667,
145 -0.4445555,
146 -0.663086601049,
147 1.451297044490,
148 -0.887998041597,
149 0.234528941479,
150 -0.023185843322};
151 const parameter CFX1 = 0.667;
152 const parameter CFX2 = 0.8;
153 const parameter CFDX1 = 1.24;
154
155 const parameter D_DAMP2 = 0.361;
156 const parameter DX_DAMP4_P = 0.232;
157 const parameter DX_DAMP4_A = 0.232;
158 const parameter B1 = 0.156632;
159 const parameter B2 = 0.12083;
160 const parameter B3 = 0.5;
161 const parameter B4 = MU * MU / K1 - 0.112654;
162
163 const parameter ALPHA_GE = 20.0 / 27.0 + ETA * 5.0 / 3.0;
164
165 // Interpolation function
166 num oma = 1.0 - alpha;
167 num ief = 0.0;
168 if (IINTERP == 0) {
169 // SCAN
170 if (alpha < 1.0) {
171 ief = exp(-CFX1 * alpha / oma);
172 } else {
173 ief = -CFDX1 * exp(CFX2 / oma);
174 }
175 } else if (IINTERP == 1) {
176 // rSCAN
177 if (alpha < 1.0e-13) {
178 ief = exp(-CFX1 * alpha / oma);
179 } else if (alpha < 2.5) {
180 for (int i = 0; i < 8; i++) {
181 ief += IE_PARAMS[i] * pow(alpha, i);
182 }
183 } else {
184 ief = -CFDX1 * exp(CFX2 / oma);
185 }
186 } else {
187 printf("ERROR: Unknown IINTERP %d\n", IINTERP);
188 }
189
190 // Single orbital enhancement
191 num h0x = 1.0 + K0;
192
193 // Slowly varying enhancement
194 num del_f2 = 0.0;
195 num C2 = 0.0;
196 num h1x = 0.0;
197
198 if (IDELFX == 0) {
199 // 2nd and 4th order gradient expansion corrections for SCAN interpolation
200 num wfac = B4 * pow2(p) * exp(-B4 * p / MU);
201 num vfac = B1 * p + B2 * oma * exp(-B3 * pow2(oma));
202 num yfac = MU * p + wfac + pow2(vfac);
203
204 h1x = 1.0 + K1 - K1 / (1.0 + yfac / K1);
205 } else if (IDELFX == 1 || IDELFX == 2) {
206 // 2nd order GE corrections for rSCAN interpolation
207 for (int i = 1; i < 8; i++) {
208 del_f2 += i * IE_PARAMS[i];
209 }
210 C2 = -del_f2 * (1.0 - h0x);
211
212 num damp = exp(-pow2(p) / pow(D_DAMP2, 4));
213 h1x = 1.0 + K1 - K1 / (1.0 + p * (MU + ALPHA_GE * C2 * damp) / K1);
214 } else {
215 printf("ERROR: Unknown IDELFX %d\n", IDELFX);
216 }
217
218 // Scaling correction
219 num gx = 1.0 - exp(-A1 / pow(p, 1.0 / 4.0));
220
221 // 4th order gradient enhancement
222 num del_fx = 0.0;
223 if (IDELFX == 2) {
224 // 4th order correction for rSCAN interpolation (r4SCAN)
225 num eta_term = ETA * 3.0 / 4.0 + 2.0 / 3.0;
226
227 num del_f4 = 0.0;
228 for (int i = 1; i < 8; i++) {
229 del_f4 += i * (i - 1) * IE_PARAMS[i];
230 }
231
232 num C_aa = 73.0 / 5000.0 - 0.5 * del_f4 * (h0x - 1.0);
233 num C_pa = 511.0 / 13500 - 73.0 / 1500.0 * ETA - del_f2 * (ALPHA_GE * C2 + MU);
234 num C_pp = 146.0 / 2025.0 * pow2(eta_term) - 73.0 / 405.0 * eta_term +
235 pow2(ALPHA_GE * C2 + MU) / K1;
236
237 num order_1 = C2 * (oma - ALPHA_GE * p);
238 num t1 = order_1 + C_aa * pow2(oma) + C_pa * p * oma + C_pp * pow2(p);
239
240 num damp_4_t1 = 2.0 * pow2(alpha) / (1.0 + pow(alpha, 4));
241 num damp_4_t2 =
242 exp(-pow2(oma) / pow2(DX_DAMP4_A) - pow2(p) / pow(DX_DAMP4_P, 4));
243 num damp_4 = damp_4_t1 * damp_4_t2;
244
245 del_fx = t1 * damp_4;
246 }
247
248 num fx = (h1x + ief * (h0x - h1x) + del_fx) * gx;
249
250 return fx;
251 }
252
253 template <class num>
SCAN_C(const densvars<num> & d,const int IALPHA,const int IINTERP,const int IDELEC)254 static num SCAN_C(const densvars<num> & d,
255 const int IALPHA,
256 const int IINTERP,
257 const int IDELEC) {
258
259 const parameter CFC1 = 0.64;
260 const parameter CFC2 = 1.5;
261 const parameter CFDC1 = 0.7;
262 const parameter IE_PARAMS[8] = {1.0,
263 -0.64,
264 -0.4352,
265 -1.535685604549,
266 3.061560252175,
267 -1.915710236206,
268 0.516884468372,
269 -0.051848879792};
270
271 const parameter ETA = 1.0e-3;
272 const parameter TAU_R = 1.0e-4;
273 const parameter A_REG = 1.0e-3;
274
275 const parameter B1C = 0.0285764;
276 const parameter B2C = 0.0889;
277 const parameter B3C = 0.125541;
278
279 num rs = pow(4.0 * PI * d.n / 3.0, -(1.0 / 3.0));
280 num sqrtrs = 0.0;
281 if (abs(rs) > 1.0e-16) {
282 sqrtrs = sqrt(rs);
283 }
284
285 num ds_z = ufunc(d.zeta, 5.0 / 3.0) / 2.0;
286
287 num s = sqrt(d.gnn) / (2.0 * pow(3.0 * PI2, 1.0 / 3.0) * pow(d.n, 4.0 / 3.0));
288
289 num tueg_con = 3.0 / 10.0 * pow(3.0 * PI2, 2.0 / 3.0);
290 num tueg = 0.0;
291 if (IALPHA == 1) {
292 tueg = (tueg_con * pow(d.n, 5.0 / 3.0) + TAU_R) * ds_z;
293 } else {
294 tueg = tueg_con * pow(d.n, 5.0 / 3.0) * ds_z;
295 }
296
297 num tauw = d.gnn / (8.0 * d.n);
298
299 num alpha = 0.0;
300 if (IALPHA == 0) {
301 // alpha (SCAN)
302 if (abs(d.tau - tauw) > 1.0e-14) {
303 alpha = (d.tau - tauw) / tueg;
304 }
305
306 } else if (IALPHA == 1) {
307 // alpha' (rSCAN)
308 alpha = (d.tau - tauw) / tueg;
309 alpha = pow(alpha, 3) / (pow2(alpha) + A_REG);
310
311 } else if (IALPHA == 2) {
312 // \bar{alpha} (r2SCAN, r4SCAN)
313 if (abs(d.tau - tauw) > 1.0e-14) {
314 alpha = (d.tau - tauw) / (tueg + ETA * tauw);
315 }
316
317 } else {
318 // Unknown IALPHA choice
319 printf("ERROR: Unknown IALPHA %d\n", IALPHA);
320 }
321
322 // Interpolation function
323 num oma = 1.0 - alpha;
324 num ief = 0.0;
325 if (IINTERP == 0) {
326 // SCAN
327 if (alpha < 1.0) {
328 ief = exp(-CFC1 * alpha / oma);
329 } else {
330 ief = -CFDC1 * exp(CFC2 / oma);
331 }
332 } else if (IINTERP == 1) {
333 // rSCAN
334 if (alpha < 1.0e-13) {
335 ief = exp(-CFC1 * alpha / oma);
336 } else if (alpha < 2.5) {
337 for (int i = 0; i < 8; i++) {
338 ief += IE_PARAMS[i] * pow(alpha, i);
339 }
340 } else {
341 ief = -CFDC1 * exp(CFC2 / oma);
342 }
343 } else {
344 printf("ERROR: Unknown IINTERP %d\n", IINTERP);
345 }
346
347 num ec0 = scan_ec0(rs, s, d.zeta, B1C, B2C, B3C);
348 num ec1 = scan_ec1(rs, s, d.zeta, IE_PARAMS, ETA, B1C, B2C, B3C, IDELEC);
349
350 num eps_c = (ec1 + ief * (ec0 - ec1)) * d.n;
351
352 return eps_c;
353 }
354
355 template <class num>
scan_ec0(const num rs,const num s,const num zeta,const parameter B1C,const parameter B2C,const parameter B3C)356 static num scan_ec0(const num rs,
357 const num s,
358 const num zeta,
359 const parameter B1C,
360 const parameter B2C,
361 const parameter B3C) {
362 const parameter CHI_LD = 0.12802585262625815;
363
364 num eclda = lda_0(rs, B1C, B2C, B3C);
365
366 num dx_z = ufunc(zeta, 4.0 / 3.0) / 2.0;
367 num gc_z = (1.0 - 2.363 * (dx_z - 1.0)) * (1 - pow(zeta, 12));
368
369 num w0 = exp(-eclda / B1C) - 1.0;
370
371 num ginf = 1.0 / pow(1.0 + 4.0 * CHI_LD * s * s, 1.0 / 4.0);
372
373 num h0 = B1C * log(1.0 + w0 * (1.0 - ginf));
374
375 return (eclda + h0) * gc_z;
376 }
377
378 template <class num>
lda_0(const num rs,const parameter B1C,const parameter B2C,const parameter B3C)379 static num lda_0(const num rs,
380 const parameter B1C,
381 const parameter B2C,
382 const parameter B3C) {
383
384 return -B1C / (1.0 + B2C * sqrt(rs) + B3C * rs);
385 }
386
387 template <class num>
scan_ec1(const num rs,const num s,const num zeta,const parameter IE_PARAMS[8],const parameter ETA,const parameter B1C,const parameter B2C,const parameter B3C,const int IDELEC)388 static num scan_ec1(const num rs,
389 const num s,
390 const num zeta,
391 const parameter IE_PARAMS[8],
392 const parameter ETA,
393 const parameter B1C,
394 const parameter B2C,
395 const parameter B3C,
396 const int IDELEC) {
397 const parameter BETA_MB = 0.066725;
398 const parameter AFACTOR = 0.1;
399 const parameter BFACTOR = 0.1778;
400 const parameter GAMMA = 0.031090690869655;
401 const parameter AFIX_T = sqrt(PI / 4.0) * pow(9.0 * PI / 4.0, 1.0 / 6.0);
402 const parameter D_DAMP2 = 0.361;
403
404 num dx_z = ufunc(zeta, 4.0 / 3.0) / 2.0;
405 num gc_z = (1.0 - 2.363 * (dx_z - 1.0)) * (1 - pow(zeta, 12));
406 num phi = ufunc(zeta, 2.0 / 3.0) / 2.0;
407 num phi3 = pow(phi, 3);
408
409 num sqrtrs = sqrt(rs);
410
411 num eclda0 = lda_0(rs, B1C, B2C, B3C);
412
413 num eclsda1 = 0.0;
414 num d_eclsda1_drs = 0.0;
415 get_lsda1(rs, sqrtrs, zeta, eclsda1, d_eclsda1_drs);
416
417 num t = AFIX_T * s / (sqrtrs * phi);
418
419 num w1 = exp(-eclsda1 / (GAMMA * phi3)) - 1.0;
420
421 num beta = BETA_MB * (1.0 + AFACTOR * rs) / (1.0 + BFACTOR * rs);
422 num y = beta / (GAMMA * w1) * pow2(t);
423
424 num del_y = 0.0;
425 if (IDELEC == 0) {
426 // No correction for 2nd order GE in correlation (SCAN)
427 del_y = 0.0;
428 } else if (IDELEC == 1 || IDELEC == 2) {
429 // correcting terms for 2nd order GE with rSCAN interpolation
430 // Note that IDELEC = 1 is identical to 2
431
432 num p = s * s;
433 num ds_z = ufunc(zeta, 5.0 / 3.0) / 2.0;
434
435 num del_f2 = 0.0;
436 for (int i = 1; i < 8; i++) {
437 del_f2 += i * IE_PARAMS[i];
438 }
439
440 num eclsda0 = eclda0 * gc_z;
441 num d_eclsda0_drs = gc_z * (B3C + B2C / (2.0 * sqrtrs)) * pow2(eclda0) / B1C;
442
443 num t1 = del_f2 / (27.0 * GAMMA * ds_z * phi3 * w1);
444 num t2 = 20.0 * rs * (d_eclsda0_drs - d_eclsda1_drs);
445 num t3 = 45.0 * ETA * (eclsda0 - eclsda1);
446
447 num k = t1 * (t2 - t3);
448
449 num damp = exp(-pow2(p) / pow(D_DAMP2, 4));
450
451 del_y = k * p * damp;
452 } else {
453 printf("ERROR: Unrecognised IDELEC %d\n", IDELEC);
454 }
455
456 num g_y = 1.0 / pow(1.0 + 4.0 * (y - del_y), 1.0 / 4.0);
457
458 num h1 = GAMMA * phi3 * log(1.0 + w1 * (1.0 - g_y));
459
460 return eclsda1 + h1;
461 }
462
463 template <class num>
get_lsda1(const num rs,const num sqrtrs,const num zeta,num & eclda1,num & d_eclda1_drs)464 static void get_lsda1(const num rs,
465 const num sqrtrs,
466 const num zeta,
467 num & eclda1,
468 num & d_eclda1_drs) {
469 const parameter GAM = 0.51984209978974632953442121455650;
470 const parameter FZZ = 8.0 / (9.0 * GAM);
471 const parameter p_eu[6] = {
472 0.03109070, 0.213700, 7.59570, 3.58760, 1.63820, 0.492940};
473 const parameter p_ep[6] = {
474 0.015545350, 0.205480, 14.11890, 6.19770, 3.36620, 0.625170};
475 const parameter p_alfm[6] = {
476 0.01688690, 0.111250, 10.3570, 3.62310, 0.880260, 0.496710};
477
478 num eu = 0.0;
479 num deudrs = 0.0;
480 gcor2(p_eu, rs, sqrtrs, eu, deudrs);
481 num ep = 0.0;
482 num depdrs = 0.0;
483 gcor2(p_ep, rs, sqrtrs, ep, depdrs);
484 num alfm = 0.0;
485 num dalfmdrs = 0.0;
486 gcor2(p_alfm, rs, sqrtrs, alfm, dalfmdrs);
487
488 num z3 = pow(zeta, 3);
489 num z4 = zeta * z3;
490
491 num f = (ufunc(zeta, 4.0 / 3.0) - 2.0) / GAM;
492
493 eclda1 = eu * (1.0 - f * z4) + ep * f * z4 - alfm * f * (1.0 - z4) / FZZ;
494 d_eclda1_drs =
495 (1.0 - z4 * f) * deudrs + z4 * f * depdrs - (1.0 - z4) * f * dalfmdrs / FZZ;
496 return;
497 }
498
499 template <class num>
gcor2(const parameter P[6],const num rs,const num sqrtrs,num & GG,num & GGRS)500 static void gcor2(const parameter P[6],
501 const num rs,
502 const num sqrtrs,
503 num & GG,
504 num & GGRS) {
505 enum IDX { A, A1, B1, B2, B3, B4 };
506
507 num Q0 = -2.0 * P[A] * (1.0 + P[A1] * rs);
508 num Q0RS = -2.0 * P[A] * P[A1];
509
510 num Q1 = 2.0 * P[A] * sqrtrs *
511 (P[B1] + sqrtrs * (P[B2] + sqrtrs * (P[B3] + P[B4] * sqrtrs)));
512 num Q1RS = P[A] * (2.0 * P[B2] + P[B1] / sqrtrs + 3.0 * P[B3] * sqrtrs +
513 4.0 * P[B4] * rs);
514
515 num Q2 = log(1.0 + 1.0 / Q1);
516 num Q2RS = -Q1RS / ((1.0 + 1.0 / Q1) * pow2(Q1));
517
518 GG = Q0 * Q2;
519 GGRS = Q0 * Q2RS + Q2 * Q0RS;
520 return;
521 }
522 } // namespace SCAN_eps
523