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

1# SPDX-FileCopyrightText: 2021 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 numpy as np 

9from scipy.linalg import norm 

10 

11from .lda_x import lda_x, lda_x_spin 

12 

13 

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

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

16 

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

18 

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

20 

21 Args: 

22 n: Real-space electronic density. 

23 

24 Keyword Args: 

25 mu: Functional parameter. 

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

27 **kwargs: Throwaway arguments. 

28 

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 return ex + gex / n, np.array([vx + gvx]), np.array([0.5 * vsigmax]) 

35 

36 

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

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

39 

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

41 

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

43 

44 Args: 

45 n: Real-space electronic density. 

46 zeta: Relative spin polarization. 

47 

48 Keyword Args: 

49 mu: Functional parameter. 

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

51 **kwargs: Throwaway arguments. 

52 

53 Returns: 

54 PBE exchange energy density, potential, and vsigma. 

55 """ 

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

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

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

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

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

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

62 

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

64 

65 vsigmax = np.array([vsigma_up, np.zeros_like(ex), vsigma_dw]) 

66 return ex + 0.5 * (ex_up + ex_dw) / n, np.array([vx[0] + vx_up, vx[1] + vx_dw]), vsigmax 

67 

68 

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

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

71 

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

73 

74 Args: 

75 n: Real-space electronic density. 

76 

77 Keyword Args: 

78 mu: Functional parameter. 

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

80 **kwargs: Throwaway arguments. 

81 

82 Returns: 

83 PBE exchange energy density, potential, and vsigma. 

84 """ 

85 kappa = 0.804 

86 

87 norm_dn = norm(dn, axis=1) 

88 kf = (3 * np.pi**2 * n) ** (1 / 3) 

89 # Handle divisions by zero 

90 # divkf = 1 / kf 

91 divkf = np.divide(1, kf, out=np.zeros_like(kf), where=(kf > 0)) 

92 # Handle divisions by zero 

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

94 s = np.divide(norm_dn * divkf, 2 * n, out=np.zeros_like(n), where=(n > 0)) 

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

96 Fx = kappa - kappa / f1 

97 exunif = -3 * kf / (4 * np.pi) 

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

99 sx = exunif * Fx 

100 

101 dsdn = -4 / 3 * s 

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

103 dexunif = exunif / 3 

104 exunifdFx = exunif * dFxds 

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

106 

107 # Handle divisions by zero 

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

109 vsigmax = np.divide( 

110 exunifdFx * divkf, 2 * norm_dn, out=np.zeros_like(norm_dn), where=(norm_dn > 0) 

111 ) 

112 return sx * n, np.array([vx]), vsigmax