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

1# SPDX-FileCopyrightText: 2023 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Chachiyo GGA correlation. 

4 

5Reference: Comput. Theor. Chem. 1172, 112669. 

6""" 

7 

8import math 

9 

10from eminus import backend as xp 

11 

12from .lda_c_chachiyo_mod import chachiyo_scaling_mod as weight_function 

13 

14 

15def gga_c_chachiyo(n, dn_spin=None, **kwargs): 

16 """Chachiyo parametrization of the correlation functional (spin-paired). 

17 

18 Corresponds to the functional with the label GGA_C_CHACHIYO and ID 309 in Libxc. 

19 

20 Reference: Comput. Theor. Chem. 1172, 112669. 

21 

22 Args: 

23 n: Real-space electronic density. 

24 

25 Keyword Args: 

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

27 **kwargs: Throwaway arguments. 

28 

29 Returns: 

30 Chachiyo correlation energy density, potential and vsigma. 

31 """ 

32 h = 0.06672632 # 0.5 * 0.00847 * 16 * (3 / math.pi)**(1 / 3) 

33 

34 # ### Start lda_c_chachiyo_mod ### # 

35 a = -0.01554535 # (xp.log(2) - 1) / (2 * math.pi**2) 

36 b = 20.4562557 

37 

38 rs = (3 / (4 * math.pi * n)) ** (1 / 3) 

39 rs2 = rs**2 

40 ecinner = 1 + b / rs + b / rs2 

41 

42 ec = a * xp.log(ecinner) 

43 # ### End lda_c_chachiyo_mod ### # 

44 

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

46 t = (math.pi / 3) ** (1 / 6) / 4 * norm_dn / n ** (7 / 6) 

47 t2 = t**2 

48 gec = 1 + t2 

49 expgec = gec ** (h / ec) 

50 

51 # ### Start lda_c_chachiyo_mod ### # 

52 vc = ec + a * b * (2 + rs) / (3 * (b + b * rs + rs2)) 

53 # ### End lda_c_chachiyo_mod ### # 

54 

55 term1 = h * (1 - 1 / gec) 

56 term2 = h * xp.log(gec) * (1 - vc / ec) 

57 gvc = (vc - 7 / 3 * term1 + term2) * expgec 

58 

59 vsigmac = n * expgec * term1 / norm_dn**2 

60 return ec * expgec, xp.stack([gvc]), xp.stack([vsigmac]) 

61 

62 

63def gga_c_chachiyo_spin(n, zeta, dn_spin=None, **kwargs): 

64 """Chachiyo parametrization of the correlation functional (spin-polarized). 

65 

66 Corresponds to the functional with the label GGA_C_CHACHIYO and ID 309 in Libxc. 

67 

68 Reference: Comput. Theor. Chem. 1172, 112669. 

69 

70 Args: 

71 n: Real-space electronic density. 

72 zeta: Relative spin polarization. 

73 

74 Keyword Args: 

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

76 **kwargs: Throwaway arguments. 

77 

78 Returns: 

79 Chachiyo correlation energy density, potential and vsigma. 

80 """ 

81 h = 0.06672632 # 0.5 * 0.00847 * 16 * (3 / math.pi)**(1 / 3) 

82 

83 # ### Start lda_c_chachiyo_spin_mod ### # 

84 a0 = -0.01554535 # (xp.log(2) - 1) / (2 * math.pi**2) 

85 a1 = -0.007772675 # (xp.log(2) - 1) / (4 * math.pi**2) 

86 b0 = 20.4562557 

87 b1 = 27.4203609 

88 

89 rs = (3 / (4 * math.pi * n)) ** (1 / 3) 

90 rs2 = rs**2 

91 

92 fzeta, dfdzeta = weight_function(zeta) 

93 

94 ec0inner = 1 + b0 / rs + b0 / rs2 

95 ec1inner = 1 + b1 / rs + b1 / rs2 

96 ec0 = a0 * xp.log(ec0inner) 

97 ec1 = a1 * xp.log(ec1inner) 

98 

99 ec = ec0 + (ec1 - ec0) * fzeta 

100 # ### End lda_c_chachiyo_spin_mod ### # 

101 

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

103 t = (math.pi / 3) ** (1 / 6) / 4 * norm_dn / n ** (7 / 6) 

104 t2 = t**2 

105 gec = 1 + t2 

106 expgec = gec ** (h / ec) 

107 

108 # ### Start lda_c_chachiyo_spin_mod ### # 

109 factor = -1 / rs2 - 2 / rs**3 

110 dec0drs = a0 / ec0inner * b0 * factor 

111 dec1drs = a1 / ec1inner * b1 * factor 

112 decdrs = dec0drs + (dec1drs - dec0drs) * fzeta 

113 # prefactor = ec - rs / 3 * decdrs 

114 decdf = (ec1 - ec0) * dfdzeta 

115 

116 # vc_up = prefactor + decdf * (1 - zeta) 

117 # vc_dw = prefactor - decdf * (1 + zeta) 

118 # ### End lda_c_chachiyo_spin_mod ### # 

119 

120 dn2 = ( 

121 xp.linalg.norm(dn_spin[0], axis=1) ** 2 

122 + 2 * xp.sum(dn_spin[0] * dn_spin[1], axis=1) 

123 + xp.linalg.norm(dn_spin[1], axis=1) ** 2 

124 ) 

125 ht2divgecdn2 = (1 - 1 / gec) * h / norm_dn**2 

126 term1 = -7 / 3 * ht2divgecdn2 * dn2 

127 term2 = 1 - h * xp.log(gec) / ec 

128 prefactor = -decdrs * rs / 3 

129 gvc = ec + term1 + term2 * prefactor 

130 gvc_up = gvc + term2 * decdf * (1 - zeta) 

131 gvc_dw = gvc - term2 * decdf * (1 + zeta) 

132 

133 vsigma = n * expgec * 2 * ht2divgecdn2 

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

135 return ec * expgec, xp.stack([gvc_up, gvc_dw]) * expgec, vsigmac