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

1# SPDX-FileCopyrightText: 2023 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Perdew-Wang LDA correlation. 

4 

5Reference: Phys. Rev. B 45, 13244. 

6""" 

7 

8import math 

9 

10from eminus import backend as xp 

11 

12 

13def lda_c_pw(n, A=0.031091, a1=0.2137, b1=7.5957, b2=3.5876, b3=1.6382, b4=0.49294, **kwargs): 

14 """Perdew-Wang parametrization of the correlation functional (spin-paired). 

15 

16 Corresponds to the functional with the label LDA_C_PW and ID 12 in Libxc. 

17 

18 Reference: Phys. Rev. B 45, 13244. 

19 

20 Args: 

21 n: Real-space electronic density. 

22 

23 Keyword Args: 

24 A: Functional parameter. 

25 a1: Functional parameter. 

26 b1: Functional parameter. 

27 b2: Functional parameter. 

28 b3: Functional parameter. 

29 b4: Functional parameter. 

30 **kwargs: Throwaway arguments. 

31 

32 Returns: 

33 PW correlation energy density and potential. 

34 """ 

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

36 rs12 = xp.sqrt(rs) 

37 rs32 = rs * rs12 

38 rs2 = rs**2 

39 

40 om = 2 * A * (b1 * rs12 + b2 * rs + b3 * rs32 + b4 * rs2) 

41 olog = xp.log(1 + 1 / om) 

42 ec = -2 * A * (1 + a1 * rs) * olog 

43 

44 dom = 2 * A * (0.5 * b1 * rs12 + b2 * rs + 1.5 * b3 * rs32 + 2 * b4 * rs2) 

45 vc = -2 * A * (1 + 2 / 3 * a1 * rs) * olog - 2 / 3 * A * (1 + a1 * rs) * dom / (om * (om + 1)) 

46 return ec, xp.stack([vc]), None 

47 

48 

49def lda_c_pw_spin(n, zeta, A=(0.031091, 0.015545, 0.016887), fzeta0=1.709921, **kwargs): 

50 """Perdew-Wang parametrization of the correlation functional (spin-polarized). 

51 

52 Corresponds to the functional with the label LDA_C_PW and ID 12 in Libxc. 

53 

54 Reference: Phys. Rev. B 45, 13244. 

55 

56 Args: 

57 n: Real-space electronic density. 

58 zeta: Relative spin polarization. 

59 

60 Keyword Args: 

61 A: Functional parameters. 

62 fzeta0: Functional parameter. 

63 **kwargs: Throwaway arguments. 

64 

65 Returns: 

66 PW correlation energy density and potential. 

67 """ 

68 a1 = (0.2137, 0.20548, 0.11125) 

69 b1 = (7.5957, 14.1189, 10.357) 

70 b2 = (3.5876, 6.1977, 3.6231) 

71 b3 = (1.6382, 3.3662, 0.88026) 

72 b4 = (0.49294, 0.62517, 0.49671) 

73 

74 ec0, vc0, _ = lda_c_pw(n, A[0], a1[0], b1[0], b2[0], b3[0], b4[0]) # Unpolarized 

75 ec1, vc1, _ = lda_c_pw(n, A[1], a1[1], b1[1], b2[1], b3[1], b4[1]) # Polarized 

76 ac, dac, _ = lda_c_pw(n, A[2], a1[2], b1[2], b2[2], b3[2], b4[2]) # Spin stiffness 

77 ac = -ac # The PW spin interpolation is parametrized with -ac instead of ac 

78 vc0, vc1, dac = vc0[0], vc1[0], dac[0] # Remove spin dimension for the correct shape 

79 

80 fzeta = ((1 + zeta) ** (4 / 3) + (1 - zeta) ** (4 / 3) - 2) / (2 ** (4 / 3) - 2) 

81 zeta3 = zeta**3 

82 zeta4 = zeta3 * zeta 

83 

84 ec = ec0 + ac * fzeta * (1 - zeta4) / fzeta0 + (ec1 - ec0) * fzeta * zeta4 

85 

86 dac = -dac # Also change the sign for the derivative 

87 dfzeta = ((1 + zeta) ** (1 / 3) - (1 - zeta) ** (1 / 3)) * 4 / (3 * (2 ** (4 / 3) - 2)) 

88 factor1 = vc0 + dac * fzeta * (1 - zeta4) / fzeta0 + (vc1 - vc0) * fzeta * zeta4 

89 factor2 = ac / fzeta0 * (dfzeta * (1 - zeta4) - 4 * fzeta * zeta3) + (ec1 - ec0) * ( 

90 dfzeta * zeta4 + 4 * fzeta * zeta3 

91 ) 

92 

93 vc_up = factor1 + factor2 * (1 - zeta) 

94 vc_dw = factor1 - factor2 * (1 + zeta) 

95 return ec, xp.stack([vc_up, vc_dw]), None