Coverage for eminus/xc/utils.py: 84.12%

170 statements  

« prev     ^ index     » next       coverage.py v7.6.4, created at 2024-11-08 12:59 +0000

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Utility functions for exchange-correlation functionals.""" 

4 

5import inspect 

6import sys 

7 

8import numpy as np 

9 

10from .. import config 

11from ..logger import log 

12from ..utils import add_maybe_none 

13from .gga_c_chachiyo import gga_c_chachiyo, gga_c_chachiyo_spin 

14from .gga_c_pbe import gga_c_pbe, gga_c_pbe_spin 

15from .gga_c_pbe_sol import gga_c_pbe_sol, gga_c_pbe_sol_spin 

16from .gga_x_chachiyo import gga_x_chachiyo, gga_x_chachiyo_spin 

17from .gga_x_pbe import gga_x_pbe, gga_x_pbe_spin 

18from .gga_x_pbe_sol import gga_x_pbe_sol, gga_x_pbe_sol_spin 

19from .lda_c_chachiyo import lda_c_chachiyo, lda_c_chachiyo_spin 

20from .lda_c_chachiyo_mod import lda_c_chachiyo_mod, lda_c_chachiyo_mod_spin 

21from .lda_c_pw import lda_c_pw, lda_c_pw_spin 

22from .lda_c_pw_mod import lda_c_pw_mod, lda_c_pw_mod_spin 

23from .lda_c_vwn import lda_c_vwn, lda_c_vwn_spin 

24from .lda_x import lda_x, lda_x_spin 

25 

26 

27def get_xc(xc, n_spin, Nspin, dn_spin=None, tau=None, xc_params=None, dens_threshold=0): 

28 """Handle and get exchange-correlation functionals. 

29 

30 Args: 

31 xc: Exchange and correlation identifier. 

32 n_spin: Real-space electronic densities per spin channel. 

33 Nspin: Number of spin states. 

34 

35 Keyword Args: 

36 dn_spin: Real-space gradient of densities per spin channel. 

37 tau: Real-space kinetic energy densities per spin channel. 

38 xc_params: Exchange-correlation functional parameters. 

39 dens_threshold: Do not treat densities smaller than the threshold. 

40 

41 Returns: 

42 Exchange-correlation energy density and potential. 

43 """ 

44 if isinstance(xc, str): 

45 xc = parse_functionals(xc) 

46 f_exch, f_corr = xc 

47 if xc_params is None: 

48 xc_params = {} 

49 

50 # Only use non-zero values of the density 

51 n = np.sum(n_spin, axis=0) 

52 nz_mask = np.where(n > dens_threshold) 

53 n_nz = n[nz_mask] 

54 # Zeta is only needed for non-zero values of the density 

55 zeta_nz = get_zeta(n_spin[:, nz_mask]) 

56 # dn_spin is only needed for non-zero values of the density 

57 if dn_spin is not None: 

58 dn_spin_nz = dn_spin[:, nz_mask[0], :] 

59 else: 

60 dn_spin_nz = None 

61 

62 def handle_functional(fxc): 

63 """Calculate a given functional fxc, same for exchange and correlation.""" 

64 # Calculate with the libxc extra... 

65 if ":" in fxc: 

66 from ..extras.libxc import libxc_functional 

67 

68 fxc = fxc.split(":")[-1] 

69 exc, vxc, vsigma, vtau = libxc_functional(fxc, n_spin, Nspin, dn_spin, tau, xc_params) 

70 # ...or use an internal functional 

71 else: 

72 if Nspin == 2 and fxc != "mock_xc": 

73 fxc += "_spin" 

74 exc_nz, vxc_nz, vsigma_nz = IMPLEMENTED[fxc]( 

75 n_nz, zeta=zeta_nz, dn_spin=dn_spin_nz, Nspin=Nspin, **xc_params 

76 ) 

77 # Map the non-zero values back to the right dimension 

78 exc = np.zeros_like(n) 

79 exc[nz_mask] = exc_nz 

80 vxc = np.zeros_like(n_spin) 

81 for s in range(Nspin): 

82 vxc[s, nz_mask] = vxc_nz[s] 

83 if vsigma_nz is not None: 

84 vsigma = np.zeros((len(vsigma_nz), len(exc))) 

85 for i in range(len(vsigma)): 

86 vsigma[i, nz_mask] = vsigma_nz[i] 

87 else: 

88 vsigma = None 

89 # There are no internal meta-GGAs 

90 vtau = None 

91 return exc, vxc, vsigma, vtau 

92 

93 ex, vx, vsigmax, vtaux = handle_functional(f_exch) # Calculate the exchange part 

94 ec, vc, vsigmac, vtauc = handle_functional(f_corr) # Calculate the correlation part 

95 return ex + ec, vx + vc, add_maybe_none(vsigmax, vsigmac), add_maybe_none(vtaux, vtauc) 

96 

97 

98def get_exc(xc, n_spin, Nspin, dn_spin=None, tau=None, xc_params=None, dens_threshold=0): 

99 """Get the exchange-correlation energy density. 

