Coverage for gpaw/atom/filter.py: 61%

1# Copyright (C) 2006 CAMP 

2# Please see the accompanying LICENSE file for further information. 

3 

4from math import pi, log, sqrt, factorial as fac 

5 

6import numpy as np 

7 

8"""Fourier filtering 

9 

10This module is an implementation of this Fourier filtering scheme: 

11 

12*A general and efficient pseudopotential Fourier filtering scheme for 

13real space methods using mask functions*, Maxim Tafipolsky, Rochus 

14Schmid, J Chem Phys. 2006 May 7;124:174102. 

15 

16Only difference is that we use a gaussian for the mask function. The 

17filtering is used for the projector functions and for the zero 

18potential.""" 

19 

20# 3D-Fourier transform: 

21# 

22# / _ _ 

23# ~ ^ | _ iq.r ^ 

24# f (q) Y (q) = | dr e f (r) Y (r) 

25# l lm | l lm 

26# / 

27# 

28# Radial part: 

29# 

30# / 

31# ~ __ l | 2 

32# f (q) = 4||i | r dr j (qr) f (r) 

33# l | l l 

34# / 

35# 

36 

37# XXX use fast bessel transform !!! 

38 

39 

40class Filter: 

41 """Mask-function Fourier filter""" 

42 

43 def __init__(self, r_g, dr_g, gcut, h): 

44 """Construct filter. 

45 

46 The radial grid is defined by r(g) and dr/dg(g) (`r_g` and 

47 `dr_g`), `gcut` is the cutoff grid point, and `h` is the target 

48 grid spacing used in the calculation.""" 

49 

50 self.gcut = gcut 

51 rcut = r_g[gcut] 

52 

53 N = 200 

54 self.r_g = r_g = r_g[:gcut].copy() # will be modified later! 

55 self.dr_g = dr_g[:gcut] 

56 

57 # Matrices for Bessel transform: 

58 q1 = 5 * pi / h / N 

59 self.q_i = q_i = q1 * np.arange(N) 

60 self.c = sqrt(2 * q1 / pi) 

61 self.sinqr_ig = np.sin(q_i[:, None] * r_g) * self.c 

62 self.cosqr_ig = np.cos(q_i[:, None] * r_g) * self.c 

63 

64 # Cutoff function: 

65 qmax = pi / h 

66 alpha = 1.1 

67 qcut = qmax / alpha 

68 icut = 1 + int(qcut / q1) 

69 beta = 5 * log(10) / (alpha - 1.0)**2 

70 self.cut_i = np.ones(N) 

71 self.cut_i[icut:] = np.exp( 

72 -np.clip(beta * (q_i[icut:] / qcut - 1.0)**2, 0, 400)) 

73 # self.cut_i[icut:] = np.exp( 

74 # -np.clip(0, 400, beta * (q_i[icut:] / qcut - 1.0)**2)) 

75 

76 # Mask function: 

77 gamma = 3 * log(10) / rcut**2 

78 self.m_g = np.exp(-gamma * r_g**2) 

79 

80 # We will need to divide by these two! Remove zeros: 

81 q_i[0] = 1.0 

82 r_g[0] = 1.0 

83 

84 def filter(self, f_g, l=0): 

85 """Filter radial function. 

86 

87 The function to be filtered is:: 

88 

89 f(r) ^ 

90 ---- Y (r) 

91 r lm 

92 

93 Output is:: 

94 

95 l ^ 

96 g(r) r Y (r), 

97 lm 

98 

99 where the filtered radial part ``g(r)`` is returned.""" 

100 

101 r_g = self.r_g 

102 q_i = self.q_i 

103 

104 fdrim_g = f_g[:self.gcut] * self.dr_g / self.m_g / r_g 

105 

106 # sin(x) 

107 # j (x) = ------, 

108 # 0 x 

109 # 

110 # sin(x) cos(x) 

111 # j (x) = ------ - ------, 

112 # 1 2 x 

113 # x 

114 # 

115 # 3 1 3 

116 # j (x) = (--- - -) sin(x) - --- cos(x), 

117 # 2 3 x 2 

118 # x x 

119 # 

120 # 15 6 15 1 

121 # j (x) = (-- - ---) sin(x) - (-- - -) cos(x). 

122 # 3 4 2 3 x 

123 # x x x 

124 # 

125 

126 if l == 0: 

