Coverage for gpaw/xc/ri/spherical_hse_kernel.py: 100%

1""" 

2 

3Implementation of spherical harmonic expansion of screened Coulomb kernel. 

4Based on 

5 

6 János G Ángyán et al 2006 J. Phys. A: Math. Gen. 39 8613 

7 

8 

9""" 

10 

11 

12import numpy as np 

13from scipy.special import erfc, comb, factorial2 

14from math import factorial 

15 

16 

17def safeerfc(x): 

18 taylor = (1 

19 - 2 * x / np.pi**0.5 

20 + 2 * x**3 / (3 * np.pi**0.5) 

21 - x**5 / (5 * np.pi)**0.5) 

22 return np.where(x < 1e-4, taylor, erfc(x)) 

23 

24 

25def Dnk(n, k, Xi): 

26 # Eq. 28 

27 if k == 0: 

28 sum = 0 

29 for m in range(1, n + 1): 

30 sum += 2**(-m) * Xi**(-2 * m) / factorial2(2 * n - 2 * m + 1) 

31 return (safeerfc(Xi) 

32 + np.exp(-Xi**2) 

33 / (np.pi**0.5) * 2**(n + 1) * Xi**(2 * n + 1) * sum) 

34 # Eq. 29 

35 sum = 0 

36 for m in range(1, k + 1): 

37 sum += (comb(m - k - 1, m - 1) * 2**(k - m) * Xi**(2 * (k - m)) 

38 / factorial2(2 * n + 2 * k - 2 * m + 1)) 

39 

40 return (np.exp(-Xi**2) 

41 * 2**(n + 1) * (2 * n + 1) * Xi**(2 * n + 1) 

42 / np.pi**0.5 / factorial(k) 

43 / (2 * n + 2 * k + 1) * sum) 

44 

45 

46def Phinj(n, j, Xi, xi): 

47 # Eq. 30 

48 sum = 0 

49 for k in range(j): 

50 sum += Dnk(n, k, Xi) / (Xi**(n + 1)) * xi**(n + 2 * k) 

51 return sum 

52 

53 

54def Hn(n, Xi, xi): 

55 """ 

56 

57 Helper function (Eq. 24) 

58 

59 """ 

60 return 1 / (2 * (xi * Xi)**(n + 1)) * ((Xi**(2 * n + 1) + xi**(2 * n + 1)) * safeerfc(Xi + xi) - (Xi**(2 * n + 1) - xi**(2 * n + 1)) * safeerfc(Xi - xi)) # noqa: E501 

61 

62 

63def Fn(n, Xi, xi): 

64 """ 

65 

66 Helper function (Eq. 22). 

67 

68 It appears, that the article has a typo, because the summation 

69 starts at p=1, but correct results require to start at p=0. 

70 

71 """ 

72 prefactor = 2 / np.pi**0.5 

73 result = 0.0 

74 

75 for p in range(0, n + 1): 

76 result += (-1 / (4 * Xi * xi))**(p + 1) * factorial(n + p) / (factorial(p) * factorial(n - p)) * ((-1)**(n - p) * np.exp(-(xi + Xi)**2) - np.exp(-(xi - Xi)**2)) # noqa: E501 

77 taylor = np.exp(-Xi**2 - xi**2) * 2**(n + 1) * (3 + 2 * n + 2 * xi**2 * Xi**2) * xi**n * Xi**n / (np.pi**0.5 * factorial2(2 * n + 3)) # noqa: E501 

78 

79 return np.where((Xi * xi)**(2 * n + 1) < 1e-6, taylor, prefactor * result) 

80 

81 

82def Phi(n, mu, R, r): 

83 """ 

84 

85 The official spherical kernel expansion 

86 

87 """ 

88 Rg = np.maximum.reduce([R, r]) 

89 Rl = np.minimum.reduce([R, r]) 

90 

91 # Scaling as given by Eq. 16 and the text above. 

92 Xi = mu * Rg 

93 xi = mu * Rl 

94 

95 # Eq. 21 

96 result = Fn(n, Xi, xi) + Hn(n, Xi, xi) 

97 for m in range(1, n + 1): 

98 result += (Fn(n - m, Xi, xi) 

99 * (Xi**(2 * m) + xi**(2 * m)) / (xi * Xi)**m) 

100 

101 result = np.where(xi < [1e-3, 1e-2, 1e-1, 1e-1, 1e-1, 1e-1][n], 

102 Phinj(n, 2, Xi, xi), result) 

103 result *= mu 

104 

105 return result 

106 

107 

108""" 

109 

110Implementation of spherical harmonic expansion ends. GPAW spesific stuff 

111remains below. 

112 

113""" 

114 

115 

116class RadialHSE: 

117 def __init__(self, rgd, omega): 

118 self.rgd = rgd 

119 self.omega = omega 

120 

121 self.r1_gg = np.zeros((rgd.N, rgd.N)) 

122 self.r2_gg = np.zeros((rgd.N, rgd.N)) 

123 self.d_gg = np.zeros((rgd.N, rgd.N)) 

124 r_g = rgd.r_g.copy() 

125 r_g[0] = r_g[1] # XXX 

126 self.r1_gg[:] = r_g[None, :] 

127 self.r2_gg[:] = r_g[:, None] 

128 self.d_gg[:] = rgd.dr_g[None, :] * rgd.r_g[None, :]**2 * 4 * np.pi 

129 self.V_lgg = {} 

130 

131 def screened_coulomb(self, n_g, l): 

132 # Buffer different l-values for optimal performance 

133 if l not in self.V_lgg: 

134 kernel_gg = np.reshape(Phi(l, self.omega, self.r1_gg.ravel(), 

135 self.r2_gg.ravel()), 

136 self.d_gg.shape) / (2 * l + 1) 

137 self.V_lgg[l] = self.d_gg * kernel_gg 

138 vr_g = (self.V_lgg[l] @ n_g) * self.rgd.r_g 

139 return vr_g 

140 

141 def screened_coulomb_dv(self, n_g, l): 

142 return self.screened_coulomb(n_g, l) * self.rgd.r_g * self.rgd.dr_g