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

1# SPDX-FileCopyrightText: 2023 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 math 

9 

10import numpy as np 

11 

12from eminus import backend as xp 

13 

14from .lda_c_pw_mod import lda_c_pw_mod, lda_c_pw_mod_spin 

15 

16 

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

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

19 

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

21 

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

23 

24 Args: 

25 n: Real-space electronic density. 

26 

27 Keyword Args: 

28 beta: Functional parameter. 

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

30 **kwargs: Throwaway arguments. 

31 

32 Returns: 

33 PBE correlation energy density, potential, and vsigma. 

34 """ 

35 gamma = (1 - math.log(2)) / math.pi**2 

36 

37 pi34 = (3 / (4 * math.pi)) ** (1 / 3) 

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

39 norm_dn = xp.linalg.norm(dn_spin[0], axis=1) 

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

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

42 

43 kf = (9 / 4 * math.pi) ** (1 / 3) / rs 

44 ks = xp.sqrt(4 * kf / math.pi) 

45 divt = 2 * ks * n 

46 t = norm_dn / divt 

47 expec = xp.exp(-ec / gamma) 

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

49 t2 = t**2 

50 At2 = A * t2 

51 A2t4 = At2**2 

52 divsum = 1 + At2 + A2t4 

53 div = (1 + At2) / divsum 

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

55 gec = gamma * xp.log(nolog) 

56 

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

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

59 gvc = gec + dgec 

60 

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

62 return ec + gec, xp.stack([vc + gvc]), xp.stack([0.5 * vsigmac]) 

63 

64 

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

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

67 

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

69 

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

71 

72 Args: 

73 n: Real-space electronic density. 

74 zeta: Relative spin polarization. 

75 

76 Keyword Args: 

77 beta: Functional parameter. 

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

79 **kwargs: Throwaway arguments. 

80 

81 Returns: 

82 PBE correlation energy density, potential, and vsigma. 

83 """ 

84 gamma = (1 - math.log(2)) / math.pi**2 

85 

86 pi34 = (3 / (4 * math.pi)) ** (1 / 3) 

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

88 norm_dn = xp.linalg.norm(dn_spin[0] + dn_spin[1], axis=1) 

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

90 vc_up, vc_dw = vc 

91 

92 kf = (9 / 4 * math.pi) ** (1 / 3) / rs 

93 ks = xp.sqrt(4 * kf / math.pi) 

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

95 phi2 = phi**2 

96 phi3 = phi2 * phi 

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

98 expec = xp.exp(-ec / (gamma * phi3)) 

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

100 t2 = t**2 

101 At2 = A * t2 

102 A2t4 = At2**2 

103 divsum = 1 + At2 + A2t4 

104 div = (1 + At2) / divsum 

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

106 gec = gamma * phi3 * xp.log(nolog) 

107 

108 # Handle divisions by zero 

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

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

111 dfz = xp.nan_to_num(dfz, nan=0, posinf=0, neginf=0) 

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

113 bfpre = expec / phi3 

114 bf_up = bfpre * (vc_up - ec) 

115 bf_dw = bfpre * (vc_dw - ec) 

116 dgecpre = beta * t2 * phi3 / nolog 

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

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

119 dgeczpre = ( 

120 3 * gec / phi 

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

122 ) * dfz 

123 dgecz_up = dgeczpre * (1 - zeta) 

124 dgecz_dw = -dgeczpre * (1 + zeta) 

125 gvc_up = gec + dgec_up + dgecz_up 

126 gvc_dw = gec + dgec_dw + dgecz_dw 

127 

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

129 vsigmac = xp.stack([0.5 * vsigma, vsigma, 0.5 * vsigma]) 

130 return ec + gec, xp.stack([vc_up + gvc_up, vc_dw + gvc_dw]), vsigmac