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

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Chachiyo GGA correlation. 

4 

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

6""" 

7 

8import numpy as np 

9from scipy.linalg import norm 

10 

11from .lda_c_chachiyo_mod import chachiyo_scaling_mod as weight_function 

12 

13 

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

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

16 

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

18 

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

20 

21 Args: 

22 n: Real-space electronic density. 

23 

24 Keyword Args: 

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

26 **kwargs: Throwaway arguments. 

27 

28 Returns: 

29 Chachiyo correlation energy density, potential and vsigma. 

30 """ 

31 h = 0.06672632 # 0.5 * 0.00847 * 16 * (3 / np.pi)**(1 / 3) 

32 

33 # ### Start lda_c_chachiyo_mod ### # 

34 a = -0.01554535 # (np.log(2) - 1) / (2 * np.pi**2) 

35 b = 20.4562557 

36 

37 rs = (3 / (4 * np.pi * n)) ** (1 / 3) 

38 rs2 = rs**2 

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

40 

41 ec = a * np.log(ecinner) 

42 # ### End lda_c_chachiyo_mod ### # 

43 

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

45 t = (np.pi / 3) ** (1 / 6) / 4 * norm_dn / n ** (7 / 6) 

46 t2 = t**2 

47 gec = 1 + t2 

48 expgec = gec ** (h / ec) 

49 

50 # ### Start lda_c_chachiyo_mod ### # 

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

52 # ### End lda_c_chachiyo_mod ### # 

53 

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

55 term2 = h * np.log(gec) * (1 - vc / ec) 

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

57 

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

59 return ec * expgec, np.array([gvc]), np.array([vsigmac]) 

60 

61 

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

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

64 

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

66 

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

68 

69 Args: 

70 n: Real-space electronic density. 

71 zeta: Relative spin polarization. 

72 

73 Keyword Args: 

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

75 **kwargs: Throwaway arguments. 

76 

77 Returns: 

78 Chachiyo correlation energy density, potential and vsigma. 

79 """ 

80 h = 0.06672632 # 0.5 * 0.00847 * 16 * (3 / np.pi)**(1 / 3) 

81 

82 # ### Start lda_c_chachiyo_spin_mod ### # 

83 a0 = -0.01554535 # (np.log(2) - 1) / (2 * np.pi**2) 

84 a1 = -0.007772675 # (np.log(2) - 1) / (4 * np.pi**2) 

85 b0 = 20.4562557 

86 b1 = 27.4203609 

87 

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

89 rs2 = rs**2 

90 

91 fzeta, dfdzeta = weight_function(zeta) 

92 

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

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

95 ec0 = a0 * np.log(ec0inner) 

96 ec1 = a1 * np.log(ec1inner) 

97 

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

99 # ### End lda_c_chachiyo_spin_mod ### # 

100 

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

102 t = (np.pi / 3) ** (1 / 6) / 4 * norm_dn / n ** (7 / 6) 

103 t2 = t**2 

104 gec = 1 + t2 

105 expgec = gec ** (h / ec) 

106 

107 # ### Start lda_c_chachiyo_spin_mod ### # 

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

109 dec0drs = a0 / ec0inner * b0 * factor 

110 dec1drs = a1 / ec1inner * b1 * factor 

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

112 # prefactor = ec - rs / 3 * decdrs 

113 decdf = (ec1 - ec0) * dfdzeta 

114 

115 # vcup = prefactor + decdf * (1 - zeta) 

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

117 # ### End lda_c_chachiyo_spin_mod ### # 

118 

119 dn2 = ( 

120 norm(dn_spin[0], axis=1) ** 2 

121 + 2 * np.sum(dn_spin[0] * dn_spin[1], axis=1) 

122 + norm(dn_spin[1], axis=1) ** 2 

123 ) 

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

125 term1 = -7 / 3 * ht2divgecdn2 * dn2 

126 term2 = 1 - h * np.log(gec) / ec 

127 prefactor = -decdrs * rs / 3 

128 gvc = ec + term1 + term2 * prefactor 

129 gvcup = gvc + term2 * decdf * (1 - zeta) 

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

131 

132 vsigma = n * expgec * 2 * ht2divgecdn2 

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

134 return ec * expgec, np.array([gvcup, gvcdw]) * expgec, vsigmac