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

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

9 

10 

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

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

13 

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

15 

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

17 

18 Args: 

19 n: Real-space electronic density. 

20 

21 Keyword Args: 

22 A: Functional parameter. 

23 b: Functional parameter. 

24 c: Functional parameter. 

25 x0: Functional parameter. 

26 **kwargs: Throwaway arguments. 

27 

28 Returns: 

29 VWN correlation energy density and potential. 

30 """ 

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

32 

33 x = np.sqrt(rs) 

34 X = rs + b * x + c 

35 Q = np.sqrt(4 * c - b**2) 

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

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

38 tx = 2 * x + b 

39 tanx = np.arctan(Q / tx) 

40 

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

42 

43 tt = tx**2 + Q**2 

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

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

46 ) 

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

48 

49 

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

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

52 

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

54 

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

56 

57 Args: 

58 n: Real-space electronic density. 

59 zeta: Relative spin polarization. 

60 

61 Keyword Args: 

62 **kwargs: Throwaway arguments. 

63 

64 Returns: 

65 VWN correlation energy density and potential. 

66 """ 

67 A = (0.0310907, 0.01554535, -1 / (6 * np.pi**2)) 

68 b = (3.72744, 7.06042, 1.13107) 

69 c = (12.9352, 18.0578, 13.0045) 

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

71 

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

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

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

75 

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

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

78 zeta4 = zeta**4 

79 fzetaz4 = fzeta * zeta4 

80 De = ec1 - ec0 - ac / d2fzeta0 

81 

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

83 

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

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

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

87 

88 vcup = dec1 + (1 - zeta) * dec2 

89 vcdw = dec1 - (1 + zeta) * dec2 

90 return ec, np.array([vcup, vcdw]), None