Coverage for eminus/xc/lda_xc_ksdt.py: 100.00%
284 statements
« prev ^ index » next coverage.py v7.8.2, created at 2025-06-02 10:16 +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
11import numpy as np
13# ### Temperature dependent coefficients ###
16@dataclasses.dataclass
17class Coefficients:
18 """Coefficients class to calculate temperature/theta dependent coefficients.
20 Reference: Phys. Rev. Lett. 112, 076403.
21 """
23 theta: float #: Reduced temperature.
24 omega: float
25 e1: float
26 e2: float
27 e3: float
28 e4: float
29 e5: float
30 d1: float
31 d2: float
32 d3: float
33 d4: float
34 d5: float
35 c1: float
36 c2: float
37 c3: float
38 b1: float
39 b2: float
40 b3: float
41 b4: float
42 a1: float = 0.75
43 a2: float = 3.04363
44 a3: float = -0.09227
45 a4: float = 1.7035
46 a5: float = 8.31051
47 a6: float = 5.1105
49 @functools.cached_property
50 def a0(self):
51 """Calculate a0."""
52 return 1 / (np.pi * (4 / (9 * np.pi)) ** (1 / 3))
54 @functools.cached_property
55 def b5(self):
56 """Calculate b5."""
57 return self.b3 * np.sqrt(3 / 2) * self.omega * (4 / (9 * np.pi)) ** (-1 / 3)
59 @functools.cached_property
60 def a(self):
61 """Calculate a."""
62 with np.errstate(divide="ignore"):
63 u = self.a0 * np.tanh(
64 1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0
65 )
66 return u * _pade(self.theta, self.a1, self.a2, self.a3, self.a4, self.a5, self.a6)
68 @property
69 def dadtheta(self):
70 """Calculate da / dtheta."""
71 with np.errstate(divide="ignore"):
72 u = self.a0 * np.tanh(
73 1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0
74 )
75 du = np.divide(
76 u**2 / self.a0 - self.a0,
77 self.theta**2,
78 out=np.zeros_like(self.theta),
79 where=self.theta > 0,
80 )
81 v, dv = _dpade(self.theta, self.a1, self.a2, self.a3, self.a4, self.a5, self.a6)
82 return du * v + u * dv
84 @functools.cached_property
85 def b(self):
86 """Calculate b."""
87 with np.errstate(divide="ignore"):
88 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0)
89 return u * _pade(self.theta, self.b1, self.b2, 0, self.b3, self.b4, self.b5)
91 @property
92 def dbdtheta(self):
93 """Calculate db / dtheta."""
94 with np.errstate(divide="ignore"):
95 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0)
96 du = np.divide(
97 u**2 - 1,
98 2 * self.theta * np.sqrt(self.theta),
99 out=np.zeros_like(self.theta),
100 where=self.theta > 0,
101 )
102 v, dv = _dpade(self.theta, self.b1, self.b2, 0, self.b3, self.b4, self.b5)
103 return du * v + u * dv
105 @functools.cached_property
106 def c(self):
107 """Calculate c."""
108 with np.errstate(divide="ignore"):
109 exp = np.exp(-self.c3 / self.theta, out=np.zeros_like(self.theta), where=self.theta > 0)
110 return (self.c1 + self.c2 * exp) * self.e
112 @property
113 def dcdtheta(self):
114 """Calculate dc / dtheta."""
115 with np.errstate(divide="ignore"):
116 exp = np.exp(-self.c3 / self.theta, out=np.zeros_like(self.theta), where=self.theta > 0)
117 u = self.c1 + self.c2 * exp
118 du = np.divide(
119 self.c2 * self.c3 * exp,
120 self.theta**2,
121 out=np.zeros_like(self.theta),
122 where=self.theta > 0,
123 )
124 v, dv = self.e, self.dedtheta
125 return du * v + u * dv
127 @functools.cached_property
128 def d(self):
129 """Calculate d."""
130 with np.errstate(divide="ignore"):
131 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0)
132 return u * _pade(self.theta, self.d1, self.d2, 0, self.d3, self.d4, self.d5)
134 @property
135 def dddtheta(self):
136 """Calculate dd / dtheta."""
137 with np.errstate(divide="ignore"):
138 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0)
139 du = np.divide(
140 u**2 - 1,
141 2 * self.theta * np.sqrt(self.theta),
142 out=np.zeros_like(self.theta),
143 where=self.theta > 0,
144 )
145 v, dv = _dpade(self.theta, self.d1, self.d2, 0, self.d3, self.d4, self.d5)
146 return du * v + u * dv
148 @functools.cached_property
149 def e(self):
150 """Calculate e."""
151 with np.errstate(divide="ignore"):
152 u = np.tanh(1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0)
153 return u * _pade(self.theta, self.e1, self.e2, 0, self.e3, self.e4, self.e5)
155 @functools.cached_property
156 def dedtheta(self):
157 """Calculate de / dtheta."""
158 with np.errstate(divide="ignore"):
159 u = np.tanh(1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0)
160 du = np.divide(u**2 - 1, self.theta**2, out=np.zeros_like(self.theta), where=self.theta > 0)
161 v, dv = _dpade(self.theta, self.e1, self.e2, 0, self.e3, self.e4, self.e5)
162 return du * v + u * dv
165# ### Parameters ###
168@dataclasses.dataclass
169class Zeta0Coeffs(Coefficients):
170 """Coefficient class using the parameters for zeta=0.
172 Reference: Phys. Rev. Lett. 112, 076403.
173 """
175 omega: float = 1
176 b1: float = 0.283997
177 b2: float = 48.932154
178 b3: float = 0.370919
179 b4: float = 61.095357
180 c1: float = 0.870089
181 c2: float = 0.193077
182 c3: float = 2.414644
183 d1: float = 0.579824
184 d2: float = 94.537454
185 d3: float = 97.839603
186 d4: float = 59.939999
187 d5: float = 24.388037
188 e1: float = 0.212036
189 e2: float = 16.731249
190 e3: float = 28.485792
191 e4: float = 34.028876
192 e5: float = 17.235515
195@dataclasses.dataclass
196class Zeta1Coeffs(Coefficients):
197 """Coefficient class using the parameters for zeta=1.
199 Reference: Phys. Rev. Lett. 112, 076403.
200 """
202 omega: float = 2 ** (1 / 3)
203 b1: float = 0.329001
204 b2: float = 111.598308
205 b3: float = 0.537053
206 b4: float = 105.086663
207 c1: float = 0.84893
208 c2: float = 0.167952
209 c3: float = 0.08882
210 d1: float = 0.55133
211 d2: float = 180.213159
212 d3: float = 134.486231
213 d4: float = 103.861695
214 d5: float = 17.75071
215 e1: float = 0.153124
216 e2: float = 19.543945
217 e3: float = 43.400337
218 e4: float = 120.255145
219 e5: float = 15.662836
222@dataclasses.dataclass
223class PhiParams:
224 """Parameter class holding the spin-interpolation function parameters.
226 Reference: Phys. Rev. Lett. 112, 076403.
227 """
229 g1: float = 2 / 3
230 g2: float = -0.0139261
231 g3: float = 0.183208
232 lambda1: float = 1.064009
233 lambda2: float = 0.572565
236# ### Functional implementation ###
239def lda_xc_ksdt(
240 n, T=0, zeta0_coeffs=Zeta0Coeffs, zeta1_coeffs=Zeta1Coeffs, phi_params=PhiParams, **kwargs
241):
242 """KSDT exchange-correlation functional (spin-polarized).
244 Similar to the functional with the label LDA_XC_KSDT and ID 259 in Libxc, but the
245 implementations differ: https://gitlab.com/libxc/libxc/-/issues/525
246 Exchange and correlation cannot be separated.
248 Reference: Phys. Rev. Lett. 112, 076403.
250 Args:
251 n: Real-space electronic density.
253 Keyword Args:
254 T: Temperature in Hartree.
255 zeta0_coeffs: Coefficient class using the parameters for zeta=0.
256 zeta1_coeffs: Coefficient class using the parameters for zeta=1.
257 phi_params: Parameter class holding the spin-interpolation function parameters.
258 **kwargs: Throwaway arguments.
260 Returns:
261 KSDT exchange-correlation energy density and potential.
262 """
263 kwargs["zeta"] = np.zeros_like(n)
264 exc, vxc, _ = lda_xc_ksdt_spin(
265 n,
266 T=T,
267 zeta0_coeffs=zeta0_coeffs,
268 zeta1_coeffs=zeta1_coeffs,
269 phi_params=phi_params,
270 **kwargs,
271 )
272 return exc, np.array([vxc[0]]), None
275def lda_xc_ksdt_spin(
276 n, zeta, T=0, zeta0_coeffs=Zeta0Coeffs, zeta1_coeffs=Zeta1Coeffs, phi_params=PhiParams, **kwargs
277):
278 """KSDT exchange-correlation functional (spin-polarized).
280 Similar to the functional with the label LDA_XC_KSDT and ID 259 in Libxc, but the
281 implementations differ: https://gitlab.com/libxc/libxc/-/issues/525
282 Exchange and correlation cannot be separated.
284 Reference: Phys. Rev. Lett. 112, 076403.
286 Args:
287 n: Real-space electronic density.
288 zeta: Relative spin polarization.
290 Keyword Args:
291 T: Temperature in Hartree.
292 zeta0_coeffs: Coefficient class using the parameters for zeta=0.
293 zeta1_coeffs: Coefficient class using the parameters for zeta=1.
294 phi_params: Parameter class holding the spin-interpolation function parameters.
295 **kwargs: Throwaway arguments.
297 Returns:
298 KSDT exchange-correlation energy density and potential.
299 """
300 # Calculate properties
301 n_up = (1 + zeta) * n / 2
302 rs = (3 / (4 * np.pi * n)) ** (1 / 3)
303 theta = _get_theta(T, n, zeta)
304 theta0 = _get_theta0(theta, zeta)
305 theta1 = _get_theta1(theta, zeta)
307 # Initialize parameters
308 # We need to calculate coefficients for specific theta
309 # Create a coefficients object for each theta with its corresponding parameters
310 phi_params = phi_params()
311 zeta0theta = zeta0_coeffs(theta)
312 zeta0theta0 = zeta0_coeffs(theta0)
313 zeta1theta = zeta1_coeffs(theta)
314 zeta1theta1 = zeta1_coeffs(theta1)
316 # Calculate fxc
317 fxc0 = _get_fxc_zeta(rs, zeta0theta0)
318 fxc1 = _get_fxc_zeta(rs, zeta1theta1)
319 phi = _get_phi(rs, theta0, zeta, phi_params)
320 fxc = fxc0 + (fxc1 - fxc0) * phi
322 # Generic derivatives
323 drsdn = -(6 ** (1 / 3)) * (1 / n) ** (1 / 3) / (6 * np.pi ** (1 / 3) * n)
324 dzetadn_up = -zeta / n**2 + 1 / n
325 dzetadn_dw = -zeta / n**2 - 1 / n
327 # fxc derivatives
328 dfxc0drs = _get_dfxc_zetadrs(rs, zeta0theta)
329 dfxc1drs = _get_dfxc_zetadrs(rs, zeta1theta)
330 dfxc0dtheta0 = _get_dfxc_zetadtheta(rs, zeta0theta0)
331 dfxc1dtheta1 = _get_dfxc_zetadtheta(rs, zeta1theta1)
333 # phi derivatives
334 dphidrs = _get_dphidrs(rs, theta0, zeta, phi_params)
335 dphidtheta = _get_dphidtheta(rs, theta0, zeta, phi_params)
336 dphidzeta = _get_dphidzeta(rs, theta0, zeta, phi_params)
338 # theta derivatives
339 dthetadn_up = _get_dthetadn_up(T, n_up)
340 dtheta0dtheta = _get_dtheta0dtheta(zeta)
341 dtheta0dzeta = _get_dtheta0dzeta(theta0, zeta)
342 dtheta1dtheta0 = _get_dtheta1dtheta0()
344 # Calculate vxc_up (using dndn_up=1)
345 dfxc0a = dfxc0drs * drsdn # * dndn_up = 1
346 dfxc1a = dfxc1drs * drsdn # * dndn_up = 1
347 dfxc0b_up = dfxc0dtheta0 * (dtheta0dtheta * dthetadn_up + dtheta0dzeta * dzetadn_up)
348 dfxc1b_up = (
349 dfxc1dtheta1 * dtheta1dtheta0 * (dtheta0dtheta * dthetadn_up + dtheta0dzeta * dzetadn_up)
350 )
351 dphi_up = dphidtheta * dthetadn_up + dphidzeta * dzetadn_up + dphidrs * drsdn # * dndn_up = 1
352 vxc_up = (
353 dfxc0a
354 + dfxc0b_up
355 - phi * (dfxc0a + dfxc0b_up - dfxc1a - dfxc1b_up)
356 - (fxc0 - fxc1) * dphi_up
357 )
359 # Calculate vxc_dw (using dndn_dw=1 and dthetadn_dw=0)
360 # dfxc0a and dfxc1a are identical for vxc_up and vxc_dw
361 dfxc0b_dw = dfxc0dtheta0 * dtheta0dzeta * dzetadn_dw
362 # dfxc0b += dfxc0dtheta0 * dtheta0dtheta * dthetadn_dw
363 dfxc1b_dw = dfxc1dtheta1 * dtheta1dtheta0 * dtheta0dzeta * dzetadn_dw
364 # dfxc1b += dfxc1dtheta1 * dtheta1dtheta0 * dtheta0dtheta * dthetadn_dw
365 dphi_dw = dphidzeta * dzetadn_dw + dphidrs * drsdn # * dndn_dw = 1
366 # dfxc1b += dphidtheta * dthetadn_dw
367 vxc_dw = (
368 dfxc0a
369 + dfxc0b_dw
370 - phi * (dfxc0a + dfxc0b_dw - dfxc1a - dfxc1b_dw)
371 - (fxc0 - fxc1) * dphi_dw
372 )
374 return fxc, fxc + np.array([vxc_up, vxc_dw]) * n, None
377# ### Pade approximation and derivative ###
380def _pade(x, n1, n2, n3, n4, d1, d2):
381 """Pade approximation.
383 Not the general case but as needed for this functional.
384 """
385 x2, x4 = x**2, x**4
386 num = n1 + n2 * x2 + n3 * x**3 + n4 * x4
387 denom = 1 + d1 * x2 + d2 * x4
388 return num / denom
391def _dpade(x, n1, n2, n3, n4, d1, d2):
392 """Pade approximation and its derivative."""
393 x2, x3, x4 = x**2, x**3, x**4
394 num = n1 + n2 * x2 + n3 * x3 + n4 * x4
395 denom = 1 + d1 * x2 + d2 * x4
396 dnum = 2 * n2 * x + 3 * n3 * x2 + 4 * n4 * x3
397 ddenom = 2 * d1 * x + 4 * d2 * x3
398 # df = (a'b - ab') / b^2
399 return num / denom, (dnum * denom - num * ddenom) / denom**2
402# ### theta and derivatives ###
405def _get_theta(T, n, zeta):
406 """Calculate the reduced temperature theta.
408 Reference: Phys. Rev. Lett. 119, 135001.
409 Only mentioned in the arXiv version: https://arxiv.org/abs/1703.08074
410 """
411 n_up = (1 + zeta) * n / 2
412 T_fermi = (6 * np.pi**2 * n_up) ** (2 / 3) / 2
413 return T / T_fermi
416def _get_dthetadn_up(T, n_up):
417 """Calculate dtheta / dn_up."""
418 return -4 / (3 * 6 ** (2 / 3)) * T / (np.pi ** (4 / 3) * n_up ** (5 / 3))
421def _get_theta0(theta, zeta):
422 """Calculate theta0.
424 Reference: Phys. Rev. Lett. 119, 135001.
425 """
426 return theta * (1 + zeta) ** (2 / 3)
429def _get_dtheta0dtheta(zeta):
430 """Calculate dtheta0 / dtheta."""
431 return (1 + zeta) ** (2 / 3)
434def _get_dtheta0dzeta(theta, zeta):
435 """Calculate dtheta0 / dzeta."""
436 return 2 * theta / (3 * (1 + zeta) ** (1 / 3))
439def _get_theta1(theta, zeta):
440 """Calculate theta1.
442 Reference: Phys. Rev. Lett. 119, 135001.
444 It is not explicitly mentioned but introduced as used in Eq. (5).
445 """
446 theta0 = _get_theta0(theta, zeta)
447 return theta0 * 2 ** (-2 / 3)
450def _get_dtheta1dtheta0():
451 """Calculate dtheta1 / dtheta0."""
452 return 2 ** (-2 / 3)
455# ### fxc_zeta and derivatives ###
458def _get_fxc_zeta(rs, p):
459 """Calculate the Pade formula f_xc^zeta.
461 Reference: Phys. Rev. Lett. 112, 076403.
462 """
463 num = p.a * p.omega + p.b * np.sqrt(rs) + p.c * rs
464 denom = 1 + p.d * np.sqrt(rs) + p.e * rs
465 return -1 / rs * num / denom
468def _get_dfxc_zetadrs(rs, p):
469 """Calculate dfxc_zeta / drs."""
470 num = p.a * p.omega + p.b * np.sqrt(rs) + p.c * rs
471 denom = 1 + p.d * np.sqrt(rs) + p.e * rs
472 tmp1 = (p.d / (2 * np.sqrt(rs)) + p.e) * num / denom**2
473 tmp2 = (p.b / (2 * np.sqrt(rs)) + p.c) / denom
474 return (tmp1 - tmp2) / rs + num / denom / rs**2
477def _get_dfxc_zetadtheta(rs, p):
478 """Calculate dfxc_zeta / dzeta."""
479 num = p.a * p.omega + p.b * np.sqrt(rs) + p.c * rs
480 denom = 1 + p.d * np.sqrt(rs) + p.e * rs
481 tmp1 = (p.dddtheta * np.sqrt(rs) + p.dedtheta * rs) * num / denom**2
482 tmp2 = (p.dadtheta * p.omega + p.dbdtheta * np.sqrt(rs) + p.dcdtheta * rs) / denom
483 return (tmp1 - tmp2) / rs
486# ### phi and derivatives ###
489def _get_phi(rs, theta, zeta, phi_params):
490 """Calculate the interpolation function phi.
492 Reference: Phys. Rev. Lett. 112, 076403.
493 """
494 alpha = _get_alpha(rs, theta, phi_params)
495 return ((1 + zeta) ** alpha + (1 - zeta) ** alpha - 2) / (2**alpha - 2)
498def _get_dphidrs(rs, theta, zeta, phi_params):
499 """Calculate dphi / drs."""
500 alpha = _get_alpha(rs, theta, phi_params)
501 dalphadrs = _get_dalphadrs(rs, theta, phi_params)
502 tmp1 = (1 - zeta) ** alpha * np.log(1 - zeta, out=np.zeros_like(zeta), where=1 - zeta > 0)
503 tmp2 = (1 + zeta) ** alpha * np.log(1 + zeta)
504 duv = (tmp1 + tmp2) * (2**alpha - 2)
505 udv = ((1 - zeta) ** alpha + (1 + zeta) ** alpha - 2) * 2**alpha * np.log(2)
506 vv = (2**alpha - 2) ** 2
507 return (duv - udv) * dalphadrs / vv
510def _get_dphidtheta(rs, theta, zeta, phi_params):
511 """Calculate dphi / dtheta."""
512 alpha = _get_alpha(rs, theta, phi_params)
513 dalphadtheta = _get_dalphadtheta(rs, theta, phi_params)
514 tmp1 = (1 - zeta) ** alpha * np.log(1 - zeta, out=np.zeros_like(zeta), where=1 - zeta > 0)
515 tmp2 = (1 + zeta) ** alpha * np.log(1 + zeta)
516 duv = (tmp1 + tmp2) * (2**alpha - 2)
517 udv = ((1 - zeta) ** alpha + (1 + zeta) ** alpha - 2) * 2**alpha * np.log(2)
518 vv = (2**alpha - 2) ** 2
519 return (duv - udv) * dalphadtheta / vv
522def _get_dphidzeta(rs, theta, zeta, phi_params):
523 """Calculate dphi / dzeta."""
524 alpha = _get_alpha(rs, theta, phi_params)
525 tmp1 = alpha * (1 + zeta) ** alpha / (1 + zeta)
526 tmp2 = np.divide(
527 alpha * (1 - zeta) ** alpha, 1 - zeta, out=np.zeros_like(zeta), where=1 - zeta > 0
528 )
529 return (tmp1 - tmp2) / (2**alpha - 2)
532# ### alpha and derivatives ###
535def _get_alpha(rs, theta, phi_params):
536 """Calculate alpha.
538 Reference: Phys. Rev. Lett. 112, 076403.
539 """
540 g = _get_g(rs, phi_params)
541 lamda = _get_lambda(rs, theta, phi_params)
542 return 2 - g * np.exp(-theta * lamda)
545def _get_dalphadrs(rs, theta, phi_params):
546 """Calculate dalpha / drs."""
547 g = _get_g(rs, phi_params)
548 lamda = _get_lambda(rs, theta, phi_params)
549 dgdrs = _get_dgdrs(rs, phi_params)
550 dlambdadrs = _get_dlambdadrs(rs, theta, phi_params)
551 return -dgdrs * np.exp(-theta * lamda) + dlambdadrs * theta * g * np.exp(-theta * lamda)
554def _get_dalphadtheta(rs, theta, phi_params):
555 """Calculate dalpha / dtheta."""
556 g = _get_g(rs, phi_params)
557 lamda = _get_lambda(rs, theta, phi_params)
558 dlambdadtheta = _get_dlambdadtheta(rs, phi_params)
559 return (dlambdadtheta * theta + lamda) * g * np.exp(-theta * lamda)
562# ### g and derivative ###
565def _get_g(rs, phi_params):
566 """Calculate g.
568 Reference: Phys. Rev. Lett. 112, 076403.
569 """
570 return (phi_params.g1 + phi_params.g2 * rs) / (1 + phi_params.g3 * rs)
573def _get_dgdrs(rs, phi_params):
574 """Calculate dg / drs."""
575 return (
576 phi_params.g2 / (phi_params.g3 * rs + 1)
577 - phi_params.g3 * (phi_params.g2 * rs + phi_params.g1) / (phi_params.g3 * rs + 1) ** 2
578 )
581# ### lambda and derivatives ###
584def _get_lambda(rs, theta, phi_params):
585 """Calculate lambda.
587 Reference: Phys. Rev. Lett. 112, 076403.
588 """
589 return phi_params.lambda1 + phi_params.lambda2 * theta * rs ** (1 / 2)
592def _get_dlambdadrs(rs, theta, phi_params):
593 """Calculate dlambda / drs."""
594 return phi_params.lambda2 * theta / (2 * np.sqrt(rs))
597def _get_dlambdadtheta(rs, phi_params):
598 """Calculate dlambda / dtheta."""
599 return phi_params.lambda2 * np.sqrt(rs)