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 

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

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

41 divt = 2 * ks * n 

42 t = norm_dn / divt 

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

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

45 t2 = t**2 

46 At2 = A * t2 

47 A2t4 = At2**2 

48 divsum = 1 + At2 + A2t4 

49 div = (1 + At2) / divsum 

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

51 gec = gamma * np.log(nolog) 

52 

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

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

55 gvc = gec + dgec 

56 

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

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

59 

60 

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

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

63 

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

65 

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

67 

68 Args: 

69 n: Real-space electronic density. 

70 zeta: Relative spin polarization. 

71 

72 Keyword Args: 

73 beta: Functional parameter. 

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

75 **kwargs: Throwaway arguments. 

76 

77 Returns: 

78 PBE correlation energy density, potential, and vsigma. 

79 """ 

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

81 

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

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

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

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

86 vcup, vcdw = vc 

87 

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

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

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

91 phi2 = phi**2 

92 phi3 = phi2 * phi 

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

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

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

96 t2 = t**2 

97 At2 = A * t2 

98 A2t4 = At2**2 

99 divsum = 1 + At2 + A2t4 

100 div = (1 + At2) / divsum 

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

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

103 

104 # Handle divisions by zero 

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

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

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

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

109 bfpre = expec / phi3 

110 bfup = bfpre * (vcup - ec) 

111 bfdw = bfpre * (vcdw - ec) 

112 dgecpre = beta * t2 * phi3 / nolog 

113 dgecup = dgecpre * (-7 / 3 * div - factor * (A * bfup / beta - 7 / 3)) 

114 dgecdw = dgecpre * (-7 / 3 * div - factor * (A * bfdw / beta - 7 / 3)) 

115 dgeczpre = ( 

116 3 * gec / phi 

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

118 ) * dfz 

119 dgeczup = dgeczpre * (1 - zeta) 

120 dgeczdw = -dgeczpre * (1 + zeta) 

121 gvcup = gec + dgecup + dgeczup 

122 gvcdw = gec + dgecdw + dgeczdw 

123 

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

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

126 return ec + gec, np.array([vcup + gvcup, vcdw + gvcdw]), vsigmac