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

1# SPDX-FileCopyrightText: 2023 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Vosko-Wilk-Nusair LDA correlation. 

4 

5Reference: Phys. Rev. B 22, 3812. 

6""" 

7 

8import math 

9 

10from eminus import backend as xp 

11 

12 

13def lda_c_vwn(n, A=0.0310907, b=3.72744, c=12.9352, x0=-0.10498, **kwargs): 

14 """Vosko-Wilk-Nusair parametrization of the correlation functional (spin-paired). 

15 

16 Corresponds to the functional with the label LDA_C_VWN and ID 7 in Libxc. 

17 

18 Reference: Phys. Rev. B 22, 3812. 

19 

20 Args: 

21 n: Real-space electronic density. 

22 

23 Keyword Args: 

24 A: Functional parameter. 

25 b: Functional parameter. 

26 c: Functional parameter. 

27 x0: Functional parameter. 

28 **kwargs: Throwaway arguments. 

29 

30 Returns: 

31 VWN correlation energy density and potential. 

32 """ 

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

34 

35 x = xp.sqrt(rs) 

36 X = rs + b * x + c 

37 Q = math.sqrt(4 * c - b**2) 

38 fx0 = b * x0 / (x0**2 + b * x0 + c) 

39 f3 = 2 * (2 * x0 + b) / Q 

40 tx = 2 * x + b 

41 tanx = xp.arctan(Q / tx) 

42 

43 ec = A * (xp.log(rs / X) + 2 * b / Q * tanx - fx0 * (xp.log((x - x0) ** 2 / X) + f3 * tanx)) 

44 

45 tt = tx**2 + Q**2 

46 vc = ec - x * A / 6 * ( 

47 2 / x - tx / X - 4 * b / tt - fx0 * (2 / (x - x0) - tx / X - 4 * (2 * x0 + b) / tt) 

48 ) 

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

50 

51 

52def lda_c_vwn_spin(n, zeta, **kwargs): 

53 """Vosko-Wilk-Nusair parametrization of the correlation functional (spin-polarized). 

54 

55 Corresponds to the functional with the label LDA_C_VWN and ID 7 in Libxc. 

56 

57 Reference: Phys. Rev. B 22, 3812. 

58 

59 Args: 

60 n: Real-space electronic density. 

61 zeta: Relative spin polarization. 

62 

63 Keyword Args: 

64 **kwargs: Throwaway arguments. 

65 

66 Returns: 

67 VWN correlation energy density and potential. 

68 """ 

69 A = (0.0310907, 0.01554535, -1 / (6 * math.pi**2)) 

70 b = (3.72744, 7.06042, 1.13107) 

71 c = (12.9352, 18.0578, 13.0045) 

72 x0 = (-0.10498, -0.325, -0.0047584) 

73 

74 ec0, vc0, _ = lda_c_vwn(n, A[0], b[0], c[0], x0[0]) # Paramagnetic 

75 ec1, vc1, _ = lda_c_vwn(n, A[1], b[1], c[1], x0[1]) # Ferromagnetic 

76 ac, dac, _ = lda_c_vwn(n, A[2], b[2], c[2], x0[2]) # Spin stiffness 

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

78 

79 d2fzeta0 = 4 / 9 / (2 ** (1 / 3) - 1) 

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

81 zeta4 = zeta**4 

82 fzetaz4 = fzeta * zeta4 

83 De = ec1 - ec0 - ac / d2fzeta0 

84 

85 ec = ec0 + ac / d2fzeta0 * fzeta + De * fzetaz4 

86 

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

88 dec1 = vc0 + dac / d2fzeta0 * fzeta + (vc1 - vc0 - dac / d2fzeta0) * fzetaz4 

89 dec2 = ac / d2fzeta0 * dfzeta + De * (4 * zeta**3 * fzeta + zeta4 * dfzeta) 

90 

91 vc_up = dec1 + (1 - zeta) * dec2 

92 vc_dw = dec1 - (1 + zeta) * dec2 

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