Coverage for eminus/xc/gga_x_pbe.py: 100.00%
35 statements
« prev ^ index » next coverage.py v7.6.10, created at 2025-01-16 10:16 +0000
1# SPDX-FileCopyrightText: 2021 The eminus developers
2# SPDX-License-Identifier: Apache-2.0
3"""Perdew-Burke-Ernzerhof GGA exchange.
5Reference: Phys. Rev. Lett. 77, 3865.
6"""
8import numpy as np
9from scipy.linalg import norm
11from .lda_x import lda_x, lda_x_spin
14def gga_x_pbe(n, mu=0.2195149727645171, dn_spin=None, **kwargs):
15 """Perdew-Burke-Ernzerhof parametrization of the exchange functional (spin-paired).
17 Corresponds to the functional with the label GGA_X_PBE and ID 101 in Libxc.
19 Reference: Phys. Rev. Lett. 77, 3865.
21 Args:
22 n: Real-space electronic density.
24 Keyword Args:
25 mu: Functional parameter.
26 dn_spin: Real-space gradient of densities per spin channel.
27 **kwargs: Throwaway arguments.
29 Returns:
30 PBE exchange energy density, potential, and vsigma.
31 """
32 ex, vx, _ = lda_x(n, **kwargs)
33 gex, gvx, vsigmax = pbe_x_base(n, mu, dn_spin[0], **kwargs)
34 vx, gvx = vx[0], gvx[0] # Remove spin dimension for the correct shape
35 return ex + gex / n, np.array([vx + gvx]), np.array([0.5 * vsigmax])
38def gga_x_pbe_spin(n, zeta, mu=0.2195149727645171, dn_spin=None, **kwargs):
39 """Perdew-Burke-Ernzerhof parametrization of the exchange functional (spin-polarized).
41 Corresponds to the functional with the label GGA_X_PBE and ID 101 in Libxc.
43 Reference: Phys. Rev. Lett. 77, 3865.
45 Args:
46 n: Real-space electronic density.
47 zeta: Relative spin polarization.
49 Keyword Args:
50 mu: Functional parameter.
51 dn_spin: Real-space gradient of densities per spin channel.
52 **kwargs: Throwaway arguments.
54 Returns:
55 PBE exchange energy density, potential, and vsigma.
56 """
57 # Use the spin-scaling relationship Exc(n_up, n_down)=(Exc(2 n_up)+Exc(2 n_down))/2
58 zeta = zeta[0] # Getting the non-zero values from zeta adds an extra dimension, remove it here
59 n_up = zeta * n + n # 2 * n_up
60 n_dw = -zeta * n + n # 2 * n_down
61 ex_up, vx_up, vsigma_up = pbe_x_base(n_up, mu, 2 * dn_spin[0], **kwargs)
62 ex_dw, vx_dw, vsigma_dw = pbe_x_base(n_dw, mu, 2 * dn_spin[1], **kwargs)
63 vx_up, vx_dw = vx_up[0], vx_dw[0] # Remove spin dimension for the correct shape
65 ex, vx, _ = lda_x_spin(n, zeta, **kwargs)
67 vsigmax = np.array([vsigma_up, np.zeros_like(ex), vsigma_dw])
68 return ex + 0.5 * (ex_up + ex_dw) / n, np.array([vx[0] + vx_up, vx[1] + vx_dw]), vsigmax
71def pbe_x_base(n, mu=0.2195149727645171, dn=None, **kwargs):
72 """Base PBE exchange functional to be used in the spin-paired and -polarized case.
74 Reference: Phys. Rev. Lett. 77, 3865.
76 Args:
77 n: Real-space electronic density.
79 Keyword Args:
80 mu: Functional parameter.
81 dn: Real-space gradient of densities per spin channel.
82 **kwargs: Throwaway arguments.
84 Returns:
85 PBE exchange energy density, potential, and vsigma.
86 """
87 kappa = 0.804
89 norm_dn = norm(dn, axis=1)
90 kf = (3 * np.pi**2 * n) ** (1 / 3)
91 # Handle divisions by zero
92 # divkf = 1 / kf
93 divkf = np.divide(1, kf, out=np.zeros_like(kf), where=kf > 0)
94 # Handle divisions by zero
95 # s = norm_dn * divkf / (2 * n)
96 s = np.divide(norm_dn * divkf, 2 * n, out=np.zeros_like(n), where=n > 0)
97 f1 = 1 + mu * s**2 / kappa
98 Fx = kappa - kappa / f1
99 exunif = -3 * kf / (4 * np.pi)
100 # In Fx a "1 + " is missing, since n * exunif is the Slater exchange that is added later
101 sx = exunif * Fx
103 dsdn = -4 / 3 * s
104 dFxds = 2 * mu * s / f1**2
105 dexunif = exunif / 3
106 exunifdFx = exunif * dFxds
107 vx = sx + dexunif * Fx + exunifdFx * dsdn # dFx/dn = dFx/ds * ds/dn
109 # Handle divisions by zero
110 # vsigmax = exunifdFx * divkf / (2 * norm_dn)
111 vsigmax = np.divide(
112 exunifdFx * divkf, 2 * norm_dn, out=np.zeros_like(norm_dn), where=norm_dn > 0
113 )
114 return sx * n, np.array([vx]), vsigmax