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

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""GDSMFB LDA exchange-correlation. 

4 

5Reference: Phys. Rev. Lett. 119, 135001. 

6""" 

7 

8import dataclasses 

9import functools 

10 

11import numpy as np 

12 

13 

14def lda_xc_gdsmfb(n, T=0, **kwargs): 

15 """GDSMFB exchange-correlation functional (spin-paired). 

16 

17 Similar to the functional with the label LDA_XC_GDSMFB and ID 577 in Libxc, but the 

18 implementations differ: https://gitlab.com/libxc/libxc/-/issues/525 

19 Exchange and correlation cannot be separated. 

20 

21 Reference: Phys. Rev. Lett. 119, 135001. 

22 

23 Args: 

24 n: Real-space electronic density. 

25 

26 Keyword Args: 

27 T: Temperature in Hartree. 

28 **kwargs: Throwaway arguments. 

29 

30 Returns: 

31 GDSMFB exchange-correlation energy density and potential. 

32 """ 

33 kwargs["zeta"] = np.zeros_like(n) 

34 exc, vxc, _ = lda_xc_gdsmfb_spin(n, T=T, **kwargs) 

35 return exc, np.array([vxc[0]]), None 

36 

37 

38def lda_xc_gdsmfb_spin(n, zeta, T=0, **kwargs): 

39 """GDSMFB exchange-correlation functional (spin-polarized). 

40 

41 Similar to the functional with the label LDA_XC_GDSMFB and ID 577 in Libxc, but the 

42 implementations differ: https://gitlab.com/libxc/libxc/-/issues/525 

43 Exchange and correlation cannot be separated. 

44 

45 Reference: Phys. Rev. Lett. 119, 135001. 

46 

47 Args: 

48 n: Real-space electronic density. 

49 zeta: Relative spin polarization. 

50 

51 Keyword Args: 

52 T: Temperature in Hartree. 

53 **kwargs: Throwaway arguments. 

54 

55 Returns: 

56 GDSMFB exchange-correlation energy density and potential. 

