Coverage for eminus/xc/lda_xc_ksdt.py: 100.00%
291 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-10-21 12:19 +0000
1# SPDX-FileCopyrightText: 2024 The eminus developers
2# SPDX-License-Identifier: Apache-2.0
3"""KSDT LDA exchange-correlation.
5Reference: Phys. Rev. Lett. 112, 076403.
6"""
8import dataclasses
9import functools
10import math
12import numpy as np
14from eminus import backend as xp
16# ### Temperature dependent coefficients ###
19@dataclasses.dataclass
20class Coefficients:
21 """Coefficients class to calculate temperature/theta dependent coefficients.
23 Reference: Phys. Rev. Lett. 112, 076403.
24 """
26 theta: float #: Reduced temperature.
27 omega: float
28 e1: float
29 e2: float
30 e3: float
31 e4: float
32 e5: float
33 d1: float
34 d2: float
35 d3: float
36 d4: float
37 d5: float
38 c1: float
39 c2: float
40 c3: float
41 b1: float
42 b2: float
43 b3: float
44 b4: float
45 a1: float = 0.75
46 a2: float = 3.04363
47 a3: float = -0.09227
48 a4: float = 1.7035
49 a5: float = 8.31051
50 a6: float = 5.1105
52 @functools.cached_property
53 def a0(self):
54 """Calculate a0."""
55 return 1 / (math.pi * (4 / (9 * math.pi)) ** (1 / 3))
57 @functools.cached_property
58 def b5(self):
59 """Calculate b5."""
60 return self.b3 * math.sqrt(3 / 2) * self.omega * (4 / (9 * math.pi)) ** (-1 / 3)
62 @functools.cached_property
63 def a(self):
64 """Calculate a."""
65 with np.errstate(divide="ignore"):
66 u = self.a0 * xp.where(self.theta > 0, xp.tanh(1 / self.theta), 1)
67 return u * _pade(self.theta, self.a1, self.a2, self.a3, self.a4, self.a5, self.a6)
69 @property
70 def dadtheta(self):
71 """Calculate da / dtheta."""
72 with np.errstate(divide="ignore", invalid="ignore"):
73 u = self.a0 * xp.where(self.theta > 0, xp.tanh(1 / self.theta), 1)
74 du = xp.where(self.theta > 0, (u**2 / self.a0 - self.a0) / self.theta**2, 0)
75 v, dv = _dpade(self.theta, self.a1, self.a2, self.a3, self.a4, self.a5, self.a6)
76 return du * v + u * dv
78 @functools.cached_property
79 def b(self):
80 """Calculate b."""
81 with np.errstate(divide="ignore"):
82 u = xp.where(self.theta > 0, xp.tanh(1 / xp.sqrt(self.theta)), 1)
83 return u * _pade(self.theta, self.b1, self.b2, 0, self.b3, self.b4, self.b5)
85 @property
86 def dbdtheta(self):
87 """Calculate db / dtheta."""
88 with np.errstate(divide="ignore", invalid="ignore"):
89 u = xp.where(self.theta > 0, xp.tanh(1 / xp.sqrt(self.theta)), 1)
90 du = xp.where(self.theta > 0, (u**2 - 1) / (2 * self.theta * xp.sqrt(self.theta)), 0)
91 v, dv = _dpade(self.theta, self.b1, self.b2, 0, self.b3, self.b4, self.b5)
92 return du * v + u * dv
94 @functools.cached_property
95 def c(self):
96 """Calculate c."""
97 with np.errstate(divide="ignore"):
98 exp = xp.where(self.theta > 0, xp.exp(-self.c3 / self.theta), 0)
99 return (self.c1 + self.c2 * exp) * self.e
101 @property
102 def dcdtheta(self):
103 """Calculate dc / dtheta."""
104 with np.errstate(divide="ignore"):
105 exp = xp.where(self.theta > 0, xp.exp(-self.c3 / self.theta), 0)
106 u = self.c1 + self.c2 * exp
107 with np.errstate(invalid="ignore"):
108 du = xp.where(self.theta > 0, self.c2 * self.c3 * exp / self.theta**2, 0)
109 v, dv = self.e, self.dedtheta
110 return du * v + u * dv
112 @functools.cached_property
113 def d(self):
114 """Calculate d."""
115 with np.errstate(divide="ignore"):
116 u = xp.where(self.theta > 0, xp.tanh(1 / xp.sqrt(self.theta)), 1)
117 return u * _pade(self.theta, self.d1, self.d2, 0, self.d3, self.d4, self.d5)
119 @property
120 def dddtheta(self):
121 """Calculate dd / dtheta."""
122 with np.errstate(divide="ignore", invalid="ignore"):
123 u = xp.where(self.theta > 0, xp.tanh(1 / xp.sqrt(self.theta)), 1)
124 du = xp.where(self.theta > 0, (u**2 - 1) / (2 * self.theta * xp.sqrt(self.theta)), 0)
125 v, dv = _dpade(self.theta, self.d1, self.d2, 0, self.d3, self.d4, self.d5)
126 return du * v + u * dv
128 @functools.cached_property
129 def e(self):
130 """Calculate e."""
131 with np.errstate(divide="ignore"):
132 u = xp.where(self.theta > 0, xp.tanh(1 / self.theta), 1)
133 return u * _pade(self.theta, self.e1, self.e2, 0, self.e3, self.e4, self.e5)
135 @functools.cached_property
136 def dedtheta(self):
137 """Calculate de / dtheta."""
138 with np.errstate(divide="ignore", invalid="ignore"):
139 u = xp.where(self.theta > 0, xp.tanh(1 / self.theta), 1)
140 du = xp.where(self.theta > 0, (u**2 - 1) / self.theta**2, 0)
141 v, dv = _dpade(self.theta, self.e1, self.e2, 0, self.e3, self.e4, self.e5)
142 return du * v + u * dv
145# ### Parameters ###
148@dataclasses.dataclass
149class Zeta0Coeffs(Coefficients):
150 """Coefficient class using the parameters for zeta=0.
152 Reference: Phys. Rev. Lett. 112, 076403.
153 """
155 omega: float = 1
156 b1: float = 0.283997
157 b2: float = 48.932154
158 b3: float = 0.370919
159 b4: float = 61.095357
160 c1: float = 0.870089
161 c2: float = 0.193077
162 c3: float = 2.414644
163 d1: float = 0.579824
164 d2: float = 94.537454
165 d3: float = 97.839603
166 d4: float = 59.939999
167 d5: float = 24.388037
168 e1: float = 0.212036
169 e2: float = 16.731249
170 e3: float = 28.485792
171 e4: float = 34.028876
172 e5: float = 17.235515
175@dataclasses.dataclass
176class Zeta1Coeffs(Coefficients):
177 """Coefficient class using the parameters for zeta=1.
179 Reference: Phys. Rev. Lett. 112, 076403.
180 """
182 omega: float = 2 ** (1 / 3)
183 b1: float = 0.329001
184 b2: float = 111.598308
185 b3: float = 0.537053
186 b4: float = 105.086663
187 c1: float = 0.84893
188 c2: float = 0.167952
189 c3: float = 0.08882
190 d1: float = 0.55133
191 d2: float = 180.213159
192 d3: float = 134.486231
193 d4: float = 103.861695
194 d5: float = 17.75071
195 e1: float = 0.153124
196 e2: float = 19.543945
197 e3: float = 43.400337
198 e4: float = 120.255145
199 e5: float = 15.662836
202@dataclasses.dataclass
203class PhiParams:
204 """Parameter class holding the spin-interpolation function parameters.
206 Reference: Phys. Rev. Lett. 112, 076403.
207 """
209 g1: float = 2 / 3
210 g2: float = -0.0139261
211 g3: float = 0.183208
212 lambda1: float = 1.064009
213 lambda2: float = 0.572565
216# ### Functional implementation ###
219def lda_xc_ksdt(
220 n, T=0, zeta0_coeffs=Zeta0Coeffs, zeta1_coeffs=Zeta1Coeffs, phi_params=PhiParams, **kwargs
221):
222 """KSDT exchange-correlation functional (spin-polarized).
224 Similar to the functional with the label LDA_XC_KSDT and ID 259 in Libxc, but the
225 implementations differ: https://gitlab.com/libxc/libxc/-/issues/525
226 Exchange and correlation cannot be separated.
228 Reference: Phys. Rev. Lett. 112, 076403.
230 Args:
231 n: Real-space electronic density.
233 Keyword Args:
234 T: Temperature in Hartree.
235 zeta0_coeffs: Coefficient class using the parameters for zeta=0.
236 zeta1_coeffs: Coefficient class using the parameters for zeta=1.
237 phi_params: Parameter class holding the spin-interpolation function parameters.
238 **kwargs: Throwaway arguments.
240 Returns:
241 KSDT exchange-correlation energy density and potential.
242 """
243 kwargs["zeta"] = xp.zeros_like(n)
244 exc, vxc, _ = lda_xc_ksdt_spin(
245 n,
246 T=T,
247 zeta0_coeffs=zeta0_coeffs,
248 zeta1_coeffs=zeta1_coeffs,
249 phi_params=phi_params,
250 **kwargs,
251 )
252 return exc, xp.stack([vxc[0]]), None
255def lda_xc_ksdt_spin(
256 n, zeta, T=0, zeta0_coeffs=Zeta0Coeffs, zeta1_coeffs=Zeta1Coeffs, phi_params=PhiParams, **kwargs
257):
258 """KSDT exchange-correlation functional (spin-polarized).
260 Similar to the functional with the label LDA_XC_KSDT and ID 259 in Libxc, but the
261 implementations differ: https://gitlab.com/libxc/libxc/-/issues/525
262 Exchange and correlation cannot be separated.
264 Reference: Phys. Rev. Lett. 112, 076403.
266 Args:
267 n: Real-space electronic density.
268 zeta: Relative spin polarization.
270 Keyword Args:
271 T: Temperature in Hartree.
272 zeta0_coeffs: Coefficient class using the parameters for zeta=0.
273 zeta1_coeffs: Coefficient class using the parameters for zeta=1.
274 phi_params: Parameter class holding the spin-interpolation function parameters.
275 **kwargs: Throwaway arguments.
277 Returns:
278 KSDT exchange-correlation energy density and potential.
279 """
280 # Calculate properties
281 n_up = (1 + zeta) * n / 2
282 rs = (3 / (4 * math.pi * n)) ** (1 / 3)
283 theta = _get_theta(T, n, zeta)
284 theta0 = _get_theta0(theta, zeta)
285 theta1 = _get_theta1(theta, zeta)
287 # Initialize parameters
288 # We need to calculate coefficients for specific theta
289 # Create a coefficients object for each theta with its corresponding parameters
290 phi_params = phi_params()
291 zeta0theta = zeta0_coeffs(theta)
292 zeta0theta0 = zeta0_coeffs(theta0)
293 zeta1theta = zeta1_coeffs(theta)
294 zeta1theta1 = zeta1_coeffs(theta1)
296 # Calculate fxc
297 fxc0 = _get_fxc_zeta(rs, zeta0theta0)
298 fxc1 = _get_fxc_zeta(rs, zeta1theta1)
299 phi = _get_phi(rs, theta0, zeta, phi_params)
300 fxc = fxc0 + (fxc1 - fxc0) * phi
302 # Generic derivatives
303 drsdn = -(6 ** (1 / 3)) * (1 / n) ** (1 / 3) / (6 * math.pi ** (1 / 3) * n)
304 dzetadn_up = -zeta / n**2 + 1 / n
305 dzetadn_dw = -zeta / n**2 - 1 / n
307 # fxc derivatives
308 dfxc0drs = _get_dfxc_zetadrs(rs, zeta0theta)
309 dfxc1drs = _get_dfxc_zetadrs(rs, zeta1theta)
310 dfxc0dtheta0 = _get_dfxc_zetadtheta(rs, zeta0theta0)
311 dfxc1dtheta1 = _get_dfxc_zetadtheta(rs, zeta1theta1)
313 # phi derivatives
314 dphidrs = _get_dphidrs(rs, theta0, zeta, phi_params)
315 dphidtheta = _get_dphidtheta(rs, theta0, zeta, phi_params)
316 dphidzeta = _get_dphidzeta(rs, theta0, zeta, phi_params)
318 # theta derivatives
319 dthetadn_up = _get_dthetadn_up(T, n_up)
320 dtheta0dtheta = _get_dtheta0dtheta(zeta)
321 dtheta0dzeta = _get_dtheta0dzeta(theta0, zeta)
322 dtheta1dtheta0 = _get_dtheta1dtheta0()
324 # Calculate vxc_up (using dndn_up=1)
325 dfxc0a = dfxc0drs * drsdn # * dndn_up = 1
326 dfxc1a = dfxc1drs * drsdn # * dndn_up = 1
327 dfxc0b_up = dfxc0dtheta0 * (dtheta0dtheta * dthetadn_up + dtheta0dzeta * dzetadn_up)
328 dfxc1b_up = (
329 dfxc1dtheta1 * dtheta1dtheta0 * (dtheta0dtheta * dthetadn_up + dtheta0dzeta * dzetadn_up)
330 )
331 dphi_up = dphidtheta * dthetadn_up + dphidzeta * dzetadn_up + dphidrs * drsdn # * dndn_up = 1
332 vxc_up = (
333 dfxc0a
334 + dfxc0b_up
335 - phi * (dfxc0a + dfxc0b_up - dfxc1a - dfxc1b_up)
336 - (fxc0 - fxc1) * dphi_up
337 )
339 # Calculate vxc_dw (using dndn_dw=1 and dthetadn_dw=0)
340 # dfxc0a and dfxc1a are identical for vxc_up and vxc_dw
341 dfxc0b_dw = dfxc0dtheta0 * dtheta0dzeta * dzetadn_dw
342 # dfxc0b += dfxc0dtheta0 * dtheta0dtheta * dthetadn_dw
343 dfxc1b_dw = dfxc1dtheta1 * dtheta1dtheta0 * dtheta0dzeta * dzetadn_dw
344 # dfxc1b += dfxc1dtheta1 * dtheta1dtheta0 * dtheta0dtheta * dthetadn_dw
345 dphi_dw = dphidzeta * dzetadn_dw + dphidrs * drsdn # * dndn_dw = 1
346 # dfxc1b += dphidtheta * dthetadn_dw
347 vxc_dw = (
348 dfxc0a
349 + dfxc0b_dw
350 - phi * (dfxc0a + dfxc0b_dw - dfxc1a - dfxc1b_dw)
351 - (fxc0 - fxc1) * dphi_dw
352 )
354 return fxc, fxc + xp.stack([vxc_up, vxc_dw]) * n, None
357# ### Pade approximation and derivative ###
360def _pade(x, n1, n2, n3, n4, d1, d2):
361 """Pade approximation.
363 Not the general case but as needed for this functional.
364 """
365 x2, x4 = x**2, x**4
366 num = n1 + n2 * x2 + n3 * x**3 + n4 * x4
367 denom = 1 + d1 * x2 + d2 * x4
368 return num / denom
371def _dpade(x, n1, n2, n3, n4, d1, d2):
372 """Pade approximation and its derivative."""
373 x2, x3, x4 = x**2, x**3, x**4
374 num = n1 + n2 * x2 + n3 * x3 + n4 * x4
375 denom = 1 + d1 * x2 + d2 * x4
376 dnum = 2 * n2 * x + 3 * n3 * x2 + 4 * n4 * x3
377 ddenom = 2 * d1 * x + 4 * d2 * x3
378 # df = (a'b - ab') / b^2
379 return num / denom, (dnum * denom - num * ddenom) / denom**2
382# ### theta and derivatives ###
385def _get_theta(T, n, zeta):
386 """Calculate the reduced temperature theta.
388 Reference: Phys. Rev. Lett. 119, 135001.
389 Only mentioned in the arXiv version: https://arxiv.org/abs/1703.08074
390 """
391 n_up = (1 + zeta) * n / 2
392 T_fermi = (6 * math.pi**2 * n_up) ** (2 / 3) / 2
393 return T / T_fermi
396def _get_dthetadn_up(T, n_up):
397 """Calculate dtheta / dn_up."""
398 return -4 / (3 * 6 ** (2 / 3)) * T / (math.pi ** (4 / 3) * n_up ** (5 / 3))
401def _get_theta0(theta, zeta):
402 """Calculate theta0.
404 Reference: Phys. Rev. Lett. 119, 135001.
405 """
406 return theta * (1 + zeta) ** (2 / 3)
409def _get_dtheta0dtheta(zeta):
410 """Calculate dtheta0 / dtheta."""
411 return (1 + zeta) ** (2 / 3)
414def _get_dtheta0dzeta(theta, zeta):
415 """Calculate dtheta0 / dzeta."""
416 return 2 * theta / (3 * (1 + zeta) ** (1 / 3))
419def _get_theta1(theta, zeta):
420 """Calculate theta1.
422 Reference: Phys. Rev. Lett. 119, 135001.
424 It is not explicitly mentioned but introduced as used in Eq. (5).
425 """
426 theta0 = _get_theta0(theta, zeta)
427 return theta0 * 2 ** (-2 / 3)
430def _get_dtheta1dtheta0():
431 """Calculate dtheta1 / dtheta0."""
432 return 2 ** (-2 / 3)
435# ### fxc_zeta and derivatives ###
438def _get_fxc_zeta(rs, p):
439 """Calculate the Pade formula f_xc^zeta.
441 Reference: Phys. Rev. Lett. 112, 076403.
442 """
443 p.theta = p.theta
444 num = p.a * p.omega + p.b * xp.sqrt(rs) + p.c * rs
445 denom = 1 + p.d * xp.sqrt(rs) + p.e * rs
446 return -1 / rs * num / denom
449def _get_dfxc_zetadrs(rs, p):
450 """Calculate dfxc_zeta / drs."""
451 num = p.a * p.omega + p.b * xp.sqrt(rs) + p.c * rs
452 denom = 1 + p.d * xp.sqrt(rs) + p.e * rs
453 tmp1 = (p.d / (2 * xp.sqrt(rs)) + p.e) * num / denom**2
454 tmp2 = (p.b / (2 * xp.sqrt(rs)) + p.c) / denom
455 return (tmp1 - tmp2) / rs + num / denom / rs**2
458def _get_dfxc_zetadtheta(rs, p):
459 """Calculate dfxc_zeta / dzeta."""
460 num = p.a * p.omega + p.b * xp.sqrt(rs) + p.c * rs
461 denom = 1 + p.d * xp.sqrt(rs) + p.e * rs
462 tmp1 = (p.dddtheta * xp.sqrt(rs) + p.dedtheta * rs) * num / denom**2
463 tmp2 = (p.dadtheta * p.omega + p.dbdtheta * xp.sqrt(rs) + p.dcdtheta * rs) / denom
464 return (tmp1 - tmp2) / rs
467# ### phi and derivatives ###
470def _get_phi(rs, theta, zeta, phi_params):
471 """Calculate the interpolation function phi.
473 Reference: Phys. Rev. Lett. 112, 076403.
474 """
475 alpha = _get_alpha(rs, theta, phi_params)
476 return ((1 + zeta) ** alpha + (1 - zeta) ** alpha - 2) / (2**alpha - 2)
479def _get_dphidrs(rs, theta, zeta, phi_params):
480 """Calculate dphi / drs."""
481 alpha = _get_alpha(rs, theta, phi_params)
482 dalphadrs = _get_dalphadrs(rs, theta, phi_params)
483 with np.errstate(divide="ignore"):
484 tmp1 = (1 - zeta) ** alpha * xp.where(1 - zeta > 0, xp.log(1 - zeta), 0)
485 tmp2 = (1 + zeta) ** alpha * xp.log(1 + zeta)
486 duv = (tmp1 + tmp2) * (2**alpha - 2)
487 udv = ((1 - zeta) ** alpha + (1 + zeta) ** alpha - 2) * 2**alpha * math.log(2)
488 vv = (2**alpha - 2) ** 2
489 return (duv - udv) * dalphadrs / vv
492def _get_dphidtheta(rs, theta, zeta, phi_params):
493 """Calculate dphi / dtheta."""
494 alpha = _get_alpha(rs, theta, phi_params)
495 dalphadtheta = _get_dalphadtheta(rs, theta, phi_params)
496 with np.errstate(divide="ignore"):
497 tmp1 = (1 - zeta) ** alpha * xp.where(1 - zeta > 0, xp.log(1 - zeta), 0)
498 tmp2 = (1 + zeta) ** alpha * xp.log(1 + zeta)
499 duv = (tmp1 + tmp2) * (2**alpha - 2)
500 udv = ((1 - zeta) ** alpha + (1 + zeta) ** alpha - 2) * 2**alpha * math.log(2)
501 vv = (2**alpha - 2) ** 2
502 return (duv - udv) * dalphadtheta / vv
505def _get_dphidzeta(rs, theta, zeta, phi_params):
506 """Calculate dphi / dzeta."""
507 alpha = _get_alpha(rs, theta, phi_params)
508 tmp1 = alpha * (1 + zeta) ** alpha / (1 + zeta)
509 with np.errstate(invalid="ignore"):
510 tmp2 = xp.where(1 - zeta > 0, alpha * (1 - zeta) ** alpha / (1 - zeta), 0)
511 return (tmp1 - tmp2) / (2**alpha - 2)
514# ### alpha and derivatives ###
517def _get_alpha(rs, theta, phi_params):
518 """Calculate alpha.
520 Reference: Phys. Rev. Lett. 112, 076403.
521 """
522 g = _get_g(rs, phi_params)
523 lamda = _get_lambda(rs, theta, phi_params)
524 return 2 - g * xp.exp(-theta * lamda)
527def _get_dalphadrs(rs, theta, phi_params):
528 """Calculate dalpha / drs."""
529 g = _get_g(rs, phi_params)
530 lamda = _get_lambda(rs, theta, phi_params)
531 dgdrs = _get_dgdrs(rs, phi_params)
532 dlambdadrs = _get_dlambdadrs(rs, theta, phi_params)
533 return -dgdrs * xp.exp(-theta * lamda) + dlambdadrs * theta * g * xp.exp(-theta * lamda)
536def _get_dalphadtheta(rs, theta, phi_params):
537 """Calculate dalpha / dtheta."""
538 g = _get_g(rs, phi_params)
539 lamda = _get_lambda(rs, theta, phi_params)
540 dlambdadtheta = _get_dlambdadtheta(rs, phi_params)
541 return (dlambdadtheta * theta + lamda) * g * xp.exp(-theta * lamda)
544# ### g and derivative ###
547def _get_g(rs, phi_params):
548 """Calculate g.
550 Reference: Phys. Rev. Lett. 112, 076403.
551 """
552 return (phi_params.g1 + phi_params.g2 * rs) / (1 + phi_params.g3 * rs)
555def _get_dgdrs(rs, phi_params):
556 """Calculate dg / drs."""
557 return (
558 phi_params.g2 / (phi_params.g3 * rs + 1)
559 - phi_params.g3 * (phi_params.g2 * rs + phi_params.g1) / (phi_params.g3 * rs + 1) ** 2
560 )
563# ### lambda and derivatives ###
566def _get_lambda(rs, theta, phi_params):
567 """Calculate lambda.
569 Reference: Phys. Rev. Lett. 112, 076403.
570 """
571 return phi_params.lambda1 + phi_params.lambda2 * theta * rs ** (1 / 2)
574def _get_dlambdadrs(rs, theta, phi_params):
575 """Calculate dlambda / drs."""
576 return phi_params.lambda2 * theta / (2 * xp.sqrt(rs))
579def _get_dlambdadtheta(rs, phi_params):
580 """Calculate dlambda / dtheta."""
581 return phi_params.lambda2 * xp.sqrt(rs)