100 

101 This is a convenience function to interface :func:`~eminus.xc.utils.get_xc`. 

102 

103 Args: 

104 xc: Exchange and correlation identifier. 

105 n_spin: Real-space electronic densities per spin channel. 

106 Nspin: Number of spin states. 

107 

108 Keyword Args: 

109 dn_spin: Real-space gradient of densities per spin channel. 

110 tau: Real-space kinetic energy densities per spin channel. 

111 xc_params: Exchange-correlation functional parameters. 

112 dens_threshold: Do not treat densities smaller than the threshold. 

113 

114 Returns: 

115 Exchange-correlation energy potential. 

116 """ 

117 exc, _, _, _ = get_xc(xc, n_spin, Nspin, dn_spin, tau, xc_params, dens_threshold) 

118 return exc 

119 

120 

121def get_vxc(xc, n_spin, Nspin, dn_spin=None, tau=None, xc_params=None, dens_threshold=0): 

122 """Get the exchange-correlation potential. 

123 

124 This is a convenience function to interface :func:`~eminus.xc.utils.get_xc`. 

125 

126 Args: 

127 xc: Exchange and correlation identifier. 

128 n_spin: Real-space electronic densities per spin channel. 

129 Nspin: Number of spin states. 

130 

131 Keyword Args: 

132 dn_spin: Real-space gradient of densities per spin channel. 

133 tau: Real-space kinetic energy densities per spin channel. 

134 xc_params: Exchange-correlation functional parameters. 

135 dens_threshold: Do not treat densities smaller than the threshold. 

136 

137 Returns: 

138 Exchange-correlation energy density. 

139 """ 

140 _, vxc, vsigma, vtau = get_xc(xc, n_spin, Nspin, dn_spin, tau, xc_params, dens_threshold) 

141 return vxc, vsigma, vtau 

142 

143 

144def parse_functionals(xc): 

145 """Parse exchange-correlation functional strings to the internal format. 

146 

147 Args: 

148 xc: Exchange and correlation identifier, separated by a comma. 

149 

150 Returns: 

151 Exchange and correlation string. 

152 """ 

153 # Check for combined aliases 

154 try: 

155 # Remove underscores when looking up in the dictionary 

156 xc_ = xc.replace("_", "") 

157 xc = ALIAS[xc_] 

158 except KeyError: 

159 pass 

160 

161 # Parse functionals 

162 functionals = [] 

163 for f in xc.split(","): 

164 if ":" in f or f in IMPLEMENTED: 

165 f_xc = f 

166 elif not f: 

167 f_xc = "mock_xc" 

168 else: 

169 try: 

170 # Remove underscores when looking up in the dictionary 

171 f_ = f.replace("_", "") 

172 f_xc = XC_MAP[f_] 

173 except KeyError: 

174 log.exception(f'No functional found for "{f}".') 

175 raise 

176 functionals.append(f_xc) 

177 

178 # If only one or no functional has been parsed append with mock functionals 

179 functionals.extend("mock_xc" for _ in range(2 - len(functionals))) 

180 return functionals 

181 

182 

183def parse_xc_type(xc): 

184 """Parse functional strings to identify the corresponding functional type. 

185 

186 Args: 

187 xc: Exchange and correlation identifier, separated by a comma. 

188 

189 Returns: 

190 Functional type. 