57 """ 

58 # Calculate properties 

59 n_up = (1 + zeta) * n / 2 

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

61 theta = _get_theta(T, n, zeta) 

62 theta0 = _get_theta0(theta, zeta) 

63 theta1 = _get_theta1(theta, zeta) 

64 

65 # Initialize parameters 

66 # We need to calculate coefficients for specific theta 

67 # Create a coefficients object for each theta with its corresponding parameters 

68 phi_params = PhiParams() 

69 zeta0theta = Zeta0Coeffs(theta) 

70 zeta0theta0 = Zeta0Coeffs(theta0) 

71 zeta1theta = Zeta1Coeffs(theta) 

72 zeta1theta1 = Zeta1Coeffs(theta1) 

73 

74 # Calculate fxc 

75 fxc0 = _get_fxc_zeta(rs, zeta0theta0) 

76 fxc1 = _get_fxc_zeta(rs, zeta1theta1) 

77 phi = _get_phi(rs, theta0, zeta, phi_params) 

78 fxc = fxc0 + (fxc1 - fxc0) * phi 

79 

80 # Generic derivatives 

81 drsdn = -(6 ** (1 / 3)) * (1 / n) ** (1 / 3) / (6 * np.pi ** (1 / 3) * n) 

82 dzetadn_up = -zeta / n**2 + 1 / n 

83 dzetadn_dw = -zeta / n**2 - 1 / n 

84 

85 # fxc derivatives 

86 dfxc0drs = _get_dfxc_zetadrs(rs, zeta0theta) 

87 dfxc1drs = _get_dfxc_zetadrs(rs, zeta1theta) 

88 dfxc0dtheta0 = _get_dfxc_zetadtheta(rs, zeta0theta0) 

89 dfxc1dtheta1 = _get_dfxc_zetadtheta(rs, zeta1theta1) 

90 

91 # phi derivatives 

92 dphidrs = _get_dphidrs(rs, theta0, zeta, phi_params) 

93 dphidtheta = _get_dphidtheta(rs, theta0, zeta, phi_params) 

94 dphidzeta = _get_dphidzeta(rs, theta0, zeta, phi_params) 

95 

96 # theta derivatives 

97 dthetadn_up = _get_dthetadn_up(T, n_up) 

98 dtheta0dtheta = _get_dtheta0dtheta(zeta) 

99 dtheta0dzeta = _get_dtheta0dzeta(theta0, zeta) 

100 dtheta1dtheta0 = _get_dtheta1dtheta0() 

101 

102 # Calculate vxc_up (using dndn_up=1) 

103 dfxc0a = dfxc0drs * drsdn # * dndn_up = 1 

104 dfxc1a = dfxc1drs * drsdn # * dndn_up = 1 

105 dfxc0b_up = dfxc0dtheta0 * (dtheta0dtheta * dthetadn_up + dtheta0dzeta * dzetadn_up) 

106 dfxc1b_up = ( 

107 dfxc1dtheta1 * dtheta1dtheta0 * (dtheta0dtheta * dthetadn_up + dtheta0dzeta * dzetadn_up) 

108 ) 

109 dphi_up = dphidtheta * dthetadn_up + dphidzeta * dzetadn_up + dphidrs * drsdn # * dndn_up = 1 

110 vxc_up = ( 

111 dfxc0a 

112 + dfxc0b_up 

113 - phi * (dfxc0a + dfxc0b_up - dfxc1a - dfxc1b_up) 

114 - (fxc0 - fxc1) * dphi_up 

115 ) 

116 

117 # Calculate vxc_dw (using dndn_dw=1 and dthetadn_dw=0) 

118 # dfxc0a and dfxc1a are identical for vxc_up and vxc_dw 

119 dfxc0b_dw = dfxc0dtheta0 * dtheta0dzeta * dzetadn_dw 

120 # dfxc0b += dfxc0dtheta0 * dtheta0dtheta * dthetadn_dw 

121 dfxc1b_dw = dfxc1dtheta1 * dtheta1dtheta0 * dtheta0dzeta * dzetadn_dw 

122 # dfxc1b += dfxc1dtheta1 * dtheta1dtheta0 * dtheta0dtheta * dthetadn_dw 

123 dphi_dw = dphidzeta * dzetadn_dw + dphidrs * drsdn # * dndn_dw = 1 

124 # dfxc1b += dphidtheta * dthetadn_dw 

125 vxc_dw = ( 

126 dfxc0a 

127 + dfxc0b_dw 

128 - phi * (dfxc0a + dfxc0b_dw - dfxc1a - dfxc1b_dw) 

129 - (fxc0 - fxc1) * dphi_dw 

130 ) 

131 

132 return fxc, fxc + np.array([vxc_up, vxc_dw]) * n, None 

133 

134 

135# ### Temperature dependent coefficients ### 

136 

137 

138@dataclasses.dataclass 

139class Coefficients: 

140 """Coefficients class to calculate temperature/theta dependent coefficients. 

141 

142 Reference: Phys. Rev. Lett. 119, 135001. 

