Coverage for gpaw/xc/lda.py: 99%

1from math import sqrt, pi 

2 

3import numpy as np 

4 

5from gpaw.new import trace 

6from gpaw.sphere.lebedev import Y_nL, weight_n 

7from gpaw.xc.functional import XCFunctional 

8 

9 

10def stress_integral(gd, a1_g, a2_g=None): 

11 return gd.integrate(a1_g, a2_g, global_integral=False) 

12 

13 

14def stress_lda_term(gd, e_g, n_sg, v_sg): 

15 if e_g is None: 

16 P = 0 

17 else: 

18 P = stress_integral(gd, e_g) 

19 for v_g, n_g in zip(v_sg, n_sg): 

20 P -= stress_integral(gd, v_g, n_g) 

21 return P * np.eye(3) 

22 

23 

24class LDARadialExpansion: 

25 def __init__(self, rcalc, collinear=True): 

26 self.rcalc = rcalc 

27 self.collinear = collinear 

28 

29 def __call__(self, rgd, D_sLq, n_qg, nc0_sg): 

30 n_sLg = np.dot(D_sLq, n_qg) 

31 if self.collinear: 

32 n_sLg[:, 0] += nc0_sg 

33 else: 

34 n_sLg[0, 0] += 4 * nc0_sg[0] 

35 

36 dEdD_sqL = np.zeros_like(np.transpose(D_sLq, (0, 2, 1))) 

37 

38 Lmax = n_sLg.shape[1] 

39 E = 0.0 

40 for n, Y_L in enumerate(Y_nL[:, :Lmax]): 

41 w = weight_n[n] 

42 

43 e_g, dedn_sg = self.rcalc(rgd, n_sLg, Y_L) 

44 dEdD_sqL += np.dot(rgd.dv_g * dedn_sg, 

45 n_qg.T)[:, :, np.newaxis] * (w * Y_L) 

46 E += w * rgd.integrate(e_g) 

47 return E, dEdD_sqL 

48 

49 

50@trace 

51def calculate_paw_correction(expansion, 

52 setup, D_sp, dEdD_sp=None, 

53 addcoredensity=True, a=None): 

54 xcc = setup.xc_correction 

55 if xcc is None: 

56 return 0.0 

57 

58 rgd = xcc.rgd 

59 nspins = len(D_sp) 

60 

61 if addcoredensity: 

62 nc0_sg = rgd.empty(nspins) 

63 nct0_sg = rgd.empty(nspins) 

64 nc0_sg[:] = sqrt(4 * pi) / nspins * xcc.nc_g 

65 nct0_sg[:] = sqrt(4 * pi) / nspins * xcc.nct_g 

66 if xcc.nc_corehole_g is not None and nspins == 2: 

67 nc0_sg[0] -= 0.5 * sqrt(4 * pi) * xcc.nc_corehole_g 

68 nc0_sg[1] += 0.5 * sqrt(4 * pi) * xcc.nc_corehole_g 

69 else: 

70 nc0_sg = 0 

71 nct0_sg = 0 

72 

73 D_sLq = np.inner(D_sp, xcc.B_pqL.T) 

74 

75 e, dEdD_sqL = expansion(rgd, D_sLq, xcc.n_qg, nc0_sg) 

76 et, dEtdD_sqL = expansion(rgd, D_sLq, xcc.nt_qg, nct0_sg) 

77 

78 if dEdD_sp is not None: 

79 dEdD_sp += np.inner((dEdD_sqL - dEtdD_sqL).reshape((nspins, -1)), 

80 xcc.B_pqL.reshape((len(xcc.B_pqL), -1))) 

81 

82 if addcoredensity: 

83 return e - et - xcc.e_xc0 

84 else: 

85 return e - et 

86 

87 

88class LDARadialCalculator: 

89 def __init__(self, kernel): 

90 self.kernel = kernel 

91 

92 def __call__(self, rgd, n_sLg, Y_L): 

93 nspins = len(n_sLg) 

94 n_sg = np.dot(Y_L, n_sLg) 

95 e_g = rgd.empty() 

96 dedn_sg = rgd.zeros(nspins) 

97 self.kernel.calculate(e_g, n_sg, dedn_sg) 

98 return e_g, dedn_sg 

99 

100 

101class LDA(XCFunctional): 

102 def __init__(self, kernel): 

103 self.kernel = kernel 

104 XCFunctional.__init__(self, kernel.name, kernel.type) 

105 

106 def calculate_impl(self, gd, n_sg, v_sg, e_g): 

107 self.kernel.calculate(e_g, n_sg, v_sg) 

108 

109 def calculate_paw_correction(self, setup, D_sp, dEdD_sp=None, 

110 addcoredensity=True, a=None): 

111 from gpaw.xc.noncollinear import NonCollinearLDAKernel 

112 collinear = not isinstance(self.kernel, NonCollinearLDAKernel) 

113 rcalc = LDARadialCalculator(self.kernel) 

114 expansion = LDARadialExpansion(rcalc, collinear) 

115 return calculate_paw_correction(expansion, 

116 setup, D_sp, dEdD_sp, 

117 addcoredensity, a) 

118 

119 def calculate_radial(self, rgd, n_sLg, Y_L): 

120 rcalc = LDARadialCalculator(self.kernel) 

121 return rcalc(rgd, n_sLg, Y_L) 

122 

123 def calculate_spherical(self, rgd, n_sg, v_sg, e_g=None): 

124 if e_g is None: 

125 e_g = rgd.empty() 

126 rcalc = LDARadialCalculator(self.kernel) 

127 e_g[:], dedn_sg = rcalc(rgd, n_sg[:, np.newaxis], [1.0]) 

