Coverage for eminus/xc/utils.py: 84.39%
173 statements
« prev ^ index » next coverage.py v7.8.2, created at 2025-06-02 10:16 +0000
1# SPDX-FileCopyrightText: 2023 The eminus developers
2# SPDX-License-Identifier: Apache-2.0
3"""Utility functions for exchange-correlation functionals."""
5import inspect
6import sys
8import numpy as np
10from eminus import config
11from eminus.logger import log
12from eminus.utils import add_maybe_none
14from .gga_c_chachiyo import gga_c_chachiyo, gga_c_chachiyo_spin
15from .gga_c_pbe import gga_c_pbe, gga_c_pbe_spin
16from .gga_c_pbe_sol import gga_c_pbe_sol, gga_c_pbe_sol_spin
17from .gga_x_chachiyo import gga_x_chachiyo, gga_x_chachiyo_spin
18from .gga_x_pbe import gga_x_pbe, gga_x_pbe_spin
19from .gga_x_pbe_sol import gga_x_pbe_sol, gga_x_pbe_sol_spin
20from .lda_c_chachiyo import lda_c_chachiyo, lda_c_chachiyo_spin
21from .lda_c_chachiyo_mod import lda_c_chachiyo_mod, lda_c_chachiyo_mod_spin
22from .lda_c_pw import lda_c_pw, lda_c_pw_spin
23from .lda_c_pw_mod import lda_c_pw_mod, lda_c_pw_mod_spin
24from .lda_c_vwn import lda_c_vwn, lda_c_vwn_spin
25from .lda_x import lda_x, lda_x_spin
26from .lda_xc_corr_ksdt import lda_xc_corr_ksdt
27from .lda_xc_gdsmfb import lda_xc_gdsmfb, lda_xc_gdsmfb_spin
28from .lda_xc_ksdt import lda_xc_ksdt, lda_xc_ksdt_spin
31def get_xc(xc, n_spin, Nspin, dn_spin=None, tau=None, xc_params=None, dens_threshold=0):
32 """Handle and get exchange-correlation functionals.
34 Args:
35 xc: Exchange and correlation identifier.
36 n_spin: Real-space electronic densities per spin channel.
37 Nspin: Number of spin states.
39 Keyword Args:
40 dn_spin: Real-space gradient of densities per spin channel.
41 tau: Real-space kinetic energy densities per spin channel.
42 xc_params: Exchange-correlation functional parameters.
43 dens_threshold: Do not treat densities smaller than the threshold.
45 Returns:
46 Exchange-correlation energy density and potential.
47 """
48 if isinstance(xc, str):
49 xc = parse_functionals(xc)
50 f_exch, f_corr = xc
51 if xc_params is None:
52 xc_params = {}
54 # Only use non-zero values of the density
55 n = np.sum(n_spin, axis=0)
56 nz_mask = np.where(n > dens_threshold)
57 n_nz = n[nz_mask]
58 # Zeta is only needed for non-zero values of the density
59 zeta_nz = get_zeta(n_spin[:, nz_mask])
60 # dn_spin is only needed for non-zero values of the density
61 if dn_spin is not None:
62 dn_spin_nz = dn_spin[:, nz_mask[0], :]
63 else:
64 dn_spin_nz = None
66 def handle_functional(fxc):
67 """Calculate a given functional fxc, same for exchange and correlation."""
68 # Calculate with the libxc extra...
69 if ":" in fxc:
70 from eminus.extras.libxc import libxc_functional
72 fxc = fxc.split(":")[-1]
73 exc, vxc, vsigma, vtau = libxc_functional(fxc, n_spin, Nspin, dn_spin, tau, xc_params)
74 # ...or use an internal functional
75 else:
76 if Nspin == 2 and fxc != "mock_xc":
77 fxc += "_spin"
78 exc_nz, vxc_nz, vsigma_nz = IMPLEMENTED[fxc](
79 n_nz, zeta=zeta_nz, dn_spin=dn_spin_nz, Nspin=Nspin, **xc_params
80 )
81 # Map the non-zero values back to the right dimension
82 exc = np.zeros_like(n)
83 exc[nz_mask] = exc_nz
84 vxc = np.zeros_like(n_spin)
85 for s in range(Nspin):
86 vxc[s, nz_mask] = vxc_nz[s]
87 if vsigma_nz is not None:
88 vsigma = np.zeros((len(vsigma_nz), len(exc)))
89 for i in range(len(vsigma)):
90 vsigma[i, nz_mask] = vsigma_nz[i]
91 else:
92 vsigma = None
93 # There are no internal meta-GGAs
94 vtau = None
95 return exc, vxc, vsigma, vtau
97 ex, vx, vsigmax, vtaux = handle_functional(f_exch) # Calculate the exchange part
98 ec, vc, vsigmac, vtauc = handle_functional(f_corr) # Calculate the correlation part
99 return ex + ec, vx + vc, add_maybe_none(vsigmax, vsigmac), add_maybe_none(vtaux, vtauc)
102def get_exc(xc, n_spin, Nspin, dn_spin=None, tau=None, xc_params=None, dens_threshold=0):
103 """Get the exchange-correlation energy density.
105 This is a convenience function to interface :func:`~eminus.xc.utils.get_xc`.
107 Args:
108 xc: Exchange and correlation identifier.
109 n_spin: Real-space electronic densities per spin channel.
110 Nspin: Number of spin states.
112 Keyword Args:
113 dn_spin: Real-space gradient of densities per spin channel.
114 tau: Real-space kinetic energy densities per spin channel.
115 xc_params: Exchange-correlation functional parameters.
116 dens_threshold: Do not treat densities smaller than the threshold.
118 Returns:
119 Exchange-correlation energy potential.
120 """
121 exc, _, _, _ = get_xc(xc, n_spin, Nspin, dn_spin, tau, xc_params, dens_threshold)
122 return exc
125def get_vxc(xc, n_spin, Nspin, dn_spin=None, tau=None, xc_params=None, dens_threshold=0):
126 """Get the exchange-correlation potential.
128 This is a convenience function to interface :func:`~eminus.xc.utils.get_xc`.
130 Args:
131 xc: Exchange and correlation identifier.
132 n_spin: Real-space electronic densities per spin channel.
133 Nspin: Number of spin states.
135 Keyword Args:
136 dn_spin: Real-space gradient of densities per spin channel.
137 tau: Real-space kinetic energy densities per spin channel.
138 xc_params: Exchange-correlation functional parameters.
139 dens_threshold: Do not treat densities smaller than the threshold.
141 Returns:
142 Exchange-correlation energy density.
143 """
144 _, vxc, vsigma, vtau = get_xc(xc, n_spin, Nspin, dn_spin, tau, xc_params, dens_threshold)
145 return vxc, vsigma, vtau
148def parse_functionals(xc):
149 """Parse exchange-correlation functional strings to the internal format.
151 Args:
152 xc: Exchange and correlation identifier, separated by a comma.
154 Returns:
155 Exchange and correlation string.
156 """
157 # Check for combined aliases
158 try:
159 # Remove underscores when looking up in the dictionary
160 xc_ = xc.replace("_", "")
161 xc = ALIAS[xc_]
162 except KeyError:
163 pass
165 # Parse functionals
166 functionals = []
167 for f in xc.split(","):
168 if ":" in f or f in IMPLEMENTED:
169 f_xc = f
170 elif not f:
171 f_xc = "mock_xc"
172 else:
173 try:
174 # Remove underscores when looking up in the dictionary
175 f_ = f.replace("_", "")
176 f_xc = XC_MAP[f_]
177 except KeyError:
178 log.exception(f'No functional found for "{f}".')
179 raise
180 functionals.append(f_xc)
182 # If only one or no functional has been parsed append with mock functionals
183 functionals.extend("mock_xc" for _ in range(2 - len(functionals)))
184 return functionals
187def parse_xc_type(xc):
188 """Parse functional strings to identify the corresponding functional type.
190 Args:
191 xc: Exchange and correlation identifier, separated by a comma.
193 Returns:
194 Functional type.
195 """
196 xc_type = []
197 for func in xc:
198 if ":" in func:
199 xc_id = func.split(":")[-1]
200 # Try to parse the functional using pylibxc at first
201 try:
202 family = parse_xc_libxc(xc_id)
203 # Otherwise parse it with PySCF
204 except (ImportError, AssertionError):
205 family = parse_xc_pyscf(xc_id)
207 if family == 1:
208 xc_type.append("lda")
209 elif family == 2:
210 xc_type.append("gga")
211 elif family == 4:
212 xc_type.append("meta-gga")
213 else:
214 msg = "Unsupported functional family."
215 raise NotImplementedError(msg)
216 # Fall back to internal xc functionals
217 elif "gga" in func:
218 xc_type.append("gga")
219 else:
220 xc_type.append("lda")
222 # When mixing functional types use the higher level of theory
223 if xc_type[0] != xc_type[1]:
224 log.warning("Detected mixing of different functional types.")
225 if "meta-gga" in xc_type:
226 return "meta-gga"
227 return "gga"
228 return xc_type[0]
231def parse_xc_libxc(xc_id):
232 """Parse functional type by its ID using pylibxc.
234 Args:
235 xc_id: Functional ID or identifier.
237 Returns:
238 Functional type.
239 """
240 if not config.use_pylibxc:
241 raise AssertionError
242 import pylibxc
244 if not xc_id.isdigit():
245 xc_id = pylibxc.util.xc_functional_get_number(xc_id)
247 func = pylibxc.LibXCFunctional(int(xc_id), 1)
248 if func._needs_laplacian:
249 msg = "meta-GGAs that need a laplacian are not supported."
250 raise NotImplementedError(msg)
251 return func.get_family()
254def parse_xc_pyscf(xc_id):
255 """Parse functional type by its ID using PySCF.
257 Args:
258 xc_id: Functional ID or identifier.
260 Returns:
261 Functional type.
262 """
263 from pyscf.dft.libxc import is_gga, is_lda, is_meta_gga, needs_laplacian, XC_CODES
265 if not xc_id.isdigit():
266 xc_id = XC_CODES[xc_id.upper()]
268 if needs_laplacian(int(xc_id)):
269 msg = "meta-GGAs that need a laplacian are not supported."
270 raise NotImplementedError(msg)
271 # Use the same values as in parse_xc_libxc
272 if is_lda(xc_id):
273 return 1
274 if is_gga(xc_id):
275 return 2
276 if is_meta_gga(xc_id):
277 return 4
278 return -1
281def get_xc_defaults(xc):
282 """Get the default parameters and values for a given set of functionals.
284 Args:
285 xc: Exchange and correlation identifier, separated by a comma.
287 Returns:
288 Default parameters and values.
289 """
290 if isinstance(xc, str):
291 xc = parse_functionals(xc)
293 # Names of special kewyword arguments that should not be used in xc_params
294 SPECIAL_NAMES = ["dn_spin", "Nspin"]
296 params = {}
297 for func in xc:
298 # If pylibxc functionals are used determine the default values from it
299 if ":" in func:
300 # This only works for pylibxc, not with PySCF
301 if not config.use_pylibxc or "pylibxc" not in sys.modules:
302 msg = "ext_params only work with pylibxc as the libxc backend, not with pyscf."
303 raise NotImplementedError(msg)
304 from pylibxc import LibXCFunctional
306 fxc = func.split(":")[-1]
307 try:
308 f_xc = LibXCFunctional(int(fxc), 1)
309 except ValueError:
310 f_xc = LibXCFunctional(fxc, 1)
311 fxc_params = dict(zip(f_xc.get_ext_param_names(), f_xc.get_ext_param_default_values()))
313 # Analyze the signature for implemented functionals
314 if func in IMPLEMENTED:
315 sig = inspect.signature(IMPLEMENTED[func])
316 fxc_params = {
317 param.name: param.default
318 for param in sig.parameters.values()
319 if param.default is not inspect.Parameter.empty
320 }
322 # Remove special names from the parsed parameters
323 for special in SPECIAL_NAMES:
324 if special in fxc_params:
325 del fxc_params[special]
327 # Append all parameters, warn if a parameter has been used before
328 for name in fxc_params:
329 if name in params:
330 log.warning(
331 f'External parameter "{name}" is used in the exchange AND correlation part. '
332 "It will be passed to both functionals if used in xc_params."
333 )
334 params[name] = fxc_params[name]
335 return params
338def get_zeta(n_spin):
339 """Calculate the relative spin polarization.
341 Args:
342 n_spin: Real-space electronic densities per spin channel.
344 Returns:
345 Relative spin polarization.
346 """
347 # If only one spin is given return an array of ones as if the density only is in one channel
348 if len(n_spin) == 1:
349 return np.ones_like(n_spin[0])
350 return (n_spin[0] - n_spin[1]) / (n_spin[0] + n_spin[1])
353def mock_xc(n, Nspin=1, **kwargs):
354 """Mock exchange-correlation functional with no effect (spin-paired).
356 Args:
357 n: Real-space electronic density.
358 Nspin: Number of spin states.
360 Keyword Args:
361 **kwargs: Throwaway arguments.
363 Returns:
364 Mock exchange-correlation energy density and potential.
365 """
366 zeros = np.zeros_like(n)
367 return zeros, np.array([zeros] * Nspin), None
370#: Map functional names with their respective implementation.
371IMPLEMENTED = {
372 i.__name__: i
373 for i in (
374 mock_xc,
375 gga_c_chachiyo,
376 gga_c_chachiyo_spin,
377 gga_c_pbe,
378 gga_c_pbe_spin,
379 gga_c_pbe_sol,
380 gga_c_pbe_sol_spin,
381 gga_x_chachiyo,
382 gga_x_chachiyo_spin,
383 gga_x_pbe,
384 gga_x_pbe_spin,
385 gga_x_pbe_sol,
386 gga_x_pbe_sol_spin,
387 lda_x,
388 lda_x_spin,
389 lda_c_pw,
390 lda_c_pw_spin,
391 lda_c_pw_mod,
392 lda_c_pw_mod_spin,
393 lda_c_vwn,
394 lda_c_vwn_spin,
395 lda_c_chachiyo,
396 lda_c_chachiyo_spin,
397 lda_c_chachiyo_mod,
398 lda_c_chachiyo_mod_spin,
399 lda_xc_corr_ksdt,
400 lda_xc_gdsmfb,
401 lda_xc_gdsmfb_spin,
402 lda_xc_ksdt,
403 lda_xc_ksdt_spin,
404 )
405}
407#: Map functional shorthands to the actual functional names.
408XC_MAP = {
409 # lda_x
410 "1": "lda_x",
411 "s": "lda_x",
412 "lda": "lda_x",
413 "slater": "lda_x",
414 # lda_c_vwn
415 "7": "lda_c_vwn",
416 "vwn": "lda_c_vwn",
417 "vwn5": "lda_c_vwn",
418 # lda_c_pw
419 "12": "lda_c_pw",
420 "pw": "lda_c_pw",
421 "pw92": "lda_c_pw",
422 # lda_c_pw_mod
423 "13": "lda_c_pw_mod",
424 "pwmod": "lda_c_pw_mod",
425 "pw92mod": "lda_c_pw_mod",
426 # gga_x_pbe
427 "101": "gga_x_pbe",
428 "pbex": "gga_x_pbe",
429 # gga_x_pbe_sol
430 "116": "gga_x_pbe_sol",
431 "pbesolx": "gga_x_pbe_sol",
432 # gga_c_pbe
433 "130": "gga_c_pbe",
434 "pbec": "gga_c_pbe",
435 # gga_c_pbe_sol
436 "133": "gga_c_pbe_sol",
437 "pbesolc": "gga_c_pbe_sol",
438 # lda_xc_ksdt
439 "259": "lda_xc_ksdt",
440 "ksdt": "lda_xc_ksdt",
441 # lda_c_chachiyo
442 "287": "lda_c_chachiyo",
443 "chachiyo": "lda_c_chachiyo",
444 # gga_x_chachiyo
445 "298": "gga_x_chachiyo",
446 "chachiyox": "gga_x_chachiyo",
447 # lda_c_chachiyo_mod
448 "307": "lda_c_chachiyo_mod",
449 "chachiyomod": "lda_c_chachiyo_mod",
450 # gga_c_chachiyo
451 "309": "gga_c_chachiyo",
452 "chachiyoc": "gga_c_chachiyo",
453 # lda_xc_corr_ksdt
454 "318": "lda_xc_corr_ksdt",
455 "corrksdt": "lda_xc_corr_ksdt",
456 # lda_xc_gdsmfb
457 "577": "lda_xc_gdsmfb",
458 "gdsmfb": "lda_xc_gdsmfb",
459}
461#: Dictionary of common functional aliases.
462ALIAS = {
463 "svwn": "slater,vwn5",
464 "svwn5": "slater,vwn5",
465 "spw92": "slater,pw92mod",
466 "pbe": "pbex,pbec",
467 "pbesol": "pbesolx,pbesolc",
468 "chachiyo": "chachiyox,chachiyoc",
469}