143 """ 

144 

145 theta: float #: Reduced temperature. 

146 a0: float = 0.610887 

147 a1: float = 0.75 

148 a2: float = 3.04363 

149 a3: float = -0.09227 

150 a4: float = 1.7035 

151 a5: float = 8.31051 

152 a6: float = 5.1105 

153 

154 @functools.cached_property 

155 def b5(self): 

156 """Calculate b5.""" 

157 return self.b3 * np.sqrt(3 / 2) * self.omega * (4 / (9 * np.pi)) ** (-1 / 3) 

158 

159 @functools.cached_property 

160 def a(self): 

161 """Calculate a.""" 

162 with np.errstate(divide="ignore"): 

163 u = self.a0 * np.tanh( 

164 1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0 

165 ) 

166 return u * _pade(self.theta, self.a1, self.a2, self.a3, self.a4, self.a5, self.a6) 

167 

168 @property 

169 def dadtheta(self): 

170 """Calculate da / dtheta.""" 

171 with np.errstate(divide="ignore"): 

172 u = self.a0 * np.tanh( 

173 1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0 

174 ) 

175 du = np.divide( 

176 u**2 / self.a0 - self.a0, 

177 self.theta**2, 

178 out=np.zeros_like(self.theta), 

179 where=self.theta > 0, 

180 ) 

181 v, dv = _dpade(self.theta, self.a1, self.a2, self.a3, self.a4, self.a5, self.a6) 

182 return du * v + u * dv 

183 

184 @functools.cached_property 

185 def b(self): 

186 """Calculate b.""" 

187 with np.errstate(divide="ignore"): 

188 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0) 

189 return u * _pade(self.theta, self.b1, self.b2, 0, self.b3, self.b4, self.b5) 

190 

191 @property 

192 def dbdtheta(self): 

193 """Calculate db / dtheta.""" 

194 with np.errstate(divide="ignore"): 

195 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0) 

196 du = np.divide( 

197 u**2 - 1, 

198 2 * self.theta * np.sqrt(self.theta), 

199 out=np.zeros_like(self.theta), 

200 where=self.theta > 0, 

201 ) 

202 v, dv = _dpade(self.theta, self.b1, self.b2, 0, self.b3, self.b4, self.b5) 

203 return du * v + u * dv 

204 

205 @functools.cached_property 

206 def c(self): 

207 """Calculate c.""" 

208 with np.errstate(divide="ignore"): 

209 exp = np.exp(-1 / self.theta, out=np.zeros_like(self.theta), where=self.theta > 0) 

210 return (self.c1 + self.c2 * exp) * self.e 

211 

212 @property 

213 def dcdtheta(self): 

214 """Calculate dc / dtheta.""" 

215 with np.errstate(divide="ignore"): 

216 exp = np.exp(-1 / self.theta, out=np.zeros_like(self.theta), where=self.theta > 0) 

217 u = self.c1 + self.c2 * exp 

218 du = np.divide( 

219 self.c2 * exp, self.theta**2, out=np.zeros_like(self.theta), where=self.theta > 0 

220 ) 

221 v, dv = self.e, self.dedtheta 

222 return du * v + u * dv 

223 

224 @functools.cached_property 

225 def d(self): 

226 """Calculate d.""" 

227 with np.errstate(divide="ignore"): 

228 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0) 

229 return u * _pade(self.theta, self.d1, self.d2, 0, self.d3, self.d4, self.d5) 

230 

231 @property 

232 def dddtheta(self): 

233 """Calculate dd / dtheta.""" 

234 with np.errstate(divide="ignore"): 

235 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0) 

236 du = np.divide( 

237 u**2 - 1, 

238 2 * self.theta * np.sqrt(self.theta), 

239 out=np.zeros_like(self.theta), 

240 where=self.theta > 0, 

241 ) 

242 v, dv = _dpade(self.theta, self.d1, self.d2, 0, self.d3, self.d4, self.d5) 

243 return du * v + u * dv 

244 

245 @functools.cached_property 

246 def e(self): 

247 """Calculate e.""" 

248 with np.errstate(divide="ignore"): 

249 u = np.tanh(1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0) 

250 return u * _pade(self.theta, self.e1, self.e2, 0, self.e3, self.e4, self.e5) 

251 

252 @functools.cached_property 

253 def dedtheta(self): 

254 """Calculate de / dtheta.""" 

255 with np.errstate(divide="ignore"): 

256 u = np.tanh(1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0) 

257 du = np.divide(u**2 - 1, self.theta**2, out=np.zeros_like(self.theta), where=self.theta > 0) 

258 v, dv = _dpade(self.theta, self.e1, self.e2, 0, self.e3, self.e4, self.e5) 

259 return du * v + u * dv 

260 

261 

262# ### Parameters ### 

263 

264 

265@dataclasses.dataclass 

266class Zeta0Coeffs(Coefficients): 

267 """Coefficient class using the parameters for zeta=0. 

268 

269 Reference: Phys. Rev. Lett. 119, 135001. 