127 fq_i = np.dot(self.sinqr_ig, fdrim_g * r_g) * self.cut_i 

128 fr_g = np.dot(fq_i, self.sinqr_ig) 

129 elif l == 1: 

130 fq_i = np.dot(self.sinqr_ig, fdrim_g) / q_i 

131 fq_i -= np.dot(self.cosqr_ig, r_g * fdrim_g) 

132 fq_i[0] = 0.0 

133 fq_i *= self.cut_i 

134 fr_g = np.dot(fq_i / q_i, self.sinqr_ig) / r_g 

135 fr_g -= np.dot(fq_i, self.cosqr_ig) 

136 elif l == 2: 

137 fq_i = 3 * np.dot(self.sinqr_ig, fdrim_g / r_g) / q_i**2 

138 fq_i -= np.dot(self.sinqr_ig, fdrim_g * r_g) 

139 fq_i -= 3 * np.dot(self.cosqr_ig, fdrim_g) / q_i 

140 fq_i[0] = 0.0 

141 fq_i *= self.cut_i 

142 fr_g = 3 * np.dot(fq_i / q_i**2, self.sinqr_ig) / r_g**2 

143 fr_g -= np.dot(fq_i, self.sinqr_ig) 

144 fr_g -= 3 * np.dot(fq_i / q_i, self.cosqr_ig) / r_g 

145 elif l == 3: # XXX This should be tested 

146 fq_i = 15 * np.dot(self.sinqr_ig, fdrim_g / r_g**2) / q_i**3 

147 fq_i -= 6 * np.dot(self.sinqr_ig, fdrim_g) / q_i 

148 fq_i -= 15 * np.dot(self.cosqr_ig, fdrim_g / r_g) / q_i**2 

149 fq_i += np.dot(self.cosqr_ig, r_g * fdrim_g) 

150 fq_i[0] = 0.0 

151 fq_i *= self.cut_i 

152 fr_g = 15 * np.dot(fq_i / q_i**3, self.sinqr_ig) / r_g**3 

153 fr_g -= 6 * np.dot(fq_i / q_i, self.sinqr_ig) / r_g 

154 fr_g -= 15 * np.dot(fq_i / q_i**2, self.cosqr_ig) / r_g**2 

155 fr_g += np.dot(fq_i, self.cosqr_ig) 

156 

157 else: 

158 raise NotImplementedError 

159 

160 a_g = np.zeros(len(f_g)) 

161 a_g[:self.gcut] = fr_g * self.m_g / r_g**(l + 1) 

162 

163 # n 

164 # 2 n! n 

165 # j (x) = --------- x for x << 1. 

166 # n (2n + 1)! 

167 # 

168 # This formula is used for finding the value of 

169 # 

170 # -l 

171 # f(r) r for r -> 0 

172 # 

173 c = 2.0**l * fac(l) / fac(2 * l + 1) * self.c 

174 a_g[0] = np.dot(fq_i, q_i**(l + 1)) * c 

175 

176 return a_g 

177 

178 

179if __name__ == '__main__': 

180 rc = 1.1 

181 gamma = 1.95 

182 rc2 = rc * gamma 

183 M = 300 

184 beta = 0.3 

185 gcut = 1 + int(M * rc / (beta + rc)) 

186 g_g = np.arange(M) 

187 r_g = beta * g_g / (M - g_g) 

188 drdg_g = beta * M / (M - g_g)**2 

189 

190 x_g = r_g / rc 

191 p_g = 1 - x_g**2 * (3 - 2 * x_g) 

192 p_g[gcut:] = 0.0 

193 # p_g = np.exp(-np.clip(5.0 * r_g**2, 0, 400)) 

194 

195 h = 0.4 

196 f = Filter(r_g, drdg_g, rc2, h) 

197 pf0_g = f.filter(p_g) 

198 pf1_g = f.filter(p_g * r_g**1, 1) 

199 pf2_g = f.filter(p_g * r_g**2, 2) 

200 

201 if 0: 

202 for i in range(200): 

203 print(5 * pi / h * i / 200, pf0_g[i], pf1_g[i], pf2_g[i]) 

204 if 1: 

205 for r, p, pf0, pf1, pf2 in zip(r_g, p_g, pf0_g, pf1_g, pf2_g): 

206 print(r, p, pf0, pf1, pf2) 

207 if r > rc2: 

208 break