128 v_sg[:] = dedn_sg 

129 return rgd.integrate(e_g) 

130 

131 def stress_tensor_contribution(self, n_sg, skip_sum=False): 

132 nspins = len(n_sg) 

133 v_sg = self.gd.zeros(nspins) 

134 e_g = self.gd.empty() 

135 self.calculate_impl(self.gd, n_sg, v_sg, e_g) 

136 stress_vv = stress_lda_term(self.gd, e_g, n_sg, v_sg) 

137 if not skip_sum: 

138 self.gd.comm.sum(stress_vv) 

139 return stress_vv 

140 

141 

142class PurePythonLDAKernel: 

143 def __init__(self): 

144 self.name = 'LDA' 

145 self.type = 'LDA' 

146 

147 def calculate(self, e_g, n_sg, dedn_sg, 

148 sigma_xg=None, dedsigma_xg=None, 

149 tau_sg=None, dedtau_sg=None): 

150 

151 e_g[:] = 0. 

152 if len(n_sg) == 1: 

153 n = n_sg[0] 

154 n[n < 1e-20] = 1e-40 

155 

156 # exchange 

157 lda_x(0, e_g, n, dedn_sg[0]) 

158 # correlation 

159 lda_c(0, e_g, n, dedn_sg[0], 0) 

160 

161 else: 

162 na = 2. * n_sg[0] 

163 na[na < 1e-20] = 1e-40 

164 nb = 2. * n_sg[1] 

165 nb[nb < 1e-20] = 1e-40 

166 n = 0.5 * (na + nb) 

167 zeta = 0.5 * (na - nb) / n 

168 

169 # exchange 

170 lda_x(1, e_g, na, dedn_sg[0]) 

171 lda_x(1, e_g, nb, dedn_sg[1]) 

172 # correlation 

173 lda_c(1, e_g, n, dedn_sg, zeta) 

174 

175 

176def lda_x(spin, e, n, v): 

177 assert spin in [0, 1] 

178 C0I, C1, CC1, CC2, IF2 = lda_constants() 

179 

180 rs = (C0I / n) ** (1 / 3.) 

181 ex = C1 / rs 

182 dexdrs = -ex / rs 

183 if spin == 0: 

184 e[:] += n * ex 

185 else: 

186 e[:] += 0.5 * n * ex 

187 v += ex - rs * dexdrs / 3. 

188 

189 

190def lda_c(spin, e, n, v, zeta): 

191 assert spin in [0, 1] 

192 C0I, C1, CC1, CC2, IF2 = lda_constants() 

193 

194 rs = (C0I / n) ** (1 / 3.) 

195 ec, decdrs_0 = G(rs ** 0.5, 

196 0.031091, 0.21370, 7.5957, 3.5876, 1.6382, 0.49294) 

197 

198 if spin == 0: 

199 e[:] += n * ec 

200 v += ec - rs * decdrs_0 / 3. 

201 else: 

202 e1, decdrs_1 = G(rs ** 0.5, 

203 0.015545, 0.20548, 14.1189, 6.1977, 3.3662, 0.62517) 

204 alpha, dalphadrs = G(rs ** 0.5, 

205 0.016887, 0.11125, 10.357, 3.6231, 0.88026, 

206 0.49671) 

207 alpha *= -1. 

208 dalphadrs *= -1. 

209 zp = 1.0 + zeta 

210 zm = 1.0 - zeta 

211 xp = zp ** (1 / 3.) 

212 xm = zm ** (1 / 3.) 

213 f = CC1 * (zp * xp + zm * xm - 2.0) 

214 f1 = CC2 * (xp - xm) 

215 zeta3 = zeta * zeta * zeta 

216 zeta4 = zeta * zeta * zeta * zeta 

217 x = 1.0 - zeta4 

218 decdrs = (decdrs_0 * (1.0 - f * zeta4) + 

219 decdrs_1 * f * zeta4 + 

220 dalphadrs * f * x * IF2) 

221 decdzeta = (4.0 * zeta3 * f * (e1 - ec - alpha * IF2) + 

222 f1 * (zeta4 * e1 - zeta4 * ec + x * alpha * IF2)) 

223 ec += alpha * IF2 * f * x + (e1 - ec) * f * zeta4 

224 e[:] += n * ec 

225 v[0] += ec - rs * decdrs / 3.0 - (zeta - 1.0) * decdzeta 

226 v[1] += ec - rs * decdrs / 3.0 - (zeta + 1.0) * decdzeta 

227 

228 

229def G(rtrs, gamma, alpha1, beta1, beta2, beta3, beta4): 

230 Q0 = -2.0 * gamma * (1.0 + alpha1 * rtrs * rtrs) 

231 Q1 = 2.0 * gamma * rtrs * (beta1 + 

232 rtrs * (beta2 + 

233 rtrs * (beta3 + 

234 rtrs * beta4))) 

235 G1 = Q0 * np.log(1.0 + 1.0 / Q1) 

236 dQ1drs = gamma * (beta1 / rtrs + 2.0 * beta2 + 

237 rtrs * (3.0 * beta3 + 4.0 * beta4 * rtrs)) 

238 dGdrs = -2.0 * gamma * alpha1 * G1 / Q0 - Q0 * dQ1drs / (Q1 * (Q1 + 1.0)) 

239 return G1, dGdrs 

240 

241 

242def lda_constants(): 

243 C0I = 0.238732414637843 

244 C1 = -0.45816529328314287 

245 CC1 = 1.9236610509315362 

246 CC2 = 2.5648814012420482 

247 IF2 = 0.58482236226346462 

248 return C0I, C1, CC1, CC2, IF2