270 """ 

271 

272 omega: float = 1 

273 b1: float = 0.3436902 

274 b2: float = 7.82159531356 

275 b3: float = 0.300483986662 

276 b4: float = 15.8443467125 

277 c1: float = 0.8759442 

278 c2: float = -0.230130843551 

279 d1: float = 0.72700876 

280 d2: float = 2.38264734144 

281 d3: float = 0.30221237251 

282 d4: float = 4.39347718395 

283 d5: float = 0.729951339845 

284 e1: float = 0.25388214 

285 e2: float = 0.815795138599 

286 e3: float = 0.0646844410481 

287 e4: float = 15.0984620477 

288 e5: float = 0.230761357474 

289 

290 

291@dataclasses.dataclass 

292class Zeta1Coeffs(Coefficients): 

293 """Coefficient class using the parameters for zeta=1. 

294 

295 Reference: Phys. Rev. Lett. 119, 135001. 

296 """ 

297 

298 omega: float = 2 ** (1 / 3) 

299 b1: float = 0.84987704 

300 b2: float = 3.04033012073 

301 b3: float = 0.0775730131248 

302 b4: float = 7.57703592489 

303 c1: float = 0.91126873 

304 c2: float = -0.0307957123308 

305 d1: float = 1.48658718 

306 d2: float = 4.92684905511 

307 d3: float = 0.0849387225179 

308 d4: float = 8.3269821188 

309 d5: float = 0.218864952126 

310 e1: float = 0.27454097 

311 e2: float = 0.400994856555 

312 e3: float = 2.88773194962 

313 e4: float = 6.33499237092 

314 e5: float = 24.823008753 

315 

316 

317@dataclasses.dataclass 

318class PhiParams: 

319 """Parameter class holding the spin-interpolation function parameters. 

320 

321 Reference: Phys. Rev. Lett. 119, 135001. 

322 """ 

323 

324 # Sign of parameters is different from the supplemental material 

325 h1: float = 3.18747258 

326 h2: float = 7.74662802 

327 lambda1: float = 1.85909536 

328 lambda2: float = 0 

329 

330 

331# ### Pade approximation and derivative ### 

332 

333 

334def _pade(x, n1, n2, n3, n4, d1, d2): 

335 """Pade approximation. 

336 

337 Not the general case but as needed for this functional. 

338 """ 

339 x2, x4 = x**2, x**4 

340 num = n1 + n2 * x2 + n3 * x**3 + n4 * x4 

341 denom = 1 + d1 * x2 + d2 * x4 

342 return num / denom 

343 

344 

345def _dpade(x, n1, n2, n3, n4, d1, d2): 

346 """Pade approximation and its derivative.""" 

347 x2, x3, x4 = x**2, x**3, x**4 

348 num = n1 + n2 * x2 + n3 * x3 + n4 * x4 

349 denom = 1 + d1 * x2 + d2 * x4 

350 dnum = 2 * n2 * x + 3 * n3 * x2 + 4 * n4 * x3 

351 ddenom = 2 * d1 * x + 4 * d2 * x3 

352 # df = (a'b - ab') / b^2 

353 return num / denom, (dnum * denom - num * ddenom) / denom**2 

354 

355 

356# ### theta and derivatives ### 

357 

358 

359def _get_theta(T, n, zeta): 

360 """Calculate the reduced temperature theta. 

361 

362 Reference: Phys. Rev. Lett. 119, 135001. 

363 Only mentioned in the arXiv version: https://arxiv.org/abs/1703.08074 

364 """ 

365 n_up = (1 + zeta) * n / 2 

366 T_fermi = (6 * np.pi**2 * n_up) ** (2 / 3) / 2 

367 return T / T_fermi 

368 

369 

370def _get_dthetadn_up(T, n_up): 

371 """Calculate dtheta / dn_up.""" 

372 return -4 / (3 * 6 ** (2 / 3)) * T / (np.pi ** (4 / 3) * n_up ** (5 / 3)) 

373 

374 

375def _get_theta0(theta, zeta): 

376 """Calculate theta0. 

377 

378 Reference: Phys. Rev. Lett. 119, 135001. 