191 """ 

192 xc_type = [] 

193 for func in xc: 

194 if ":" in func: 

195 xc_id = func.split(":")[-1] 

196 # Try to parse the functional using pylibxc at first 

197 try: 

198 family = parse_xc_libxc(xc_id) 

199 # Otherwise parse it with PySCF 

200 except (ImportError, AssertionError): 

201 family = parse_xc_pyscf(xc_id) 

202 

203 if family == 1: 

204 xc_type.append("lda") 

205 elif family == 2: 

206 xc_type.append("gga") 

207 elif family == 4: 

208 xc_type.append("meta-gga") 

209 else: 

210 msg = "Unsupported functional family." 

211 raise NotImplementedError(msg) 

212 # Fall back to internal xc functionals 

213 elif "gga" in func: 

214 xc_type.append("gga") 

215 else: 

216 xc_type.append("lda") 

217 

218 # When mixing functional types use the higher level of theory 

219 if xc_type[0] != xc_type[1]: 

220 log.warning("Detected mixing of different functional types.") 

221 if "meta-gga" in xc_type: 

222 return "meta-gga" 

223 return "gga" 

224 return xc_type[0] 

225 

226 

227def parse_xc_libxc(xc_id): 

228 """Parse functional type by its ID using pylibxc. 

229 

230 Args: 

231 xc_id: Functional ID or identifier. 

232 

233 Returns: 

234 Functional type. 

235 """ 

236 if not config.use_pylibxc: 

237 raise AssertionError 

238 import pylibxc 

239 

240 if not xc_id.isdigit(): 

241 xc_id = pylibxc.util.xc_functional_get_number(xc_id) 

242 

243 func = pylibxc.LibXCFunctional(int(xc_id), 1) 

244 if func._needs_laplacian: 

245 msg = "meta-GGAs that need a laplacian are not supported." 

246 raise NotImplementedError(msg) 

247 return func.get_family() 

248 

249 

250def parse_xc_pyscf(xc_id): 

251 """Parse functional type by its ID using PySCF. 

252 

253 Args: 

254 xc_id: Functional ID or identifier. 

255 

256 Returns: 

257 Functional type. 

258 """ 

259 from pyscf.dft.libxc import is_gga, is_lda, is_meta_gga, needs_laplacian, XC_CODES 

260 

261 if not xc_id.isdigit(): 

262 xc_id = XC_CODES[xc_id.upper()] 

263 

264 if needs_laplacian(int(xc_id)): 

265 msg = "meta-GGAs that need a laplacian are not supported." 

266 raise NotImplementedError(msg) 

267 # Use the same values as in parse_xc_libxc 

268 if is_lda(xc_id): 

269 return 1 

270 if is_gga(xc_id): 

271 return 2 

272 if is_meta_gga(xc_id): 

273 return 4 

274 return -1 

275 

276 

277def get_xc_defaults(xc): 

278 """Get the default parameters and values for a given set of functionals. 

279 

280 Args: 

281 xc: Exchange and correlation identifier, separated by a comma. 

282 

283 Returns: 

284 Default parameters and values. 

285 """ 

286 if isinstance(xc, str): 

287 xc = parse_functionals(xc) 

288 

289 # Names of special kewyword arguments that should not be used in xc_params 

290 SPECIAL_NAMES = ["dn_spin", "Nspin"] 

291 

292 params = {} 

293 for func in xc: 

294 # If pylibxc functionals are used determine the default values from it 

295 if ":" in func: 

296 # This only works for pylibxc, not with PySCF 

297 if not config.use_pylibxc or "pylibxc" not in sys.modules: 

298 msg = "ext_params only work with pylibxc as the libxc backend, not with pyscf." 

299 raise NotImplementedError(msg) 

300 from pylibxc import LibXCFunctional 

301 

302 fxc = func.split(":")[-1] 

303 try: 

304 f_xc = LibXCFunctional(int(fxc), 1) 

305 except ValueError: 

306 f_xc = LibXCFunctional(fxc, 1) 

307 fxc_params = dict(zip(f_xc.get_ext_param_names(), f_xc.get_ext_param_default_values())) 

308 

309 # Analyze the signature for implemented functionals 

310 if func in IMPLEMENTED: 

311 sig = inspect.signature(IMPLEMENTED[func]) 

312 fxc_params = { 

313 param.name: param.default 

314 for param in sig.parameters.values() 

315 if param.default is not inspect.Parameter.empty 

316 } 

317 

318 # Remove special names from the parsed parameters 

319 for special in SPECIAL_NAMES: 

320 if special in fxc_params: 

321 del fxc_params[special] 

322 

323 # Append all parameters, warn if a parameter has been used before 

324 for name in fxc_params: 

325 if name in params: 

326 log.warning( 

327 f'External parameter "{name}" is used in the exchange AND correlation part. ' 

328 "It will be passed to both functionals if used in xc_params." 

329 ) 

330 params[name] = fxc_params[name] 

331 return params 

332 

333 

334def get_zeta(n_spin): 

335 """Calculate the relative spin polarization. 

