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

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

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

4 

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

6""" 

7 

8import numpy as np 

9from scipy.linalg import norm 

10 

11from .lda_c_pw_mod import lda_c_pw_mod, lda_c_pw_mod_spin 

12 

13 

14def gga_c_pbe(n, beta=0.06672455060314922, dn_spin=None, **kwargs): 

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

16 

17 Corresponds to the functional with the label GGA_C_PBE and ID 130 in Libxc. 

18 

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

20 

21 Args: 

22 n: Real-space electronic density. 

23 

24 Keyword Args: 

25 beta: Functional parameter. 

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

27 **kwargs: Throwaway arguments. 

28 

29 Returns: 

30 PBE correlation energy density, potential, and vsigma. 

31 """ 

32 gamma = (1 - np.log(2)) / np.pi**2 

33 

34 pi34 = (3 / (4 * np.pi)) ** (1 / 3) 

35 rs = pi34 * n ** (-1 / 3) 

36 norm_dn = norm(dn_spin[0], axis=1) 

37 ec, vc, _ = lda_c_pw_mod(n, **kwargs) 

38 vc = vc[0] # Remove spin dimension for the correct shape 

39 

40 kf = (9 / 4 * np.pi) ** (1 / 3) / rs 

41 ks = np.sqrt(4 * kf / np.pi) 

42 divt = 2 * ks * n 

43 t = norm_dn / divt 

44 expec = np.exp(-ec / gamma) 

45 A = beta / (gamma * (expec - 1)) 

46 t2 = t**2 

47 At2 = A * t2 

48 A2t4 = At2**2 

49 divsum = 1 + At2 + A2t4 

50 div = (1 + At2) / divsum 

51 nolog = 1 + beta / gamma * t2 * div 

52 gec = gamma * np.log(nolog) 

53 

54 factor = A2t4 * (2 + At2) / divsum**2 

55 dgec = beta * t2 / nolog * (-7 / 3 * div - factor * (A * expec * (vc - ec) / beta - 7 / 3)) 

56 gvc = gec + dgec 

57 

58 vsigmac = beta / (divt * ks) * (div - factor) / nolog 

59 return ec + gec, np.array([vc + gvc]), np.array([0.5 * vsigmac]) 

60 

61 

62def gga_c_pbe_spin(n, zeta, beta=0.06672455060314922, dn_spin=None, **kwargs): 

63 """Perdew-Burke-Ernzerhof parametrization of the correlation functional (spin-polarized). 

64 

65 Corresponds to the functional with the label GGA_C_PBE and ID 130 in Libxc. 

66 

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

68 

69 Args: 

70 n: Real-space electronic density. 

71 zeta: Relative spin polarization. 

72 

73 Keyword Args: 

74 beta: Functional parameter. 

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

76 **kwargs: Throwaway arguments. 

77 

78 Returns: 

79 PBE correlation energy density, potential, and vsigma. 

80 """ 

81 gamma = (1 - np.log(2)) / np.pi**2 

82 

83 pi34 = (3 / (4 * np.pi)) ** (1 / 3) 

84 rs = pi34 * n ** (-1 / 3) 

85 norm_dn = norm(dn_spin[0] + dn_spin[1], axis=1) 

86 ec, vc, _ = lda_c_pw_mod_spin(n, zeta, **kwargs) 

87 vc_up, vc_dw = vc 

88 

89 kf = (9 / 4 * np.pi) ** (1 / 3) / rs 

90 ks = np.sqrt(4 * kf / np.pi) 

91 phi = ((1 + zeta) ** (2 / 3) + (1 - zeta) ** (2 / 3)) / 2 

92 phi2 = phi**2 

93 phi3 = phi2 * phi 

94 t = norm_dn / (2 * phi * ks * n) 

95 expec = np.exp(-ec / (gamma * phi3)) 

96 A = beta / (gamma * (expec - 1)) 

97 t2 = t**2 

98 At2 = A * t2 

99 A2t4 = At2**2 

100 divsum = 1 + At2 + A2t4 

101 div = (1 + At2) / divsum 

102 nolog = 1 + beta / gamma * t2 * div 

103 gec = gamma * phi3 * np.log(nolog) 

104 

105 # Handle divisions by zero 

106 with np.errstate(divide="ignore", invalid="ignore"): 

107 dfz = ((1 + zeta) ** (-1 / 3) - (1 - zeta) ** (-1 / 3)) / 3 

108 dfz = np.nan_to_num(dfz, nan=0, posinf=0, neginf=0) 

109 factor = A2t4 * (2 + At2) / divsum**2 

110 bfpre = expec / phi3 

111 bf_up = bfpre * (vc_up - ec) 

112 bf_dw = bfpre * (vc_dw - ec) 

113 dgecpre = beta * t2 * phi3 / nolog 

114 dgec_up = dgecpre * (-7 / 3 * div - factor * (A * bf_up / beta - 7 / 3)) 

115 dgec_dw = dgecpre * (-7 / 3 * div - factor * (A * bf_dw / beta - 7 / 3)) 

116 dgeczpre = ( 

117 3 * gec / phi 

118 - beta * t2 * phi2 / nolog * (2 * div - factor * (3 * A * expec * ec / phi3 / beta + 2)) 

119 ) * dfz 

120 dgecz_up = dgeczpre * (1 - zeta) 

121 dgecz_dw = -dgeczpre * (1 + zeta) 

122 gvc_up = gec + dgec_up + dgecz_up 

123 gvc_dw = gec + dgec_dw + dgecz_dw 

124 

125 vsigma = beta * phi / (2 * ks * ks * n) * (div - factor) / nolog 

126 vsigmac = np.array([0.5 * vsigma, vsigma, 0.5 * vsigma]) 

127 return ec + gec, np.array([vc_up + gvc_up, vc_dw + gvc_dw]), vsigmac