379 """ 

380 return theta * (1 + zeta) ** (2 / 3) 

381 

382 

383def _get_dtheta0dtheta(zeta): 

384 """Calculate dtheta0 / dtheta.""" 

385 return (1 + zeta) ** (2 / 3) 

386 

387 

388def _get_dtheta0dzeta(theta, zeta): 

389 """Calculate dtheta0 / dzeta.""" 

390 return 2 * theta / (3 * (1 + zeta) ** (1 / 3)) 

391 

392 

393def _get_theta1(theta, zeta): 

394 """Calculate theta1. 

395 

396 Reference: Phys. Rev. Lett. 119, 135001. 

397 

398 It is not explicitly mentioned but introduced as used in Eq. (5). 

399 """ 

400 theta0 = _get_theta0(theta, zeta) 

401 return theta0 * 2 ** (-2 / 3) 

402 

403 

404def _get_dtheta1dtheta0(): 

405 """Calculate dtheta1 / dtheta0.""" 

406 return 2 ** (-2 / 3) 

407 

408 

409# ### fxc_zeta and derivatives ### 

410 

411 

412def _get_fxc_zeta(rs, p): 

413 """Calculate the Pade formula f_xc^zeta. 

414 

415 Reference: Phys. Rev. Lett. 119, 135001. 

416 """ 

417 num = p.a * p.omega + p.b * np.sqrt(rs) + p.c * rs 

418 denom = 1 + p.d * np.sqrt(rs) + p.e * rs 

419 return -1 / rs * num / denom 

420 

421 

422def _get_dfxc_zetadrs(rs, p): 

423 """Calculate dfxc_zeta / drs.""" 

424 num = p.a * p.omega + p.b * np.sqrt(rs) + p.c * rs 

425 denom = 1 + p.d * np.sqrt(rs) + p.e * rs 

426 tmp1 = (p.d / (2 * np.sqrt(rs)) + p.e) * num / denom**2 

427 tmp2 = (p.b / (2 * np.sqrt(rs)) + p.c) / denom 

428 return (tmp1 - tmp2) / rs + num / denom / rs**2 

429 

430 

431def _get_dfxc_zetadtheta(rs, p): 

432 """Calculate dfxc_zeta / dzeta.""" 

433 num = p.a * p.omega + p.b * np.sqrt(rs) + p.c * rs 

434 denom = 1 + p.d * np.sqrt(rs) + p.e * rs 

435 tmp1 = (p.dddtheta * np.sqrt(rs) + p.dedtheta * rs) * num / denom**2 

436 tmp2 = (p.dadtheta * p.omega + p.dbdtheta * np.sqrt(rs) + p.dcdtheta * rs) / denom 

437 return (tmp1 - tmp2) / rs 

438 

439 

440# ### phi and derivatives ### 

441 

442 

443def _get_phi(rs, theta, zeta, phi_params): 

444 """Calculate the interpolation function phi. 

445 

446 Reference: Phys. Rev. Lett. 119, 135001. 

