Coverage for eminus/xc/gga_x_pbe.py: 100.00%

1# SPDX-FileCopyrightText: 2023 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Perdew-Burke-Ernzerhof GGA exchange. 

4 

5Reference: Phys. Rev. Lett. 77, 3865. 

6""" 

7 

8import math 

9 

10import numpy as np 

11 

12from eminus import backend as xp 

13 

14from .lda_x import lda_x, lda_x_spin 

15 

16 

17def gga_x_pbe(n, mu=0.2195149727645171, dn_spin=None, **kwargs): 

18 """Perdew-Burke-Ernzerhof parametrization of the exchange functional (spin-paired). 

19 

20 Corresponds to the functional with the label GGA_X_PBE and ID 101 in Libxc. 

21 

22 Reference: Phys. Rev. Lett. 77, 3865. 

23 

24 Args: 

25 n: Real-space electronic density. 

26 

27 Keyword Args: 

28 mu: Functional parameter. 

29 dn_spin: Real-space gradient of densities per spin channel. 

30 **kwargs: Throwaway arguments. 

31 

32 Returns: 

33 PBE exchange energy density, potential, and vsigma. 

34 """ 

35 ex, vx, _ = lda_x(n, **kwargs) 

36 gex, gvx, vsigmax = pbe_x_base(n, mu, dn_spin[0], **kwargs) 

37 vx, gvx = vx[0], gvx[0] # Remove spin dimension for the correct shape 

38 return ex + gex / n, xp.stack([vx + gvx]), xp.stack([0.5 * vsigmax]) 

39 

40 

41def gga_x_pbe_spin(n, zeta, mu=0.2195149727645171, dn_spin=None, **kwargs): 

42 """Perdew-Burke-Ernzerhof parametrization of the exchange functional (spin-polarized). 

43 

44 Corresponds to the functional with the label GGA_X_PBE and ID 101 in Libxc. 

45 

46 Reference: Phys. Rev. Lett. 77, 3865. 

47 

48 Args: 

49 n: Real-space electronic density. 

50 zeta: Relative spin polarization. 

51 

52 Keyword Args: 

53 mu: Functional parameter. 

54 dn_spin: Real-space gradient of densities per spin channel. 

55 **kwargs: Throwaway arguments. 

56 

57 Returns: 

58 PBE exchange energy density, potential, and vsigma. 

59 """ 

60 # Use the spin-scaling relationship Exc(n_up, n_down)=(Exc(2 n_up)+Exc(2 n_down))/2 

61 zeta = zeta[0] # Getting the non-zero values from zeta adds an extra dimension, remove it here 

62 n_up = zeta * n + n # 2 * n_up 

63 n_dw = -zeta * n + n # 2 * n_down 

64 ex_up, vx_up, vsigma_up = pbe_x_base(n_up, mu, 2 * dn_spin[0], **kwargs) 

65 ex_dw, vx_dw, vsigma_dw = pbe_x_base(n_dw, mu, 2 * dn_spin[1], **kwargs) 

66 vx_up, vx_dw = vx_up[0], vx_dw[0] # Remove spin dimension for the correct shape 

67 

68 ex, vx, _ = lda_x_spin(n, zeta, **kwargs) 

69 

70 vsigmax = xp.stack([vsigma_up, xp.zeros_like(ex), vsigma_dw]) 

71 return ex + 0.5 * (ex_up + ex_dw) / n, xp.stack([vx[0] + vx_up, vx[1] + vx_dw]), vsigmax 

72 

73 

74def pbe_x_base(n, mu=0.2195149727645171, dn=None, **kwargs): 

75 """Base PBE exchange functional to be used in the spin-paired and -polarized case. 

76 

77 Reference: Phys. Rev. Lett. 77, 3865. 

78 

79 Args: 

80 n: Real-space electronic density. 

81 

82 Keyword Args: 

83 mu: Functional parameter. 

84 dn: Real-space gradient of densities per spin channel. 

85 **kwargs: Throwaway arguments. 

86 

87 Returns: 

88 PBE exchange energy density, potential, and vsigma. 

89 """ 

90 kappa = 0.804 

91 

92 norm_dn = xp.linalg.norm(dn, axis=1) 

93 kf = (3 * math.pi**2 * n) ** (1 / 3) 

94 # Handle divisions by zero 

95 # divkf = 1 / kf 

96 with np.errstate(divide="ignore"): 

97 divkf = xp.where(kf > 0, 1 / kf, 0) 

98 # Handle divisions by zero 

99 # s = norm_dn * divkf / (2 * n) 

100 with np.errstate(invalid="ignore"): 

101 s = xp.where(n > 0, norm_dn * divkf / (2 * n), 0) 

102 f1 = 1 + mu * s**2 / kappa 

103 Fx = kappa - kappa / f1 

104 exunif = -3 * kf / (4 * math.pi) 

105 # In Fx a "1 + " is missing, since n * exunif is the Slater exchange that is added later 

106 sx = exunif * Fx 

107 

108 dsdn = -4 / 3 * s 

109 dFxds = 2 * mu * s / f1**2 

110 dexunif = exunif / 3 

111 exunifdFx = exunif * dFxds 

112 vx = sx + dexunif * Fx + exunifdFx * dsdn # dFx/dn = dFx/ds * ds/dn 

113 

114 # Handle divisions by zero 

115 # vsigmax = exunifdFx * divkf / (2 * norm_dn) 

116 with np.errstate(invalid="ignore"): 

117 vsigmax = xp.where(norm_dn > 0, exunifdFx * divkf / (2 * norm_dn), 0) 

118 

119 return sx * n, xp.stack([vx]), vsigmax