336 

337 Args: 

338 n_spin: Real-space electronic densities per spin channel. 

339 

340 Returns: 

341 Relative spin polarization. 

342 """ 

343 # If only one spin is given return an array of ones as if the density only is in one channel 

344 if len(n_spin) == 1: 

345 return np.ones_like(n_spin[0]) 

346 return (n_spin[0] - n_spin[1]) / (n_spin[0] + n_spin[1]) 

347 

348 

349def mock_xc(n, Nspin=1, **kwargs): 

350 """Mock exchange-correlation functional with no effect (spin-paired). 

351 

352 Args: 

353 n: Real-space electronic density. 

354 Nspin: Number of spin states. 

355 

356 Keyword Args: 

357 **kwargs: Throwaway arguments. 

358 

359 Returns: 

360 Mock exchange-correlation energy density and potential. 

361 """ 

362 zeros = np.zeros_like(n) 

363 return zeros, np.array([zeros] * Nspin), None 

364 

365 

366#: Map functional names with their respective implementation. 

367IMPLEMENTED = { 

368 i.__name__: i 

369 for i in ( 

370 mock_xc, 

371 gga_c_chachiyo, 

372 gga_c_chachiyo_spin, 

373 gga_c_pbe, 

374 gga_c_pbe_spin, 

375 gga_c_pbe_sol, 

376 gga_c_pbe_sol_spin, 

377 gga_x_chachiyo, 

378 gga_x_chachiyo_spin, 

379 gga_x_pbe, 

380 gga_x_pbe_spin, 

381 gga_x_pbe_sol, 

382 gga_x_pbe_sol_spin, 

383 lda_x, 

384 lda_x_spin, 

385 lda_c_pw, 

386 lda_c_pw_spin, 

387 lda_c_pw_mod, 

388 lda_c_pw_mod_spin, 

389 lda_c_vwn, 

390 lda_c_vwn_spin, 

391 lda_c_chachiyo, 

392 lda_c_chachiyo_spin, 

393 lda_c_chachiyo_mod, 

394 lda_c_chachiyo_mod_spin, 

395 ) 

396} 

397 

398#: Map functional shorthands to the actual functional names. 

399XC_MAP = { 

400 # lda_x 

401 "1": "lda_x", 

402 "s": "lda_x", 

403 "lda": "lda_x", 

404 "slater": "lda_x", 

405 # lda_c_vwn 

406 "7": "lda_c_vwn", 

407 "vwn": "lda_c_vwn", 

408 "vwn5": "lda_c_vwn", 

409 # lda_c_pw 

410 "12": "lda_c_pw", 

411 "pw": "lda_c_pw", 

412 "pw92": "lda_c_pw", 

413 # lda_c_pw_mod 

414 "13": "lda_c_pw_mod", 

415 "pwmod": "lda_c_pw_mod", 

416 "pw92mod": "lda_c_pw_mod", 

417 # gga_x_pbe 

418 "101": "gga_x_pbe", 

419 "pbex": "gga_x_pbe", 

420 # gga_x_pbe_sol 

421 "116": "gga_x_pbe_sol", 

422 "pbesolx": "gga_x_pbe_sol", 

423 # gga_c_pbe 

424 "130": "gga_c_pbe", 

425 "pbec": "gga_c_pbe", 

426 # gga_c_pbe_sol 

427 "133": "gga_c_pbe_sol", 

428 "pbesolc": "gga_c_pbe_sol", 

429 # lda_c_chachiyo 

430 "287": "lda_c_chachiyo", 

431 "chachiyo": "lda_c_chachiyo", 

432 # gga_x_chachiyo 

433 "298": "gga_x_chachiyo", 

434 "chachiyox": "gga_x_chachiyo", 

435 # lda_c_chachiyo_mod 

436 "307": "lda_c_chachiyo_mod", 

437 "chachiyomod": "lda_c_chachiyo_mod", 

438 # gga_c_chachiyo 

439 "309": "gga_c_chachiyo", 

440 "chachiyoc": "gga_c_chachiyo", 

441} 

442 

443#: Dictionary of common functional aliases. 

444ALIAS = { 

445 "svwn": "slater,vwn5", 

446 "svwn5": "slater,vwn5", 

447 "spw92": "slater,pw92mod", 

448 "pbe": "pbex,pbec", 

449 "pbesol": "pbesolx,pbesolc", 

450 "chachiyo": "chachiyox,chachiyoc", 

451}