447 """ 

448 alpha = _get_alpha(rs, theta, phi_params) 

449 return ((1 + zeta) ** alpha + (1 - zeta) ** alpha - 2) / (2**alpha - 2) 

450 

451 

452def _get_dphidrs(rs, theta, zeta, phi_params): 

453 """Calculate dphi / drs.""" 

454 alpha = _get_alpha(rs, theta, phi_params) 

455 dalphadrs = _get_dalphadrs(rs, theta, phi_params) 

456 tmp1 = (1 - zeta) ** alpha * np.log(1 - zeta, out=np.zeros_like(zeta), where=1 - zeta > 0) 

457 tmp2 = (1 + zeta) ** alpha * np.log(1 + zeta) 

458 duv = (tmp1 + tmp2) * (2**alpha - 2) 

459 udv = ((1 - zeta) ** alpha + (1 + zeta) ** alpha - 2) * 2**alpha * np.log(2) 

460 vv = (2**alpha - 2) ** 2 

461 return (duv - udv) * dalphadrs / vv 

462 

463 

464def _get_dphidtheta(rs, theta, zeta, phi_params): 

465 """Calculate dphi / dtheta.""" 

466 alpha = _get_alpha(rs, theta, phi_params) 

467 dalphadtheta = _get_dalphadtheta(rs, theta, phi_params) 

468 tmp1 = (1 - zeta) ** alpha * np.log(1 - zeta, out=np.zeros_like(zeta), where=1 - zeta > 0) 

469 tmp2 = (1 + zeta) ** alpha * np.log(1 + zeta) 

470 duv = (tmp1 + tmp2) * (2**alpha - 2) 

471 udv = ((1 - zeta) ** alpha + (1 + zeta) ** alpha - 2) * 2**alpha * np.log(2) 

472 vv = (2**alpha - 2) ** 2 

473 return (duv - udv) * dalphadtheta / vv 

474 

475 

476def _get_dphidzeta(rs, theta, zeta, phi_params): 

477 """Calculate dphi / dzeta.""" 

478 alpha = _get_alpha(rs, theta, phi_params) 

479 tmp1 = alpha * (1 + zeta) ** alpha / (1 + zeta) 

480 tmp2 = np.divide( 

481 alpha * (1 - zeta) ** alpha, 1 - zeta, out=np.zeros_like(zeta), where=1 - zeta > 0 

482 ) 

483 return (tmp1 - tmp2) / (2**alpha - 2) 

484 

485 

486# ### alpha and derivatives ### 

487 

488 

489def _get_alpha(rs, theta, phi_params): 

490 """Calculate alpha. 

491 

492 Reference: Phys. Rev. Lett. 119, 135001. 

493 """ 

494 h = _get_h(rs, phi_params) 

495 lamda = _get_lambda(rs, theta, phi_params) 

496 return 2 - h * np.exp(-theta * lamda) 

497 

498 

499def _get_dalphadrs(rs, theta, phi_params): 

500 """Calculate dalpha / drs.""" 

501 h = _get_h(rs, phi_params) 

502 lamda = _get_lambda(rs, theta, phi_params) 

503 dhdrs = _get_dhdrs(rs, phi_params) 

504 dlambdadrs = _get_dlambdadrs(rs, theta, phi_params) 

505 return -dhdrs * np.exp(-theta * lamda) + dlambdadrs * theta * h * np.exp(-theta * lamda) 

506 

507 

508def _get_dalphadtheta(rs, theta, phi_params): 

509 """Calculate dalpha / dtheta.""" 

510 h = _get_h(rs, phi_params) 

511 lamda = _get_lambda(rs, theta, phi_params) 

512 dlambdadtheta = _get_dlambdadtheta(rs, phi_params) 

513 return (dlambdadtheta * theta + lamda) * h * np.exp(-theta * lamda) 

514 

515 

516# ### h and derivative ### 

517 

518 

519def _get_h(rs, phi_params): 

520 """Calculate h. 

521 

522 Reference: Phys. Rev. Lett. 119, 135001. 

523 """ 

524 return (2 / 3 + phi_params.h1 * rs) / (1 + phi_params.h2 * rs) 

525 

526 

527def _get_dhdrs(rs, phi_params): 

528 """Calculate dh / drs.""" 

529 return ( 

530 phi_params.h1 / (phi_params.h2 * rs + 1) 

531 - phi_params.h2 * (phi_params.h1 * rs + 2 / 3) / (phi_params.h2 * rs + 1) ** 2 

532 ) 

533 

534 

535# ### lambda and derivatives ### 

536 

537 

538def _get_lambda(rs, theta, phi_params): 

539 """Calculate lambda. 

540 

541 Reference: Phys. Rev. Lett. 119, 135001. 

542 """ 

543 return phi_params.lambda1 + phi_params.lambda2 * theta * rs ** (1 / 2) 

544 

545 

546def _get_dlambdadrs(rs, theta, phi_params): 

547 """Calculate dlambda / drs.""" 

548 return phi_params.lambda2 * theta / (2 * np.sqrt(rs)) 

549 

550 

551def _get_dlambdadtheta(rs, phi_params): 

552 """Calculate dlambda / dtheta.""" 

553 return phi_params.lambda2 * np.sqrt(rs)