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

1# SPDX-FileCopyrightText: 2021 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 numpy as np 

9 

10 

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

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

13 

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

15 

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

17 

18 Args: 

19 n: Real-space electronic density. 

20 

21 Keyword Args: 

22 A: Functional parameter. 

23 a1: Functional parameter. 

24 b1: Functional parameter. 

25 b2: Functional parameter. 

26 b3: Functional parameter. 

27 b4: Functional parameter. 

28 **kwargs: Throwaway arguments. 

29 

30 Returns: 

31 PW correlation energy density and potential. 

32 """ 

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

34 rs12 = np.sqrt(rs) 

35 rs32 = rs * rs12 

36 rs2 = rs**2 

37 

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

39 olog = np.log(1 + 1 / om) 

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

41 

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

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

44 return ec, np.array([vc]), None 

45 

46 

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

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

49 

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

51 

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

53 

54 Args: 

55 n: Real-space electronic density. 

56 zeta: Relative spin polarization. 

57 

58 Keyword Args: 

59 A: Functional parameters. 

60 fzeta0: Functional parameter. 

61 **kwargs: Throwaway arguments. 

62 

63 Returns: 

64 PW correlation energy density and potential. 

65 """ 

66 a1 = (0.2137, 0.20548, 0.11125) 

67 b1 = (7.5957, 14.1189, 10.357) 

68 b2 = (3.5876, 6.1977, 3.6231) 

69 b3 = (1.6382, 3.3662, 0.88026) 

70 b4 = (0.49294, 0.62517, 0.49671) 

71 

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

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

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

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

76 

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

78 zeta3 = zeta**3 

79 zeta4 = zeta3 * zeta 

80 

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

82 

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

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

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

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

87 dfzeta * zeta4 + 4 * fzeta * zeta3 

88 ) 

89 

90 vc0p = factor1 + factor2 * (1 - zeta) 

91 vcdw = factor1 - factor2 * (1 + zeta) 

92 return ec, np.array([vc0p, vcdw]), None